function [iR,iR2,dpsi,dpsi2,dpsi2oR2,dpsidV,F,epsilon,Rgeom,kappa,delta,Phib,tgridin,rhodatain] = LIUQEfile2RAPTOR(config)
% Get the LIUQE output from a file which store LY LIUQE output. 
% LIUQE had to be run with 'itert' > 1 option to activated the geometrical
% post processing similar to what done in RT liuqe.

fileout = load(config.equi.filenames);

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

% Need to find where rhotor is used as well as the other 2 outputs
[Phib,rhodatain,~] = calcPhi(fileout.L.P.pq.^2',iota,psib,psia);

% numel(rho_index) = numel(L.P.pq.^2) - 1 since the first value of
% L.P.pq.^2 is 0 and the flux surface is not defined.
% I will put the first values equal to the last ones
rq = fileout.LY.rq; % ([r], rho_index, time_index) 
zq = fileout.LY.zq; % ([z], rho_index, time_index) 

np  = size(rq,1); % number of point for the flux surface
nrho = size(rq,2); % number of rho points
ntime = size(rq,3); % number of time points

R_IN = zeros(nrho,ntime);
R_OUT = zeros(nrho,ntime);
R_IN_MID = zeros(nrho,ntime);
R_OUT_MID = zeros(nrho,ntime);
R_TOP = zeros(nrho,ntime);
R_BOT = zeros(nrho,ntime);
Z_TOP =  zeros(nrho,ntime);
Z_BOT =  zeros(nrho,ntime);
DELTA_BOT =  zeros(nrho,ntime);
DELTA_TOP =  zeros(nrho,ntime);
A_MINOR_MID = zeros(nrho,ntime);
R_GEOM_MID = zeros(nrho,ntime);

% Computed from LCFS R,Z points epsilon(rho,t)
epsilon = ones(size(iR2)); % inverse aspect ratio
delta = ones(size(iR2)); % triangularity
kappa = ones(size(iR2)); % elongation
Rgeom = ones(size(iR2)); % geometrical major radius

for tt=1:ntime
    for rr=1:nrho
        
        R_IN(rr,tt) = min(rq(:,rr,tt));
        R_OUT(rr,tt) = max(rq(:,rr,tt));
        R_IN_MID(rr,tt) = R_IN(rr,tt); % I should get the value at the midplane but this is a good approx
        R_OUT_MID(rr,tt) = R_OUT(rr,tt); % I should get the value at the midplane but this is a good approx
        [Z_TOP(rr,tt), index_top] = max(zq(:,rr,tt));
        [Z_BOT(rr,tt), index_bottom]= min(zq(:,rr,tt));
        R_TOP(rr,tt) = rq(8,rr,tt); % Since the discretization of the point is not good also this are noisy. I take a fix index for the maximum which is wrong.
        R_BOT(rr,tt) = rq(27,rr,tt); % Since the discretization of the point is not good also this are noisy. I take a fix index for the maximum.
        A_MINOR_MID(rr,tt) = (R_OUT_MID(rr,tt)-R_IN_MID(rr,tt))/2;
        R_GEOM_MID(rr,tt) = (R_OUT_MID(rr,tt)+R_IN_MID(rr,tt))/2;
        
         % Warning this has to be changed. Due to few point in the LCFS in rq and zq the r(z_max) is jumping between 2 points making R_TOP noisy.
                                    % Since this quantities are used just
                                    % for the evaluation of the trap
                                    % particle and the triangularity
                                    % doesn't vary much I will take one
                                    % from good flux surface and keep it.
        
        DELTA_TOP(rr,tt) = (R_IN(rr,tt)+R_OUT(rr,tt)-2*R_TOP(rr,tt))/(R_OUT(rr,tt)-R_IN(rr,tt));
        DELTA_BOT(rr,tt) = (R_IN(rr,tt)+R_OUT(rr,tt)-2*R_BOT(rr,tt))/(R_OUT(rr,tt)-R_IN(rr,tt));
        
        % Output requested by RAPTOR
        kappa(rr+1,tt) = (Z_TOP(rr,tt)- Z_BOT(rr,tt))/(R_OUT(rr,tt)-R_IN(rr,tt));
        delta(rr+1,tt) = DELTA_TOP(rr,tt);%(DELTA_TOP(rr,tt)+DELTA_BOT(rr,tt))/2; % Since the discretization is bad I will take just delta top
        Rgeom(rr+1,tt) = R_GEOM_MID(rr,tt);
        epsilon(rr+1,tt) = A_MINOR_MID(rr,tt)/R_GEOM_MID(rr,tt);
    end
    kappa(1,tt) = kappa(2,tt);
    delta(1,tt) = delta(2,tt);
    Rgeom(1,tt) = Rgeom(2,tt);
    epsilon(1,tt) = epsilon(2,tt);
    
end

tgridin = fileout.LY.t;

% Put some smoothing
for  rr=1:nrho+1
    kappa(rr,:) = interpos(tgridin, kappa(rr,:),tgridin, -10);
    delta(rr,:) = interpos(tgridin, delta(rr,:) , tgridin, -10);
    Rgeom(rr,:) = interpos(tgridin, Rgeom(rr,:), tgridin, -10);
    epsilon(rr,:) = interpos(tgridin, epsilon(rr,:), tgridin, -10);
end

end