Contents
seq=mr.Sequence();
fov=224e-3; Nx=112; Ny=112;
thickness=2e-3;
Nslices=3;
bFactor=1000;
TE=100e-3;
pe_enable=1;
ro_os=1;
readoutTime=6.3e-4;
partFourierFactor=0.75;
tRFex=3e-3;
tRFref=3e-3;
lims = mr.opts('MaxGrad',38,'GradUnit','mT/m',...
'MaxSlew',180,'SlewUnit','T/m/s',...
'rfRingdownTime', 10e-6, 'rfDeadtime', 100e-6);
B0=2.89;
sat_ppm=-3.45;
sat_freq=sat_ppm*1e-6*B0*lims.gamma;
rf_fs = mr.makeGaussPulse(110*pi/180,'system',lims,'Duration',8e-3,...
'bandwidth',abs(sat_freq),'freqOffset',sat_freq);
gz_fs = mr.makeTrapezoid('z',lims,'delay',mr.calcDuration(rf_fs),'Area',1/1e-4);
[rf, gz, gzReph] = mr.makeSincPulse(pi/2,'system',lims,'Duration',tRFex,...
'SliceThickness',thickness,'apodization',0.5,'timeBwProduct',4);
[rf180, gz180] = mr.makeSincPulse(pi,'system',lims,'Duration',tRFref,...
'SliceThickness',thickness,'apodization',0.5,'timeBwProduct',4,'PhaseOffset',pi/2,'use','refocusing');
[~, gzr_t, gzr_a]=mr.makeExtendedTrapezoidArea('z',gz180.amplitude,0,-gzReph.area+0.5*gz180.amplitude*gz180.fallTime,lims);
gz180n=mr.makeExtendedTrapezoid('z','system',lims,'times',[0 gz180.riseTime gz180.riseTime+gz180.flatTime+gzr_t]+gz180.delay, 'amplitudes', [0 gz180.amplitude gzr_a]);
trig=mr.makeDigitalOutputPulse('osc0','duration', 100e-6);
deltak=1/fov;
kWidth = Nx*deltak;
blip_dur = ceil(2*sqrt(deltak/lims.maxSlew)/10e-6/2)*10e-6*2;
gy = mr.makeTrapezoid('y',lims,'Area',-deltak,'Duration',blip_dur);
extra_area=blip_dur/2*blip_dur/2*lims.maxSlew;
gx = mr.makeTrapezoid('x',lims,'Area',kWidth+extra_area,'duration',readoutTime+blip_dur);
actual_area=gx.area-gx.amplitude/gx.riseTime*blip_dur/2*blip_dur/2/2-gx.amplitude/gx.fallTime*blip_dur/2*blip_dur/2/2;
gx.amplitude=gx.amplitude/actual_area*kWidth;
gx.area = gx.amplitude*(gx.flatTime + gx.riseTime/2 + gx.fallTime/2);
gx.flatArea = gx.amplitude*gx.flatTime;
adcDwellNyquist=deltak/gx.amplitude/ro_os;
adcDwell=floor(adcDwellNyquist*1e7)*1e-7;
adcSamples=floor(readoutTime/adcDwell/4)*4;
adc = mr.makeAdc(adcSamples,'Dwell',adcDwell,'Delay',blip_dur/2);
time_to_center=adc.dwell*((adcSamples-1)/2+0.5);
adc.delay=round((gx.riseTime+gx.flatTime/2-time_to_center)*1e6)*1e-6;
gy_parts = mr.splitGradientAt(gy, blip_dur/2, lims);
[gy_blipup, gy_blipdown]=mr.align('right',gy_parts(1),'left',gy_parts(2),gx);
gy_blipdownup=mr.addGradients({gy_blipdown, gy_blipup}, lims);
gy_blipup.waveform=gy_blipup.waveform*pe_enable;
gy_blipdown.waveform=gy_blipdown.waveform*pe_enable;
gy_blipdownup.waveform=gy_blipdownup.waveform*pe_enable;
Ny_pre=round(partFourierFactor*Ny/2-1);
Ny_post=round(Ny/2+1);
Ny_meas=Ny_pre+Ny_post;
gxPre = mr.makeTrapezoid('x',lims,'Area',-gx.area/2);
gyPre = mr.makeTrapezoid('y',lims,'Area',Ny_pre*deltak);
[gxPre,gyPre]=mr.align('right',gxPre,'left',gyPre);
gyPre = mr.makeTrapezoid('y',lims,'Area',gyPre.area,'Duration',mr.calcDuration(gxPre,gyPre));
gyPre.amplitude=gyPre.amplitude*pe_enable;
durationToCenter = (Ny_pre+0.5)*mr.calcDuration(gx);
rfCenterInclDelay=rf.delay + mr.calcRfCenter(rf);
rf180centerInclDelay=rf180.delay + mr.calcRfCenter(rf180);
delayTE1=ceil((TE/2 - mr.calcDuration(rf,gz) + rfCenterInclDelay - rf180centerInclDelay)/lims.gradRasterTime)*lims.gradRasterTime;
delayTE2tmp=ceil((TE/2 - mr.calcDuration(rf180,gz180n) + rf180centerInclDelay - durationToCenter)/lims.gradRasterTime)*lims.gradRasterTime;
assert(delayTE1>=0);
gxPre.delay=0;
gyPre.delay=0;
delayTE2=delayTE2tmp-mr.calcDuration(gxPre,gyPre);
[gxPre,gyPre]=mr.align('right',gxPre,'left',gyPre);
assert(delayTE2>=0);
small_delta=delayTE2-ceil(lims.maxGrad/lims.maxSlew/lims.gradRasterTime)*lims.gradRasterTime;
big_delta=delayTE1+mr.calcDuration(rf180,gz180n);
g=sqrt(bFactor*1e6/bFactCalc(1,small_delta,big_delta));
gr=ceil(g/lims.maxSlew/lims.gradRasterTime)*lims.gradRasterTime;
gDiff=mr.makeTrapezoid('z','amplitude',g,'riseTime',gr,'flatTime',small_delta-gr,'system',lims);
assert(mr.calcDuration(gDiff)<=delayTE1);
assert(mr.calcDuration(gDiff)<=delayTE2);
for s=1:Nslices
seq.addBlock(rf_fs,gz_fs);
rf.freqOffset=gz.amplitude*thickness*(s-1-(Nslices-1)/2);
rf180.freqOffset=gz180.amplitude*thickness*(s-1-(Nslices-1)/2);
seq.addBlock(rf,gz,trig);
seq.addBlock(mr.makeDelay(delayTE1),gDiff);
rf180.freqOffset=gz180.amplitude*thickness*(s-1-(Nslices-1)/2);
seq.addBlock(rf180,gz180n);
seq.addBlock(mr.makeDelay(delayTE2),gDiff);
seq.addBlock(gxPre,gyPre);
for i=1:Ny_meas
if i==1
seq.addBlock(gx,gy_blipup,adc);
elseif i==Ny_meas
seq.addBlock(gx,gy_blipdown,adc);
else
seq.addBlock(gx,gy_blipdownup,adc);
end
gx.amplitude = -gx.amplitude;
end
end
check whether the timing of the sequence is correct
[ok, error_report]=seq.checkTiming;
if (ok)
fprintf('Timing check passed successfully\n');
else
fprintf('Timing check failed! Error listing follows:\n');
fprintf([error_report{:}]);
fprintf('\n');
end
Timing check passed successfully
do some visualizations
seq.plot();
[ktraj_adc, ktraj, t_excitation, t_refocusing, t_adc] = seq.calculateKspace();
time_axis=(1:(size(ktraj,2)))*lims.gradRasterTime;
figure; plot(time_axis, ktraj');
hold on; plot(t_adc,ktraj_adc(1,:),'.');
figure; plot(ktraj(1,:),ktraj(2,:),'b');
axis('equal');
hold on;plot(ktraj_adc(1,:),ktraj_adc(2,:),'r.');
prepare the sequence output for the scanner
seq.setDefinition('FOV', [fov fov thickness*Nslices]);
seq.setDefinition('Name', 'epi');
seq.write('epidiff_rs.seq');
very optional slow step, but useful for testing during development e.g. for the real TE, TR or for staying within slewrate limits
rep = seq.testReport;
fprintf([rep{:}]);
function b=bFactCalc(g, delta, DELTA)
sigma=1;
kappa_minus_lambda=1/3-1/2;
b= (2*pi * g * delta * sigma)^2 * (DELTA + 2*kappa_minus_lambda*delta);
end
Number of blocks: 312
Number of events:
RF: 9
Gx: 297
Gy: 297
Gz: 15
ADC: 294
Delay: 6
Sequence duration: 0.465330s
TE: 0.100003s
TR: 0.012020s
Flip angle: 90.00°
Flip angle: 110.00°
Flip angle: 180.00°
Unique k-space positions (a.k.a. columns, rows, etc): 124
Unique k-space positions (a.k.a. columns, rows, etc): 98
Dimensions: 2
Spatial resolution: 2.02 mm
Spatial resolution: 2.00 mm
Repetitions/slices/contrasts: 3 range: [3 3]
Cartesian encoding trajectory detected
Block timing check passed successfully
Max. Gradient: 1298246 Hz/m == 30.49 mT/m
Max. Gradient: 610119 Hz/m == 14.33 mT/m
Max. Gradient: 1615509 Hz/m == 37.94 mT/m
Max. Slew Rate: 7.40741e+09 Hz/m/s == 173.98 T/m/s
Max. Slew Rate: 7.62649e+09 Hz/m/s == 179.13 T/m/s
Max. Slew Rate: 7.59687e+09 Hz/m/s == 178.43 T/m/s
Max. Absolute Gradient: 1615509 Hz/m == 37.94 mT/m
Max. Absolute Slew Rate: 1.03638e+10 Hz/m/s == 243.42 T/m/s