function [qlkNNinputs,dqlkNNinputs_dx,GBfac,dGBfac_dx] = RAPTORtoQLKNN4Dkin(stap,geop,model)

% Get physical quantities from state
te      =   stap.te; dte_dx  =   stap.dte_dx;
ti      =   stap.ti; dti_dx  =   stap.dti_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;

R0 = geop.Rgeom(end); % LCFS R0 for R0/LTi etc

% % transform gradients from d/drhotor to d/dr to be the same as CRONOS
[~,drhodr] = drho_dr(geop);
[tipr,dtipr_dx] = ddrho_to_ddr(tip,dtip_dx,drhodr);

%NN trained for dimensionless R/LTi = -R*grad(Ti)/Ti , Ti/Te , q, magnetic shear.
%Not some as function input, thus, must initialize input vector

rlti = -R0 .* (tipr) ./ ti; 
tite = ti./te ;

drlti_dx = bsxfun(@times,rlti./tipr,dtipr_dx) - bsxfun(@times,rlti./ti,dti_dx);
dtite_dx = bsxfun(@times,tite./ti,dti_dx) - bsxfun(@times,tite./te,dte_dx);

qlkNNinputs = [q,shear,tite,rlti];
dqlkNNinputs_dx = cat(3,dq_dx,dshear_dx,dtite_dx,drlti_dx);

%% GB normalization
R0 = geop.Rgeom(end); % LCFS R0 for R0/LTi etc
a0 = geop.epsilon(end).*geop.Rgeom(end); %LCFS a0 for gB units
B0 = model.equi.B0;
Amain = model.atom.Ai;
qel = 1.6e-19;

%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 = sqrt(Amain*1.67e-27)/(qel^2*B0^2*a0);
GBfac  = teSI32.*chifac; % for chie,chii,de

dGBfac_dte = 3/2*qel*teSI12 * chifac;
dGBfac_dx = bsxfun(@times,dGBfac_dte,dte_dx);

% array of GB factors per output
GBfac     = [GBfac,GBfac,GBfac,GBfac/R0];
dGBfac_dx = cat(3,dGBfac_dx,dGBfac_dx,dGBfac_dx,dGBfac_dx/R0);
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

