function [out,exitflag] = RAPTOR_out(varargin)
% extract RAPTOR raw simulation results to post-processed output structure
%
% Documentation of the output structure is given in this script as
% well. Where documentation is missing, see the eval_XX function
% documentation.
%
% All units are in SI, temperatures in eV
%
% [out,exitflag] = RAPTOR_out(simres,model,params); % all data for a full simulation
% timeslice = RAPTOR_out(simres,model,params,it); % timeslice for full simulation
% timeslice = RAPTOR_out(xk,gk,vk,xdot,uk,it,stap,geop,trap,model,params);
%
% exitflag = 0 if everything is ok, else exitflag=1;

if isstruct(varargin{1})
  % called for whole simulation
  if nargin == 3
    [simres,model,params] = deal(varargin{:});
    outmode = 'fullsim';
  elseif nargin==4
    [simres,model,params,itt] = deal(varargin{:});
    outmode = 'simtimeslice';
  end
else
  % called for one time step
  assert(numel(varargin)==11,'wrong number of inputs')
  [xk,gk,vk,xdot,uk,itt,stap,geop,trap,model,params] = deal(varargin{:});
  outmode = 'onetimeslice';
end

%%
switch outmode
  case 'fullsim'
    exitflag = max(simres.exitflags);
    % check if simulation result makes sense
    
    t_ok = 1:find(all(~isnan(simres.X)),1,'last'); % ok time indices up to first NaN
    
    X = simres.X(:,t_ok);
    Xdot = simres.Xdot(:,t_ok);
    G = simres.G(:,t_ok);
    V = simres.V(:,t_ok);
    U = simres.U(:,t_ok);
    traps = simres.traps;
    staps = simres.staps;
    geops = simres.geops;
    itt = t_ok;
    
  case 'simtimeslice'
    X = simres.X(:,itt);
    Xdot = simres.Xdot(:,itt);
    G = simres.G(:,itt);
    V = simres.V(:,itt);
    U = simres.U(:,itt);
    traps = simres.traps{itt};
    staps = simres.staps{itt};
    geops = simres.geops{itt};
    
    t_ok = 1;
    
  case 'onetimeslice'
    % single time point
    X = xk;
    G = gk;
    V = vk;
    U = uk;
    Xdot = xdot;
    traps = {trap};
    staps = {stap};
    geops = {geop};
    t_ok = 1; % only time index
end
%%

%% Grids
out.it = itt;
if ~strcmp(outmode,'onetimeslice') && isfield(simres,'time')
  out.time = simres.time(itt);
else
  out.time = params.tgrid(itt);
end
out.ntime = numel(out.it);

out.rho       =  model.rgrid.rho;
out.rhogauss  =  model.rgrid.rhogauss;
out.nrho      =  numel(out.rho);
out.nrhogauss =  numel(out.rhogauss);

% input signals
out.U = U;

%% Physical output quantities
% all profile quantities are on the out.rho rho-toroidal grid
% unless they are called *gauss, then they are on the out.rhogauss grid

% Psi related quantities

% poloidal flux profile
out.psi = eval_psi(X,[],[],model,false); % values of psi on rho grid

% normalized poloidal flux profile
out.psiN = bsxfun(@times,bsxfun(@minus,out.psi,out.psi(1,:)),1./(out.psi(end,:)-out.psi(1,:)));

% normalized square root poloidal flux
out.rhopol = sqrt(abs(out.psiN));

% dpsi/drho
out.psip = eval_psip(X,[],[],model,false);
% 'loop voltage profile' d/dt (psi)
out.upl = eval_psi(Xdot,[],[],model,false); % loop voltage
out.dupldrho = eval_psip(Xdot,[],[],model,false);

% steady-stateness factor. Zero = steady-state. High=transient
% (from F. Felici PPCF2012 paper)
out.f_ss = sum(out.dupldrho.^2); % f_ss (steady state factor)

% geometry factors: see doc/RAPTOR_equation_summary/RAPTOR_equilibrium_coupling.pdf
out.H  = eval_H([],G,[],model,false); % dpsi/dV
out.g0 = eval_g0([],G,[],model,false);
out.g1 = eval_g1([],G,[],model,false);
out.g2 = eval_g2([],G,[],model,false);
out.g3 = eval_g3([],G,[],model,false);
out.g23or = eval_g23or([],G,[],model,false);
out.Vp = eval_Vp([],G,[],model,false);
out.F  = eval_F([],G,[],model,false);
out.Phib  = eval_Phib([],G,[],model,false);
out.kappa = eval_kappa([],G,[],model,false);
out.epsilon = eval_epsilon([],G,[],model,false);

%%
for it=1:numel(t_ok);
  trap = traps{it};
  stap = staps{it};
  geop = geops{it};
  
  if  numel(t_ok) >= 1 && exist('simres','var')
    xk = simres.X(:,it);
    gk = simres.G(:,it);
    vk = simres.V(:,it);
  end % else, for t_ok=1, xk,gk,vk already defined
  
  out.rhob(it) = sqrt(out.Phib(it)/(pi*model.equi.B0));
  
  % Te related quantities
  % electron temperature profile
  out.te(:,it) = eval_te(xk,gk,vk,model,false); % electron temperature
  out.te0(:,it) = out.te(1,it); % central electron temperature
  out.tep(:,it) = eval_tep(xk,gk,vk,model,false); % dTe/drhotor
  
  % other kinetic profiles
  out.ne(:,it) = eval_ne(xk,gk,vk,model,false); % electron density
  out.nep(:,it) = eval_nep(xk,gk,vk,model,false); % dne/drhotor
  out.ti(:,it) = eval_ti(xk,gk,vk,model,false); % ion temperature
  out.tip(:,it) = eval_tip(xk,gk,vk,model,false); % dTi/drhotor
  
  out.ni(:,it) = eval_ni(xk,gk,vk,model,false); % ion density
  out.nip(:,it) = eval_nip(xk,gk,vk,model,false); % ion density gradient
  
  out.n1(:,it) = eval_nx(1,xk,gk,vk,model,false,false); % minority species No1 density
  out.n2(:,it) = eval_nx(2,xk,gk,vk,model,false,false); % minority species No2 density
  out.n3(:,it) = eval_nx(3,xk,gk,vk,model,false,false); % minority species No2 density
  out.nitot(:,it) = out.ni(:,it) + out.n1(:,it) + out.n2(:,it) + out.n3(:,it);
  
  out.ze(:,it) = eval_ze(xk,gk,vk,model,false,false); % effective Z
  
  out.vt(:,it)  = eval_y(xk,vk,model.vt,model.rgrid,false);
  out.vtp(:,it) = eval_yp(xk,vk,model.vt,model.rgrid,false);
  
  % pressures
  out.presse(:,it) = 1.6e-19*(out.ne(:,it).*out.te(:,it)); % electron thermal pressure
  out.pressi(:,it) = 1.6e-19*(out.ni(:,it).*out.ti(:,it)); % ion thermal pressure
  out.press1(:,it) = 1.6e-19*(out.n1(:,it).*out.ti(:,it)); % impurity 1 thermal pressure
  out.press2(:,it) = 1.6e-19*(out.n2(:,it).*out.ti(:,it)); % impurity 2 thermal pressure
  out.press3(:,it) = 1.6e-19*(out.n3(:,it).*out.ti(:,it)); % impurity 3 thermal pressure
  out.p(:,it) =  out.presse(:,it) + out.pressi(:,it) + out.press1(:,it) + out.press2(:,it) + out.press3(:,it); % total (thermal) pressure profile
  
  % magnetic profiles
  out.iota(:,it) = eval_iota(xk,gk,vk,model,false); % iota (=1/q)
  out.q(:,it)    = eval_q(xk,gk,vk,model,false);
  out.shear(:,it)= eval_shear(xk,gk,vk,model,false); % rho/q dq/drho with rho=rhotorN
  
  out.qgauss(:,it)    = stap.q;
  out.psigauss(:,it)   = stap.psi;
  out.sheargauss(:,it)= stap.shear;
  out.uplgauss(:,it) =  stap.upl;
  out.psiNgauss(:,it) = bsxfun(@times,bsxfun(@minus,out.psigauss(:,it),out.psi(1,it)),1./(out.psi(end,it)-out.psi(1,it)));
  out.psipgauss(:,it) = stap.psip;
  
  out.Ip(:,it)   = eval_Ipl(xk,gk,vk,model,false); % Plasma current (rho)
  out.jtor(:,it) = eval_jphioR(xk,gk,vk,model,false); % jtor =DEF= R0 * <jphi/R>
  %out.jtorgauss(:,it) = eval_jphioR(xk,gk,vk,model, true); % jtor =DEF= R0 * <jphi/R>
  out.jpar(:,it) = eval_jdotB(xk,gk,vk,model,false)/model.equi.B0; % j|| =DEF= <j dot B> /B0
  
  out.q0(it)     = 1./out.iota(1,it); % q at axis
  out.q95(:,it)  = 1./interp1(model.rgrid.rho,out.iota(:,it),interp1(out.psiN(:,it),model.rgrid.rho,0.95)); % q at 95% normalized flux
  out.qmin(it)   = min(out.q(:,it)); % q min
  
  out.jbsgauss(:,it)      = trap.jbsB/model.equi.B0;% <jbs dot B> /B0 (bootstrap)
  out.jecgauss(:,it)      = trap.jec;% <jec dot B> /B0 (EC)
  out.jnbgauss(:,it)      = trap.jnb;% <jnb dot B> /B0 (NB)
  out.jlhgauss(:,it)      = trap.jlh;% <jlh dot B> /B0 (LH)
  out.jicgauss(:,it)      = trap.jic;% <jic dot B> /B0 (IC)
  out.jauxgauss(:,it)     = trap.jauxB/model.equi.B0; % <jaux dot B> /B0 (aux = EC+NB+LH+IC)
  out.signeogauss(:,it)   = trap.signeo; % neoclassical conductivity sigmal_parallel
  out.jpargauss(:,it)     = eval_jdotB(xk,gk,vk,model,true)/model.equi.B0; % see above
  % ohmic current = sigma_parallel*dpsi/dt * rho/(dV/drho)* 2*Phib / B0
  out.johgauss(:,it)      = out.signeogauss(:,it).*out.uplgauss(:,it).*out.rhogauss./eval_Vp([],G(:,it),[],model,true)*2*eval_Phib([],G(:,it),[],model,false)/model.equi.B0;
  out.jregauss(:,it)           = trap.jreB/model.equi.B0;
  
  % same but linearly interpolated to rho. % BEWARE THESE MAY BE INACCURATE!!
  out.jbs(:,it)    = interp1(out.rhogauss,out.jbsgauss(:,it),out.rho,'linear','extrap');
  out.jaux(:,it)   = interp1(out.rhogauss,out.jauxgauss(:,it),out.rho,'linear','extrap');
  out.joh(:,it)    = interp1(out.rhogauss,out.johgauss(:,it),out.rho,'linear','extrap');
  out.jec(:,it)    = interp1(out.rhogauss,out.jecgauss(:,it),out.rho,'linear','extrap');
  out.jnb(:,it)    = interp1(out.rhogauss,out.jnbgauss(:,it),out.rho,'linear','extrap');
  out.jlh(:,it)    = interp1(out.rhogauss,out.jlhgauss(:,it),out.rho,'linear','extrap');
  out.jic(:,it)    = interp1(out.rhogauss,out.jicgauss(:,it),out.rho,'linear','extrap');
  out.jre(:,it)    = interp1(out.rhogauss,out.jregauss(:,it),out.rho,'linear','extrap');
  
  out.signeo(:,it) = interp1(out.rhogauss,trap.signeo,out.rho,'linear','extrap');
  
  % power densities [W/m^3]
  out.pohgauss(:,it)       = trap.poh;     % ohmic power density
  out.pecgauss(:,it)       = trap.pec;     % EC
  out.pnbegauss(:,it)      = trap.pnbe;    % NB power to electrons
  out.pnbigauss(:,it)      = trap.pnbi;    % NB power to electrons
  out.plhgauss(:,it)       = trap.plh;     % LH power
  out.picegauss(:,it)      = trap.pice;    % IC power to electrons
  out.picigauss(:,it)      = trap.pici;    % IC power to electrons
  out.pauxegauss(:,it)     = trap.pauxe;    % total auxiliary power EC+NB+LH+IC
  out.pauxigauss(:,it)     = trap.pauxi;    % total auxiliary power EC+NB+LH+IC
  out.palphaegauss(:,it)   = trap.palphae; % Alpha power to electrons
  out.palphaigauss(:,it)   = trap.palphai; % Alpha power to ions
  out.pbremgauss(:,it)     = trap.pbrem;   % Bremsstrahlung loss
  out.pradgauss(:,it)      = trap.prad;    % Radiation power loss
  out.peigauss(:,it)       = trap.pei;     % electron->ion equipartition power
  
  out.chiegauss(:,it)      = trap.chie;    % Electron thermal conductivity
  out.chiigauss(:,it)      = trap.chii;    % Electron thermal conductivity
  out.dnegauss(:,it)      = trap.dne;    % Electron particle diffusivity
  out.vnegauss(:,it)      = trap.vne;    % Electron particle pinch
  % powers on normal grid % BEWARE! may be inaccurate
  out.poh(:,it)       = interp1(out.rhogauss,trap.poh,out.rho,'linear','extrap');
  out.pec(:,it)       = interp1(out.rhogauss,trap.pec,out.rho,'linear','extrap');
  out.pnbe(:,it)      = interp1(out.rhogauss,trap.pnbe,out.rho,'linear','extrap');
  out.pnbi(:,it)      = interp1(out.rhogauss,trap.pnbi,out.rho,'linear','extrap');
  out.plh(:,it)       = interp1(out.rhogauss,trap.plh,out.rho,'linear','extrap');
  out.pice(:,it)      = interp1(out.rhogauss,trap.pice,out.rho,'linear','extrap');
  out.pici(:,it)      = interp1(out.rhogauss,trap.pici,out.rho,'linear','extrap');
  out.pauxe(:,it)     = interp1(out.rhogauss,trap.pauxe,out.rho,'linear','extrap');
  out.pauxi(:,it)     = interp1(out.rhogauss,trap.pauxi,out.rho,'linear','extrap');
  out.palphae(:,it)   = interp1(out.rhogauss,trap.palphae,out.rho,'linear','extrap');
  out.palphai(:,it)   = interp1(out.rhogauss,trap.palphai,out.rho,'linear','extrap');
  out.pbrem(:,it)     = interp1(out.rhogauss,trap.pbrem,out.rho,'linear','extrap');
  out.prad(:,it)      = interp1(out.rhogauss,trap.prad,out.rho,'linear','extrap');
  out.pei(:,it)       = interp1(out.rhogauss,trap.pei,out.rho,'linear','extrap');
  out.chie(:,it)      = interp1(out.rhogauss,trap.chie,out.rho,'linear','extrap');
  out.chii(:,it)      = interp1(out.rhogauss,trap.chii,out.rho,'linear','extrap');
  out.dne(:,it)       = interp1(out.rhogauss,trap.dne,out.rho,'linear','extrap');
  out.vne(:,it)       = interp1(out.rhogauss,trap.vne,out.rho,'linear','extrap');
  out.sne(:,it)       = interp1(out.rhogauss,trap.sne,out.rho,'linear','extrap');
  % sum powers in/out of electron species
  out.pein(:,it)      = out.poh(:,it) + out.pauxe(:,it) + ...
    out.palphae(:,it);
  out.peout(:,it)     = out.pbrem(:,it) + out.prad(:,it) + out.pei(:,it);
  out.penet(:,it)     = out.pein(:,it) - out.peout(:,it);
  % ion
  out.piin(:,it)      = out.pauxi(:,it) + out.palphai(:,it) + out.pei(:,it);
  out.piout(:,it)     = zeros(model.rgrid.nrho,1);
  out.pinet(:,it)     = out.piin(:,it) - out.piout(:,it);
  
  % integral quantities
  out.tegauss(:,it) = stap.te;
  out.negauss(:,it) = stap.ne;
  out.tigauss(:,it) = stap.ti;
  out.nigauss(:,it) = stap.ni;
  out.n1gauss(:,it) = stap.n1;
  out.n2gauss(:,it) = stap.n2;
  out.n3gauss(:,it) = stap.n3;
  
  out.We(:,it)   = int_Vtot(3/2*1.6e-19*(out.tegauss(:,it).*out.negauss(:,it)),geop,model); % electron thermal energy profile
  out.Wi(:,it)   = int_Vtot(3/2*1.6e-19*(out.tigauss(:,it).*out.nigauss(:,it)),geop,model); % ion thermal energy profile
  out.Wi1(:,it)   = int_Vtot(3/2*1.6e-19*(out.tigauss(:,it).*out.n1gauss(:,it)),geop,model); % impurity 1 thermal energy profile
  out.Wi2(:,it)   = int_Vtot(3/2*1.6e-19*(out.tigauss(:,it).*out.n2gauss(:,it)),geop,model); % impurity 2 thermal energy profile
  out.Wi3(:,it)   = int_Vtot(3/2*1.6e-19*(out.tigauss(:,it).*out.n3gauss(:,it)),geop,model); % impurity 3 thermal energy profile
  
  out.Wth(:,it) = out.We(:,it) + (out.Wi(:,it) + out.Wi1(:,it) + out.Wi2(:,it) + out.Wi3(:,it));
  out.Wpol(:,it) = int_Wpoltot(out.psipgauss(:,it).^2,geop,model); % Wi (magnetic energy)
  
  % surface integrated currents: edge value = total current
  out.Ibs(:,it)  = int_SjdotB(trap.jbsB,gk,model); % surface-integrated BS current profile
  out.Iec(:,it)  = int_SjdotB(trap.jec*model.equi.B0,gk,model); % surface-integrated EC-driven current profile
  out.Ilh(:,it)  = int_SjdotB(trap.jlh*model.equi.B0,gk,model); % surface-integrated LH-driven current profile
  out.Iic(:,it)  = int_SjdotB(trap.jic*model.equi.B0,gk,model); % surface-integrated IC-driven current profile
  out.Inb(:,it)  = int_SjdotB(trap.jnb*model.equi.B0,gk,model); % surface-integrated NB-driven current profile
  out.Iaux(:,it) = int_SjdotB(trap.jauxB,gk,model); % surface-integrated auxiliary current profile
  out.Ire(:,it)  = int_SjdotB(trap.jreB,gk,model); % surface-integrated auxiliary current profile
  
  % volume integrated power densities. edge value = total current
  out.Poh(:,it)       =  int_Vtot(out.pohgauss(:,it),geop,model); % OH power
  out.Pauxe(:,it)     =  int_Vtot(out.pauxegauss(:,it),geop,model); % AUX power to electrons
  out.Pauxi(:,it)     =  int_Vtot(out.pauxigauss(:,it),geop,model); % AUX power to ions
  out.Pauxtot(:,it)   =  out.Pauxe(:,it) + out.Pauxi(:,it);
  out.Palphatot(:,it) =  int_Vtot(out.palphaegauss(:,it)+out.palphaigauss(:,it),geop,model); % total alpha power
  out.Palphai(:,it)   =  int_Vtot(out.palphaigauss(:,it),geop,model); % alpha power to ions
  out.Palphae(:,it)   =  int_Vtot(out.palphaegauss(:,it),geop,model); % alpha power to electrons
  out.Pbrem(:,it)     =  int_Vtot(out.pbremgauss(:,it),geop,model); % bremsstrahlung
  out.Prad(:,it)      =  int_Vtot(out.pradgauss(:,it),geop,model); % radation
  out.Pcycl(:,it)     =  0; % cyclotron (placeholder)
  out.Pei(:,it)       =  int_Vtot(out.peigauss(:,it),geop,model); % ei equipartition
  out.Pec(:,it)       =  int_Vtot(out.pecgauss(:,it),geop,model); % EC power
  out.Pnbe(:,it)      =  int_Vtot(out.pnbegauss(:,it),geop,model); % NB power to electrons
  out.Pnbi(:,it)      =  int_Vtot(out.pnbigauss(:,it),geop,model); % NB power to ions
  out.Pice(:,it)      =  int_Vtot(out.picegauss(:,it),geop,model); % IC power to electrons
  out.Pici(:,it)      =  int_Vtot(out.picigauss(:,it),geop,model); % IC power to ions
  out.Volume(:,it)    =  int_V(1,gk,model); % Volume
  
  % greenwald density fraction
  out.ne20_vol_av(:,it)   =   int_Vtot(out.negauss(:,it),geop,model)./ out.Volume(end,it)/1e20; % Volume averaged plasmadenstity in 10^20 m^{-3}
  out.ne20_line_av(:,it)  =   int_Ltot(out.negauss(:,it),model)/1e20; % Volume averaged plasmadenstity in 10^20 m^{-3}
  out.neGR(:,it)      =   eval_neGR(xk,gk,vk,model); % greenwald density at this current
  out.fne_gr(1,it)    =   out.ne20_vol_av(:,it)/out.neGR(:,it);
  out.nel(:,it)       =   sum(bsxfun(@times,model.rgrid.wgauss,out.negauss(:,it)));
  
  % other profiles for equilibrium coupling
  out.pprime(:,it)    =   eval_pprime(xk,[],vk,model,false); % p' = dpress/dpsi
  out.ttprime(:,it)   =   eval_ttprime(xk,gk,vk,model,false); % tt' = t*dt/dpsi with t=R*Bphi
  %out.rbphi(:,it)     =   eval_rbphi(xk,gk,vk,model,false); % t=RBphi
  
  % some betas
  out.beta(:,it)      = eval_beta( xk,gk,vk,model); % thermal beta w.r.t. B0
  out.betaN(:,it)     = eval_betaN(xk,gk,vk,model); % normalized beta
  out.betapol(:,it)   = eval_betap(xk,gk,vk,model); % poloidal beta
  
  %% NTM width
  % NTM island widths if present in the state
  % % w_mn: width of NTM with given n and m number
  % % wdot_mn: time derivative of NTM with given n and m number
  
  for iw = 1:model.ntm.n_ntm
    % island name (e.g. w_21)
    wname = sprintf('w_%d%d',params.ntm.m(iw),params.ntm.n(iw));
    wdotname = sprintf('wdot_%d%d',params.ntm.m(iw),params.ntm.n(iw));
    
    if strcmp(model.ntm.method,'state')
      out.(wname)(it) = staps{it}.wntms(iw);
      out.(wdotname)(it) = staps{it}.wdotntms(iw);
    else
      out.(wname)(it) = V(model.ntm.vind(iw),it);
      if it==1
        out.(wdotname) = 0;
      else
        out.(wdotname) = diff(V(model.ntm.vind(iw),[it-1,it]),[],2)./diff(out.time);
      end
    end
  end
  
  %% rational rho locations
  iotagauss = stap.iota;
  rhogauss = model.rgrid.rhogauss;
  
  % out.rhoqn_m % rho_toroidal of rational surface
  for n=[1,2]
    for m = [1,2,3]
      imn = find(iotagauss>(n/m),1,'last'); % from out to in
      strr = sprintf('rhoq%d_%d',m,n);
      %             if isempty(imn)
      if isempty(imn) || imn == numel(iotagauss)
        out.(strr)(it) = NaN;
      else
        out.(strr)(it) = rhogauss(imn) + ...
          (n/m-iotagauss(imn))*diff(rhogauss(imn+(0:1)))./diff(iotagauss(imn+(0:1)));
      end
    end
  end
end

out.Li  = 2*out.Wpol(end,:)./(out.Ip(end,:).^2); % internal inductance in Henry
% normalized internal inductance per unit length
% (ITER definition, e.g. Jackson NF2008)
out.li3 = 1e7/(pi*model.equi.R0)*out.Wpol(end,:)./(out.Ip(end,:).^2);

out.pav = 2/3*out.Wth ./ out.Volume(end,:); % volume averaged pressure

out.Ini = out.Ibs + out.Iaux; % fraction of non-inductive current
out.Ioh = out.Ip - out.Ini - out.Ire; % surface-averaged Ohmic current profile

% non-inductive fractions
out.fbs = out.Ibs(end,:)./out.Ip(end,:); % BS current fraction
out.faux = out.Iaux(end,:)./out.Ip(end,:); % auxiliary (driven) current fraction
out.fni = out.Ini(end,:)./out.Ip(end,:);   % non-inductive current fraction
out.foh = out.Ioh(end,:)./out.Ip(end,:);  % ohmic current fraction
out.fre = out.Ire(end,:)./out.Ip(end,:); % fraction of runaway currents

% LH threshold power
out.PLH = eval_PLHthr(X,G,V,model);

% integrated powers

if params.echcd.active
  if ~strcmp(model.echcd.name,'hcd_external') && ~strcmp(model.echcd.name,'hcd_manual')
    iupow = model.echcd.uind(model.echcd.ind_power);
    out.Pecin = sum(U(iupow,:),1); % total NB input power
  else
    out.Pecin = []; % total input power can not be obtained for external hcd module
  end
else
  out.Pecin = zeros(1,numel(t_ok));
end

if params.nbhcd.active
  if ~strcmp(model.nbhcd.name,'hcd_external') && ~strcmp(model.nbhcd.name,'hcd_manual')
    iupow = model.nbhcd.uind(model.nbhcd.ind_power);
    out.Pnbin = sum(U(iupow,:),1); % total NB input power
  else
    out.Pnbin = []; % total input NB power can not be obtained for external hcd module
  end
else
  out.Pnbin = zeros(1,numel(t_ok));
end

if params.ichcd.active
  if ~strcmp(model.ichcd.name,'hcd_external') && ~strcmp(model.ichcd.name,'hcd_manual')
    iupow = model.ichcd.uind(model.ichcd.ind_power);
    out.Picin = sum(U(iupow,:),1); % total IC input power
  else
    out.Pecin = []; % total input power can not be obtained for external hcd module
  end
else
  out.Picin = zeros(1,numel(t_ok));
end

% sum integrated powers
out.Pein  = out.Poh + out.Pauxe + out.Palphae; % total electron input power
out.Peout = out.Pbrem + out.Prad + out.Pei;   % total electron loss power
out.Penet = out.Pein-out.Peout; % net power (in-out)

out.Piin  = out.Pauxi + out.Palphai + out.Pei; % total electron input power
out.Piout = zeros(size(out.Piin));   % total electron loss power
out.Pinet = out.Pein-out.Peout; % net power (in-out)


% Fusion gain
out.Pfusion = 5*out.Palphatot(end,:); % total fusion power (incl neutrons)
out.Palpha  = out.Palphatot(end,:); % total alpha power
out.Pheat = out.Pauxe(end,:) + out.Pauxi(end,:) + out.Poh(end,:); % total heating power
out.Q = out.Pfusion./out.Pheat; % fusion gain

% total input power
out.Ptotin = out.Poh(end,:) + out.Pnbi(end,:) + out.Pnbe(end,:)  + ...
  + out.Pec(end,:) + out.Pice(end,:) + out.Pici(end,:) + out.Palphatot(end,:);  % total input power to plasma

% Confinement time
if strcmp(outmode,'fullsim') 
  % time derivative of thermal energy
  out.dWtdt = diff(out.Wth)./diff(out.time); % simple finite difference derivative
  out.dWtdt = [out.dWtdt out.dWtdt(end)]; % append to get right size
  
  %     % CRONOS and Citrin2010 definition
  %     out.Ploss = out.Ptotin - out.Pbrem(end,:) - out.Pcycl(end,:) - 1/3*out.Prad(end,:) -  out.dWtdt;
  %     out.tauE  = (out.Wth./(out.Ptotin - out.Pbrem(end,:) -  out.dWtdt)); % Energy confinement time
  
  % Power loss
  out.Ploss = out.Ptotin - out.dWtdt;
  out.tauE = out.Wth./out.Ploss; % Energy confinement time
%   out.tauEe = out.We./out.Ploss; % Electron energy confinement time
%   out.tauEi = out.Wi./out.Ploss; % Ion energy confinement time
  
  % check if epsilon and kappa are known to calculate tauth and H_98
  if isfield(model.equi,'epsilon') && isfield(model.equi,'kappa')
    % H98
    ne19=out.ne(1,:)/1e19;
    epsilon = out.epsilon(end,:);
    kappa   = out.kappa(end,:);
    % energy confinement time for H98 scaling
    out.tauE_H98 = (5.62e-2*(out.Ip(end,:)/1e6).^0.93.*model.equi.B0^0.15.*ne19.^0.41).*(out.Ploss/1e6).^-0.69.*...
      model.equi.R0^1.97.*epsilon.^0.58.*kappa.^0.78.*model.atom.Ai.^0.19;
    out.H98 = out.tauE./out.tauE_H98;
    out.He98 = out.H98.*(1-out.Wi./out.Wth); 
    out.Hi98 = out.H98.*(1-out.We./out.Wth); 
  end
end

%% ECH deposition locations
if params.echcd.active && isfield(params.echcd,'rdep')
  n_units = numel(model.echcd.n_units);
  out.rhodep = zeros(n_units,out.ntime);
  idirect = find(params.echcd.rdep ~= -1);
  if isempty(idirect)
    % rho_dep idixes
    iinput  = model.echcd.uind(n_units + (1:n_units));
    out.rhodep = U(iinput,:);
  else
    iinput  = model.echcd.uind(n_units + (1:sum(~idirect)));
    out.rhodep(idirect,:) = params.echcd.rdep(idirect)'*ones(1,out.ntime);
    out.rhodep(~idirect,:) = U(iinput,:);
  end
else
  out.rhodep = [];
end
% EC deposition location on rhotor
%%

% Coupling to equilibrium
%        out.jdotB = out.jpar*model.equi.B0;
%        out.Ipar  = out.jdotB./out.F./out.g3./model.equi.R0; % I|| (CHEASE)

% MHD
% sawtooth detector

out.ST = [false,(diff(out.te0)./diff(out.time)<-2e5)]; % 'sawtooth detector'
out.STtimes = out.time(out.ST); % times of sawtooth crashes
out.STper = [0,diff(out.STtimes)]; % sawtooth period

if strcmp(outmode,'fullsim') && isfield(simres,'df_dx')
  % Jacobians
  % continuous time
  out.df_dx = simres.df_dx;
  out.df_dxdot = simres.df_dxdot;
  % discrete time
  out.df_dxk = simres.df_dxk;
  out.df_dxkm1 = simres.df_dxkm1;
  % output
  out.df_du = simres.df_du;
end

if numel(out.time)==1
  B0 = model.equi.B0;
  R0 = model.equi.R0;
  a = model.equi.epsilon*R0;
  Ai = model.atom.Ai;
  tau_E = out.Wth/out.Ptotin;
  out.tau_E = tau_E;
  % ELMy H mode 98y,2 scaling law
  tau_scl = (5.621e-2*(out.Ip(end)/1e6).^0.93.*B0^0.15.*(out.ne20_vol_av*10)^0.41).*(out.Ptotin/1e6).^-0.69.*...
             R0^1.97.*out.epsilon(end).^0.58.*out.kappa(end).^0.78.*Ai.^0.19; 
  out.H98y2 = tau_E/tau_scl; 
  % power crossing separatrix
  out.Psep = out.Pheat+out.Palphae(end)+out.Palphai(end)-out.Prad-out.Pbrem;

  Wthiprime = int_Vtot(3/2*1.6e-19*(stap.tip.*stap.ni+stap.ti.*stap.nip),geop,model);
  Wtheprime = int_Vtot(3/2*1.6e-19*(stap.tep.*stap.ne+stap.te.*stap.nep),geop,model); 
  Wthprime = Wtheprime + Wthiprime;
  beta_constant = 2/3 ./ (model.equi.B0.^2 / (2*4e-7*pi) .* out.Volume(end));
  betaprime = beta_constant.*Wthprime;
  out.alpha = -R0 .* betaprime .* out.q.^2;

  out.RLTi = - (R0/a) * (out.tip./out.ti);
  out.RLTe = - (R0/a) * (out.tep./out.te);
  out.RLne = - (R0/a) * (out.nep./out.ne);
end

s_sim_3o2 = zeros(1,numel(params.tgrid));
r_sim_3o2 = zeros(1,numel(params.tgrid));
d_sim_3o2 = zeros(1,numel(params.tgrid));
rflip = flip(out.rhogauss);
for tii = 1:numel(out.time)
  rii3o2 = find(flip(out.qgauss(:,tii))<=1.5,1,'first');
  riirev = find(flip(out.sheargauss(:,tii))<=0,1,'first');
  if isempty(riirev)
    d_sim_3o2(tii) = abs(out.qgauss(1,tii)-1.5);
  else
    qflip = flip(out.qgauss(:,tii));
    d_sim_3o2(tii) = abs(qflip(riirev)-1.5);
  end
  if isempty(rii3o2)
    s_sim_3o2(tii) = NaN;
    r_sim_3o2(tii) = NaN;
  else
      shearflip = flip(out.sheargauss(:,tii));
      s_sim_3o2(tii) = shearflip(rii3o2);
      r_sim_3o2(tii) = rflip(rii3o2);
  end
end
out.s_sim_3o2 = s_sim_3o2;
out.r_sim_3o2 = r_sim_3o2;
out.d_sim_3o2 = d_sim_3o2;

return

