function [image,kData,Mss] = EPI(FOV,res,pd,T1,T2,TR,TE,alpha,TI,freqMap,BW,SNR)

%% Usage:  [image,kData] = EPI(FOV,res,pd,T1,T2,TR,TE,alpha,TI,freqMap,BW,SNR)
%%     Where:   FOV is field of view (in meters)
%%              res is resolution (isotropic, read and phase)
%%              pd is a proton density map (mm resolution)
%%              T1 is a T1 decay map (in milliseconds)
%%              T2 is a T2 decay map(in milliseconds)
%%              TR is the repetition time (milliseconds)
%%              TE is the echo time (milliseconds)
%%              alpha is the (optional) flip angle (degrees) - if not
%%                  given, ERnst is calculated
%%              TI is the (optional) inversion time (milliseconds)
%%              freqMap is an (optional) map of local B0 (shim imperfections, in Hz)
%%              BW is the (optional) bandwidth (in Hz) - 1/dwellTime
%%              SNR is the (optional) thermal SNR (added in k-space)
%%      The returned image is generated by simulating the imaging process,
%%      allowing for T2 decay and frequency offset errors after excitation.
%%      T1 contrast is calculated for steady-state, and inversion recovery
%%      if TI is given - returned as Mss.  
%%      Image is complex.  Simulated kData is also
%%      returned - image is simply the inverse FT of kData.

if ~exist('alpha','var') | isempty(alpha)
    avgT1 = mean(T1(find(T1)));
    alpha = acos(exp(-TR/avgT1))*180/pi; end

if ~exist('BW','var') | isempty(BW)
    BW = 100000; end

if ~exist('freqMap','var') | isempty(freqMap)
    freqMap = zeros(size(pd)); end

if size(pd,1) ~= size(pd,2)
    error('Programmer was too lazy.  Need pd size to be square.')
end

subplot(2,3,1); imagesc(pd); axis image; colorbar; title('pd')
subplot(2,3,2); imagesc(T1); axis image; colorbar; title('T1')
subplot(2,3,3); imagesc(T2); axis image; colorbar; title('T2')

% BW is given in Hz, and is 1/dwellTime (true BW, not Siemens Hz/pixel ...)
dwellTime = 1/BW; % this is the time it takes to read each data point
TE = TE/1000;
TR = TR/1000; % convert to seconds
T2 = T2/1000;
T1 = T1/1000;
alpha = alpha*pi/180; % convert from degrees to radians

% Calculate steady-state longitudinal magnetization, using TR, T1 and pd
mask = T1 > 0;
Mss = zeros(size(T1));
Mss(find(mask)) = (1-exp(-TR./T1(find(mask))))./(1-cos(alpha)*exp(-TR./T1(find(mask))));
Mss = Mss./max(Mss(:)); Mss(find(isnan(Mss))) = 0;
Mss = pd.*Mss;

% Do inversion recovery, if TI exists
warning off;
if exist('TI','var') & ~isempty(TI)
    TI = TI/1000;
    Mss = (Mss.*(1 - 2*exp(-TI./T1)));
    Mss(find(isnan(Mss))) = 0;
end

% Set gradients
Gro = BW/FOV; % sampling rate and field of view determine required gradient strength,
  % which is therefore calculated in units of Hz/m
check_Gro = Gro/42.58e3; % Gradients are actually in mT/m (G/cm)
  % conversion uses gyromagnetic ratio for the proton, which is 42.58 MHz/T
  % ... which is 42.58e-3 MHz/mT which is 42.58e3 Hz/mT
  % this conversion is checking to see if the requested read gradient is
  % too strong.  If so, one has to lower bandwidth (sample more slowly).
if check_Gro > 40 % slightly arbitrary limit, 40 mT/m
    error('Requested Gro too large; lower bandwidth or increase field of view.')
end
% disp(['G_ro = ' num2str(check_Gro,2) ' mT/m.'])
% disp(['dwell time (1/BW) = ' num2str(1e6*dwellTime) ' microseconds.'])
% for simplicity in the following calculations, convert the gradient from
% T/m to Hz/m
Gpe = Gro; % keep things square for simplicity

% Spit out some parameters to think about k-space
%
% Definition of k:  integral(gamma * G * t)
%  Letting G be a constant (ignoring ramps ...)
%  k = gamma*G*timeSinceExcitation;
%
% We calculated Gro as gamma*G, so k information is simple
% delta_k = Gro*dwellTime; %
% disp(['Delta k:  ' num2str(delta_k,3) ' m-1.'])
% disp(['  (Check that FOV = 1/delta_k)'])
% 
% k_max = Gro*dwellTime*res/2;
% disp(['k_max:  ' num2str(k_max,3) ' m-1.'])
% disp(['  (Check that resolution = 1/(2*k_max) )'])

Nro = res; Npe = res;
xMin = size(pd,1)/2 - 1; xMax = size(pd,1)/2;
xcoord = -xMin:xMax; % assuming mm resolution in input
[X Y] = meshgrid(xcoord*1e-3,xcoord*1e-3); % assume square

% Calculate timing info for EPI image  
sampleTimes = 0:dwellTime:Nro*dwellTime; % for one read line
minTE = Nro/2*dwellTime*Npe/2; % approximately
if TE < minTE
    error(['Requested TE is shorter than minimum (' num2str(1000*minTE,3) 'ms) allowed by bandwidth and resolution.'])
end
deadTime = TE - minTE;
T2decay = zeros(size(pd));

minTR = TE + minTE;
if TR < minTR
    error(['Requested TR is shorter than minimum (' num2str(1000*minTR,3) 'ms) allowed by TE, bandwidth and resolution.'])
end
readPolarity = ones([1 Npe]); readPolarity(2:2:end) = -1;
phasePolarity = 1; % can reverse up/down

% Define thermal noise relative to DC signal (from fully relaxed sample ...
% physiological noise is also defined relative to proton density ...
if exist('SNR')
    fudgeFactor = 50e-5/sqrt(dwellTime);
    thermalNoise = fudgeFactor*mean(pd(:))/sqrt(2)*complex(randn([Nro Npe]),randn([Nro Npe]));
    disp('adding uniform thermal noise')
else
    thermalNoise = zeros([Nro Npe]);
end

% "Acquire" EPI image
kData = zeros(res);
w=waitbar(0,'Stepping through PE  ...');
phase0 = zeros(size(pd)); % phase right after excitation
% Apply RF pulse
% Moving this out of PE loop is most significant change from FLASH to EPI
timeSinceExcitation = deadTime; % new RF pulse every repetition
Mt = Mss*sin(alpha*pi/180);
% Apply initial read and phase dephasing gradients
phaseDephased = phase0 - readPolarity(1)*Gro*dwellTime*X*Nro/2 - ...
    phasePolarity*Gpe*dwellTime*Y*Npe/2;
for iPE = 1:Npe
    waitbar(iPE/Npe,w)
    for it = 1:Nro
        phase = phaseDephased + readPolarity(iPE)*Gro*X*sampleTimes(it);
        % add in the evolution of spins due to off-resonance effects
        timeSinceExcitation = timeSinceExcitation + dwellTime; % should be = deadtime + sampleTimes(it)
        evolvedPhase = freqMap*timeSinceExcitation; % in Hz
        % calculate T2 weighting
        T2decay(find(T2)) = exp(-timeSinceExcitation./T2(find(T2)));
        % Calculate phase and magnitude of signal in sample
        weighted = Mt.*T2decay.*exp(2*pi*i*(phase + evolvedPhase));
        % Signal in coil is complex sum of signal in sample
        % Noise in coil is constant 
        if phasePolarity > 0
            if readPolarity(iPE) > 0
                kData(Npe + 1 - iPE,Nro + 1 - it) = sum(weighted(:)) + thermalNoise(it,iPE);
            else
                kData(Npe + 1 - iPE,it) = sum(weighted(:)) + thermalNoise(it,iPE);
            end
        else
            if readPolarity(iPE) > 0
                kData(iPE,Nro + 1 - it) = sum(weighted(:)) + thermalNoise(it,iPE);
            else
                kData(iPE,it) = sum(weighted(:)) + thermalNoise(it,iPE);
            end
        end
        % Display phase patten on a diagonal across k-space
        if iPE == it
            subplot(2,3,4); imagesc(real(weighted)); axis image; colorbar; drawnow
        end
    end
    % Apply phase blip for next line
    phaseDephased = phase + phasePolarity*dwellTime*Gpe*Y;
end    ; close(w) 

% We should recover the image by doing the inverse
% Fourier transform of the kData(t):
kData(find(isnan(kData))) = 0; % NaN's blow everything ...
image = fftshift(ifft2(fftshift(kData)));

subplot(2,3,4); imagesc(imresize(pd,res/size(pd,1))); axis image off; title('brain at sampled resolution')
subplot(2,3,5); imagesc(freqMap); axis image off; colorbar; title('off resonance map')
subplot(2,3,6); imagesc(abs(image)); axis image off; title(['T_R_O = ' num2str(1000*timeSinceExcitation,2) 'ms'])
colormap(gray)