function [iR,iR2,dpsi,dpsi2,dpsi2oR2,dpsidV,F,epsilon,Rgeom,kappa,delta,Phib,tgridin,rhotor,config] = CHEASE2SI(config)
% function [iR,iR2,dpsi,dpsi2,dpsi2oR2,dpsidV,F,epsilon,Rgeom,kappa,delta,Phib,tgridin,rhotor] = CHEASE2SI(config)
% Loads CHEASE data from .mat file
% and returns geometric quantities, time and radial grids.

% CHEASE file name(s)
cheasefile = config.equi.filenames;
if ~isempty(config.envopts.equil_path)
    cheasedatapath = config.envopts.equil_path;
else
    cheasedatapath = fullfile(config.envopts.RAPTOR_path,'chease_equils');
end
if iscell(cheasefile)
    nbt = size(cheasefile,2);
else
    nbt = 1;
end
% time grid for input data
tgridin = config.equi.times;

for ii=1:nbt    
    % Check if cheasefile is a cell array
    assert(iscell(cheasefile),'chease file must be cell');
    cheasefilepath = fullfile(cheasedatapath,cheasefile{ii});
     
    if config.debug.verbosity>0
        fprintf('loading CHEASE equilibrium from %s\n',cheasefilepath);
    end

    load(cheasefilepath); % contains CHEASE data for equilibrium

    if ~exist('chease_data','var') && exist('cheasedata','var')
        chease_data = cheasedata;
    end
    
    % Magnetic field: constant or time varying parameter
    if isempty(config.equi.B0)
        B0(ii) = chease_data.eqdsk.b0;
        if ii==1
            config.equi.B0 = B0(1);
        end
    else
        B0(ii) = config.equi.B0;
    end

    if config.debug.verbosity>0
        fprintf('B0 = %2.2f\n',B0(ii));
    end

    if isfield(chease_data,'profiledata');
        b = chease_data.profiledata;
    else
        b = chease_data;
    end

    % experimental data from eq_params used for scaling   
    R0 = config.equi.R0;

    %
    %rhoe = b.rho_tor_eq_sqrt_phi_over_pi_over_b0.data(end) * R0;
    %rhotorN = b.rho_tor_norm.data ;

    iR2_ch = b.av_1_over_rpow2_av.data /(R0^2);
    iR_ch = b.av_1_over_r_avprofile.data/(R0);
    dpsi_ch = b.av_grad_psi_av.data* B0(ii)*R0;
    dpsi2_ch = b.av_grad_psi_pow2_av.data * B0(ii)^2*R0^2;
    dpsi2oR2_ch = b.av_bppow2_av.data * B0(ii)^2;
    dpsidV_ch = 1./b.dvdpsi.data * (B0(ii)/R0);
    
    %C1_ch = b.int_rdlp_over_grad_psi_eq_int_jdchi.data.* R0/B0;
    %C2_ch = b.av_1_over_rpow2_av.data.* b.int_rdlp_over_grad_psi_eq_int_jdchi.data./(R0*B0);
    %C3_ch = b.av_bppow2_av.data .* b.int_rdlp_over_grad_psi_eq_int_jdchi.data.* R0*B0;
    %C4_ch = b.av_grad_psi_pow2_av.data.* b.int_rdlp_over_grad_psi_eq_int_jdchi.data.* R0^3*B0;
    F_ch = b.t_eq_rbphi.data* B0(ii)*R0;

    % [TODO: Is there an easier way to get Phib directly from chease? 
    % I guess there should be]
    
    %q0 = b.qprofile.data(1);
    Psi_cheaseSi = b.psichease_eq_psi_over_2pi.data* R0^2* B0(ii); 
    s = b.s_mesh.data;
    C2_ch = b.av_1_over_rpow2_av.data.* b.int_rdlp_over_grad_psi_eq_int_jdchi.data./(R0*B0(ii));
    y = F_ch.*C2_ch;
    psi_edge = Psi_cheaseSi(end);
    psi_axis = Psi_cheaseSi(1);
    Phi_ch = cumtrapz(s.*s.*(psi_edge-psi_axis),y);

    %correct value for C1 and C2 at the first point
    %C2_ch(1) = 2*pi* q0/F_ch(1);
    %C1_ch(1) = C2_ch(1)* R0^2;

    Rgeom_ch = R0*b.rgeomprofile.data;
    epsilon_ch = b.a_over_rgeom.data;
    delta_up = b.delta_upper.data;
    delta_b = b.delta_bottom.data;
    delta_ch = (delta_up + delta_b)*0.5;
    kappa_ch = b.elongation.data;

    % decimate if too large
    nCHEASE = numel(iR_ch);
%     if nCHEASE>20;
%         % reduce
%         ndec = fix(nCHEASE/20);
%         % also eliminate first point
%     end
%    ndec = 1; %[TODO: decide how to fix this]
%    ij = unique([2:ndec:nCHEASE,nCHEASE]);
     ij = 2:numel(iR_ch);

    % prepare output profiles
    if nbt==1
        iR2 = iR2_ch(ij);
        iR  = iR_ch(ij);
        dpsi = dpsi_ch(ij);
        dpsi2 = dpsi2_ch(ij);
        dpsi2oR2 = dpsi2oR2_ch(ij);
        dpsidV = dpsidV_ch(ij);
        F = F_ch(ij);
        Rgeom = Rgeom_ch(ij);
        epsilon = epsilon_ch(ij);
        kappa = kappa_ch(ij);
        delta = delta_ch(ij);
        Phib = Phi_ch(end);
        rhotor = sqrt(Phi_ch(ij)/Phi_ch(end));
    else
        iR2_ca{ii} = iR2_ch(ij);
        iR_ca{ii}  = iR_ch(ij);
        dpsi_ca{ii} = dpsi_ch(ij);
        dpsi2_ca{ii} = dpsi2_ch(ij);
        dpsi2oR2_ca{ii} = dpsi2oR2_ch(ij);
        dpsidV_ca{ii} = dpsidV_ch(ij);
        F_ca{ii} = F_ch(ij);
        Rgeom_ca{ii} = Rgeom_ch(ij);
        kappa_ca{ii} = kappa_ch(ij);
        epsilon_ca{ii} = epsilon_ch(ij);
        delta_ca{ii} = delta_ch(ij);
        Phib(ii) = Phi_ch(end);
        rhotor_ca{ii} = sqrt(Phi_ch(ij)/Phi_ch(end));    
    end
end
if nbt>1
    % unified radial grid for input data
    [rhotor,iR,iR2,dpsi,dpsi2,dpsi2oR2,dpsidV,F,epsilon,Rgeom,kappa,delta] = ...
        unified_radial_grid(rhotor_ca,iR_ca,iR2_ca,dpsi_ca,dpsi2_ca,dpsi2oR2_ca,dpsidV_ca,F_ca,epsilon_ca,Rgeom_ca,kappa_ca,delta_ca);
end
return
