classdef espic2dhdf5
    %UNTITLED Summary of this class goes here
    %   Detailed explanation goes here
    
    properties
        filename
        timestamp
        info
        t0d
        t1d
        t2d
        tpart
        it0
        it1
        it2
        
        %% Physical constants
        vlight=299792458;
        qe=1.60217662E-19;
        me=9.109383E-31;
        eps_0=8.85418781762E-12;
        kb=1.38064852E-23;
        
        %% Run parameters
        dt
        nrun
        nlres
        nlsave
        nlclassical
        nlPhis
        nz
        nnr
        lz
        qsim
        msim
        nplasma
        potinn
        potout
        B0
        Rcurv
        width
        n0
        temp
        femorder
        ngauss
        plasmadim
        radii
        H0
        P0
        
        
        %% Frequencies
        omepe
        omece
        
        %% Normalizations
        tnorm
        rnorm
        bnorm
        enorm
        phinorm
        vnorm
        
        %% Grid data
        rgrid
        zgrid
        dz
        dr
        Vcell
        
        %% Magnetic field
        Br
        Bz
        Athet
        B
        
        %% Energies
        epot
        ekin
        etot
        etot0
        eerr
        
        %% 2D time data
        N
        fluidUR
        fluidUZ
        fluidUTHET
        pot
        Er
        Ez
        CellVol
        Presstens
        
        %% Splines
        knotsr
        knotsz
        
        %% Particle parameters
        nbparts
        H
        P
        partepot
        R
        Z
        Rindex
        Zindex
        partindex
        VR
        VZ
        VTHET
        THET
        
    end
    
    methods
        function file=file(obj)
            file=obj.filename;
        end
        
        function obj = espic2dhdf5(filename,readparts)
            filedata=dir(filename);
            obj.filename=filename;
            if (isempty(filedata))
                error("File: ""%s"" doesn't exist",filename)
            end
            obj.timestamp=filedata.date;
            if nargin==1
                readparts=true;
            end
            %obj.info=h5info(filename);
            
            obj.dt = h5readatt(obj.filename,'/data/input.00/','dt');
            obj.nrun = h5readatt(obj.filename,'/data/input.00/','nrun');
            obj.nlres = strcmp(h5readatt(obj.filename,'/data/input.00/','nlres'),'y');
            obj.nlsave = strcmp(h5readatt(obj.filename,'/data/input.00/','nlsave'),'y');
            obj.nlclassical =strcmp(h5readatt(obj.filename,'/data/input.00/','nlclassical'),'y');
            obj.nlPhis =strcmp(h5readatt(obj.filename,'/data/input.00/','nlPhis'),'y');
            obj.nz = h5readatt(obj.filename,'/data/input.00/','nz');
            obj.nnr = h5read(obj.filename,'/data/input.00/nnr');
            obj.lz = h5read(obj.filename,'/data/input.00/lz');
            obj.qsim = h5readatt(obj.filename,'/data/input.00/','qsim');
            obj.msim = h5readatt(obj.filename,'/data/input.00/','msim');
            obj.nplasma = h5readatt(obj.filename,'/data/input.00/','nplasma');
            obj.potinn = h5readatt(obj.filename,'/data/input.00/','potinn');
            obj.potout = h5readatt(obj.filename,'/data/input.00/','potout');
            obj.B0 = h5readatt(obj.filename,'/data/input.00/','B0');
            obj.Rcurv = h5readatt(obj.filename,'/data/input.00/','Rcurv');
            obj.width = h5readatt(obj.filename,'/data/input.00/','width');
            obj.n0 = h5readatt(obj.filename,'/data/input.00/','n0');
            obj.temp = h5readatt(obj.filename,'/data/input.00/','temp');
            try
                obj.it0 = h5readatt(obj.filename,'/data/input.00/','it0d');
                obj.it1 = h5readatt(obj.filename,'/data/input.00/','it2d');
                obj.it2 = h5readatt(obj.filename,'/data/input.00/','itparts');
            catch
                obj.it0 = h5readatt(obj.filename,'/data/input.00/','it0');
                obj.it1 = h5readatt(obj.filename,'/data/input.00/','it1');
                obj.it1 = h5readatt(obj.filename,'/data/input.00/','it2');
            end
            
            obj.omepe=sqrt(abs(obj.n0)*obj.qe^2/(obj.me*obj.eps_0));
            obj.omece=obj.qe*obj.B0/obj.me;
            
            obj.nbparts= h5read(obj.filename, '/data/var0d/nbparts');
            try
                obj.H0 = h5read(obj.filename,'/data/input.00/H0');
                obj.P0 = h5read(obj.filename,'/data/input.00/P0');
            catch
                obj.H0=3.2e-14;
                obj.P0=8.66e-25;
            end
            
            % Normalizations
            obj.tnorm=1/obj.omepe;
            obj.rnorm=obj.vlight*obj.tnorm;
            obj.bnorm=obj.B0;
            obj.enorm=obj.vlight*obj.bnorm;
            obj.phinorm=obj.enorm*obj.rnorm;
            obj.vnorm=obj.vlight;
            
            % Grid data
            obj.rgrid= h5read(obj.filename, '/data/var1d/rgrid')*obj.rnorm;
            obj.zgrid= h5read(obj.filename, '/data/var1d/zgrid')*obj.rnorm;
            obj.dz=(obj.zgrid(end)-obj.zgrid(1))/double(obj.nz);
            rid=1;
            for i=1:length(obj.nnr)
                obj.dr(i)=(obj.rgrid(sum(obj.nnr(1:i))+1)-obj.rgrid(rid))/double(obj.nnr(i));
                rid=rid+obj.nnr(i);
            end
            obj.Vcell=(obj.zgrid(2:end)-obj.zgrid(1:end-1))*((obj.rgrid(2:end).^2-obj.rgrid(1:end-1).^2)*pi)';
            
            
            Br = h5read(obj.filename,'/data/fields/Br')*obj.bnorm;
            obj.Br= reshape(Br,length(obj.zgrid),length(obj.rgrid));
            Bz = h5read(obj.filename,'/data/fields/Bz')*obj.bnorm;
            obj.Bz= reshape(Bz,length(obj.zgrid),length(obj.rgrid));
            try
                Atheta = h5read(obj.filename,'/data/fields/Athet')*obj.bnorm;
                obj.Athet= reshape(Atheta,length(obj.zgrid),length(obj.rgrid));
            catch
            end
            obj.B=sqrt(obj.Bz.^2+obj.Br.^2);
            clear Br Bz
            try
                obj.t0d=h5read(obj.filename,'/data/var0d/time');
            catch
                obj.t0d=obj.dt.*double(0:length(obj.epot)-1);
            end
            
            
            obj.femorder = h5read(obj.filename,'/data/input.00/femorder');
            obj.ngauss = h5read(obj.filename,'/data/input.00/ngauss');
            obj.plasmadim = h5read(obj.filename,'/data/input.00/plasmadim');
            obj.radii = h5read(obj.filename,'/data/input.00/radii');
            
            obj.epot = h5read(obj.filename,'/data/var0d/epot');
            obj.ekin = h5read(obj.filename,'/data/var0d/ekin');
            obj.etot = h5read(obj.filename,'/data/var0d/etot');
            try
                obj.etot0 = h5read(filename,'/data/var0d/etot0');
                obj.eerr = obj.etot-obj.etot0;
            catch
                obj.eerr = obj.etot-obj.etot(2);
            end
            
            obj.pot=gridquantity(obj.filename,'/data/fields/pot',sum(obj.nnr)+1, obj.nz+1,obj.phinorm);
            obj.Er=gridquantity(obj.filename,'/data/fields/Er',sum(obj.nnr)+1, obj.nz+1,obj.enorm);
            obj.Ez=gridquantity(obj.filename,'/data/fields/Ez',sum(obj.nnr)+1, obj.nz+1,obj.enorm);
            
            try
                obj.t2d = h5read(obj.filename,'/data/fields/time');
            catch
                info=h5info(obj.filename,'/data/fields/partdensity');
                obj.t2d=obj.dt*(0:info.objspace.Size(2)-1)*double(obj.it1);
            end
            
            try
                info=h5info(obj.filename,'/data/fields/moments');
                obj.femorder = h5read(obj.filename,'/data/input.00/femorder');
                kr=obj.femorder(2)+1;
                obj.knotsr=augknt(obj.rgrid,kr);
                
                kz=obj.femorder(1)+1;
                obj.knotsz=augknt(obj.zgrid,kz);
                try
                    obj.CellVol= reshape(h5read(obj.filename,'/data/fields/volume'),length(obj.knotsz)-kz,length(obj.knotsr)-kr);
                    obj.CellVol=permute(obj.CellVol,[2,1,3])*obj.rnorm^3;
                catch
                    zvol=fnder(spmak(obj.knotsz,ones(1,length(obj.knotsz)-kz)), -1 );
                    rvol=fnder(spmak(obj.knotsr,2*pi*[obj.rgrid' 2*obj.rgrid(end)-obj.rgrid(end-1)]), -1 );
                    ZVol=diff(fnval(zvol,obj.knotsz));
                    RVol=diff(fnval(rvol,obj.knotsr));
                    obj.CellVol=RVol(3:end-1)*ZVol(3:end-1)';
                    obj.CellVol=padarray(obj.CellVol,[1,1],'replicate','post');
                end
                
                obj.N=splinedensity(obj.filename, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.CellVol, abs(obj.qsim/obj.qe), 1);
                obj.fluidUR=splinevelocity(obj.filename, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.vnorm, 2);
                obj.fluidUTHET=splinevelocity(obj.filename, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.vnorm, 3);
                obj.fluidUZ=splinevelocity(obj.filename, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.vnorm, 4);
                obj.Presstens=splinepressure(obj.filename, '/data/fields/moments', obj.knotsr, obj.knotsz, obj.femorder, obj.CellVol, obj.vnorm^2*obj.msim);
            catch
                obj.CellVol=(obj.zgrid(2:end)-obj.zgrid(1:end-1))*((obj.rgrid(2:end).^2-obj.rgrid(1:end-1).^2)*pi)';
                obj.CellVol=obj.CellVol';
                obj.N=griddensity(obj.filename, '/data/fields/partdensity', sum(obj.nnr)+1, obj.nz+1, obj.CellVol, abs(obj.qsim/obj.qe), true);
                obj.fluidUR=gridquantity(obj.filename, '/data/fields/fluidur', sum(obj.nnr)+1, obj.nz+1, obj.vnorm, true);
                obj.fluidUTHET=gridquantity(obj.filename, '/data/fields/fluiduthet', sum(obj.nnr)+1, obj.nz+1, obj.vnorm, true);
                obj.fluidUZ=gridquantity(obj.filename, '/data/fields/fluiduz', sum(obj.nnr)+1, obj.nz+1, obj.vnorm, true);
            end
            
            if(readparts)
            
            obj.R = h5read(filename,'/data/part/R')*obj.rnorm;
            obj.Z = h5read(filename,'/data/part/Z')*obj.rnorm;
            try
                obj.THET = h5read(filename,'/data/part/THET');
            catch
                clear obj.THET
            end
            try
                obj.Rindex=h5read(filename,'/data/part/Rindex');
                obj.Zindex=h5read(filename,'/data/part/Zindex');
            catch
                clear obj.Rindex obj.Zindex
            end
            UR = h5read(filename,'/data/part/UR');
            UZ = h5read(filename,'/data/part/UZ');
            UTHET= h5read(filename,'/data/part/UTHET');
            if(obj.nlclassical)
                obj.VR = UR*obj.vnorm;
                clear UR
                obj.VZ = UZ*obj.vnorm;
                clear UZ
                obj.VTHET = UTHET*obj.vnorm;
                clear UTHET
            else
                GAMM = sqrt(1+UR.^2+UZ.^2+UTHET.^2);
                obj.VR = UR./GAMM*obj.vnorm;
                clear UR
                obj.VZ = UZ./GAMM*obj.vnorm;
                clear UZ
                obj.VTHET = UTHET./GAMM*obj.vnorm;
                clear UTHET
                clear GAMM
            end
            %obj.VPERP = sqrt(obj.VR.*obj.VR+obj.VTHET.*obj.VTHET);
            obj.partepot = sign(obj.qsim)*h5read(filename,'/data/part/pot')*obj.phinorm*obj.qe;
            obj.H=0.5*obj.me*(obj.VR.^2+obj.VTHET.^2+obj.VZ.^2)+obj.partepot;
            %clear partepot
            obj.P=obj.R.*(obj.VTHET*obj.me+sign(obj.qsim)*obj.qe*Athet(obj.R,obj.Z));
            try
                obj.partindex = h5read(filename,'/data/part/partindex');
                partindex=obj.partindex;
                [~,Indices]=sort(partindex,'descend');
                for i=1:size(Indices,2)
                    obj.H(:,i)=obj.H(Indices(:,i),i);
                    obj.P(:,i)=obj.P(Indices(:,i),i);
                    obj.R(:,i)=obj.R(Indices(:,i),i);
                    obj.Z(:,i)=obj.Z(Indices(:,i),i);
                    obj.Rindex(:,i)=obj.Rindex(Indices(:,i),i);
                    obj.Zindex(:,i)=obj.Zindex(Indices(:,i),i);
                    obj.VR(:,i)=obj.VR(Indices(:,i),i);
                    obj.VZ(:,i)=obj.VZ(Indices(:,i),i);
                    obj.VTHET(:,i)=obj.VTHET(Indices(:,i),i);
                    obj.THET(:,i)=obj.THET(Indices(:,i),i);
                end
            catch
            end
            clear partindex;
            end 
            
            try
                obj.tpart = h5read(filename,'/data/part/time');
            catch
                obj.tpart=obj.dt*(0:size(obj.R,2)-1)*double(obj.it2);
            end
            
            function Atheta=Athet(R,Z)
                Atheta=0.5*obj.B0*(R-obj.width/pi*(obj.Rcurv-1)/(obj.Rcurv+1)...
                .*besseli(1,2*pi*R/obj.width).*cos(2*pi*Z/obj.width));
            end
        end
        
        function changed=ischanged(obj)
            %ischanged Check if the file has been changed and some data
            %must be reloaded
            try
                filedata=dir(obj.filename);
                checkedtimestamp=filedata.date;
                if (max(checkedtimestamp > obj.timestamp) )
                    changed=true;
                    return
                end
                changed=false;
                return
            catch
                changed=true;
                return
            end
        end
        
        function plotEnergy(obj)
            tmin=2;
            tmax=length(obj.ekin);
            [~, name, ~] = fileparts(obj.file);
            f=figure('Name', sprintf('%s Energy',name));
            subplot(2,1,1)
            plot(obj.t0d(tmin:tmax),obj.ekin(tmin:tmax),'o-',...
                 obj.t0d(tmin:tmax),obj.epot(tmin:tmax),'d-',...
                 obj.t0d(tmin:tmax),obj.etot(tmin:tmax),'h-',...
                 obj.t0d(tmin:tmax),obj.eerr(tmin:tmax),'x--')
            legend('ekin', 'epot', 'etot','eerr')
            xlabel('Time [s]')
            ylabel('Energies [J]')
            grid on
            xlimits=xlim();
            
            subplot(2,1,2)
            try
                semilogy(obj.t0d(tmin:tmax),abs(obj.eerr(tmin:tmax)./obj.etot0(tmin:tmax)),'h-')
            catch
                semilogy(obj.t0d(tmin:tmax),abs(obj.eerr(tmin:tmax)/obj.etot(2)),'h-')
            end
            xlabel('t [s]')
            ylabel('E_{err}/E_{tot}')
            xlim(xlimits)
            grid on
            try
                yyaxis right
                plot(obj.t0d(tmin:tmax),abs(obj.nbparts(tmin:tmax)./obj.nbparts(1)*100),'d--')
                ylabel('Nparts %')
                %ylim([0 110])
            catch
            end
            obj.savegraph(f,sprintf('%sEnergy',name));
        end
        
        function savegraph(obj, figure, name, papsize)
            figure.PaperUnits='centimeters';
            if (nargin < 4)
                papsize=[14 16];
            end
            figure.PaperSize=papsize;
            print(figure,name,'-dpdf','-fillpage')
            savefig(figure,name)
        end
    end
end