function  [g,rhotornorm] = LY2RAPTOR(psiN,LY,geom,gprev)


[iR,iR2,dpsi,dpsi2,dpsi2oR2,dpsidV,F,iota,psib,psia] = LY2geom(LY);


[Phib,rhotornorm,~] = calcPhi(psiN,iota,psib,psia);


% dummy, substitute later
epsilon = ones(size(iR2));
delta = ones(size(iR2));
kappa = ones(size(iR2));
Rgeom = ones(size(iR2));

g = geom2RAPTOR(rhotornorm,iR,iR2,dpsi,dpsi2,dpsi2oR2,dpsidV,F,epsilon,Rgeom,kappa,delta,Phib,geom);

% substitute default epsilon,Rgeom,delta while they are not given by RT-LIUQE
ii = [geom.ind_epsilon, geom.ind_Rgeom, geom.ind_kappa, geom.ind_delta];
g(ii) = gprev(ii);


% 
% keyboard
% 
% psi = psiN*(LY.FB- LY.FA)+ LY.FA;
% dVdrhotor = 4*pi*Phib*rhotornorm./(F.*LY.Q2Q);
% 
% [dVdpsi_interp, ~, ~, V_psi] = interpos(psi, -1./LY.Q1Q, psi, -0.001);
% [dVdrho_interp, ~, ~, V_rhotor] = interpos(rhotornorm, dVdrhotor, rhotornorm, -0.001);
% 
% figure
% plot(psi, -1./LY.Q1Q,'.',  psi,dVdpsi_interp, '-')
% figure
% plot(rhotornorm, dVdrhotor,'.',  rhotornorm,dVdrho_interp, '-')
% figure
% plot(psiN, V_psi, psiN, V_rhotor) 



% 
%  keyboard
% Ip = (4*pi*4*pi*pi*1e-7)^-1*LY.Q3Q.*LY.Q1Q.^-1;
% psi = psiN*(psib-psia)+psia;
% [~, dpsidrho]  = interpos(rhotor, psi, -0.01, [1 0],[0 0]);  
% [~, dIpdpsi] = interpos(psi, Ip);
% [~, ~,~, Volume] = interpos(psi, 1./dpsidV);
% jtor2 = dIpdpsi.*dpsidV;
% [~,jtor1] = interpos(Volume, Ip);
% gradrho2R2 = dpsi2oR2./dpsidV.^2./dpsidrho.^2;
% plot(psiN, jtor1, psiN, jtor2)
% 
% figure
% plot(psiN, dpsidV.^2, '*')
% plot(psiN, dpsidrho)
% plot(psiN, gradrho2R2, '*')

% Debug
% psi = psiN*(psib-psia)+psia;
% [~, ~, psipp] = interpos(rhotor, psi, -0.01);
% figure
% plot(psiN, psipp);
% figure(400)
% hold on 
% plot(g(geom.ind_Vp(1:end)))


return