function [astra_profiles,chiegauss,pechgauss,jecgauss,pnbe,kappa,B0,te_BC] = get_ASTRA_profiles(config,model)
% function [astra_profiles,chiegauss,pechgauss,jecgauss,pnbe,kappa,B0,te_BC] = get_ASTRA_profiles(config,model)
% Get ASTRA parameters: astra_profiles contains profiles for ne, ni, ti, ze,
% chiegauss: electron heat conductivity on RAPTOR gauss grid,
% pechgauss: ECH power density RAPTOR gauss grid,
% jecgauss: EC CD density RAPTOR gauss grid,
% pnbe: NBI power density from ASTRA data.
% kappa: plasma boundary elongation
% B0: magnetic field
% te_BC: boundary conditions for Te.
%
% Load in data from ASTRA file
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
astrafile_fullpath = fullfile(astradatapath,astrafile{1});
% contains ASTRA data for equilibrium
  
astra = load(astrafile_fullpath);
if isfield(astra,'TCV_astra')
    astra = astra.TCV_astra;
end

if strfind(config.equi.filenames{1},'RES')  
    data.rho = astra.data.rho;
    data.ne = astra.data.ne;
    data.ni = astra.data.ni;
    data.ti = astra.data.ti;
    data.zef = astra.data.zeffr;
    %data.zef = repmat(astra.zef',size(data.rho,1),1);
    data.he = astra.data.he;
    data.peecr = astra.data.peecr;
    data.cdr = astra.data.cdr;
    data.pnb = astra.data.pbeam;
    data.ipl = astra.data.ipl';   
    data.t = astra.data.t';
    data.te = astra.data.te;
    data.qe = astra.data.qe;
    data.g11 = astra.data.g11;
    data.elon = astra.data.elon;
    data.btor = astra.data.btor';    
else
    data.rho = astra.out.RHO;
    data.ne = astra.out.NE;
    data.ni = astra.out.NI;
    data.ti = astra.out.TI;
    data.zef = astra.out.ZEF;
    data.he = astra.out.HE;
    data.peecr = astra.out.PEECR;
    data.cdr = astra.out.CD;
    data.pnb = astra.out.PBEAM;
    data.ipl = astra.out.IPL;
    data.t = astra.out.T;
    data.te = astra.out.TE;
    data.qe = astra.out.QE;
    data.g11 = astra.out.G11;
    data.elon = astra.out.ELON;
    data.btor = astra.out.BTOR;
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

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

% Get plasma parameters
ne_astra = 1e19*data.ne(:,ik);
ni_astra = 1e19*data.ni(:,ik);
ti_astra = 1e03*data.ti(:,ik);
ze_astra = data.zef(:,ik);
chie_astra = data.he(:,ik);
pech_astra = 1e06*data.peecr(:,ik);
jec_astra = 1e06*data.cdr(:,ik);
pnb_astra = 1e06*data.pnb(:,ik);
Ip_astra = 1e06*data.ipl(:,ik);
te_BC_astra = 1e03*data.te(end,ik);
kappa_astra = data.elon(end,ik);
B0_astra = data.btor(:,ik);

% normalized input radial grid (normalized toroidal rho)
rho = bsxfun(@rdivide,data.rho,data.rho(end,:));
% unified radial grid for input data
[rhou,neu,niu,tiu,zeu,chieu,pechu,jecu,pnbu] = unified_radial_grid(rho,ne_astra,ni_astra,ti_astra,ze_astra,...
    chie_astra,pech_astra,jec_astra,pnb_astra);

% interpolate to RAPTOR time grid if necessary  
if isempty(config.equi.times)
    [rhoin,astra_profiles.ne,astra_profiles.ni,astra_profiles.ti,astra_profiles.ze,chie,...
        pech,jec,pnb,astra_profiles.Ip,te_BC,kappa,B0] = interp2RAPTOR_tgrid(data.t(ik),config.grid.tgrid,...
        rhou,neu,niu,tiu,zeu,chieu,pechu,jecu,pnbu,Ip_astra,te_BC_astra,kappa_astra,B0_astra);
    astra_profiles.rho = rhoin(:,end);
else
    astra_profiles.rho = rhou;
    astra_profiles.ne = neu;
    astra_profiles.ni = niu;
    astra_profiles.ti = tiu;
    astra_profiles.ze = zeu;
    astra_profiles.Ip = Ip_astra;
    chie = chieu;
    pech = pechu;
    jec = jecu;
    pnb = pnbu;
    te_BC = te_BC_astra;
    kappa = kappa_astra;
    B0 = B0_astra;
end

% Recalculate chie and ECH power density pec on RAPTOR radial grid.
nbt = size(chie,2);
for ii=1:nbt
    chiegauss(:,ii) = spline(astra_profiles.rho,chie(:,ii),model.rgrid.rhogauss);
    pechgauss(:,ii) = spline(astra_profiles.rho,pech(:,ii),model.rgrid.rhogauss);
    jecgauss(:,ii) = spline(astra_profiles.rho,jec(:,ii),model.rgrid.rhogauss);
    pnbe(:,ii) = spline(astra_profiles.rho,pnb(:,ii),model.rgrid.rho);
end

return