% 2D RF-spoiled gradient echo (GRE), 128x128, single slice.
% The classic pulseq demo: slice-selective sinc excitation, phase encoding,
% RF spoiling, and reconstruction labels via seq.autoLabel.
% Adapted from pulseq's demoSeq/writeGradientEcho.m.

% set system limits
sys = mr.opts('MaxGrad', 22, 'GradUnit', 'mT/m', ...
    'MaxSlew', 120, 'SlewUnit', 'T/m/s', ...
    'rfRingdownTime', 20e-6, 'rfDeadTime', 100e-6, 'adcDeadTime', 10e-6);

seq = mr.Sequence(sys);         % Create a new sequence object
fov = 256e-3; Nx = 128; Ny = Nx; % Define FOV and resolution
alpha = 10;                     % flip angle
sliceThickness = 3e-3;          % slice
TR = 12e-3;                     % repetition time TR
TE = 5e-3;                      % echo time TE

rfSpoilingInc = 84;             % RF spoiling increment
roDuration = 3.2e-3;            % ADC duration

% autoLabel adds on-scanner reconstruction labels (LIN/PAR) by recomputing
% the full k-space trajectory. It is NOT needed to generate or view the
% sequence and dominates the runtime (~25 s here vs ~2 s without), so it is
% off by default. Set to true to include the labels.
doAutoLabel = false;

% Create alpha-degree slice selection pulse and gradient
[rf, gz] = mr.makeSincPulse(alpha*pi/180, sys, 'Duration', 3e-3, ...
    'SliceThickness', sliceThickness, 'apodization', 0.42, 'timeBwProduct', 4, ...
    'use', 'excitation');

% Define other gradients and ADC events
deltak = 1/fov;
gx = mr.makeTrapezoid('x', sys, 'FlatArea', Nx*deltak, 'FlatTime', roDuration);
adc = mr.makeAdc(Nx, sys, 'Duration', gx.flatTime, 'Delay', gx.riseTime);
gxPre = mr.makeTrapezoid('x', sys, 'Area', -gx.area/2, 'Duration', 1e-3);
gzReph = mr.makeTrapezoid('z', sys, 'Area', -gz.area/2, 'Duration', 1e-3);
phaseAreas = ((0:Ny-1)-Ny/2)*deltak;
gyPre = mr.makeTrapezoid('y', sys, 'Area', max(abs(phaseAreas)), 'Duration', mr.calcDuration(gxPre));
peScales = phaseAreas/gyPre.area;

% gradient spoiling
gxSpoil = mr.makeTrapezoid('x', sys, 'Area', 2*Nx*deltak);
gzSpoil = mr.makeTrapezoid('z', sys, 'Area', 4/sliceThickness);

% Calculate timing
delayTE = ceil((TE - mr.calcDuration(gxPre) - gz.fallTime - gz.flatTime/2 ...
    - mr.calcDuration(gx)/2)/seq.gradRasterTime)*seq.gradRasterTime;
delayTR = ceil((TR - mr.calcDuration(gz) - mr.calcDuration(gxPre) ...
    - mr.calcDuration(gx) - delayTE)/seq.gradRasterTime)*seq.gradRasterTime;
assert(all(delayTE >= 0));
assert(all(delayTR >= mr.calcDuration(gxSpoil, gzSpoil)));

rf_phase = 0;
rf_inc = 0;

% Loop over phase encodes and define sequence blocks
for i = 1:Ny
    rf.phaseOffset = rf_phase/180*pi;
    adc.phaseOffset = rf_phase/180*pi;
    rf_inc = mod(rf_inc + rfSpoilingInc, 360.0);
    rf_phase = mod(rf_phase + rf_inc, 360.0);
    %
    seq.addBlock(rf, gz);
    seq.addBlock(gxPre, mr.scaleGrad(gyPre, peScales(i)), gzReph);
    seq.addBlock(mr.makeDelay(delayTE));
    seq.addBlock(gx, adc);
    seq.addBlock(mr.makeDelay(delayTR), gxSpoil, mr.scaleGrad(gyPre, -peScales(i)), gzSpoil)
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

% add data labels to make image reconstruction on the scanner possible
if doAutoLabel
    seq.autoLabel('mirrorFourier', true, 'sortSlices', 'descending');
end

% prepare sequence export
seq.setDefinition('FOV', [fov fov sliceThickness]);
seq.setDefinition('Name', 'gre');

seq.write('gre.seq')       % Write to pulseq file