function [iR,iR2,dpsi,dpsi2,dpsi2oR2,dpsidV,F,epsilon,Rgeom,kappa,delta,Phib,tgridin,rhodata] = ASTRA_to_RAPTOR(config)
% function [iR,iR2,dpsi,dpsi2,dpsi2oR2,dpsidV,F,epsilon,Rgeom,delta,Phib,tgridin,rhodata] = ASTRA_to_RAPTOR(config)
% Loads ASTRA data from .mat file
% and returns geometric quantities, time and radial grids.
%% Load in data from ASTRA
% full path to ASTRA
assert(numel(config.equi.filenames) == 1, 'RAPTOR:SingleAstraEqFile', 'There must 1 single ASTRA file for equilibrium')

astrafile = config.equi.filenames;
if ~isempty(config.envopts.equil_path)	
   astradatapath = config.envopts.equil_path;
else
    astradatapath = fullfile(config.envopts.RAPTOR_path,'equils', 'ASTRA');
end
% Path to ASTRA data
astrafilepath = fullfile(astradatapath,astrafile{1});
% contains ASTRA data for equilibrium
astra = load(astrafilepath);

% Check data file structure
if isfield(astra,'TCV_astra')
    astra = astra.TCV_astra;
end
%
if strfind(config.equi.filenames{1},'RES')
    % for EF data files
    data.btor = astra.data.btor';
    data.rho = astra.data.rho;
    data.vr = astra.data.vr;
    data.ipol = astra.data.ipol;
    data.g11 = astra.data.g11;
    data.g22 = astra.data.g22;
    data.g33 = astra.data.g33;
    data.grrho = astra.out.gradro;
    data.sqeps = astra.data.sqeps;
    data.rtor = astra.data.rtor;
    data.shif = astra.data.shif;
    data.kappa = astra.data.elon;
    data.tria = astra.data.tria;
    data.t = astra.data.t'; 
    data.iota = astra.data.mu;
else    
    data.btor = astra.out.BTOR;
    data.rho = astra.out.RHO;
    data.vr = astra.out.VR;
    data.ipol = astra.out.IPOL;
    data.g11 = astra.out.G11;
    data.g22 = astra.out.G22;
    data.g33 = astra.out.G33;
    data.grrho = astra.out.GRADRO;
    data.sqeps = astra.out.SQEPS;
    data.rtor = astra.out.RTOR;
    data.shif = astra.out.SHIF;
    data.kappa = astra.out.ELON;
    data.tria = astra.out.TRIA;
    data.t = astra.out.T;    
    data.iota = astra.out.MU;
end

if isempty(config.equi.times)
    indt = [1:1:numel(data.t)];
else
     it = config.equi.times;
        if it==inf
            indt = numel(data.t);
        else           
            diff = abs(data.t - it);
            indt = find(diff==min(diff));
        end
end

%
% experimental data used for scaling
% Main radius
R0 = config.equi.R0;
% Magnetic field: constant or time-varying parameter
B0 = data.btor(:,indt);
%
%% interpolate data, to match rhotorN grid of ASTRA
% Normalized radial grid (rho toroidal normalized)
rhotorN = bsxfun(@rdivide,data.rho(:,indt),data.rho(end,indt));
rhotorN(1,:)= 0.0;
%
F_a = bsxfun(@times,data.ipol(:,indt),B0)*R0;
epsilon_a = data.sqeps(:,indt).*data.sqeps(:,indt);
Rgeom_a = data.rtor(1)+data.shif(:,indt);
kappa_a = data.kappa(:,indt);
delta_a = data.tria(:,indt);
Phi_a = pi* bsxfun(@times,data.rho(:,indt).* data.rho(:,indt),B0);

%%% NOTE on definition g1,g2 and g3 ASTRA and RAPTOR.
% RAPTOR needs g1, g2 and g3 as defined in the paper of E. Fable (Nucl fusion 2013) 
% In ASTRA g11 = V'*G1, g22 = R0*G2/J, g33 = G3 
% Source : http://www.ipp.mpg.de/~git/astra/section.php?sec=36&str=G2
% In the data set of Emiliano Fable b_flattop.g33 is defined as <1/(R^2) >, in the
% data set generated at TCV b_flattop.g33 is defined as <R0^2/R^2>
% In order to calculate g33_ASTRA (<1/R^2>) a check is made to
% automatically identify which definition of b_flattop.g33 was used;

% In the definition <R0^2/R^2> for b_flattop.g33 the array should contain a
% value close to one (around the magnetic axis)
if any ((data.g33(data.g33(:,indt) > 0.95)) < 1.05 )
    iR2_a = (data.g33(:,indt))/(R0^2); % = <1/R^2>, TCV data use this convention 
else
    iR2_a = (data.g33(:,indt)); % = <1/R^2>, ASDEX data use this convention 
    fprintf(' The <1/(R^2)> definition was used for g33 \n');
end

% dpsi/dV: goes to inf near axis
dpsidV_a = bsxfun(@times,2.0*pi*data.iota(:,indt).*data.rho(:,indt)./data.vr(:,indt),B0);
coeff = data.vr(2,indt)./data.rho(2,indt);
dpsidV_a(1,:) = 2.0*pi*B0.*data.iota(1,indt)./coeff;
% <|grad psi|>
dpsi_a = data.vr(:,indt).* data.grrho(:,indt).*dpsidV_a;
% <|grad psi|^2>
dpsi2_a = data.vr(:,indt).* data.g11(:,indt).*dpsidV_a.^2;
% <(|gard psi|/R)^2>
dpsi2oR2_a = (4.0*(pi^2) * data.ipol(:,indt)/R0).* data.vr(:,indt).* data.g22(:,indt).*dpsidV_a.^2;
% <1/R>
iR_a = sqrt(iR2_a);

% decimate if too large: radial grid
nASTRA = size(F_a,1);
if nASTRA>20;
    % reduce
    ndec = fix(nASTRA/20);
    % also eliminate first point
    ij = unique([2:ndec:nASTRA,nASTRA]);
end

% decimate if too large: time grid
if numel(indt)>1
    ndt = (config.grid.tgrid(2)-config.grid.tgrid(1))/(data.t(2)-data.t(1));
    nASTRAt = size(F_a,2);
    if ndt>10
        % also eliminate first point
        ik = unique([2:10:nASTRAt,nASTRAt]);
    else
        ik = [1:1:nASTRAt];
    end
else
    ik = 1;
end

iR = iR_a(ij,ik);
iR2 = iR2_a(ij,ik);
dpsi = dpsi_a(ij,ik);
dpsi2 = dpsi2_a(ij,ik);
dpsi2oR2 = dpsi2oR2_a(ij,ik);
dpsidV = dpsidV_a(ij,ik);
F = F_a(ij,ik);
epsilon = epsilon_a(ij,ik);
Rgeom = Rgeom_a(ij,ik);
delta = delta_a(ij,ik);
kappa = kappa_a(ij,ik);
Phib = Phi_a(end,ik);
rhodata = rhotorN(ij,ik);
tgridin = data.t(ik);

%% checks
% % % Check volume:
% V_check = cumtrapz(rhotor,b_flattop.vr);
% V = b_flattop.volum;
% figure;plot(rhotorN,V_check,rhotorN,V,'r')
% legend('RAPTOR', 'ASTRA')
% title('Volume')
% 
% % % Check G3:
% G3_check = 4*pi^2*rhotor*R0./(eq_profiles.J.*eq_profiles.Vprime);
% G3 = (eq_profiles.g3).*R0^2;
% figure;plot(G3,'r');hold on;plot(G3_check)
% title('G3')
% 
% % % Check of the q profile of ASTRA and RAPTOR
% iota_raptor = 1./(2*b_flattop.PHI*rhotorN(2:end)).*diff(eq_profiles.psi)./diff(rhotorN);
% figure;
% plot(rhotorN(2:end),1./iota_raptor,rhotorN(1:end),eq_profiles.q,'r')
% legend('RAPTOR', 'ASTRA')
% title('q profile')
% 
% % % Check of the Ip profile
% F = eq_profiles.J*R0*B0;
% mu0 = 4e-7*pi;
% Ip_raptor = ((F(2:end))./(8*pi^3 * mu0 * b_flattop.PHI * 2 * rhotorN(2:end))) .*eq_profiles.g2(2:end) .*eq_profiles.g3(2:end) .*diff(eq_profiles.psi)./diff(rhotorN);
% Ip_raptor2 = (eq_profiles.G2(2:end)./mu0).*diff(eq_profiles.psi)./diff(rhotor);
% figure;
% plot(rhotorN(2:end),Ip_raptor,rhotorN(1:end-1),1e6*b_flattop.itot(1:end-1) ,'r')
% legend('RAPTOR', 'ASTRA','location','best')
% title('I_p')

 return
