% Echo-planar imaging (EPI), 64x64, 3 slices.
% A basic educational EPI without ramp sampling: one excitation per slice,
% then the whole k-space plane in a single train of alternating readouts
% with phase blips. Adapted from pulseq's demoSeq/writeEpi.m.

seq = mr.Sequence();            % Create a new sequence object
fov = 220e-3; Nx = 64; Ny = 64; % Define FOV and resolution
thickness = 3e-3;               % slice thickness
Nslices = 3;

% Set system limits
lims = mr.opts('MaxGrad', 32, 'GradUnit', 'mT/m', ...
               'MaxSlew', 130, 'SlewUnit', 'T/m/s', ...
               'rfRingdownTime', 30e-6, 'rfDeadTime', 100e-6);

% Create 90 degree slice selection pulse and gradient
[rf, gz] = mr.makeSincPulse(pi/2, 'system', lims, 'Duration', 3e-3, ...
    'SliceThickness', thickness, 'apodization', 0.5, 'timeBwProduct', 4, ...
    'use', 'excitation');

% Define other gradients and ADC events
deltak = 1/fov;
kWidth = Nx*deltak;
dwellTime = 4e-6;
readoutTime = Nx*dwellTime;
flatTime = ceil(readoutTime*1e5)*1e-5; % round-up to the gradient raster
gx = mr.makeTrapezoid('x', lims, 'Amplitude', kWidth/readoutTime, 'FlatTime', flatTime);
adc = mr.makeAdc(Nx, 'Duration', readoutTime, ...
    'Delay', gx.riseTime + flatTime/2 - (readoutTime - dwellTime)/2);

% Pre-phasing gradients
preTime = 8e-4;
gxPre = mr.makeTrapezoid('x', lims, 'Area', -gx.area/2, 'Duration', preTime);
gzReph = mr.makeTrapezoid('z', lims, 'Area', -gz.area/2, 'Duration', preTime);
gyPre = mr.makeTrapezoid('y', lims, 'Area', -Ny/2*deltak, 'Duration', preTime);

% Phase blip in shortest possible time
dur = ceil(2*sqrt(deltak/lims.maxSlew)/10e-6)*10e-6;
gy = mr.makeTrapezoid('y', lims, 'Area', deltak, 'Duration', dur);

% Define sequence blocks
for s = 1:Nslices
    rf.freqOffset = gz.amplitude*thickness*(s-1-(Nslices-1)/2);
    seq.addBlock(rf, gz);
    seq.addBlock(gxPre, gyPre, gzReph);
    for i = 1:Ny
        seq.addBlock(gx, adc);          % Read one line of k-space
        seq.addBlock(gy);               % Phase blip
        gx.amplitude = -gx.amplitude;   % Reverse polarity of read gradient
    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

seq.setDefinition('FOV', [fov fov thickness]);
seq.setDefinition('Name', 'epi');

seq.write('epi.seq');