function Data = load_espic2d( file, Data, type, partstime)
% Old function used to load simulations data from espic2d HDF5 results
% This function does not handle b-splines data apart from the density and load the full file to memory
% This can be problematic for big files
if(nargin==1)
    type='all';
else
    filedata=dir(file);
    if (isempty(filedata))
        error("File: ""%s"" doesn't exist",file)
    end
    timestamp=filedata.date;
    try
        if (strcmp(Data.file,file) && ~max(timestamp > Data.timestamp) && (strcmp(Data.type,type)||strcmp(Data.type,'all')) )
            return
        end
    catch
    end
    if(nargin<3)
        type='all';
    end
    if(nargin<4)
        partstime=[0 1];
    end
end
Data.type=type;
Data.file=file;
% Physical constants
Data.vlight=299792458;
Data.qe=1.60217662E-19;
Data.me=9.109383E-31;
Data.eps_0=8.85418781762E-12;
Data.kb=1.38064852E-23;

%Data.job_time = h5readatt(file,'/data/input.00/','job_time');
%Data.extra_time = h5readatt(file,'/data/input.00/','extra_time');
Data.dt = h5readatt(file,'/data/input.00/','dt');
Data.tmax= h5readatt(file,'/data/input.00/','tmax');
Data.nrun = h5readatt(file,'/data/input.00/','nrun');
Data.nlres = strcmp(h5readatt(file,'/data/input.00/','nlres'),'y');
Data.nlsave = strcmp(h5readatt(file,'/data/input.00/','nlsave'),'y');
Data.nlclassical =strcmp(h5readatt(file,'/data/input.00/','nlclassical'),'y');
Data.nlPhis =strcmp(h5readatt(file,'/data/input.00/','nlPhis'),'y');
Data.nz = h5readatt(file,'/data/input.00/','nz');
Data.nnr = h5read(file,'/data/input.00/nnr');
Data.femorder = h5read(file,'/data/input.00/femorder');
Data.lz = h5read(file,'/data/input.00/lz');
Data.lr = h5read(file,'/data/input.00/radii');
Data.qsim = h5readatt(file,'/data/input.00/','qsim');
Data.msim = h5readatt(file,'/data/input.00/','msim');
Data.nplasma = h5readatt(file,'/data/input.00/','nplasma');
Data.potinn = h5readatt(file,'/data/input.00/','potinn');
Data.potout = h5readatt(file,'/data/input.00/','potout');
Data.B0 = h5readatt(file,'/data/input.00/','B0');
Data.Rcurv = h5readatt(file,'/data/input.00/','Rcurv');
Data.width = h5readatt(file,'/data/input.00/','width');
Data.n0 = h5readatt(file,'/data/input.00/','n0');
Data.temp = h5readatt(file,'/data/input.00/','temp');
try
    Data.it0 = h5readatt(file,'/data/input.00/','it0d');
    Data.it1 = h5readatt(file,'/data/input.00/','it2d');
    Data.it2 = h5readatt(file,'/data/input.00/','itparts');
catch
    Data.it0 = h5readatt(file,'/data/input.00/','it0');
    Data.it1 = h5readatt(file,'/data/input.00/','it1');
    Data.it1 = h5readatt(file,'/data/input.00/','it2');
end

Data.omepe=sqrt(abs(Data.n0)*Data.qe^2/(Data.me*Data.eps_0));
Data.omece=Data.qe*Data.B0/Data.me;


% Normalizations
Data.tnorm=1/Data.omepe;
Data.rnorm=Data.vlight*Data.tnorm;
Data.bnorm=Data.B0;
Data.enorm=Data.vlight*Data.bnorm;
Data.phinorm=Data.enorm*Data.rnorm;
Data.vnorm=Data.vlight;

% Grid data
Data.rgrid= h5read(file, '/data/var1d/rgrid')*Data.rnorm;
Data.zgrid= h5read(file, '/data/var1d/zgrid')*Data.rnorm;
Data.dz=(Data.zgrid(end)-Data.zgrid(1))/double(Data.nz);
Data.dr=(Data.lr(2:end)-Data.lr(1:end-1))./double(Data.nnr);


Br = h5read(file,'/data/fields/Br')*Data.bnorm;
Data.Br= reshape(Br,length(Data.zgrid),length(Data.rgrid));
Bz = h5read(file,'/data/fields/Bz')*Data.bnorm;
Data.Bz= reshape(Bz,length(Data.zgrid),length(Data.rgrid));
try
    Atheta = h5read(file,'/data/fields/Athet')*Data.bnorm;
    Data.Athet= reshape(Atheta,length(Data.zgrid),length(Data.rgrid));
catch
end
Data.B=sqrt(Data.Bz.^2+Data.Br.^2);
clear Br Bz


Data.femorder = h5read(file,'/data/input.00/femorder');
Data.ngauss = h5read(file,'/data/input.00/ngauss');
Data.plasmadim = h5read(file,'/data/input.00/plasmadim');
Data.radii = h5read(file,'/data/input.00/radii');

Data.epot = h5read(file,'/data/var0d/epot');
Data.ekin = h5read(file,'/data/var0d/ekin');
Data.etot = h5read(file,'/data/var0d/etot');
try
    Data.etot0 = h5read(file,'/data/var0d/etot0');
    Data.eerr = Data.etot-Data.etot0;
catch
    Data.eerr = Data.etot-Data.etot(2);
end

if(strcmp(type,'all') || strcmp(type,'fields'))
    pot = h5read(file,'/data/fields/pot')*Data.phinorm;
    Data.pot= reshape(pot,length(Data.zgrid),length(Data.rgrid),size(pot,2));
    clear pot
    Data.pot=permute(Data.pot,[2,1,3]);
    Er = h5read(file,'/data/fields/Er')*Data.enorm;
    Data.Er= reshape(Er,length(Data.zgrid),length(Data.rgrid),size(Er,2));
    clear Er
    Data.Er=permute(Data.Er,[2,1,3]);
    Ez = h5read(file,'/data/fields/Ez')*Data.enorm;
    Data.Ez= reshape(Ez,length(Data.zgrid),length(Data.rgrid),size(Ez,2));
    clear Ez
    Data.Ez=permute(Data.Ez,[2,1,3]);
    try
        dens = h5read(file,'/data/fields/moments',[1 1 1],[1,Inf,Inf]);
        kr=Data.femorder(2)+1;
        Data.knotsr=augknt(Data.rgrid,kr);
        kz=Data.femorder(1)+1;
        Data.knotsz=augknt(Data.zgrid,kz);
        
        zvol=fnder(spmak(Data.knotsz,ones(1,length(Data.zgrid))), -1 );
        rvol=fnder(spmak(Data.knotsr,2*pi*Data.rgrid'), -1 );
        ZVol=diff(fnval(zvol,Data.knotsz));
        RVol=diff(fnval(rvol,Data.knotsr));
        Data.V=ZVol(2:end-1)*RVol(2:end-1)';
        Data.invV=1./Data.V';
        Data.invV(isinf(Data.invV))=0;
        
        dens = h5read(file,'/data/fields/moments',[1 1 1],[1,Inf,Inf]);
        Data.dens=reshape(dens,Data.nz+Data.femorder(1),sum(Data.nnr)+Data.femorder(2),[]);
        Data.dens=permute(Data.dens,[2,1,3]);
        clear dens
        
        fluidvr = h5read(file,'/data/fields/moments',[2 1 1],[1,Inf,Inf]);
        Data.fluidvr=reshape(fluidvr,Data.nz+Data.femorder(1),sum(Data.nnr)+Data.femorder(2),[]);
        clear fluidvr
        Data.fluidvr=permute(Data.fluidvr,[2,1,3])*Data.vnorm./Data.dens;
        Data.fluidvr(isinf(Data.fluidvr))=0;
        
        fluidvthet = h5read(file,'/data/fields/moments',[3 1 1],[1,Inf,Inf]);
        Data.fluidvthet=reshape(fluidvthet,Data.nz+Data.femorder(1),sum(Data.nnr)+Data.femorder(2),[]);
        Data.fluidvthet=permute(Data.fluidvthet,[2,1,3])*Data.vnorm./Data.dens;
        Data.fluidvthet(isinf(Data.fluidvthet))=0;
        clear fluidvthet
        
        fluidvz = h5read(file,'/data/fields/moments',[4 1 1],[1,Inf,Inf]);
        Data.fluidvz=reshape(fluidvz,Data.nz+Data.femorder(1),sum(Data.nnr)+Data.femorder(2),[]);
        Data.fluidvz=permute(Data.fluidvz,[2,1,3])*Data.vnorm./Data.dens;
        Data.fluidvz(isinf(Data.fluivz))=0;
        clear fluidvz
        
        vrvr = h5read(file,'/data/fields/moments',[5 1 1],[1,Inf,Inf]);
        Data.fluidvz=reshape(vrvr,Data.nz+Data.femorder(1),sum(Data.nnr)+Data.femorder(2),[]);
        Data.fluidvz=permute(Data.fluidvz,[2,1,3]);
        
        %Data.Presstens=zeroes(6,size(Data.zgrid,1),size(Data.rgrid,1), size(moments,3));
        %Data.Presstens(1,:,:,:)=reshape(moments(5,:,:),Data.nz+Data.femorder(1),sum(Data.nnr)+Data.femorder(2),[])
        %Data.Presstens=Data.me*Data.vlight^2*reshape(h5read(file,'/data/fields/presstens'),6,length(Data.zgrid)-1,length(Data.rgrid)-1,[]);
        %Data.Presstens=permute(Data.Presstens,[1,3,2,4]);
        clear moments;
        Data.N=zeros(length(Data.rgrid),length(Data.zgrid),size(Data.dens,3));
        for i=1:size(Data.dens,3)
            Data.dens(:,:,i)=Data.dens(:,:,i).*Data.invV*abs(Data.qsim/Data.qe);
            Data.N(:,:,i)=fnval(spmak({Data.knotsr,Data.knotsz},Data.dens(:,:,i)),{Data.rgrid,Data.zgrid});
        end
        
        %Data.moments=reshape(Data.moments,10,length(Data.zgrid),length(Data.rgrid),[]);
        %Data.moments=permute(Data.moments,[1,3,2,4]);
    catch
        N = h5read(file,'/data/fields/partdensity');
        Data.N=reshape(N,length(Data.zgrid)-1,length(Data.rgrid)-1,[]);
        Data.N=permute(Data.N,[2,1,3]);
        Data.Npartscell=Data.N;
        Vcell=(Data.zgrid(2:end)-Data.zgrid(1:end-1))*((Data.rgrid(2:end).^2-Data.rgrid(1:end-1).^2)*pi)';
        for i=1:size(Data.N,3)
            Data.N(:,:,i)=abs(Data.qsim/Data.qe)*Data.N(:,:,i)./Vcell';
        end
        Data.N=padarray(Data.N,[1,1],0,'post');
        Data.Npartscell=padarray(Data.Npartscell,[1,1],0,'post');
        clear Vcell
        clear N
        Data.fluidUR=reshape(h5read(file,'/data/fields/fluidur'),length(Data.zgrid)-1,length(Data.rgrid)-1,[]);
        Data.fluidUR=permute(Data.fluidUR,[2,1,3])*Data.vnorm;
        Data.fluidUR=padarray(Data.fluidUR,[1,1],0,'post');
        Data.fluidUZ=reshape(h5read(file,'/data/fields/fluiduz'),length(Data.zgrid)-1,length(Data.rgrid)-1,[]);
        Data.fluidUZ=permute(Data.fluidUZ,[2,1,3])*Data.vnorm;
        Data.fluidUZ=padarray(Data.fluidUZ,[1,1],0,'post');
        Data.fluidUTHET=reshape(h5read(file,'/data/fields/fluiduthet'),length(Data.zgrid)-1,length(Data.rgrid)-1,[]);
        Data.fluidUTHET=permute(Data.fluidUTHET,[2,1,3])*Data.vnorm;
        Data.fluidUTHET=padarray(Data.fluidUTHET,[1,1],0,'post');
        
        Data.Presstens=Data.me*Data.vlight^2*reshape(h5read(file,'/data/fields/presstens'),6,length(Data.zgrid)-1,length(Data.rgrid)-1,[]);
        Data.Presstens=permute(Data.Presstens,[1,3,2,4]);
        Data.Presstens=padarray(Data.Presstens,[0,1,1,0],0,'post');
        for i=1:6
            Data.Presstens(i,:,:,:)=squeeze(Data.Presstens(i,:,:,:)).*Data.N;
        end
    end
    try
        Data.t2d = h5read(file,'/data/fields/time');
    catch
        Data.t2d=Data.dt*(0:size(Data.N,3)-1)*double(Data.it1);
    end
    if(strcmp(Data.type,'parts'))
        Data.type='all';
    end
end

if(strcmp(type,'all') || strcmp(type,'parts'))
    Data.R = h5read(file,'/data/part/R')*Data.rnorm;
    partsbegin=floor(partstime(1)*size(Data.R,2))+1;
    partsend=floor(partstime(2)*size(Data.R,2));
    Data.R=Data.R(:,partsbegin:partsend);
    Data.Z = h5read(file,'/data/part/Z',[1,partsbegin],[size(Data.R,1) partsend])*Data.rnorm;
    try
        Data.THET = h5read(file,'/data/part/THET',[1,partsbegin],[size(Data.R,1) partsend]);
    catch
        clear Data.THET
    end
    try
        Data.Rindex=h5read(file,'/data/part/Rindex',[1,partsbegin],[size(Data.R,1) partsend]);
        Data.Zindex=h5read(file,'/data/part/Zindex',[1,partsbegin],[size(Data.R,1) partsend]);
    catch
        clear Data.Rindex Data.Zindex
    end
    UR = h5read(file,'/data/part/UR',[1,partsbegin],[size(Data.R,1) partsend]);
    UZ = h5read(file,'/data/part/UZ',[1,partsbegin],[size(Data.R,1) partsend]);
    UTHET= h5read(file,'/data/part/UTHET',[1,partsbegin],[size(Data.R,1) partsend]);
    if(Data.nlclassical)
        Data.VR = UR*Data.vnorm;
        clear UR
        Data.VZ = UZ*Data.vnorm;
        clear UZ
        Data.VTHET = UTHET*Data.vnorm;
        clear UTHET
    else
        GAMM = sqrt(1+UR.^2+UZ.^2+UTHET.^2);
        Data.VR = UR./GAMM*Data.vnorm;
        clear UR
        Data.VZ = UZ./GAMM*Data.vnorm;
        clear UZ
        Data.VTHET = UTHET./GAMM*Data.vnorm;
        clear UTHET
        clear GAMM
    end
    %Data.VPERP = sqrt(Data.VR.*Data.VR+Data.VTHET.*Data.VTHET);
    partepot = sign(Data.qsim)*h5read(file,'/data/part/pot',[1,partsbegin],[size(Data.R,1) partsend])*Data.phinorm*Data.qe;
    Data.H=0.5*Data.me*(Data.VR.^2+Data.VTHET.^2+Data.VZ.^2)+partepot;
    clear partepot
    Data.P=Data.R.*(Data.VTHET*Data.me+sign(Data.qsim)*Data.qe*Athet(Data.R,Data.Z));
    try
        partindex = h5read(file,'/data/part/partindex',[1,partsbegin],[size(Data.R,1) partsend]);
        [~,Indices]=sort(partindex);
        for i=1:size(Indices,2)
            Data.H(:,i)=Data.H(Indices(:,i),i);
            Data.P(:,i)=Data.P(Indices(:,i),i);
            Data.R(:,i)=Data.R(Indices(:,i),i);
            Data.Z(:,i)=Data.Z(Indices(:,i),i);
            Data.Rindex(:,i)=Data.Rindex(Indices(:,i),i);
            Data.Zindex(:,i)=Data.Zindex(Indices(:,i),i);
            Data.VR(:,i)=Data.VR(Indices(:,i),i);
            Data.VZ(:,i)=Data.VZ(Indices(:,i),i);
            Data.VTHET(:,i)=Data.VTHET(Indices(:,i),i);
            Data.THET(:,i)=Data.THET(Indices(:,i),i);
        end
    catch
    end
    clear partindex;
    
    try
        Data.tpart = h5read(file,'/data/part/time');
    catch
        Data.tpart=Data.dt*(0:size(Data.R,2)-1)*double(Data.it2);
    end
    if(strcmp(Data.type,'fields'))
        Data.type='all';
    end
end

Data.t0d=Data.dt.*double(0:length(Data.epot)-1);


filedata=dir(file);
Data.timestamp=filedata.date;
disp('M reloaded')

    function Atheta=Athet(R,Z)
        Atheta=0.5*Data.B0*(R-Data.width/pi*(Data.Rcurv-1)/(Data.Rcurv+1)...
            .*besseli(1,2*pi*R/Data.width).*cos(2*pi*Z/Data.width));
    end

end