1% a basic UTE-like sequence
2% achieves "TE" below 100 us
4% set system limits
5sys = mr.opts('MaxGrad', 28, 'GradUnit', 'mT/m', ...
6 'MaxSlew', 170, 'SlewUnit', 'T/m/s', 'rfRingdownTime', 0e-6, ...
7 'rfDeadTime', 100e-6, 'adcDeadTime', 10e-6);
9seq=mr.Sequence(sys); % Create a new sequence object
10fov=240e-3; Nx=240; % Define FOV and resolution
11alpha=10; % flip angle
12sliceThickness=5e-3; % slice
13TR=20e-3; % TR
14%Nr=round(Nx*pi/2); % number of radial spokes
15Nr=Nx*2; % number of radial spokes
16Ndummy=20; % number of dummy scans
17delta= 2* pi / Nr; % angular increment; try golden angle pi*(3-5^0.5) or 0.5 of it
18rf_duration=0.5e-3; % duration of the excitation pulse
19ro_duration=0.720e-3; % read-out time: controls RO bandwidth and T2-blurring
20ro_os=2; % oversampling
21minRF_to_ADC_time=70e-6; % the parameter wich defines TE together with ro_discard
22ro_discard=0; % dummy ADC samples to discard (due to ADC filter)
23ro_spoil=1; % extend RO to achieve spoiling
25% more in-depth parameters
26% RF spoiling increment = 84° for smoother transient decay, https://doi.org/10.1002/mrm.1910350216, 169° for diffusion independent rf spoiling in steady-state https://doi.org/10.1371/journal.pone.0324455
27rfSpoilingInc=84; % RF spoiling increment
29%% Create alpha-degree slice selection pulse and gradient
30[rf, gz] = mr.makeSincPulse(alpha*pi/180,'Duration',rf_duration,...
31 'SliceThickness',sliceThickness,'apodization',0.5,'timeBwProduct',2,...
32 'centerpos',1,'system',sys,'use','excitation');
34% resample the RF pulse to the ramp
35gza=[0 1 1 0];
36gzt=cumsum([0 gz.riseTime gz.flatTime gz.fallTime]);
37gzas_0=interp1(gzt+gz.delay,gza,rf.t+rf.delay);
38rft_1=[sys.rfRasterTime:sys.rfRasterTime:rf_duration+0.5*gz.fallTime];
39gzas_1=interp1(gzt+gz.delay,gza,rft_1+rf.delay+gz.fallTime*0.5);
40gzas_1(~isfinite(gzas_1))=0; % we are getting a NaN sometimes
41kzs_0=cumsum(gzas_0);
42kzs_1=cumsum(gzas_1);
43kzs_0=kzs_0-max(kzs_0);
44kzs_1=kzs_1-max(kzs_1);
45rfs_1=diff([0 interp1(kzs_0,cumsum(rf.signal),kzs_1)]);
46rf.t=rft_1(1:end-1); % remove the 0 at the end
47rf.signal=rfs_1(1:end-1); % remove the 0 at the end
48rf.shape_dur=length(rf.signal)*sys.rfRasterTime;
49gz.flatTime=ceil((gz.flatTime-gz.fallTime*0.5)/sys.gradRasterTime)*sys.gradRasterTime;
50rf.delay=mr.calcDuration(rf,gz)-rf.shape_dur; % fix the possible time shift due to the rounding-up step above
51rf.center=rf.t(end); % or rf.shape_dur?
53% Align RO assymmetry to ADC samples
54Nxo=round(ro_os*Nx);
55% Define other gradients and ADC events
56deltak=1/fov/2;
57ro_area=Nx*deltak;
58gx = mr.makeTrapezoid('x','FlatArea',ro_area,'FlatTime',ro_duration,'system',sys);
59adc_dur=floor(gx.flatTime/Nxo*1e7)*1e-7*Nxo; % round down dwell time to 100ns (Siemens ADC raster)
60adc = mr.makeAdc(Nxo,'Duration',adc_dur,'system',sys);
62% ro-spoiling
63gx.flatTime=gx.flatTime*ro_spoil;
65% Calculate timing
66%ceil((TE - mr.calcDuration(gxPre) - gz.fallTime - gz.flatTime/2 ...
67% - mr.calcDuration(gx)/2)/seq.gradRasterTime)*seq.gradRasterTime;
68% calculate actual achieved TE
69TE = ceil((minRF_to_ADC_time + adc.dwell*ro_discard)/seq.gradRasterTime)*seq.gradRasterTime;
70delayTR=ceil((TR - mr.calcDuration(gz) ...
71 - mr.calcDuration(gx) - TE)/seq.gradRasterTime)*seq.gradRasterTime;
73fprintf('TE= %d us; delay in TR:= %d us\n', round(TE*1e6), floor(delayTR*1e6));
75% set up timing
76gx.delay=mr.calcDuration(gz)+TE;
77adc.delay=floor((gx.delay-adc.dwell*0.5-adc.dwell*ro_discard)/sys.gradRasterTime)*sys.gradRasterTime; % take into accout 0.5 samples ADC shift
79rf_phase=0;
80rf_inc=0;
82if Ndummy>0
83 seq.addBlock(mr.makeLabel('SET','ONCE', 1)); % label the few following scans as preparing scans
84end
86for i=(-Ndummy+1):Nr
87 if Ndummy>0 && i==1
88 seq.addBlock(mr.makeLabel('SET','ONCE', 0)); % remove the preparing scan label
89 end
90 for c=1:2
91 rf.phaseOffset=rf_phase/180*pi;
92 adc.phaseOffset=rf_phase/180*pi;
93 rf_inc=mod(rf_inc+rfSpoilingInc, 360.0);
94 rf_phase=mod(rf_phase+rf_inc, 360.0);
95 % UTE: alternate GZ
96 gz.amplitude=-gz.amplitude;
97 %
98 phi=delta*(i-1);
99 grc=gx; grs=gx; grc.amplitude=gx.amplitude*cos(phi); grs.amplitude=gx.amplitude*sin(phi); grs.channel='y';
100 if (i>0)
101 seq.addBlock(rf,gz,grc,grs,adc);
102 else
103 seq.addBlock(rf,gz,grc,grs);
104 end
105 seq.addBlock(delayTR);
106 end
107end
110%% check whether the timing of the sequence is correct
111[ok, error_report]=seq.checkTiming;
113if (ok)
114 fprintf('Timing check passed successfully\n');
115else
116 fprintf('Timing check failed! Error listing follows:\n');
117 fprintf([error_report{:}]);
118 fprintf('\n');
119end
121%% export
122seq.setDefinition('FOV', [fov fov sliceThickness]);
123seq.setDefinition('Name', 'ute_rs');
125seq.write('ute_rs.seq'); % Write to pulseq file