function [NNinputs,dNNinputs_dx,kfac,dkfac_dx,QLKNNparams] = RAPTORtoQLKNN10D(stap,geop,trap,model,QLKNNparams)

% QLKNN inputs:
% R/LTi = -dTi/dr * R0/Ti
% R/LTe = -dTe/dr * R0/Te
% R/Lne = -dne/dr * R0/ne
% Ti/Te
% q = q-profile
% s = magnetic shear
% -- alternatively s-.5*alpha (where alpha = dbeta/dr * q^2 * R0)
% -- gives good match with full QLK for high shear hybrids
% x = r / a   (where a is the minor radius and r is the local coordinate in the midplane averaged radius coordinate r = (Rmax-Rmin)/2 , where Rmax and Rmin correspond to the local flux surface and =r_LCFS)
% nustar (defined in PDF)
% Zeff (defined in PDF)
% gammaE_GB =abs( dvtor/dr) *epsilon/q / (c_sou/a)    where c_sou=sqrt(Te/mi) and epsilon = x * a/R0 
%
% Order into QLKNN: [Zeff, Ati, Ate, An, q, smag, x, Ti_Te, ln(Nustar), gammaE_GB or Machtor]

%% some physics constants
mp = 1.67e-27; %proton mass
sqrtmi = sqrt(model.atom.Ai*mp);
qel = 1.6e-19;
cref = sqrt(qel * 10^3 / mp);

%% Get physical quantities from state
te      =   stap.te; dte_dx  =   stap.dte_dx;
ti      =   stap.ti; dti_dx  =   stap.dti_dx;
tep     =   stap.tep; dtep_dx =   stap.dtep_dx;
tip     =   stap.tip; dtip_dx =   stap.dtip_dx;
q       = stap.q; dq_dx   = stap.dq_dx;
shear   = stap.shear; dshear_dx = stap.dshear_dx;

ne = stap.ne; dne_dx = stap.dne_dx;
nep = stap.nep; dnep_dx = stap.dnep_dx;

ze = stap.ze; dze_dx = stap.dze_dx;
vtp = stap.vtp; dvtp_dx = stap.dvtp_dx;

%% 
if isfield(QLKNNparams, 'shear_minus_alpha') && QLKNNparams.shear_minus_alpha
  ni = stap.ni; dni_dx = stap.dni_dx;
  nip = stap.nip; dnip_dx = stap.dnip_dx;
  
  piprime = 1.6e-19*(tip.*ni+ti.*nip);
  dpiprime_dx = 1.6e-19*(...
        + bsxfun(@times,ni,dtip_dx) ...
        + bsxfun(@times,tip,dni_dx) ...
        + bsxfun(@times,ti,dnip_dx) ...
        + bsxfun(@times,nip,dti_dx) );
  peprime = 1.6e-19*(tep.*ne+te.*nep); 
  dpeprime_dx = 1.6e-19*(...
          bsxfun(@times,ne,dtep_dx) ...
        + bsxfun(@times,tep,dne_dx) ...
        + bsxfun(@times,te,dnep_dx) ...
        + bsxfun(@times,nep,dte_dx)  );
      
  pprime = peprime + piprime;
  dpprime_dx = dpeprime_dx + dpiprime_dx;
      
  beta_constant = (2*4e-7*pi)/(model.equi.B0.^2);
  betaprime = beta_constant.*pprime;
  dbetaprime_dx = bsxfun(@times,beta_constant,dpprime_dx);
end

%% Geometry issues
[r,drhodr] = drho_dr(geop); % rho to r conversion, r = epsilon(psi)*Rgeom(psi)
if isfield(QLKNNparams, 'force_Rmin')
  aLCFS = QLKNNparams.force_Rmin;
else
  aLCFS = r(end); % effective minor radius for space normalization
end
epsilon = geop.epsilon;
epsilonNN = 1/3; % epsilon used to train 10D NN (see documentation)
xin = r/aLCFS * epsilon(end)/epsilonNN; % normalized x
if isfield(QLKNNparams, 'force_R0')
  R0 = QLKNNparams.force_R0;
else
  R0 = geop.Rgeom(end); % LCFS R0 for R0/LTi etc
end

%% Calculate QLKNN inputs and derivatives
% % transform gradients from d/drhotor to d/dr for QuaLiKiz

% gradients are given w.r.t. rhotorN, need w.r.t. r = (Rmax-Rmin)/2
% since this is the way it's done in CRONOS/Qualikiz as well
[tepr,dtepr_dx] = ddrho_to_ddr(tep,dtep_dx,drhodr);
[tipr,dtipr_dx] = ddrho_to_ddr(tip,dtip_dx,drhodr);
[nepr,dnepr_dx] = ddrho_to_ddr(nep,dnep_dx,drhodr);
[vtpr,dvtpr_dx] = ddrho_to_ddr(vtp,dvtp_dx,drhodr);
[shearr,dshearr_dx] = ddrho_to_ddr(shear,dshear_dx,drhodr);
if isfield(QLKNNparams, 'shear_minus_alpha') && QLKNNparams.shear_minus_alpha
  [betaprimer,dbetaprimer_dx] = ddrho_to_ddr(betaprime,dbetaprime_dx,drhodr);
  alphar = -R0 .* betaprimer .* q.^2;
  dalphar_dx = -R0 .* ( bsxfun(@times, q.^2, dbetaprimer_dx) + bsxfun(@times, 2.*q.*betaprimer, dq_dx) );
end
% sRAP = rhotor/q dq/drhotor
% sQLK = r/q dq/dr
% sQLK = rhotor/r * (sRAP) * dr/drhotor

rlti = -R0 .* (tipr) ./ ti; 
rlte = -R0 .* (tepr) ./ te; 
rlne = -R0 .* (nepr) ./ ne; 
tite = ti./te ;
shearr = shearr .* r ./ model.rgrid.rhogauss;

drlti_dx = bsxfun(@times,rlti./tipr,dtipr_dx) - bsxfun(@times,rlti./ti,dti_dx);
drlte_dx = bsxfun(@times,rlte./tepr,dtepr_dx) - bsxfun(@times,rlte./te,dte_dx);
drlne_dx = bsxfun(@times,rlne./nepr,dnepr_dx) - bsxfun(@times,rlne./ne,dne_dx);
dtite_dx = bsxfun(@times,tite./ti,dti_dx) - bsxfun(@times,tite./te,dte_dx);
dshearr_dx = bsxfun(@times, dshearr_dx, r ./ model.rgrid.rhogauss);

% x = r/a;
dxin_dx = zeros(size(dtite_dx)); % geometry does not depend on x

% ln(nuestar)
% Calculation of nuestar and partial derivatives
 [nuestar,dnuestar_dne,dnuestar_dte,dnuestar_dze,dnuestar_dq] = ...
   nue_star(ne,te,ze,q,R0,epsilon);
 
dnuestar_dx = bsxfun(@times,dnuestar_dne,dne_dx) + ...
  bsxfun(@times,dnuestar_dte,dte_dx) + ...
    bsxfun(@times,dnuestar_dze,dze_dx) + ...
    bsxfun(@times,dnuestar_dq,dq_dx);
  
log10nuestar = log10(nuestar);
dlnnuestar_dx = bsxfun(@times,1./nuestar/log(10),dnuestar_dx);
if isfield(QLKNNparams, 'coll_proxy') && QLKNNparams.coll_proxy
  log10nuestar = log10nuestar./4;
  dlnnuestar_dx = dlnnuestar_dx./4;
end

% GammaE_GB
gamE = vtpr .* r ./ q * R0 / cref;
dgamE_dx = bsxfun(@times, R0 / cref .* r.^2 ./ q, dvtpr_dx) - ...
  bsxfun(@times, gamE ./ q, dq_dx);

% Victor rule normalization parameters
QLKNNparams.norms.A1 = model.atom.Ai;
QLKNNparams.norms.a = aLCFS;
QLKNNparams.norms.R0 = R0;

if QLKNNparams.em_stab == 1
  % Calculate total pressure. Taken from eval_press.m
  niall = stap.ni+stap.n1+stap.n2+stap.n3;
  dniall_dx = stap.dni_dx + stap.dn1_dx + stap.dn2_dx + stap.dn3_dx;

  press = 1.602e-19*(te.*ne + ti.*niall);
  dpress_dx = 1.602e-19*(...
      bsxfun(@times,ne,dte_dx) + bsxfun(@times,niall,dti_dx) + ...
      bsxfun(@times,te,dne_dx) + bsxfun(@times,ti,dniall_dx));

  % rlti of all ions should be adjusted, but QLKNN10D just sees the main ion
  % Note the (3/2) difference with JETTO, J/m3 vs Pa
  drlti_dx = bsxfun(@times, (3/2) * 1./(press + trap.pfast).^2, ...
    (bsxfun(@times, trap.pfast .* rlti, dpress_dx) + ...
      bsxfun(@times, press .* (trap.pfast + press), drlti_dx)));
    
  rlti = (3/2) .* rlti .* (press ./ (press + trap.pfast));
end

%% Assemble inputs
% Order into QLKNN: [Zeff, Ati, Ate, An, q, smag, x, Ti_Te, ln(Nustar),gammaE_GB, Te]
if isfield(QLKNNparams, 'force_q')
    q_orig = q;
    dq_dx_orig = dq_dx;
    q = QLKNNparams.force_q;
    dq_dx = zeros(size(dq_dx));
end
if isfield(QLKNNparams, 'force_s')
    s_orig = shearr;
    dshearr_dx_orig = dshearr_dx;
    shearr = QLKNNparams.force_s;
    dshearr_dx = zeros(size(dshearr_dx));
end
if isfield(QLKNNparams, 'shear_minus_alpha') && QLKNNparams.shear_minus_alpha
  alphfac = .5; % from tuning with full QLK for hybrid
  shearr = shearr - alphfac.*alphar;
  dshearr_dx = bsxfun(@minus,dshearr_dx, alphfac.*dalphar_dx) ;
end

NNinputs = [ze,rlti,rlte,rlne,q,shearr,xin,tite,log10nuestar,gamE,te];
dNNinputs_dx = cat(3,dze_dx,drlti_dx,drlte_dx,drlne_dx,...
  dq_dx,dshearr_dx,dxin_dx,dtite_dx,dlnnuestar_dx,dgamE_dx,dte_dx);

%% GB normalization for outputs
B0 = model.equi.B0;
%Rescale temperatures to SI units
teSI  = te*qel; 

%Define gyroBohm normalization. NN outputs are in GB normalized units
%chiGB = sqrt(mi)/(a*e^2*B^2)*Te^1.5

teSI12 = sqrt(teSI);
teSI32 = teSI12.*teSI;
chifac = sqrtmi/(qel^2*B0^2)/aLCFS;
chiGB  = teSI32.*chifac; % for chie,chii,de
dchiGB_dte = 3/2*qel*teSI12 * chifac;
dchiGB_dx = bsxfun(@times,dchiGB_dte,dte_dx);

%% Specific conversion factors per output (from documentation table)
irlte = 1./rlte;
irlti = 1./rlti;
irlne = 1./rlne;

kchie = R0/aLCFS .* chiGB .* irlte;
kchii = R0/aLCFS .* chiGB .* irlti;
kDeff   = R0/aLCFS .* chiGB .* irlne;
kD   = chiGB;
kV  = chiGB / aLCFS;

dkchie_dx = R0/aLCFS .* bsxfun(@times,irlte,dchiGB_dx) - bsxfun(@times,kchie.*irlte,drlte_dx);
dkchii_dx = R0/aLCFS .* bsxfun(@times,irlti,dchiGB_dx) - bsxfun(@times,kchii.*irlti,drlti_dx);
dkDeff_dx = R0/aLCFS .* bsxfun(@times,irlne,dchiGB_dx) - bsxfun(@times,kDeff.*irlne,drlne_dx);
dkD_dx = dchiGB_dx;
dkV_dx = dchiGB_dx / aLCFS;

%% matrix of k factors mapping NN outputs to transport coefficients (see documentation)
kfac = [kchie,kchii,kDeff,kD,kV];
dkfac_dx = cat(3,dkchie_dx,dkchii_dx,dkDeff_dx,dkD_dx,dkV_dx);
end

function [ypr,dypr_dx] = ddrho_to_ddr(yp,dyp_dx,drhodr)
% apply drho_dr to gradients and gradient derivatives 
ypr = yp .* drhodr;
dypr_dx = bsxfun(@times,drhodr,dyp_dx); 
end

function nustar = calc_nustar_qlk(ne_SI,te_ev,zeff,q,a,R0,x,eps)
qel = 1.6e-19;
me = 9.10938356e-31;

ne = ne_SI / 10^19;
te = te_ev / 10^3;

Lam = 15.2 - 0.5 .* log(0.1 .* ne) + log(te);
nu = 914.7 .* zeff .* ne .* Lam .* te.^-1.5;
bounce_xin = q .* R0 ./ ((x .* a ./ R0).^1.5 .* sqrt(10^3 .* te .* qel ./ me));
bounce = q .* R0 ./ ((eps).^1.5 .* sqrt(10^3 .* te .* qel ./ me));
nustar = nu .* bounce;
end
