function RAPTOR_timelapse(simres,model,params,dit)
% RAPTOR_timelapse(simres,model,params,dit)
% dit: difference of it between curves
%%
set(gcf,'position',[800 400 900 600]);
collims = [0 0 1; 1 0 0; 0 0 0]; % colors

hax = multiaxes(3,4);
iPsi = model.psi.ind;
iTe = model.te.ind;
tgrid = params.grid.tgrid;
ntgrid = numel(tgrid);
colors = interp1([1,ntgrid/2,ntgrid],collims,1:ntgrid);

X = simres.X;
Xdot = simres.Xdot;
U = simres.U;

Psi = model.psi.Lam*X(iPsi,:);
%jpar = model.out.Cjpargauss*X(iPsi,:);
jpar = NaN*ones(numel(model.profiles.rhogaussN),size(X,2));
q = 1./(model.out.Ciota*X(iPsi,:));
Te = model.te.LamTe*X(iTe,:);
Upl = model.psi.Lam*Xdot(iPsi,:);
 
nt = sum(~isnan(X(1,:)));

jbs = NaN*jpar; jaux  = jbs; chie = jbs; T = jbs; %init
for it = 1:numel(simres.add.srcs)
    if isfield(simres.add.srcs(it),'jbs');
        jbs(:,it) = simres.add.srcs(it).jbs;
    else
        jbs(:,it) = NaN;
    end
    if isfield(simres.add.srcs(it),'jaux');
        jaux(:,it) = simres.add.srcs(it).jaux;
    else
        jaux(:,it) = NaN;
    end
    if isfield(simres.add.srcs(it),'chie');
        chie(:,it) = simres.add.srcs(it).chie;
    else
        chie(:,it) = NaN;
    end
    if isfield(simres.add.srcs(it),'T')
        T(:,it) = simres.add.srcs(it).T;
    else
        T(:,it) = NaN;
    end
end

joh = jpar-jbs-jaux;

Ibs = model.profiles.SIpar(end,:)*jbs;
Ioh = model.profiles.SIpar(end,:)*joh;
Ip =    model.out.CIp(end,:)*X(iPsi,:);

rtorN = model.profiles.rtorN;
rgaussN = model.profiles.rhogaussN;

for it = 1:dit:nt
    
    iax = 1;
    plot(hax(iax),rtorN,Psi(:,it),'color',colors(it,:)); hold(hax(iax),'on');
    labelff(hax(iax),'\psi','y'); labelff(hax(iax),'\rho','x');
    
    iax = 2;
    plot(hax(iax),rgaussN,jbs(:,it)/1e6,'color',colors(it,:)); hold(hax(iax),'on');
    labelff(hax(iax),'j_{BS} [MA/m^2]','y'); labelff(hax(iax),'\rho','x');

    iax = 7;
    plot(hax(iax),rtorN,q(:,it),'color',colors(it,:)); hold(hax(iax),'on');
    labelff(hax(iax),'q','y'); labelff(hax(iax),'\rho','x');
    set(hax(iax),'ylim',[0 10])
    
    iax = 10;
    plot(hax(iax),rgaussN,T(:,it),'color',colors(it,:)); hold(hax(iax),'on');
    labelff(hax(iax),'T = RB_\phi [Tm]','y'); labelff(hax(iax),'\rho','x');
    axis(hax(iax),'tight');
    
    iax = 11;
    plot(hax(iax),rgaussN,chie(:,it),'color',colors(it,:)); hold(hax(iax),'on');
    labelff(hax(iax),'\chi_e [m^2 s]','y'); labelff(hax(iax),'\rho','x');
    set(hax(iax),'ylim',[0 20]);
    
    iax = 8;
    plot(hax(iax),rtorN,Te(:,it)/1e3,'color',colors(it,:)); hold(hax(iax),'on');
    labelff(hax(iax),'Te [keV]','y'); labelff(hax(iax),'\rho','x');
    
    iax = 4;
    plot(hax(iax),rtorN,Upl(:,it),'color',colors(it,:)); hold(hax(iax),'on');
    labelff(hax(iax),'Upl [V]','y'); labelff(hax(iax),'\rho','x');
    axis(hax(iax),'tight');

    iax = 5;
    plot(hax(iax),rgaussN,joh(:,it)/1e6,'color',colors(it,:)); hold(hax(iax),'on');
    labelff(hax(iax),'j_{OH} [MA/m^2]','y'); labelff(hax(iax),'\rho','x');
    axis(hax(iax),'tight');

end

% time traces
gy = 0.7*[1 1 1];
iax = 12 ;
plot(hax(iax),tgrid,Te(1,:)/1e3,'color',gy);
hold(hax(iax),'on');
labelff(hax(iax),'Te0 [keV]','y'); labelff(hax(iax),'t[s]','x');
axis(hax(iax),'tight');


iax = 9;
plot(hax(iax),tgrid,1e-3*Ioh,'color',gy); hold(hax(iax),'on');
plot(hax(iax),tgrid,1e-3*Ibs,'color',gy);
plot(hax(iax),tgrid,1e-3*Ip,'color',gy);
axis(hax(iax),'tight');

iax = 6 ;
plot(hax(iax),tgrid,U(2,:)/1e6,'color',gy);
hold(hax(iax),'on');
labelff(hax(iax),'Paux [MW]','y'); labelff(hax(iax),'t[s]','x');
axis(hax(iax),'tight');

iax=3;
switch params.tran.psi.BCtype
    case 'Psioh'
        Voh = diff(U(1,:))./diff(tgrid);    
        lab = 'V_{oh}/10 [V],\Psi_{oh}[Wb]';
        
        plot(hax(iax),tgrid(1:end-1),Voh/10,'color',gy); hold(hax(iax),'on');
        plot(hax(iax),tgrid,U(1,:),'color',gy);
        axis(hax(iax),'tight');

        hold(hax(iax),'on');
        labelff(hax(iax),lab,'y');
        
    case 'Ip'
end


labelff(hax(iax),'t[s]','x');

set(hax([3,6,9]),'xlim',[tgrid(1),tgrid(end)]);

for it=1:nt-1
%     plot(hax(iax),tgrid([it,it+1]),1e-3*Ioh(:,[it,it+1]),'color',colors(it,:));
%     plot(hax(iax),tgrid([it,it+1]),1e-3*Ibs(:,[it,it+1]),'color',colors(it,:));
%     plot(hax(iax),tgrid([it,it+1]),1e-3*Ip(:,[it,it+1]),'color',colors(it,:));
%     
if rem(it-1,dit)==0
    iax = 9;
    hold(hax(iax),'on');
    plot(hax(iax),tgrid(it),1e-3*Ip(:,it),'color',colors(it,:),'marker','o','markersize',6);
    plot(hax(iax),tgrid(it),1e-3*Ibs(:,it),'color',colors(it,:),'marker','s','markersize',6);
    plot(hax(iax),tgrid(it),1e-3*Ioh(:,it),'color',colors(it,:),'marker','^','markersize',6);
end
    hl=labelff(hax(iax),'Ip (o), Ibs (s), Ioh(^) [kA]','y'); set(hl,'interpreter','none')
    labelff(hax(iax),'t[s]','x');
    
    end






%%