function varargout = chi_MS(varargin)
% Heat diffusivity for electrons (and ions if Ti=='state')
% Use analytical transport model from Kim, Merle, Sauter, Goodman
% PPCF 2016(58) 055002
%
% chi_e(i) = q_e(i)/(Vp*G1*n_e(i)*T_e(i))*f(rho_inv-rho/wrho_inv)...
%         [mu_Te(i)/(T_e(i)*rho_edge)*f([rho_ped-rho]/wrho_ped)+...
%         lambda_Te(i)/rho_edge*f([rho-rho_ped]/wrho_ped)]^-1+...
%         chi_ST*f([rho-rho_inv]/wrho_inv)
%         f(x) = 1/(1+exp(x))
% + add neoclassical transport cneo to avoid chi_e(i) values close to zero.
% Returns heat conductivity and derivatives. All outputs are on grid [rhogauss]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% global muTetrace

if nargin == 0
  varargout{1} = 'chi_e'; % type;
  varargout{2} = [];
  % no config for this module
  return
elseif nargin==2
  %% setup step
    mm.name = mfilename;
    mm.type = 'chi_e';

    % check h-mode module compatibility
    RAPTORmodel = varargin{2};
    assert(strcmp(RAPTORmodel.hmode.modeltype,'none'),...
      'RAPTOR:chi_MS:hmode','H-mode module set to ''%s''! chi_MS transport module already handles H-mode pedestal so it can only be used in combination with hmode module: ''none''',RAPTORmodel.hmode.modeltype)
    
    %% DEFAULT PARAMETERS
    module_params=struct(...
        'muTe'    , 300   ,... % typical muTe (gradient in pedestal region)
        'muTi'    , 300   ,... % typical muTi (gradient in pedestal region)
        'frp'     , 2.    ,... % ratio Te/Ti at rho_ped (~0.8)
        'prescribed', false ,... % false - use controlled muTe, true - use prescribed muTe
        'lambdaTe', [3.2 2.3]   ,... % typical lambdaTe values for L-, H-modes (scalelength in core region)
        'lambdaTi', [3.0 2.5]   ,... % typical lambdaTi values for L-, H-modes (scalelength in core region)
        'rhoped'  , [0.8 0.9]   ,... % typical normalised rho pedestal, defines core/edge region (for L-, H-modes)
        'wrhoped' , 0.02  ,... % width of the transition between core and edge regions 
        'rhoinv'  , 0.15   ,... % typical normalised sawthooth inversion radius, defines center/core region
        'wrhoinv' , 0.02  ,... % width of the transition between center and core regions
        'chiST'   , 5.0   ,... % flattering effect of sawtooth
        'cneo'    , 0.15 ,... % neoclassical transport
        'hh'      , [0.35 0.8] ,... % reference H-factor: tauE_e/tau_scaling or tauE/tau_scaling on RAPTOR time grid or L/H-mode
        'gp'      , 1e+02 ,... % default value for proportional gain
        'gi'      , 1e+03 ,... % default value for integration gain
        'fb_on'  , true ,...    % feedback: true - on, false - off
        'ff_on'  , true ,...    % feedforward: true - on, false - off
        'muTeff_pred', -1, ... % predictive feed-forward; -1 - calculated
        'rhodep'  , [], ... % optional parameter to take into account He reduction in case of off-axis ECRH
        'sawtooth', false, ... % enable internal chiMS sawtooth model
        's1crit'  , 1,... % critical shear for sawtooth model
        'implicit', true  ,... % implicit method
        'check'   , false ... % option to check gradients etc
        );
      
    varargout{1} = mm;
    varargout{2} = module_params;
    return %empty call, probably to get default structures
elseif nargin==6 % change this depending on number of inputs to module
    stap = varargin{1};
    geop = varargin{2};
    trap = varargin{3};
    it = varargin{4};    
    model = varargin{5};
    cp = varargin{6};
else
    error('must call with 0 or 6 inputs');
end

persistent muTe muTi muTe_km1 HH_km1 HH_km2 HHref_km1 HHref_km2 u_ff_km1;


firstNewtonIteration = isFirstTime(it);


%% chi_e module parameters 
% Electron temperature and its derivatives over the plasma state and time
te = stap.te;
dte_dx = stap.dte_dx;
dte_dt = stap.tedot;
% Ion temperature and its derivative over time
ti = stap.ti;
dti_dx = stap.dti_dx;
dti_dt = stap.tidot;
% Electron density and its derivative over time
ne = stap.ne;
dne_dx = stap.dne_dx;
dne_dt = stap.nedot;
% Ion density and its derivative over time
ni = stap.ni;
dni_dx = stap.dni_dx;
dni_dt = stap.nidot;
% Geometry: dV/drho_norm_gauss
Vpgauss = geop.Vp;
% Geometry: g1
g1 = geop.g1;
% Aspect ratio
epsilon = geop.epsilon;
% Elongation
kappa = geop.kappa;
% Safety factor
q = stap.q;
% chie params
% Gauss radial grid, normalised
rhogauss = model.rgrid.rhogauss;
% Plasma current
Ip = stap.Ip;
% Electron charge
qelectron = 1.6e-19;
% Flag: solve for Ti => 1 else 0
flg = strcmp(model.ti.method,'state');

%% Power
% Power density for electrons
pe = trap.petot;
dpe_dx = trap.dpetot_dx;
dpe_dxdot = trap.dpetot_dxdot;
dpe_du = trap.dpetot_du;

% Input power density for electrons
pedenstot = trap.poh + trap.pec + trap.pice + trap.pnbe + trap.palphae;
% Input power density for ions
pidenstot = trap.pici + trap.pnbi + trap.palphai;
% shouldn't we substract Prad and Pbrem?
% Total input power
Ptotin = int_Vtot(pedenstot+pidenstot,geop,model);

%% chi_e module
if numel(cp.lambdaTe)==2
    lambdaTe = (1.-stap.lamH)*cp.lambdaTe(1) + stap.lamH*cp.lambdaTe(2);
else
    lambdaTe = cp.lambdaTe(it);
end
% Get current rho_ped
rhoped = (1.-stap.lamH)*cp.rhoped(1) + stap.lamH*cp.rhoped(2);
% Heat flux q_e==int Pe*Vp drho
qe = int_CompSimpsonRule(pe.*Vpgauss,model);
%qe = cumtrapz(rhogauss,pe.*Vpgauss);
% Check to avoid negative power heat flux
if any(qe<0.0)
    qe(qe<0.0) = 1.;
    warning('Negative power heat flux for electrons at it=%2.0d is replaced by 1.',it);
end

%% To solve for Ti
if flg
    % Power density for ions
    pi = trap.pitot;
    dpi_dx = trap.dpitot_dx;
    dpi_dxdot = trap.dpitot_dxdot;
    dpi_du = trap.dpitot_du;
    
    % Heat flux q_i==int Pi*Vp drho
    qi = int_CompSimpsonRule(pi.*Vpgauss,model);
    %qe = cumtrapz(rhogauss,pe.*Vpgauss);
    % Check to avoid negative power heat flux
    if any(qi<0.0)
        qi(qi<0.0) = 0.0;
        warning('Negative power heat flux for ions at it=%2.0d is replaced by 1',it);
    end
    
    % Get current lambdaTi value
    if numel(cp.lambdaTi)==2
        lambdaTi = (1.-stap.lamH)*cp.lambdaTi(1) + stap.lamH*cp.lambdaTi(2);
    else
        lambdaTi = cp.lambdaTi(it);
    end
end

% Get rho_inv
% Define rho_tor for q=1.0 surface
if isempty(cp.rhoinv)|| all(q>1)
  rhoinv = 0.2;
else
  rhoinv = find_rhoinv(q,rhogauss,stap.shear,cp);
end

%% The controller
% Find muTe for current time iteration, keep it the same during Newton
% iterations
if firstNewtonIteration   
% only for first Newton iteration
    % Thermal energy for 1) electrons+ions if Te and Ti equations are
    % solved 2) electrons if Te equation is solved
    if flg
        Wth = int_Vtot(3/2*qelectron*(te.*ne+ti.*ni),geop,model);
    else
        Wth = int_Vtot(3/2*qelectron*te.*ne,geop,model);
    end
    % Time derivative dWtot/dt
    geop_dt = geop;
    geop_dt.Vp = geop_dt.dVpdt;
    dWtotdt = int_Vtot(3/2*qelectron*(dte_dt.*ne + dti_dt.*ni),geop,model)+...
              int_Vtot(3/2*qelectron*(te.*dne_dt + ti.*dni_dt),geop,model)+...
              int_Vtot(3/2*qelectron*(te.*ne+ti.*ni),geop_dt,model);
    % Loss power
    Ploss = Ptotin - dWtotdt;
    % Avoid negative power
    if Ploss<0.0
        Ploss = 0.1*abs(Ploss);
        warning('Negative power loss at it=%2.0d is replaced by 0.1*abs(Ploss)',it);
    end
    % tau_E: 1) energy confinement time 2) electron energy confinement time 
    tau_E = Wth./Ploss;
    % tau_scaling: 
    tau_scl = (5.62e-2*(Ip(end)/1e6).^0.93.*model.equi.B0^0.15.*(ne(1)/1e19)^0.41).*(Ploss/1e6).^-0.69.*...
        model.equi.R0^1.97.*epsilon(end).^0.58.*kappa(end).^0.78.*model.atom.Ai.^0.19;
    % HH factor: 1) He+Hi 2) He 
    HH = tau_E/tau_scl;
    
    if cp.prescribed
        % use prescribed muTe
        muTe = cp.muTe(it);
    else
        % mu(ii) is calculated from mu(ii-1) and He(ii-1)
        if it==1
            if cp.muTeff_pred==-1
                muTe = cp.muTe(1);                
            else
                muTe = cp.muTeff_pred(1);
            end
            muTe_km1 = muTe;
            u_ff_km1 = muTe;
            HH_km1 = HH; 
            HH_km2 = 0.;
            HHref_km1 = cp.hh(1); 
            HHref_km2 = cp.hh(1);   

        else
            % Feedforward control based on muTe dependence on Ip and input power.
            % + Feedback PI controller
            % muTe(it) = muTe(it-1) + gi*(He_ref(it-1) - He_eff(it-1)) +
            % gp*(He_ref(it-1) - He_eff(it-1) - (He_ref(it-2) - He_eff(it-2)))
            % + muTe_ref(it) - muTe_ref(it-1)
            
            % Feedback controller                   
            if cp.fb_on
                err_km1 = HHref_km1 - HH_km1;
                err_km2 = HHref_km2 - HH_km2; 
                if (it<50)
                    gainfactor = 0.5;
                else
                    gainfactor = 1;
                end
                du_fb =  gainfactor * (cp.gi*err_km1 + cp.gp*(err_km1 - err_km2));
            else
                du_fb = 0.;
            end
            
            % Feedforward controller
            if cp.ff_on
                if cp.muTeff_pred==-1
                    % calculated feedforward
                    if strcmp(model.equi.tokamak,'TCV')
                        coeff = 5.*1.5e+03;
                    elseif strcmp(model.equi.tokamak,'AUG')
                        coeff = 3.3*1.e+03;
                    elseif strcmp(model.equi.tokamak,'JET')
                        coeff = 4.1*1.e+04;
                    end
                    
                    teb = stap.te(end); % temporary fix? [TODO] - to be looked at
                    Te0 = coeff*((1.e-6*Ip(end)).^0.93).*((1.e-6*Ptotin).^0.3).*((1.e-19*ne(1)./1.5).^(-0.6));
                    u_ff_k = (Te0*exp(-lambdaTe*(rhoped -rhoinv))- teb)./(1.0-rhoped);
                else
                    % prescribed feedforward
                    u_ff_k = cp.muTeff_pred(it);
                end
                du_ff = u_ff_k - u_ff_km1; 
                % save for the following iteration
                u_ff_km1 = u_ff_k;
            else
                du_ff = 0.;
            end
            
            % Calculate muTe
            muTe = muTe_km1 + du_ff + du_fb;        
            % Avoid too low muTe values
            if muTe<=5 % 50
                muTe = 5.;
            end
            
            % Save for the following iteration
            muTe_km1 = muTe;
            HH_km2 = HH_km1; 
            HH_km1 = HH;  
            HHref_km2 = HHref_km1;
            if numel(cp.hh)==2
                HHref_km1 = (1.-stap.lamH)*cp.hh(1) + stap.lamH*cp.hh(2);
                if ~isempty(cp.rhodep)
                    rhodep = cp.rhodep(it);
                    % Reduction of He: D Kim thesis 2015 EPFL, sec. 3.3.1
                    se = 1. - 0.25*(1. - exp(-4.*((rhodep^3 - rhoinv^3)./0.03).^2)).*myheaviside(rhodep - rhoinv);
                    HHref_km1 = se*HHref_km1;
                end
            else
                HHref_km1 = cp.hh(it);
            end            
        end
    end
%     if it == 1
%       muTetrace = zeros(1,numel(cp.hh));
%       muTetrace(it) = muTe;
% 	  elseif it == numel(cp.hh)
%       muTetrace(it) = muTe;
%       save('/home/contre/NoTivoli/aug_40405/muTetrace_cste_nel_2','muTetrace')
%     else 
%       muTetrace(it) = muTe;
%     end
end
% if firstNewtonIteration
%   plot(it, HH,'b*');hold on
% end

% f functions
f1 = 1./(1.+exp((rhoped - rhogauss)./cp.wrhoped));
f2 = 1-f1;
f3 = 1./(1.+exp((rhoinv - rhogauss)./cp.wrhoinv));
f4 = 1-f3;

% Expression in the brackets
expr_br = muTe*f1./te + lambdaTe*f2;

% Formula chi_e
chie = (Vpgauss.*qe./(qelectron*te.*g1.*ne.*expr_br)).*f3 + cp.chiST.*f4;
% Avoid negative and close-to-zero chi_e: take a max value between local
% chi_e and fixed chi_e neoclassical.
chie = max(chie,cp.cneo);

if flg
    muTi = muTe./cp.frp + (model.te.BCval - cp.frp*model.ti.BCval)/(cp.frp*(1. - rhoped));
    % Expression in the brackets
    expr_bri = muTi*f1./ti + lambdaTi*f2;

    % Formula chi_i
    chii = (Vpgauss.*qi./(qelectron*ti.*g1.*ni.*expr_bri)).*f3 + cp.chiST.*f4;
    % Avoid negative and close-to-zero chi_i: take a max value between local
    % chi_i and fixed chi_i neocalssical.
    chii = max(chii,cp.cneo);
else
    muTi = 0.;
    chii = 0.*chie;
end

varargout{1} = chie;
varargout{2} = chii;

if (nargout>2) % then compute detivatives      
    %%% Heat flux qe derivatives
    % dqe_dx
    %dqe_dx = cumtrapz(rhogauss,bsxfun(@times,dpe_dx,Vpgauss));
    dqe_dx = cumsum(bsxfun(@times,model.rgrid.wgauss.*Vpgauss,dpe_dx));
    % dqe_dxdot
    %dqe_dxdot = cumtrapz(rhogauss,bsxfun(@times,dpe_dxdot,Vpgauss));
    dqe_dxdot = cumsum(bsxfun(@times,model.rgrid.wgauss.*Vpgauss,dpe_dxdot));
    % dqe_du
    %dqe_du = cumtrapz(rhogauss,bsxfun(@times,dpe_du,Vpgauss));
    dqe_du = cumsum(bsxfun(@times,model.rgrid.wgauss.*Vpgauss,dpe_du));
    
    %%% Derivatives of the expression in the brackets
    % dbrexpr_dx
    dexprbr_dx = -bsxfun(@rdivide,muTe*bsxfun(@times,f1,dte_dx),te.^2);   
    
    %%% Heat diffusivity chi_e derivatives    
    % dchie_dx
    qegrad_dx = bsxfun(@times,dqe_dx,Vpgauss.*f3./(qelectron*te.*g1.*ne.*expr_br));
    grad_dx = - bsxfun(@rdivide,dte_dx,te)...
              - bsxfun(@rdivide,dne_dx,ne)...
              - bsxfun(@rdivide,dexprbr_dx,expr_br);
    dchie_dx = qegrad_dx + bsxfun(@times,(chie - cp.chiST.*f4),grad_dx);    
    % dchie_dxdot
    dchie_dxdot = bsxfun(@times,dqe_dxdot,Vpgauss.*f3./(qelectron*te.*g1.*ne.*expr_br));    
    % dchie_du
    dchie_du = bsxfun(@times,dqe_du,Vpgauss.*f3./(qelectron*te.*g1.*ne.*expr_br));
    
    varargout{3} = dchie_dx;
    varargout{5} = dchie_dxdot;    
    varargout{7} = dchie_du;    
    
    if flg
        %%% Heat flux qi derivatives
        % dqi_dx
        dqi_dx = cumsum(bsxfun(@times,model.rgrid.wgauss.*Vpgauss,dpi_dx));
        % dqi_dxdot
        dqi_dxdot = cumsum(bsxfun(@times,model.rgrid.wgauss.*Vpgauss,dpi_dxdot));
        % dqi_du
        dqi_du = cumsum(bsxfun(@times,model.rgrid.wgauss.*Vpgauss,dpi_du));

        %%% Derivatives of the expression in the brackets
        dexprbr_dxi = -bsxfun(@rdivide,muTi*bsxfun(@times,f1,dti_dx),ti.^2);   

        %%% Heat diffusivity chi_e derivatives    
        % dchii_dx
        qigrad_dx = bsxfun(@times,dqi_dx,Vpgauss.*f3./(qelectron*ti.*g1.*ni.*expr_bri));
        grad_dxi = - bsxfun(@rdivide,dti_dx,ti)...
                   - bsxfun(@rdivide,dni_dx,ni)...
                   - bsxfun(@rdivide,dexprbr_dxi,expr_bri);
        dchii_dx = qigrad_dx + bsxfun(@times,(chii - cp.chiST.*f4),grad_dxi);    
        % dchii_dxdot
        dchii_dxdot = bsxfun(@times,dqi_dxdot,Vpgauss.*f3./(qelectron*ti.*g1.*ni.*expr_bri));    
        % dchii_du
        dchii_du = bsxfun(@times,dqi_du,Vpgauss.*f3./(qelectron*ti.*g1.*ni.*expr_bri));
        
        varargout{4} = dchii_dx;
        varargout{6} = dchii_dxdot;
        varargout{8} = dchii_du;
    else
        varargout{4} = 0.*dchie_dx;
        varargout{6} = 0.*dchie_dxdot;
        varargout{8} = 0.*dchie_du;
    end    
    
    % optional: test gradients
    if cp.check == 1
        check_gradients(x,xdot,g,v,vdot,u,it,model,cp,chie,dchie_dx,dchie_dxdot,dchie_du)
    end
else
    varargout{3} = zeros(numel(chie),model.dims.nx);
    varargout{4} = zeros(numel(chie),model.dims.nx);
    varargout{5} = zeros(numel(chie),model.dims.nx);
    varargout{6} = zeros(numel(chie),model.dims.nx);
    varargout{7} = zeros(numel(chie),model.dims.nu);
    varargout{8} = zeros(numel(chie),model.dims.nu);
end

end

function varout = myheaviside(varin)
varout=1.;
varout(varin<0.)=0.;
end

function check_gradients(x,xdot,g,v,vdot,u,it,model,cp,chie,dchie_dx,dchie_dxdot,dchie_du)
% Perturbation level
alpha = [1e-04 1e-03 1e-02];
alpha = [alpha -alpha];
na = numel(alpha);
% Initialize
ratio_dchiedte = zeros(model.rgrid.nrhogauss,na);
ratio_dchiedpsi = zeros(model.rgrid.nrhogauss,na);
ratio_dchiedtedot = zeros(model.rgrid.nrhogauss,na);
ratio_dchiedpsidot = zeros(model.rgrid.nrhogauss,na);
ratio_dchiedIp = zeros(model.rgrid.nrhogauss,na);
ratio_dchiedu = zeros(model.rgrid.nrhogauss,na);

for ii=1:na
    % dpsi_drho perturbation (in jtor dpsi_drho)
    dpsip = alpha(ii)*ones(size(model.rgrid.rhogauss));    
    % Spline coefficients for perturbed psi
    psihat_pert = model.psi.Lampgauss\dpsip;
    % Perturbed state vector
    dx = [psihat_pert; zeros(size(psihat_pert))];
    % Get perturbed part of chie
    dchie_num = chi_MS(x + dx,xdot,g,v,vdot,u,it,model,cp) - chie;
    % Get analytical dchie from dchie_dx
    dchie_anl = dchie_dx*dx;
    % Comparison of anlytical and numerical derivatives
    dchiedp_num = dchie_num./dpsip;
    dchiedp_anl = dchie_anl./dpsip;
    % Ratio
    ratio_dchiedpsi(:,ii) = dchiedp_anl./(1e-06*dchiedp_num(end) + dchiedp_num);
    
    % Electron temperature perturbation
    dte = alpha(ii)*ones(size(model.rgrid.rhogauss));    
    % Spline coefficients for perturbed te
    tehat_pert = model.te.Lamgauss\dte;
    % Perturbed state vector
    dx = [zeros(size(tehat_pert));tehat_pert];
    % Get perturbed part of chie
    dchie_num = chi_MS(x + dx,xdot,g,v,vdot,u,it,model,cp) - chie;
    % Get analytical dchie from dchie_dx
    dchie_anl = dchie_dx*dx;
    % Comparison of anlytical and numerical derivatives
    dchiedp_num = dchie_num./dte;
    dchiedp_anl = dchie_anl./dte;  
    % Ratio
    ratio_dchiedte(:,ii) = dchiedp_anl./(1e-06*dchiedp_num(end) + dchiedp_num);

    % Psi dot perturbation
    dpsidot = alpha(ii)*ones(size(model.rgrid.rhogauss));  
    % Spline coefficients for perturbed te
    psihatdot_pert = model.psi.Lamgauss\dpsidot;
    % Perturbed state vector
    dxdot = [psihatdot_pert; zeros(size(psihatdot_pert))];
    % Get perturbed part of chie
    dchie_num = chi_MS(x,xdot+dxdot,g,v,vdot,u,it,model,cp) - chie;
    % Get analytical dchie from dchie_dxdot
    dchie_anl = dchie_dxdot*dxdot;
    % Comparison of anlytical and numerical derivatives
    dchiedp_num = dchie_num./dpsidot;
    dchiedp_anl = dchie_anl./dpsidot;
    % Ratio
    ratio_dchiedpsidot(:,ii) = dchiedp_anl./(1e-06*dchiedp_num(end) + dchiedp_num);
    
    % Electron temperature time derivative perturbation
    dtedot = alpha(ii)*ones(size(model.rgrid.rhogauss));    
    % Spline coefficients for perturbed te
    tehatdot_pert = model.te.Lamgauss\dtedot;
    % Perturbed state vector
    dxdot = [zeros(size(tehatdot_pert));tehatdot_pert];
    % Get perturbed part of chie
    dchie_num = chi_MS(x,xdot+dxdot,g,v,vdot,u,it,model,cp) - chie;
    % Get analytical dchie from dchie_dxdot
    dchie_anl = dchie_dxdot*dxdot;
    % Comparison of anlytical and numerical derivatives
    dchiedp_num = dchie_num./dtedot;
    dchiedp_anl = dchie_anl./dtedot;
    % Ratio
    ratio_dchiedtedot(:,ii) = dchiedp_anl./(1e-06 + dchiedp_num);
    
    % Actuators perturbation
    du = alpha(ii)*ones(size(u));  
    % Get perturbed part of chie
    dchie_num = chi_MS(x,xdot,g,v,vdot,u+du,it,model,cp) - chie;
    % Get analytical dchie from dchie_dxdot
    dchie_anl = dchie_du*du;
    % Comparison of anlytical and numerical derivatives
    % Plasma current
    dchiedp_num = dchie_num./du(1);  
    dchiedp_anl = dchie_anl./du(1);
    % Ratio
    ratio_dchiedIp(:,ii) = dchiedp_anl./(1e-06*dchiedp_num(end) + dchiedp_num);
    % Anothe actuator
    if du>1
    dchiedp_num = dchie_num./du(2);  
    dchiedp_anl = dchie_anl./du(2);
    ratio_dchiedu(:,ii) = dchiedp_anl./(1e-06*dchiedp_num(end) + dchiedp_num);
    end
end

figure;
subplot(2,3,1)
 plot(model.rgrid.rhogauss,ratio_dchiedpsi,'o');
 legend('alpha=1e-04','alpha=1e-03','alpha=1e-02','alpha=-1e-04','alpha=-1e-03','alpha=-1e-02');
 str = sprintf('dchie/dpsi anl/num: it=%2.0f',it);title(str);
subplot(2,3,2)
 plot(model.rgrid.rhogauss,ratio_dchiedte,'o');
 str = sprintf('dchie/dte anl/num: it=%2.0f',it);title(str);
subplot(2,3,3)
 plot(model.rgrid.rhogauss,ratio_dchiedpsidot,'o');
 str = sprintf('dchie/dpsidot anl/num: it=%2.0f',it);title(str);
subplot(2,3,4)
 plot(model.rgrid.rhogauss,ratio_dchiedtedot,'o');
 str = sprintf('dchie/dtedot anl/num: it=%2.0f',it);title(str);
subplot(2,3,5)
 plot(model.rgrid.rhogauss,ratio_dchiedIp,'o');
 str = sprintf('dchie/dIp anl/num: it=%2.0f',it);title(str);
subplot(2,3,6)
 plot(model.rgrid.rhogauss,ratio_dchiedu,'o');
 str = sprintf('dchie/du anl/num: it=%2.0f',it);title(str);
end

function firstTime = isFirstTime(it)
persistent itp
if isempty(itp)
  itp = 0; % init
end

firstTime = (itp ~= it);
itp = it; % update counter
end
