function estQALAS =  modelQALAS(R1,R2,QALASinfo,alpha,beta)
% estQALAS =  modelQALAS(R1,R2,QALASinfo,B1map)
% estQALS is a 4 dimentional matirx with the 5 differnet dynamics in the
% 4th dimention
% 
% Writen by Anders Tisell 2024-12-10

t2prep = QALASinfo.timeT2prep;

tdelay = QALASinfo.timeDelay;
Ntfe = QALASinfo.TurboFactor;
FA = QALASinfo.FA .*alpha;
IRflip = pi*beta;
EchoSpace = QALASinfo.EchoSpace;
% dymmy shot to get system in equlib.  

estQALAS = repmat(ones(size(R1),'single'),1,1,1,5);

for dymmyIDX = 1:2
    
    estQALAS(:,:,:,1) = T2prep(R1,R2,t2prep,Ntfe,FA,EchoSpace,estQALAS(:,:,:,5),tdelay);
    estQALAS(:,:,:,2) = IRprep(R1,tdelay,Ntfe,FA,IRflip,EchoSpace,estQALAS(:,:,:,1));
    for dynIDX = 3:5
        estQALAS(:,:,:,dynIDX) = nextDyn(R1,tdelay,Ntfe,FA,EchoSpace,estQALAS(:,:,:,dynIDX-1));        
    end
end

end

function S1 = T2prep(R1,R2,t2prep,N,FA,tTFE,S5,tdelay)

%S5end = (1-(cos(FA)).^N.*exp(-R1.*tread)) + S5.*((cos(FA)).^N.*exp(-R1.*tread));
%p = ones(size(R1)).*cos(FA);
p = cos(FA);
q = exp(-R1.*tTFE);

S5end = (1-q)./(1-p.*q).*(1-(p.*q).^N) + S5.*(p.*q).^N;

SpostT2prep = S5end.*exp(-t2prep.*R2);
S1 = 1 + (SpostT2prep-1).*exp(-R1.*(tdelay-t2prep));

end

function S = IRprep(R1,tdelay,N,FA,IRflip,tTFE,Sin)

%p = ones(size(R1)).*cos(FA);
p = cos(FA);
q = exp(-R1.*tTFE);
r = cos(IRflip);
SendTFE = (1-q)./(1-p.*q).*(1-(p.*q).^N) + Sin.*(p.*q).^N;

S =  1  - (1-r.*SendTFE).*exp(-R1.*tdelay);  % IR pulse

%S = 1 + (Send-1).*((cos(FA)).^N.*exp(-R1.*tread));

end

function S = nextDyn(R1,tdelay,N,FA,tTFE,Sin)

%p = ones(size(R1)).*cos(FA);
p = cos(FA);
q = exp(-R1.*tTFE);
   
    %S = 1 + (Send-1).*((cos(FA)).^N.*exp(-R1.*tread));  
    SendTFE = (1-q)./(1-p.*q).*(1-(p.*q).^N) + Sin.*(p.*q).^N;
     S =  1 + (SendTFE-1).*exp(-R1.*tdelay);
end