function  [dwdt] = mre_terms(stap,trap,geop,u,model,params,it)
% The function evaluates the mode growth rate based on the modified rutherford equation (MRE).
% Note: heatgin term added Dick Hogeweij 01/06.2016
% rDp = rsDeltap_cl + rsDeltap_bs + rsDeltap_ggj + rsDeltap_cd + rsDeltap_H
% dw/dt = rDp/pre_dwdt
% A detailed description of the
% different terms can be found in the document labelled
% "Details on the expressions used in the various terms of the modified
% Rutherford equation (MRE)" by Sauter et. al, which should be shipped with
% RAPTOR. The specific equations that are implemented in this function are
% also listed in the section Neoclassical Tearing Modes in the
% RAPTOR_equation_summary pdf-file.
%
% In this function rho is the normalized toroidal flux (this is
% $\hat{\rho}$ in the RAPTOR_equation_summary.pdf).
% The equations of the modified rutherford equation rely on rho*a, with a
% the minor radius. In this function, this value is labelled rho_bar=rho*a.
% The labels mn refer to values that are evaluated on the flux surface
% labelled by m and n.
%
% Note that the expressions in this function can be further extended.
% One obvious extension would be to make the code implicit, by calculating
% derivatives.
% Another possible extension is: (feel free to add items)
% * The marginal width for electrons w_d_e can be calculated directly from
%   equation (11) in "Details on the expressions used ..." instead of being
%   specified as input.

% version taken from mengdi 01.06.2017, added log_a_bs/ggj

%% Use global variables if debug is set
global log_rhodeltap_cl
global log_rhodeltap_ggj
global log_rhodeltap_bs
global log_rhodeltap_cd
global log_rhodeltap_h
global log_eta_H
global log_pre_dwdt
global log_dwdt
global log_j_cd_ntm
global log_j_bs_mn_up
global log_L_bs_mn
global log_rho_mn
global log_rho_mn_test
global log_Ncd
global log_a_bs log_a_ggj
if params.ntm.debug
  if it==1	% Initialize global params
    n_time=numel(params.tgrid);
    n_ntm = model.ntm.n_ntm;
    log_dwdt = zeros(n_ntm,n_time);
    log_pre_dwdt = zeros(n_ntm,n_time);
    log_rhodeltap_cl = zeros(n_ntm,n_time);
    log_rhodeltap_ggj = zeros(n_ntm,n_time);
    log_rhodeltap_bs = zeros(n_ntm,n_time);
    log_rhodeltap_cd = zeros(n_ntm,n_time);
    log_rhodeltap_h = zeros(n_ntm,n_time);
    log_eta_H = zeros(n_ntm,n_time);
    log_j_bs_mn_up = zeros(n_ntm,n_time);
    log_L_bs_mn = zeros(n_ntm,n_time);
    log_j_cd_ntm = zeros(n_ntm,n_time);
    log_rho_mn = zeros(n_ntm,n_time);
    log_rho_mn_test = zeros(n_ntm,n_time);
    log_Ncd = zeros(n_ntm,n_time);
    log_a_bs = zeros(n_ntm,n_time);
    log_a_ggj = zeros(n_ntm,n_time);
  end
end

h_fact = params.ntm.h_fact;
g_coeff = params.ntm.g_coeff;

%% Physical constants mu0, e, and minor radius a
mu0 = 4*pi*10^-7;
ee = 1.602e-19;
a_minor=model.equi.epsilon*model.equi.R0;

%% Define the minor radius, and get beta_pol, iota, q, s, Te and Grad Te on rhogauss
rhogauss  = model.rgrid.rhogauss;
iotagauss = stap.iota;
qprof     = stap.q;
shear     = stap.shear;
Te        = ee * stap.te;	% in J
den       = stap.ne;
chie      = trap.chie;
epar      = stap.upl / (2*pi*model.equi.R0);
signeo    = trap.signeo;
joh       = epar .* signeo;	% neglect neo-classical effects

%% Get the bootstrap current and neocl. conductivity which needs to be known at the islands
j_bs   = trap.jbsB/model.equi.B0;
L_bs   = trap.Lbs;
f_bse  = trap.fbse;
signeo = trap.signeo;

%% Get the deposition positions in rho_bar_dep
if strcmp(model.echcd.name,'echcd_gaussian')
  nbeams = model.echcd.n_units;
  iupow = model.echcd.uind(model.echcd.ind_power); % indices (in u) of the powers
  iurdep = model.echcd.uind(model.echcd.ind_rhodep);
  irdep = model.echcd.rdep_ext; % indices of rdep taken externally
  rho_dep(~irdep) = params.echcd.rdep(~irdep); % some deposition locations are directly passed as fixed parameters
  if any(irdep)
    rho_dep(irdep) = u(iurdep); % other deposition locations are time varying, taken from u
  end
else
  %TODO: proper error handling if not defined
end

rho_bar_dep = rho_dep*a_minor; %convert rho to r
n_ntm = model.ntm.n_ntm;
% init
dwdt = zeros(n_ntm,1);

rho_mn = zeros(n_ntm,1);
rho_bar_mn = zeros(n_ntm,1);

%% test rho_mn
rho_mn_test = zeros(n_ntm,1);
rho_bar_mn_test = zeros(n_ntm,1);
for intm=1:n_ntm
  m=params.ntm.m(intm);
  n=params.ntm.n(intm);
  for imn = model.rgrid.nrhogauss:-1:1; % search from out to in
    if iotagauss(imn)>n/m; break; end
  end
  % interpolate to find rho_mn
  rho_mn_test(intm) = rhogauss(imn) + ...
    (n/m-iotagauss(imn))*diff(rhogauss(imn+(0:1)))./diff(iotagauss(imn+(0:1)));
  % Use rho_bar_mn in meter by multiplying normalized rho with a_minor
  rho_bar_mn_test(intm)=rho_mn_test(intm)*a_minor;
end
log_rho_mn_test(:,it)          = rho_mn_test;
%% Loop over all specified modes
for intm=1:n_ntm
  m=params.ntm.m(intm);
  n=params.ntm.n(intm);
  % Get the island width and evaluate the MRE if w>0 and in the plasma
  w=stap.wntms(intm);
  % make sure variables are defined
  pre_dwdt = 0;
  rhodeltap_cl = 0;
  rhodeltap_ggj = 0;
  rhodeltap_bs = 0;
  rhodeltap_cd = 0;
  rhodeltap_h = 0;
  eta_H = 0;
  j_bs_mn_up = 0;
  L_bs_mn = 0;
  j_cd_ntm=0;
  a_bs = 0;
  a_ggj = 0;
  if (w > 0) && (any(iotagauss>=n/m))
    %% Find the position of the mode and associated rho/rho_bar
    % Search for q=m/n or iota=n/m
    for imn = model.rgrid.nrhogauss:-1:1; % search from out to in
      if iotagauss(imn)>n/m; break; end
    end
    % interpolate to find rho_mn
    rho_mn(intm) = rhogauss(imn) + ...
      (n/m-iotagauss(imn))*diff(rhogauss(imn+(0:1)))./diff(iotagauss(imn+(0:1)));
    % Use rho_bar_mn in meter by multiplying normalized rho with a_minor
    rho_bar_mn(intm)=rho_mn(intm)*a_minor;
    %% Evaluate unperturbed j_bs
    % Make sure the area in which we look at j_bs_ntm is large enough
    % but not too large, for chiefact_w=2 this values is 1.2 and there
    % is a small dependence of the value on chiefact_w=1 to 3
    fact_j_bs_up=1.2;
    % Calculate normalized width to compare with rho value
    w_lim=0.5*w/a_minor*params.ntm.chiefact_w(intm)*fact_j_bs_up;
    % find rho=rho-w and rho+w, to obtain unperturbed jbs
    for imnleft = imn:-1:1
      if rhogauss(imnleft)<rho_mn(intm)-w_lim; break; end
    end
    for imnright = imn:length(rhogauss)
      if rhogauss(imnright)>rho_mn(intm)+w_lim; break; end
    end
    % Subtract 1 such that rho_mn(i)+w_norm is in the interval
    % rhogauss(imnright+(0:1))
    imnright=imnright-1;
    % Unperturbed bootstrap current density at rational surface (interpolate from values outside island)
    % The interpolation is based on the values j_bs_l and j_bs_r which
    % are also interpolated.
    nm_rhog_l=(rho_mn(intm)-w_lim-rhogauss(imnleft))/diff(rhogauss(imnleft+(0:1)));
    nm_rhog_r=(rho_mn(intm)+w_lim-rhogauss(imnright))/diff(rhogauss(imnright+(0:1)));
    j_bs_l=j_bs(imnleft) + diff(j_bs(imnleft+(0:1))) * nm_rhog_l;
    j_bs_r=j_bs(imnright) + diff(j_bs(imnright+(0:1))) * nm_rhog_r;
    rho_l=rho_mn(intm)-w_lim;
    rho_r=rho_mn(intm)+w_lim;
    j_bs_mn_up=((rho_r-rho_mn(intm))*j_bs_l+(rho_mn(intm)-rho_l)*j_bs_r)./(rho_r-rho_l);
    %% Get values at rho_bar_mn
    % Get q, s, signeo, L_bs, and f_bse on the NTM surface
    % also get Te, den, chie and joh n the NTM surface - added dick h 02.06.2017
    mn_iota=(n/m-iotagauss(imn))/diff(iotagauss(imn+(0:1)));
    q_mn       = qprof(imn)   + diff(qprof(imn+(0:1)))    * mn_iota;
    s_mn       = shear(imn)   + diff(shear(imn+(0:1)))    * mn_iota;
    signeo_mn  = signeo(imn)  + diff(signeo(imn+(0:1)))   * mn_iota;
    L_bs_mn    = L_bs(imn)    + diff(L_bs(imn+(0:1)))     * mn_iota;
    f_bse_mn   = f_bse(imn)   + diff(f_bse(imn+(0:1)))    * mn_iota;
    Te_mn      = Te(imn)      + diff(Te(imn+(0:1)))       * mn_iota;
    den_mn     = den(imn)     + diff(den(imn+(0:1)))      * mn_iota;
    chie_mn    = chie(imn)    + diff(chie(imn+(0:1)))     * mn_iota;
    joh_mn     = joh(imn)     + diff(joh(imn+(0:1)))      * mn_iota;
    % Aspect ratio at rational surface
    eps_mn=rho_bar_mn(intm)/model.equi.R0;
    %% Calculate the prefactor before dw/dt
    % Factor before dw/dt which is tau_r/rho_bar_mn (Sauter PPCF 2010)
    % The local resistivity can be increased or decreased using resist_fact.
    pre_dwdt(intm)=mu0*rho_bar_mn(intm)*signeo_mn/(1.22*params.ntm.resis_fact(intm));
    %% Calculate classical growth rate
    % The classical growth rate could be both negative and positive and
    % this value can be set to params.ntm.rhodeltap_cl0(i). It is known
    % that for large values the classical growth rate become -m.
    % Both effects are taken into account below.
    rhodeltap0m=(params.ntm.rhodeltap_cl0(intm)+m);
    rhodeltap_cl=params.ntm.rhodeltap_cl0(intm)-params.ntm.alpha(intm)*(rhodeltap0m * ...
      w/sqrt(params.ntm.alpha(intm)^2*w^2+rhodeltap0m^2*rho_bar_mn(intm)^2+10*eps));
    %% Calculate bootstrap term
    % Prefactor for bootstrap current term
    a_bs=2*params.ntm.a2_bs(intm)*rho_bar_mn(intm)^2*mu0*q_mn/(s_mn*eps_mn*model.equi.B0)*j_bs_mn_up;
    rhodeltap_bs=a_bs*w*(f_bse_mn./(w.^2+params.ntm.w_d_e(intm).^2)+(1-f_bse_mn)./(w.^2+params.ntm.w_d_i(intm).^2));
    %% Calculate Glasser Greene Johnson term
    a_ggj=params.ntm.a_ggj(intm)*12*mu0*model.equi.R0*rho_bar_mn(intm)*eps_mn^2/(q_mn*model.equi.B0*s_mn^2)*(q_mn^2-1)*j_bs_mn_up/abs(L_bs_mn);
    rhodeltap_ggj=-a_ggj./sqrt(w.^2+0.2*params.ntm.w_d_e(intm).^2);
    %% Calculate Current Drive term
    % Loop over all ech sources to check whether they influence the mode
    rhodeltap_cd=0;
    rhodeltap_h=0;
    wstar = w/params.ntm.w_cd_ntm(intm);
    for ibeam=1:nbeams
      % If the beam is on and sufficiently close to the island
      if (u(iupow(ibeam)) > 0)
        if (abs(rho_bar_mn(intm)-rho_bar_dep(ibeam))<=w+params.ntm.w_cd_ntm(intm))
          % Get only the effect of this beam
          %u_beam=u_beamsoff;
          %u_beam(iupow(ibeam))=u(iupow(ibeam));
          %[~,jec] = echcd(x,g,v,u_beam,model,params.echcd); %% replace by direct jecperunit from trap
          jec = trap.jecperunit(:,ibeam);
          % Integrate to find the the total current and the
          % multiplication factor due to the Gaussian profile
          Icd=int_SjdotB_edge(jec*model.equi.B0,geop,model);
          expdw=int_SjdotB_edge(exp(-4*(model.rgrid.rhogauss-rho_dep(ibeam)).^2./params.echcd.wdep(ibeam).^2)*model.equi.B0,geop,model);
          % The peak current for the original ECH profile is found by dividing Icd by Expdw
          % The smaller width for island stabilization needs to be
          % taken into account and the normalized wdep from echcd
          % needs to be multiplied  by a=epsilon*R0.
          j_cd_ntm=Icd(end)/expdw(end)*(params.echcd.wdep(ibeam)*a_minor)/params.ntm.w_cd_ntm(intm);
          % Determine the inefficiency as a result of the offset
          % rho_bar_dep-rho_bar_mn.
          x_norm=abs(rho_bar_dep(ibeam)-rho_bar_mn(intm))/max(w,params.ntm.w_cd_ntm(intm));
          % The efficiency Gcd as determined by De Lazzari is
          % approximated by an equation by Olivier Sauter.
          gcd_ofset=(1+g_coeff)*(1-tanh(3.75*x_norm-1.5))/(1-tanh(-1.5)+2*x_norm^3)-...
            g_coeff*exp(-x_norm^2);
          % Determine current drive peak, with efficiency correction due to offset
          rhodeltap_cd=rhodeltap_cd+rho_bar_dep(ibeam)*j_cd_ntm*gcd_ofset;
          % heating efficiency correction due to offset (Eq.29 of De Lazzari NF 2009)
          if h_fact>0
            ggh = 3.5e-4*wstar.^4 - 8e-3*wstar.^3 + 0.07*wstar.^2 + 0.02*wstar + 0.5;
            gh_ofset =  exp(-(x_norm/ggh)^2);
            eta_H = 3 * params.ntm.w_cd_ntm(intm)^2 * joh_mn / ...
              (8 * pi * model.equi.R0* den_mn * chie_mn * Te_mn);	% Eq.26 of De Lazzari NF 2009
            % rhodeltap_h = rhodeltap_h+rho_bar_dep(ibeam)*eta_H*u(iupow(ibeam))*gh_ofset;
            rhodeltap_h = rhodeltap_h + eta_H*u(iupow(ibeam))*gh_ofset;
          end
        end
      end
    end
    % Use the auxilliary power CD efficiency as defined in (18) of Sauter PPCF2010
    rhodeltap_cd=rhodeltap_cd*0.25./(1+2/3*wstar^2);
    % Calculate the total Delta'_CD term
    rhodeltap_cd = -8*sqrt(pi)*mu0*model.equi.R0*q_mn/(s_mn*model.equi.B0)* ...
      rhodeltap_cd/params.ntm.w_cd_ntm(intm);
    % Use the auxilliary power H efficiency as defined in (28) of De Lazzari NF 2009
    if h_fact>0
      rhodeltap_h = rhodeltap_h*(0.27*wstar^2+0.39*wstar)/(wstar^2+wstar+2.5);
      % Calculate the total Delta'_H term
      % rhodeltap_h = -8*sqrt(pi)*mu0*model.equi.R0*q_mn/(s_mn*model.equi.B0)* ...
      %		rhodeltap_h/params.ntm.w_cd_ntm(intm);
      rhodeltap_h = -16*mu0*model.equi.R0*q_mn/(pi*s_mn*model.equi.B0)* ...
        rhodeltap_h/params.ntm.w_cd_ntm(intm)^2;
      rhodeltap_h = h_fact * rhodeltap_h;
    end
    %% terms output
    dwdt(intm)  = (rhodeltap_cl + rhodeltap_cd + rhodeltap_h + ...
      rhodeltap_bs + rhodeltap_ggj)/pre_dwdt(intm);
    
  end
  %% Add values to local structure
  if params.ntm.debug
    log_rhodeltap_cl(intm,it) = rhodeltap_cl;
    log_rhodeltap_ggj(intm,it)= rhodeltap_ggj;
    log_rhodeltap_bs(intm,it) = rhodeltap_bs;
    log_rhodeltap_cd(intm,it) = rhodeltap_cd;
    log_rhodeltap_h(intm,it)  = rhodeltap_h;
    log_eta_H(intm,it)        = eta_H;
    log_pre_dwdt(intm,it)     = pre_dwdt;
    log_j_bs_mn_up(intm,it)   = j_bs_mn_up;
    log_L_bs_mn(intm,it)      = L_bs_mn;
    log_j_cd_ntm(intm,it)     = j_cd_ntm;
    log_dwdt(:,it)            = dwdt;
    log_rho_mn(:,it)          = rho_mn;
    log_Ncd(intm,it)          = 0.25./(1+2/3*(w./params.ntm.w_cd_ntm(intm)).^2);
    log_a_bs(intm,it)         = a_bs;
    log_a_ggj(intm,it)        = a_ggj;
  end
end
