1% set system limits (slew rate 130 and max_grad 30 work on Prisma)
2sys = mr.opts('MaxGrad', 28, 'GradUnit', 'mT/m', ...
3 'MaxSlew', 120, 'SlewUnit', 'T/m/s', 'rfRingdownTime', 10e-6, ...
4 'rfDeadTime', 100e-6, 'adcDeadTime', 10e-6);
6seq=mr.Sequence(sys); % Create a new sequence object
7fov=240e-3; Nx=240; % Define FOV and resolution
8alpha=5; % flip angle
9sliceThickness=6e-3; % slice
10Nr=256; % number of radial spokes
11Ndummy=10; % number of dummy scans
12delta= pi / Nr; % angular increment; try golden angle pi*(3-5^0.5) or 0.5 of it
13ro_dur=1200e-6; % RO duration
14ro_os=2; % readout oversampling
15ro_spoil=0.5; % additional k-max excursion for RO spoiling
16sl_spoil=1.5; % spoil area compared to the slice thickness
18% TE & TR are as short as possible derived from the above parameters and
19% the system specs below
21% more in-depth parameters
22rfSpoilingInc=84; % RF spoiling increment
23rfTimeBwProd = 2; % time-bandwidth product for the RF pulse
25% Create alpha-degree slice selection pulse and gradient
26[rf, gz, gzReph] = mr.makeSincPulse(alpha*pi/180,'Duration',400e-6,...
27 'SliceThickness',sliceThickness,'apodization',0.5,'timeBwProduct',rfTimeBwProd,'system',sys,...
28 'use','excitation');
29% gradient spoiling in slice direction
30if sl_spoil>0
31 sp_area_needed=sl_spoil/sliceThickness*rfTimeBwProd-gz.flatArea/2;
32 %gzSpoil=mr.makeTrapezoid('z','Area',sp_area_needed,'system',sys,'Delay',gx.riseTime+gx.flatTime);
33 gzSpoil=mr.makeExtendedTrapezoidArea('z', 0, gz.amplitude, sp_area_needed, sys);
34 gz=mr.makeExtendedTrapezoid('z','times',[0, gz.flatTime,gz.flatTime+gz.fallTime],'amplitudes',[gz.amplitude,gz.amplitude,0],'system',sys);
35 gz.delay=gzSpoil.shape_dur;
36 gz=mr.addGradients({gz, gzSpoil}, 'system', sys);
37 rf.delay=max(gzSpoil.shape_dur,sys.rfDeadTime);
38 gz.delay=rf.delay-gzSpoil.shape_dur;
39end
40% join rephaser
41gzReph.delay=mr.calcDuration(gz);
42gzComb=mr.addGradients({gz, gzReph}, 'system', sys);
43[gzSpoil,gz]=mr.splitGradientAt(gzComb,rf.delay-sys.rfDeadTime);
44gz.delay=gz.delay-mr.calcDuration(gzSpoil);
45rf.delay=sys.rfDeadTime;
46clear gzComb;
48% Define other gradients and ADC events
49deltak=1/fov;
50gx = mr.makeTrapezoid('x','Amplitude',Nx*deltak/ro_dur,'FlatTime',ceil(ro_dur/sys.gradRasterTime)*sys.gradRasterTime,'system',sys);
51adc = mr.makeAdc(Nx*ro_os,'Duration',ro_dur,'Delay',gx.riseTime,'system',sys);
52gxPre = mr.makeTrapezoid('x','Area',-gx.amplitude*(ro_dur/Nx/ro_os*(Nx*ro_os/2-0.5)+0.5*gx.riseTime),'system',sys); % 0.5 is necessary to acount for the Siemens sampling in the center of the dwell periods
53%
54[gxPre,~]=mr.align('right', gxPre, 'right', gz);
55addDelay=mr.calcDuration(rf)-gxPre.delay;
56if addDelay>0
57 gxPre.delay = gxPre.delay+ceil(addDelay/sys.gradRasterTime)*sys.gradRasterTime;
58end
60% gradient spoiling in slice direction
61if ro_spoil>0
62 %ro_spoil_area=(gx.area-gx.flatArea)/2;
63 ro_add_time=ceil(((gx.area/Nx*(Nx/2+1)*ro_spoil)/gx.amplitude)/sys.gradRasterTime)*sys.gradRasterTime;
64 gx.flatTime=gx.flatTime+ro_add_time; % careful, areas stored in the object are now wrong
65end
67% join slice spoiler with the slice selection
68%if (rf.delay>mr.calcDuration()) no, this does not work to be really optimal we need a new function with start, stop and area
69% could be done with mr.makeExtendedTrapezoidArea()
71% Calculate timing
72% TODO: just calculate actual TE and TR here
74% start the sequence
75rf_phase=0;
76rf_inc=0;
77TR=0;
78seq.addBlock(gzSpoil);
79[gzSpoil,~,~]=mr.align('right',gzSpoil,'left',gx,mr.makeDelay(ceil((adc.delay+adc.dwell*adc.numSamples)/sys.gradRasterTime)*sys.gradRasterTime+gzSpoil.shape_dur));
80for i=(1-Ndummy):Nr
81 rf.phaseOffset=rf_phase/180*pi;
82 adc.phaseOffset=rf_phase/180*pi;
83 rf_inc=mod(rf_inc+rfSpoilingInc, 360.0);
84 rf_phase=mod(rf_phase+rf_inc, 360.0);
85 %
86 phi=delta*(i-1);
87 seq.addBlock(mr.rotate('z',phi,rf,gz,gxPre));
88 if (i>0)
89 seq.addBlock(mr.rotate('z',phi,gx,adc,gzSpoil));
90 else
91 seq.addBlock(mr.rotate('z',phi,gx,gzSpoil));
92 end
93 if TR<=0
94 TR=seq.duration;
95 end
96end
97% dummy slice select to ramp down to 0 Z spoil gradient
98seq.addBlock(gz);
100%% check whether the timing of the sequence is correct
101[ok, error_report]=seq.checkTiming;
103if (ok)
104 fprintf('Timing check passed successfully\n');
105else
106 fprintf('Timing check failed! Error listing follows:\n');
107 fprintf([error_report{:}]);
108 fprintf('\n');
109end
111%% plot, etc
112%return;
113%% trajectory calculation
114%[ktraj_adc, ktraj, t_excitation, t_refocusing, t_adc] = seq.calculateKspace();
116% plot k-spaces
118%
119seq.setDefinition('FOV', [fov fov sliceThickness]);
120seq.setDefinition('Name', 'gre_rad');
122seq.write('fast_gre_rad.seq') % Write to pulseq file