From 096bab2609ecbc9b73d93bad1098bf44a952157a Mon Sep 17 00:00:00 2001
From: VinithSamson11 <109979733+VinithSamson11@users.noreply.github.com>
Date: Sat, 6 Jan 2024 07:17:54 +0530
Subject: [PATCH 01/10] AMsmileih5.py version1 upload

---
 postpic/datareader/AMsmileih5.py | 403 +++++++++++++++++++++++++++++++
 1 file changed, 403 insertions(+)
 create mode 100644 postpic/datareader/AMsmileih5.py

diff --git a/postpic/datareader/AMsmileih5.py b/postpic/datareader/AMsmileih5.py
new file mode 100644
index 0000000..9128d02
--- /dev/null
+++ b/postpic/datareader/AMsmileih5.py
@@ -0,0 +1,403 @@
+from __future__ import absolute_import, division, print_function, unicode_literals
+
+from . import Dumpreader_ifc
+from . import Simulationreader_ifc
+import re
+import os.path
+import h5py
+import glob
+import numpy as np
+import helper 
+from helper_fft import fft
+
+__all__ = ['AMSmileiReader', 'SmileiSeries']
+
+
+def _generateh5indexfile(indexfile, fnames):
+    '''
+    Creates a h5 index file called indexfile containing external links to all h5
+    datasets which are in in the h5 filelist fnames. This is fast as datasets
+    will only be externally linked within the new indexfile.
+    Therefore the indexfile will also be small in size.
+
+    Will throw an error, if any of the h5 files contain a dataset under the same name.
+    '''
+    if os.path.isfile(indexfile):
+        # indexfile already exists. do not recreate
+        return
+
+    dirname = os.path.dirname(fnames[0])
+    indexfile = os.path.join(dirname, indexfile)
+
+    def visitf(key):
+        # key is a string
+        # only link if key points to a dataset. Do not link groups
+        if isinstance(hf[key], h5py._hl.dataset.Dataset):
+            ih[key] = h5py.ExternalLink(fname, key)
+
+    with h5py.File(indexfile, 'w') as ih:
+        for fname in fnames:
+            with h5py.File(fname, 'r') as hf:
+                hf.visit(visitf)
+
+
+def _getindexfile(path):
+    '''
+    Returns the name of the index file after the file has been
+    generated (File generation only if it doesnt exist)
+    '''
+    indexfile = os.path.join(path, '.postpic-smilei-index.h5')
+    if not os.path.isfile(indexfile):
+        # find all h5 files
+        h5files = glob.glob(os.path.join(path, '*.h5'))
+        print('generating index file "{}" from the following'
+              'h5 files: {}'.format(indexfile, h5files))
+        _generateh5indexfile(indexfile, h5files)
+    return indexfile
+
+class AMSmileiReader:
+
+    def __init__(self, h5file, iteration=None):
+        #super(OpenPMDreader, self).__init__(h5file)
+        
+        if os.path.isfile(h5file):
+            self._h5 = h5py.File(h5file, 'r')
+        elif os.path.isdir(h5file):
+            indexfile = _getindexfile(h5file)
+            self._h5 = h5py.File(indexfile, 'r')
+        else:
+            raise IOError('"{}" is neither a h5 file nor a directory'
+                          'containing h5 files'.format(h5file))
+
+        #<<<<<Evaluating time in both SI units and smilei units>>>>>
+        self.timestep_arr = np.array(self._h5['/data/'])  #returns available iterations
+
+        if '{:010d}'.format(iteration) not in self.timestep_arr:
+            #Finding the closest integer iteration available from the user iteration
+            if iteration is not None:
+                int_iter = min(self.timestep_arr, key=lambda x: (abs(x - int(iteration)), x))
+            else:
+                int_iter = self.timestep_arr[-1]
+        
+            self.sim_iter = '{:010d}'.format(int_iter)   #timestep format like in the simulation dump     
+        else:
+            self.sim_iter = '{:010d}'.format(iteration)
+        
+        
+        #-------------------------------------------------------------------------------------------------------------------------
+        #Directly importing data files, because the field data and particle data were dumped in two different h5 files by smilei
+        #And I am having confusion with the index file :(
+        try:
+            #field data file
+            self._dataF = h5py.File('/path/Fields0.h5', 'r')['/data/{}'.format(self.sim_iter)]
+        except:
+            pass
+
+        try:
+            #particle data file
+            self._dataP = h5py.File('/path/TrackParticlesDisordered_electron.h5', 'r')['/data/{}/particles'.format(self.sim_iter)]
+        except:
+            pass
+        #-------------------------------------------------------------------------------------------------------------------------
+        #time in SI units
+        self.real_time=(self._dataF.attrs['time']) * (self._dataF.attrs['timeUnitSI'])
+
+    def attrs(self, data):
+        return data.attrs
+
+    def __del__(self):
+        del self._dataF
+
+    # --- Level 0 methods ---
+
+    def keys(self):
+        return list(self._dataF.keys())
+
+    def __getitem__(self, key):
+        return self._dataF[key]
+    
+    # --- Level 1 methods ---
+
+    #To return the number of AM modes and Field names in the dump.
+    def getAMmodes(self):    
+        strings=np.array(self._dataF.keys())
+        max_suffix = float('-inf')
+        max_suffix_string = None
+        prefix_list=[]
+
+        for s in strings:
+            prefix, suffix = s.split('_mode_')
+            suffix_int = int(suffix)
+            if suffix_int > max_suffix:
+                max_suffix = suffix_int
+                max_suffix_string = s
+            prefix_list.append(prefix)
+
+        return [prefix_list, int(max_suffix_string[-1])+1]    #[field_names_preffix, number of modes]
+    
+    def Data(self, key, modes):
+        '''
+        should work with any key, that contains data, thus on every hdf5.Dataset,
+        but not on hdf5.Group. Will extract the data, convert it to SI and return it
+        as a numpy array of complex field data.
+        '''
+        if modes is None:
+        
+            field_array = np.array(self._dataF[key])
+            field_array_shape = field_array.shape
+            reshaped_array = field_array.reshape(field_array_shape[0], field_array_shape[1]//2, 2)
+            complex_array = reshaped_array[:, :, 0] + 1j * reshaped_array[:, :, 1]
+
+            return complex_array   #complex_array consist of field data in complex numbers
+        
+        elif modes =="all":
+            array_list=[]
+            for mode in range(self.getAMmodes[-1]):
+
+                field_name = key+"_mode_"+str(mode)
+                field_array = np.array(self._dataF[field_name])
+                field_array_shape = field_array.shape
+                reshaped_array = field_array.reshape(field_array_shape[0], field_array_shape[1]//2, 2)
+                complex_array = reshaped_array[:, :, 0] + 1j * reshaped_array[:, :, 1]
+                array_list.append(complex_array)
+        
+            mod_complex_data= np.stack(array_list,axis=0)     #Modified array of shape (Nmodes, Nx, Nr)
+            return mod_complex_data
+    
+    #-----------MODE EXPANSION METHODS-----------
+
+    @staticmethod
+    def _modeexpansion_naiv_single(complex_data, theta=0):
+    
+        F = np.zeros_like(np.real(complex_data[0]))                       
+        
+        '''
+        for m in self.modes:
+            F += self.mod_data[m]*np.exp(-1j*m*self.theta)
+        '''
+        for m in range(getAMmodes[-1]):
+            F += np.real(complex_data[m])*np.cos(m*theta)+np.imag(complex_data[m])*np.sin(m*theta)
+
+        return F 
+    
+    @staticmethod
+    def _modeexpansion_naiv(complex_data, theta=0):
+        '''
+        converts to radial data using `modeexpansion`, possibly for multiple
+        theta at once.
+        '''
+        if np.asarray(theta).shape is ():
+            # single theta
+            theta = [theta]
+        # multiple theta
+        data = np.asarray([_modeexpansion_naiv_single(complex_data, theta=t)
+                           for t in theta])
+        # switch from (theta, r, z) to (r, theta, z)
+        data = data.swapaxes(0, 1)
+        return data
+    
+    @staticmethod
+    def _modeexpansion_fft(complex_data, Ntheta=None):
+        '''
+        calculate the radialdata using an fft. This is by far the fastest
+        way to do the modeexpansion.
+        '''
+        Nm, Nx, Nr = complex_data.shape
+        Nth = (Nm+1)//2
+        if Ntheta is None or Ntheta < Nth:
+            Ntheta = Nth
+        fd = np.empty((Nr, Ntheta, Nx), dtype=np.complex128)
+
+        fd[:, 0, :].real = complex_data[0, :, :]
+        rawdatasw = np.swapaxes(complex_data, 0, 1)
+        fd[:, 1:Nth, :].real = rawdatasw[:, 1::2, :]
+        fd[:, 1:Nth, :].imag = rawdatasw[:, 2::2, :]
+
+        fd = fft.fft(fd, axis=1).real
+        return fd
+    
+    @staticmethod
+    def mode_expansion(complex_data, theta=None, Ntheta=None):
+        
+        Nm, Nr, Nz = complex_data.shape
+        if Ntheta is not None or theta is None:
+            return _modeexpansion_fft(complex_data, Ntheta=Ntheta)
+        else:
+            return _modeexpansion_naiv(complex_data, theta=theta)
+
+    #----------------------------------------------------------------------------------    
+
+    def gridoffset(self, key, axis):
+        axid = helper.axesidentify[axis]
+        if axid ==90 or axis=='r':
+            axid=axid%90
+
+        if "gridUnitSI" in self._dataF[key].attrs:
+            attrs = self._dataF[key].attrs
+        else:
+            attrs = self._dataF[key].parent.attrs
+        return attrs['gridGlobalOffset'][axid] * attrs['gridUnitSI']
+
+    def gridspacing(self, key, axis):
+        axid = helper.axesidentify[axis]
+        if axid ==90 or axis=='r':
+            axid=axid%90
+
+        if "gridUnitSI" in self._dataF[key].attrs:
+            attrs = self._dataF[key].attrs
+        else:
+            attrs = self._dataF[key].parent.attrs
+        return attrs['gridSpacing'][axid] * attrs['gridUnitSI']
+    
+    #----------------Here def gridpoints(self) is missing temprorarily-----------------
+
+    # --- Level 2 methods ---
+
+    def timestep(self):
+        return self.sim_iter
+
+    def time(self):
+        return self.real_time
+
+    def simdimensions(self):
+        '''
+        the number of spatial dimensions the simulation was using.
+        '''
+        for k in self._simgridkeys():
+            try:
+                gs = self.gridspacing(k, None)
+                return len(gs)
+            except KeyError:
+                pass
+        raise KeyError('number of simdimensions could not be retrieved for {}'.format(self))
+
+    def _keyE(self, component, **kwargs):
+        # Smilei deviates from openPMD standard: Ex instead of E/x
+        axsuffix = {0: 'x', 1: 'y', 2: 'z', 90: 'r', 91: 't'}[helper.axesidentify[component]]
+        return 'E{}'.format(axsuffix)
+
+    def _keyB(self, component, **kwargs):
+        # Smilei deviates from openPMD standard: Bx instead of B/x
+        axsuffix = {0: 'x', 1: 'y', 2: 'z', 90: 'r', 91: 't'}[helper.axesidentify[component]]
+        return 'B{}'.format(axsuffix)
+
+    def _simgridkeys(self):
+        return ['El', 'Er', 'Et',
+                'Bl', 'Br', 'Bt']
+
+    def listSpecies(self):
+        ret = list(self._dataP.keys())
+        return ret
+
+    def getSpecies(self, species, attrib):
+        """
+        Returns one of the attributes out of (x,y,z,px,py,pz,weight,ID,mass,charge) of
+        this particle species.
+        """
+        attribid = helper.attribidentify[attrib]
+        options = {9: 'particles/{}/weighting',
+                   0: 'particles/{}/position/x',
+                   1: 'particles/{}/position/y',
+                   2: 'particles/{}/position/z',
+                   3: 'particles/{}/momentum/x',
+                   4: 'particles/{}/momentum/y',
+                   5: 'particles/{}/momentum/z',
+                   10: 'particles/{}/id',
+                   11: 'particles/{}/mass',
+                   12: 'particles/{}/charge'}
+        optionsoffset = {0: 'particles/{}/positionOffset/x',
+                         1: 'particles/{}/positionOffset/y',
+                         2: 'particles/{}/positionOffset/z'}
+        key = options[attribid]
+        offsetkey = optionsoffset.get(attribid)
+        try:
+            data = self._dataP(key.format(species))
+            if offsetkey is not None:
+                data += self._dataP(offsetkey.format(species))
+            ret = np.asarray(data, dtype=np.float64)
+        except IndexError:
+            raise KeyError
+        return ret
+
+    #To get the axis array.
+    def getAxis(self, axis):
+    
+        namelist = list(self._dataF.keys())
+        Nx,Nr = self._dataF[namelist[0]].shape
+        Nr=Nr/2
+
+        if axis is None:
+            raise IOError("Invalid axis")
+        elif len(axis)>1:
+            raise IOError("Only one axis at a time")
+        
+        if axis == "x":      #If moving = True, the x axis data modifies according to the moving window.
+            x_min = self._dataF.attrs['x_moved'] * self._dataF[namelist[0]].attrs['gridUnitSI'] 
+            x_max = x_min + self.gridspacing(key=namelist[0],axis=0)*(Nx-1)
+            x_axis = np.linspace(x_min, x_max, Nx-1)
+            return x_axis   
+        elif axis == "r":
+            r_max = self.gridspacing(key=namelist[0],axis='r')*(Nr-1)
+            r_axis = np.linspace(0, r_max, Nr-1)
+            return r_axis
+        else:
+            raise IOError("Invalid axis")
+    
+    '''
+    def getderived(self):
+        
+        #return all other fields dumped, except E and B.
+        
+        ret = []
+        self['fields'].visit(ret.append)
+        ret = ['fields/{}'.format(r) for r in ret if not (r.startswith('E') or r.startswith('B'))]
+        ret = [r for r in ret if hasattr(self[r], 'value')]
+        ret.sort()
+        return ret
+    '''
+
+    def __str__(self):
+        return '<OpenPMDh5reader at "' + str(self.dumpidentifier) + '">'
+    
+
+class SmileiSeries:
+    '''
+    Reads a time series of dumps from a given directory.
+
+    Point this to the directory and you will get all dumps,
+    but possibly containing different data.
+    `simreader = SmileiSeries('path/to/simulation/')`
+
+    Alternatively point this to a single file and you will only get
+    the iterations which are present in that file:
+    `simreader = SmileiSeries('path/to/simulation/Fields0.h5')`
+    '''
+
+    def __init__(self, h5file, dumpreadercls=AMSmileiReader, **kwargs):
+        super(SmileiSeries, self).__init__(h5file, **kwargs)
+        self.dumpreadercls = dumpreadercls
+        self.h5file = h5file
+        self.path = os.path.basename(h5file)
+        if os.path.isfile(h5file):
+            with h5py.File(h5file, 'r') as h5:
+                self._dumpkeys = list(h5['data'].keys())
+        elif os.path.isdir(h5file):
+            indexfile = _getindexfile(h5file)
+            self._h5 = h5py.File(indexfile, 'r')
+            with h5py.File(indexfile, 'r') as h5:
+                self._dumpkeys = list(h5['data'].keys())
+        else:
+            raise IOError('{} does not exist.'.format(h5file))
+
+    def _getDumpreader(self, n):
+        '''
+        Do not use this method. It will be called by __getitem__.
+        Use __getitem__ instead.
+        '''
+        return self.dumpreadercls(self.h5file, self._dumpkeys[n])
+
+    def __len__(self):
+        return len(self._dumpkeys)
+
+    def __str__(self):
+        return '<SmileiSeries based on "{}">'.format(self.simidentifier)

From 32265af58a9d284054ce3ac93519b8b7e9226442 Mon Sep 17 00:00:00 2001
From: VinithSamson11 <109979733+VinithSamson11@users.noreply.github.com>
Date: Tue, 16 Jan 2024 02:19:14 +0530
Subject: [PATCH 02/10] Update smileih5.py

---
 postpic/datareader/smileih5.py | 139 +++++++++++++++++++++++++++++++--
 1 file changed, 131 insertions(+), 8 deletions(-)

diff --git a/postpic/datareader/smileih5.py b/postpic/datareader/smileih5.py
index 95c3597..fe4c0a3 100644
--- a/postpic/datareader/smileih5.py
+++ b/postpic/datareader/smileih5.py
@@ -26,6 +26,9 @@
 import os.path
 import h5py
 import glob
+import numpy as np
+from .. import helper
+from ..helper_fft import fft
 
 __all__ = ['SmileiReader', 'SmileiSeries']
 
@@ -72,7 +75,6 @@ def _getindexfile(path):
         _generateh5indexfile(indexfile, h5files)
     return indexfile
 
-
 class SmileiReader(OpenPMDreader):
     '''
     The Smilei reader can read  a single h5 file or combine the information of
@@ -114,16 +116,117 @@ def __init__(self, h5file, iteration=None):
 
         if iteration is None:
             self._iteration = int(list(self._h5['data'].keys())[0])
+        elif iteration not in [int(i) for i in list(self._h5['data'].keys())]:
+            raise IOError("Iteration {} is in valid".format(iteration))
         else:
             self._iteration = int(iteration)
+        
         self._data = self._h5['/data/{:010d}/'.format(self._iteration)]
         self.attrs = self._data.attrs
 
-# --- Level 0 methods ---
-
-# --- Level 1 methods ---
-
-# --- Level 2 methods ---
+    #To return the number of AM modes and Field names in the dump.
+    def getAMmodes(self):    
+        strings=np.array(self._data.keys())
+        mask =  np.array(["_mode_" in s for s in strings])
+        arr=strings[mask]
+        max_suffix = float('-inf')
+        max_suffix_string = None
+        prefix_list=[]
+
+        for i in arr:
+            prefix, suffix = i.split('_mode_')
+            suffix_int = int(suffix)
+            if suffix_int > max_suffix:
+                max_suffix = suffix_int
+                max_suffix_string = i
+            prefix_list.append(prefix)
+
+        availModes=[i for i in range(0,int(max_suffix_string[-1])+1)]
+        return [prefix_list, int(max_suffix_string[-1])+1, availModes]
+    
+    def getComplex(self,key,modes):
+        
+        array_list=[]
+        for mode in modes:
+
+            field_name = key+"_mode_"+str(mode)
+            field_array = np.array(self._data[field_name])
+            field_array_shape = field_array.shape
+            reshaped_array = field_array.reshape(field_array_shape[0], field_array_shape[1]//2, 2)
+            complex_array = reshaped_array[:, :, 0] + 1j * reshaped_array[:, :, 1]
+            array_list.append(complex_array)
+        
+        mod_complex_data= np.stack(array_list,axis=0)                    #Modified array of shape (Nmodes, Nx, Nr)
+        mod_complex_data = np.transpose(mod_complex_data, (0, 2, 1))     #Modified array of shape (Nmodes, Nr, Nx)
+        return mod_complex_data
+    
+    def rawData(self, key, modes):
+        '''
+        should work with any key, that contains data, thus on every hdf5.Dataset,
+        but not on hdf5.Group. Will extract the data, convert it to SI and return it
+        as a numpy array of complex field data.
+        '''
+        def validateModes(lst):
+            for i in range(1, len(lst)):
+                if lst[i] != lst[i - 1] + 1:
+                    return False
+            return True
+
+        if modes is not None:
+            if validateModes(modes)==False:
+                raise IOError('Invalid modes {} order'.format(modes))
+            
+            
+            if set(modes).issubset(set(self.getAMmodes[-1]))==False:
+                raise IOError("Mode(s) {} not available in dump".format(modes))
+            if key not in self.getAMmodes[0]:
+                raise IOError("Key {} not available in dump".format(key))
+            
+            return self.getComplex(key=key,modes=modes)
+            
+        if modes is None:
+        
+            if key not in list(self._data):
+                raise IOError("Key {} not available in dump".format(key))
+            
+            else:
+                return np.array(self._data[key]) 
+        
+        if modes =="all":
+            return self.getComplex(key=key,modes=self.getAMmodes[-1]) 
+        
+    @staticmethod
+    def _modeexpansion_naiv(rawdata, theta=0, M=None):
+        '''
+        converts to radial data using `modeexpansion`, possibly for multiple
+        theta at once.
+        '''
+        if np.asarray(theta).shape is ():
+            # single theta
+            theta = [theta]
+        # multiple theta
+        #data = np.asarray([SmileiReader._modeexpansion_naiv_single(rawdata, theta=t)
+        #                   for t in theta])
+        # switch from (theta, r, z) to (r, theta, z)
+        
+        Nt=len(theta)
+        (Nm, Nr, Nx) = rawdata.shape
+        F_total = np.zeros((Nt,Nr,Nx))                       
+        mode =[m for m in range(0,Nm)]
+
+        if M is not None:
+            if Nm ==1 and M not in SmileiReader.getAMmodes[-1]:
+                raise IOError('Mode {} not available in dump'.format(M))
+
+        for t in theta:
+            index=theta.index(t)
+            if M is None:
+                for m in mode:
+                    F_total[index] += np.real(rawdata[m])*np.cos(m*t)+np.imag(rawdata[m])*np.sin(m*t)
+            else:
+                F_total[index] += np.real(rawdata[0])*np.cos(M*t)+np.imag(rawdata[0])*np.sin(M*t)
+        #F_total = F_total.swapaxes(0, 1)
+        return F_total
 
     def _keyE(self, component, **kwargs):
         # Smilei deviates from openPMD standard: Ex instead of E/x
@@ -137,8 +240,8 @@ def _keyB(self, component, **kwargs):
 
     def _simgridkeys(self):
         # Smilei deviates from openPMD standard: Ex instead of E/x
-        return ['Ex', 'Ey', 'Ez',
-                'Bx', 'By', 'Bz']
+        return ['Ex', 'Ey', 'Ez','Er','El','Et',
+                'Bx', 'By', 'Bz','Br','Bl','Bt']
 
     def getderived(self):
         '''
@@ -153,6 +256,26 @@ def getderived(self):
     def __str__(self):
         return '<SmileiReader at "{}" at iteration {:d}>'.format(self.dumpidentifier,
                                                                  self.timestep())
+    
+    # override inherited method to count points after mode expansion
+    def gridoffset(self, key, axis):
+        axid = helper.axesidentify[axis]
+        if axid == 91:  # theta
+            return 0
+        else:
+            # r, theta, z
+            axidremap = {90: 0, 2: 1}[axid]
+            return super(SmileiReader, self).gridoffset(key, axidremap)
+
+    # override inherited method to count points after mode expansion
+    def gridspacing(self, key, axis):
+        axid = helper.axesidentify[axis]
+        if axid == 91:  # theta
+            return 2 * np.pi / self.gridpoints(key, axis)
+        else:
+            # r, theta, z
+            axidremap = {90: 0, 2: 1}[axid]
+            return super(SmileiReader, self).gridspacing(key, axidremap)
 
 
 class SmileiSeries(Simulationreader_ifc):

From d0b30a408778ba70619fb57728853a0134d463b5 Mon Sep 17 00:00:00 2001
From: VinithSamson11 <109979733+VinithSamson11@users.noreply.github.com>
Date: Sat, 20 Jan 2024 06:16:04 +0530
Subject: [PATCH 03/10] Updated the modeexpansion method and data reading
 function

---
 postpic/datareader/smileih5.py | 106 ++++++++++++++++-----------------
 1 file changed, 52 insertions(+), 54 deletions(-)

diff --git a/postpic/datareader/smileih5.py b/postpic/datareader/smileih5.py
index fe4c0a3..cb33e45 100644
--- a/postpic/datareader/smileih5.py
+++ b/postpic/datareader/smileih5.py
@@ -17,7 +17,6 @@
 # Stephan Kuschel, 2023
 # Carolin Goll, 2023
 
-
 from __future__ import absolute_import, division, print_function, unicode_literals
 
 from .openPMDh5 import OpenPMDreader
@@ -29,6 +28,7 @@
 import numpy as np
 from .. import helper
 from ..helper_fft import fft
+import re
 
 __all__ = ['SmileiReader', 'SmileiSeries']
 
@@ -126,7 +126,7 @@ def __init__(self, h5file, iteration=None):
 
     #To return the number of AM modes and Field names in the dump.
     def getAMmodes(self):    
-        strings=np.array(self._data.keys())
+        strings=np.array(self._data)
         mask =  np.array(["_mode_" in s for s in strings])
         arr=strings[mask]
         max_suffix = float('-inf')
@@ -144,57 +144,6 @@ def getAMmodes(self):
         availModes=[i for i in range(0,int(max_suffix_string[-1])+1)]
         return [prefix_list, int(max_suffix_string[-1])+1, availModes]
     
-    def getComplex(self,key,modes):
-        
-        array_list=[]
-        for mode in modes:
-
-            field_name = key+"_mode_"+str(mode)
-            field_array = np.array(self._data[field_name])
-            field_array_shape = field_array.shape
-            reshaped_array = field_array.reshape(field_array_shape[0], field_array_shape[1]//2, 2)
-            complex_array = reshaped_array[:, :, 0] + 1j * reshaped_array[:, :, 1]
-            array_list.append(complex_array)
-        
-        mod_complex_data= np.stack(array_list,axis=0)                    #Modified array of shape (Nmodes, Nx, Nr)
-        mod_complex_data = np.transpose(mod_complex_data, (0, 2, 1))     #Modified array of shape (Nmodes, Nr, Nx)
-        return mod_complex_data
-    
-    def rawData(self, key, modes):
-        '''
-        should work with any key, that contains data, thus on every hdf5.Dataset,
-        but not on hdf5.Group. Will extract the data, convert it to SI and return it
-        as a numpy array of complex field data.
-        '''
-        def validateModes(lst):
-            for i in range(1, len(lst)):
-                if lst[i] != lst[i - 1] + 1:
-                    return False
-            return True
-
-        if modes is not None:
-            if validateModes(modes)==False:
-                raise IOError('Invalid modes {} order'.format(modes))
-            
-            
-            if set(modes).issubset(set(self.getAMmodes[-1]))==False:
-                raise IOError("Mode(s) {} not available in dump".format(modes))
-            if key not in self.getAMmodes[0]:
-                raise IOError("Key {} not available in dump".format(key))
-            
-            return self.getComplex(key=key,modes=modes)
-            
-        if modes is None:
-        
-            if key not in list(self._data):
-                raise IOError("Key {} not available in dump".format(key))
-            
-            else:
-                return np.array(self._data[key]) 
-        
-        if modes =="all":
-            return self.getComplex(key=key,modes=self.getAMmodes[-1]) 
-        
     @staticmethod
     def _modeexpansion_naiv(rawdata, theta=0, M=None):
         '''
@@ -227,7 +176,47 @@ def _modeexpansion_naiv(rawdata, theta=0, M=None):
                 F_total[index] += np.real(rawdata[0])*np.cos(M*t)+np.imag(rawdata[0])*np.sin(M*t)
         #F_total = F_total.swapaxes(0, 1)
         return F_total
+    
+    def getComplex(self, key,theta=0):
+        
+        array_list=[]
+        modes=self.getAMmodes()[-1]
+        for mode in modes:
 
+            field_name = key+"_mode_"+str(mode)
+            field_array = np.array(self._data[field_name])*(self._data[field_name].attrs['unitSI'])
+            field_array_shape = field_array.shape
+            reshaped_array = field_array.reshape(field_array_shape[0], field_array_shape[1]//2, 2)
+            complex_array = reshaped_array[:, :, 0] + 1j * reshaped_array[:, :, 1]
+            array_list.append(complex_array)
+        
+        mod_complex_data= np.stack(array_list,axis=0)                    #Modified array of shape (Nmodes, Nx, Nr)
+        #mod_complex_data = np.transpose(mod_complex_data, (0, 2, 1))     #Modified array of shape (Nmodes, Nr, Nx)
+        return SmileiReader._modeexpansion_naiv(rawdata=mod_complex_data,theta=theta)
+    
+    def data(self, key, **kwargs):
+        '''
+        should work with any key, that contains data, thus on every hdf5.Dataset,
+        but not on hdf5.Group. Will extract the data, convert it to SI and return it
+        as a numpy array. Constant records will be detected and converted to
+        a numpy array containing a single value only.
+        '''
+        field_pattern=re.compile(r'^{}.*_mode_'.format(key))
+        match = [True if re.match(field_pattern, s) else False for s in list(self._data)]
+
+        if True in match:
+            return self.getComplex(key=key, **kwargs)
+        else:
+            record = self._data[key]
+        
+        if "value" in record.attrs:
+            # constant data (a single int or float)
+            ret = np.float64(record.attrs['value']) * record.attrs['unitSI']
+        else:
+            # array data
+            ret = np.float64(record[()]) * record.attrs['unitSI']
+        return ret
+        
     def _keyE(self, component, **kwargs):
         # Smilei deviates from openPMD standard: Ex instead of E/x
         axsuffix = {0: 'x', 1: 'y', 2: 'z', 90: 'r', 91: 't'}[helper.axesidentify[component]]
@@ -241,7 +230,8 @@ def _keyB(self, component, **kwargs):
     def _simgridkeys(self):
         # Smilei deviates from openPMD standard: Ex instead of E/x
         return ['Ex', 'Ey', 'Ez','Er','El','Et',
-                'Bx', 'By', 'Bz','Br','Bl','Bt']
+                'Bx', 'By', 'Bz','Br','Bl','Bt',
+                'Jx','Jy','Jz','Jr','Jl','Jt','Rho']
 
     def getderived(self):
         '''
@@ -276,6 +266,14 @@ def gridspacing(self, key, axis):
             # r, theta, z
             axidremap = {90: 0, 2: 1}[axid]
             return super(SmileiReader, self).gridspacing(key, axidremap)
+    
+    # override inherited method to count points after mode expansion
+    def gridpoints(self, key, axis):
+        axid = helper.axesidentify[axis]
+        axid =int(axid/90)
+        (Nx, Nr) = self._data[key].shape
+        Nr=Nr/2
+        return (Nx,Nr)[axid]
 
 
 class SmileiSeries(Simulationreader_ifc):

From 6bd72e01c8d82618635629726bebfea0f04f198c Mon Sep 17 00:00:00 2001
From: VinithSamson11 <109979733+VinithSamson11@users.noreply.github.com>
Date: Sat, 20 Jan 2024 06:27:54 +0530
Subject: [PATCH 04/10] Delete postpic/datareader/AMsmileih5.py

---
 postpic/datareader/AMsmileih5.py | 403 -------------------------------
 1 file changed, 403 deletions(-)
 delete mode 100644 postpic/datareader/AMsmileih5.py

diff --git a/postpic/datareader/AMsmileih5.py b/postpic/datareader/AMsmileih5.py
deleted file mode 100644
index 9128d02..0000000
--- a/postpic/datareader/AMsmileih5.py
+++ /dev/null
@@ -1,403 +0,0 @@
-from __future__ import absolute_import, division, print_function, unicode_literals
-
-from . import Dumpreader_ifc
-from . import Simulationreader_ifc
-import re
-import os.path
-import h5py
-import glob
-import numpy as np
-import helper 
-from helper_fft import fft
-
-__all__ = ['AMSmileiReader', 'SmileiSeries']
-
-
-def _generateh5indexfile(indexfile, fnames):
-    '''
-    Creates a h5 index file called indexfile containing external links to all h5
-    datasets which are in in the h5 filelist fnames. This is fast as datasets
-    will only be externally linked within the new indexfile.
-    Therefore the indexfile will also be small in size.
-
-    Will throw an error, if any of the h5 files contain a dataset under the same name.
-    '''
-    if os.path.isfile(indexfile):
-        # indexfile already exists. do not recreate
-        return
-
-    dirname = os.path.dirname(fnames[0])
-    indexfile = os.path.join(dirname, indexfile)
-
-    def visitf(key):
-        # key is a string
-        # only link if key points to a dataset. Do not link groups
-        if isinstance(hf[key], h5py._hl.dataset.Dataset):
-            ih[key] = h5py.ExternalLink(fname, key)
-
-    with h5py.File(indexfile, 'w') as ih:
-        for fname in fnames:
-            with h5py.File(fname, 'r') as hf:
-                hf.visit(visitf)
-
-
-def _getindexfile(path):
-    '''
-    Returns the name of the index file after the file has been
-    generated (File generation only if it doesnt exist)
-    '''
-    indexfile = os.path.join(path, '.postpic-smilei-index.h5')
-    if not os.path.isfile(indexfile):
-        # find all h5 files
-        h5files = glob.glob(os.path.join(path, '*.h5'))
-        print('generating index file "{}" from the following'
-              'h5 files: {}'.format(indexfile, h5files))
-        _generateh5indexfile(indexfile, h5files)
-    return indexfile
-
-class AMSmileiReader:
-
-    def __init__(self, h5file, iteration=None):
-        #super(OpenPMDreader, self).__init__(h5file)
-        
-        if os.path.isfile(h5file):
-            self._h5 = h5py.File(h5file, 'r')
-        elif os.path.isdir(h5file):
-            indexfile = _getindexfile(h5file)
-            self._h5 = h5py.File(indexfile, 'r')
-        else:
-            raise IOError('"{}" is neither a h5 file nor a directory'
-                          'containing h5 files'.format(h5file))
-
-        #<<<<<Evaluating time in both SI units and smilei units>>>>>
-        self.timestep_arr = np.array(self._h5['/data/'])  #returns available iterations
-
-        if '{:010d}'.format(iteration) not in self.timestep_arr:
-            #Finding the closest integer iteration available from the user iteration
-            if iteration is not None:
-                int_iter = min(self.timestep_arr, key=lambda x: (abs(x - int(iteration)), x))
-            else:
-                int_iter = self.timestep_arr[-1]
-        
-            self.sim_iter = '{:010d}'.format(int_iter)   #timestep format like in the simulation dump     
-        else:
-            self.sim_iter = '{:010d}'.format(iteration)
-        
-        
-        #-------------------------------------------------------------------------------------------------------------------------
-        #Directly importing data files, because the field data and particle data were dumped in two different h5 files by smilei
-        #And I am having confusion with the index file :(
-        try:
-            #field data file
-            self._dataF = h5py.File('/path/Fields0.h5', 'r')['/data/{}'.format(self.sim_iter)]
-        except:
-            pass
-
-        try:
-            #particle data file
-            self._dataP = h5py.File('/path/TrackParticlesDisordered_electron.h5', 'r')['/data/{}/particles'.format(self.sim_iter)]
-        except:
-            pass
-        #-------------------------------------------------------------------------------------------------------------------------
-        #time in SI units
-        self.real_time=(self._dataF.attrs['time']) * (self._dataF.attrs['timeUnitSI'])
-
-    def attrs(self, data):
-        return data.attrs
-
-    def __del__(self):
-        del self._dataF
-
-    # --- Level 0 methods ---
-
-    def keys(self):
-        return list(self._dataF.keys())
-
-    def __getitem__(self, key):
-        return self._dataF[key]
-    
-    # --- Level 1 methods ---
-
-    #To return the number of AM modes and Field names in the dump.
-    def getAMmodes(self):    
-        strings=np.array(self._dataF.keys())
-        max_suffix = float('-inf')
-        max_suffix_string = None
-        prefix_list=[]
-
-        for s in strings:
-            prefix, suffix = s.split('_mode_')
-            suffix_int = int(suffix)
-            if suffix_int > max_suffix:
-                max_suffix = suffix_int
-                max_suffix_string = s
-            prefix_list.append(prefix)
-
-        return [prefix_list, int(max_suffix_string[-1])+1]    #[field_names_preffix, number of modes]
-    
-    def Data(self, key, modes):
-        '''
-        should work with any key, that contains data, thus on every hdf5.Dataset,
-        but not on hdf5.Group. Will extract the data, convert it to SI and return it
-        as a numpy array of complex field data.
-        '''
-        if modes is None:
-        
-            field_array = np.array(self._dataF[key])
-            field_array_shape = field_array.shape
-            reshaped_array = field_array.reshape(field_array_shape[0], field_array_shape[1]//2, 2)
-            complex_array = reshaped_array[:, :, 0] + 1j * reshaped_array[:, :, 1]
-
-            return complex_array   #complex_array consist of field data in complex numbers
-        
-        elif modes =="all":
-            array_list=[]
-            for mode in range(self.getAMmodes[-1]):
-
-                field_name = key+"_mode_"+str(mode)
-                field_array = np.array(self._dataF[field_name])
-                field_array_shape = field_array.shape
-                reshaped_array = field_array.reshape(field_array_shape[0], field_array_shape[1]//2, 2)
-                complex_array = reshaped_array[:, :, 0] + 1j * reshaped_array[:, :, 1]
-                array_list.append(complex_array)
-        
-            mod_complex_data= np.stack(array_list,axis=0)     #Modified array of shape (Nmodes, Nx, Nr)
-            return mod_complex_data
-    
-    #-----------MODE EXPANSION METHODS-----------
-
-    @staticmethod
-    def _modeexpansion_naiv_single(complex_data, theta=0):
-    
-        F = np.zeros_like(np.real(complex_data[0]))                       
-        
-        '''
-        for m in self.modes:
-            F += self.mod_data[m]*np.exp(-1j*m*self.theta)
-        '''
-        for m in range(getAMmodes[-1]):
-            F += np.real(complex_data[m])*np.cos(m*theta)+np.imag(complex_data[m])*np.sin(m*theta)
-
-        return F 
-    
-    @staticmethod
-    def _modeexpansion_naiv(complex_data, theta=0):
-        '''
-        converts to radial data using `modeexpansion`, possibly for multiple
-        theta at once.
-        '''
-        if np.asarray(theta).shape is ():
-            # single theta
-            theta = [theta]
-        # multiple theta
-        data = np.asarray([_modeexpansion_naiv_single(complex_data, theta=t)
-                           for t in theta])
-        # switch from (theta, r, z) to (r, theta, z)
-        data = data.swapaxes(0, 1)
-        return data
-    
-    @staticmethod
-    def _modeexpansion_fft(complex_data, Ntheta=None):
-        '''
-        calculate the radialdata using an fft. This is by far the fastest
-        way to do the modeexpansion.
-        '''
-        Nm, Nx, Nr = complex_data.shape
-        Nth = (Nm+1)//2
-        if Ntheta is None or Ntheta < Nth:
-            Ntheta = Nth
-        fd = np.empty((Nr, Ntheta, Nx), dtype=np.complex128)
-
-        fd[:, 0, :].real = complex_data[0, :, :]
-        rawdatasw = np.swapaxes(complex_data, 0, 1)
-        fd[:, 1:Nth, :].real = rawdatasw[:, 1::2, :]
-        fd[:, 1:Nth, :].imag = rawdatasw[:, 2::2, :]
-
-        fd = fft.fft(fd, axis=1).real
-        return fd
-    
-    @staticmethod
-    def mode_expansion(complex_data, theta=None, Ntheta=None):
-        
-        Nm, Nr, Nz = complex_data.shape
-        if Ntheta is not None or theta is None:
-            return _modeexpansion_fft(complex_data, Ntheta=Ntheta)
-        else:
-            return _modeexpansion_naiv(complex_data, theta=theta)
-
-    #----------------------------------------------------------------------------------    
-
-    def gridoffset(self, key, axis):
-        axid = helper.axesidentify[axis]
-        if axid ==90 or axis=='r':
-            axid=axid%90
-
-        if "gridUnitSI" in self._dataF[key].attrs:
-            attrs = self._dataF[key].attrs
-        else:
-            attrs = self._dataF[key].parent.attrs
-        return attrs['gridGlobalOffset'][axid] * attrs['gridUnitSI']
-
-    def gridspacing(self, key, axis):
-        axid = helper.axesidentify[axis]
-        if axid ==90 or axis=='r':
-            axid=axid%90
-
-        if "gridUnitSI" in self._dataF[key].attrs:
-            attrs = self._dataF[key].attrs
-        else:
-            attrs = self._dataF[key].parent.attrs
-        return attrs['gridSpacing'][axid] * attrs['gridUnitSI']
-    
-    #----------------Here def gridpoints(self) is missing temprorarily-----------------
-
-    # --- Level 2 methods ---
-
-    def timestep(self):
-        return self.sim_iter
-
-    def time(self):
-        return self.real_time
-
-    def simdimensions(self):
-        '''
-        the number of spatial dimensions the simulation was using.
-        '''
-        for k in self._simgridkeys():
-            try:
-                gs = self.gridspacing(k, None)
-                return len(gs)
-            except KeyError:
-                pass
-        raise KeyError('number of simdimensions could not be retrieved for {}'.format(self))
-
-    def _keyE(self, component, **kwargs):
-        # Smilei deviates from openPMD standard: Ex instead of E/x
-        axsuffix = {0: 'x', 1: 'y', 2: 'z', 90: 'r', 91: 't'}[helper.axesidentify[component]]
-        return 'E{}'.format(axsuffix)
-
-    def _keyB(self, component, **kwargs):
-        # Smilei deviates from openPMD standard: Bx instead of B/x
-        axsuffix = {0: 'x', 1: 'y', 2: 'z', 90: 'r', 91: 't'}[helper.axesidentify[component]]
-        return 'B{}'.format(axsuffix)
-
-    def _simgridkeys(self):
-        return ['El', 'Er', 'Et',
-                'Bl', 'Br', 'Bt']
-
-    def listSpecies(self):
-        ret = list(self._dataP.keys())
-        return ret
-
-    def getSpecies(self, species, attrib):
-        """
-        Returns one of the attributes out of (x,y,z,px,py,pz,weight,ID,mass,charge) of
-        this particle species.
-        """
-        attribid = helper.attribidentify[attrib]
-        options = {9: 'particles/{}/weighting',
-                   0: 'particles/{}/position/x',
-                   1: 'particles/{}/position/y',
-                   2: 'particles/{}/position/z',
-                   3: 'particles/{}/momentum/x',
-                   4: 'particles/{}/momentum/y',
-                   5: 'particles/{}/momentum/z',
-                   10: 'particles/{}/id',
-                   11: 'particles/{}/mass',
-                   12: 'particles/{}/charge'}
-        optionsoffset = {0: 'particles/{}/positionOffset/x',
-                         1: 'particles/{}/positionOffset/y',
-                         2: 'particles/{}/positionOffset/z'}
-        key = options[attribid]
-        offsetkey = optionsoffset.get(attribid)
-        try:
-            data = self._dataP(key.format(species))
-            if offsetkey is not None:
-                data += self._dataP(offsetkey.format(species))
-            ret = np.asarray(data, dtype=np.float64)
-        except IndexError:
-            raise KeyError
-        return ret
-
-    #To get the axis array.
-    def getAxis(self, axis):
-    
-        namelist = list(self._dataF.keys())
-        Nx,Nr = self._dataF[namelist[0]].shape
-        Nr=Nr/2
-
-        if axis is None:
-            raise IOError("Invalid axis")
-        elif len(axis)>1:
-            raise IOError("Only one axis at a time")
-        
-        if axis == "x":      #If moving = True, the x axis data modifies according to the moving window.
-            x_min = self._dataF.attrs['x_moved'] * self._dataF[namelist[0]].attrs['gridUnitSI'] 
-            x_max = x_min + self.gridspacing(key=namelist[0],axis=0)*(Nx-1)
-            x_axis = np.linspace(x_min, x_max, Nx-1)
-            return x_axis   
-        elif axis == "r":
-            r_max = self.gridspacing(key=namelist[0],axis='r')*(Nr-1)
-            r_axis = np.linspace(0, r_max, Nr-1)
-            return r_axis
-        else:
-            raise IOError("Invalid axis")
-    
-    '''
-    def getderived(self):
-        
-        #return all other fields dumped, except E and B.
-        
-        ret = []
-        self['fields'].visit(ret.append)
-        ret = ['fields/{}'.format(r) for r in ret if not (r.startswith('E') or r.startswith('B'))]
-        ret = [r for r in ret if hasattr(self[r], 'value')]
-        ret.sort()
-        return ret
-    '''
-
-    def __str__(self):
-        return '<OpenPMDh5reader at "' + str(self.dumpidentifier) + '">'
-    
-
-class SmileiSeries:
-    '''
-    Reads a time series of dumps from a given directory.
-
-    Point this to the directory and you will get all dumps,
-    but possibly containing different data.
-    `simreader = SmileiSeries('path/to/simulation/')`
-
-    Alternatively point this to a single file and you will only get
-    the iterations which are present in that file:
-    `simreader = SmileiSeries('path/to/simulation/Fields0.h5')`
-    '''
-
-    def __init__(self, h5file, dumpreadercls=AMSmileiReader, **kwargs):
-        super(SmileiSeries, self).__init__(h5file, **kwargs)
-        self.dumpreadercls = dumpreadercls
-        self.h5file = h5file
-        self.path = os.path.basename(h5file)
-        if os.path.isfile(h5file):
-            with h5py.File(h5file, 'r') as h5:
-                self._dumpkeys = list(h5['data'].keys())
-        elif os.path.isdir(h5file):
-            indexfile = _getindexfile(h5file)
-            self._h5 = h5py.File(indexfile, 'r')
-            with h5py.File(indexfile, 'r') as h5:
-                self._dumpkeys = list(h5['data'].keys())
-        else:
-            raise IOError('{} does not exist.'.format(h5file))
-
-    def _getDumpreader(self, n):
-        '''
-        Do not use this method. It will be called by __getitem__.
-        Use __getitem__ instead.
-        '''
-        return self.dumpreadercls(self.h5file, self._dumpkeys[n])
-
-    def __len__(self):
-        return len(self._dumpkeys)
-
-    def __str__(self):
-        return '<SmileiSeries based on "{}">'.format(self.simidentifier)

From b638eb7e48e9a105f49c810bcd2b0c06e2730f2b Mon Sep 17 00:00:00 2001
From: VinithSamson11 <109979733+VinithSamson11@users.noreply.github.com>
Date: Sun, 21 Jan 2024 23:02:03 +0530
Subject: [PATCH 05/10] updated functions gridpoints, gridspacing, gridoffset

---
 postpic/datareader/smileih5.py | 47 ++++++++++++++++++++++------------
 1 file changed, 30 insertions(+), 17 deletions(-)

diff --git a/postpic/datareader/smileih5.py b/postpic/datareader/smileih5.py
index cb33e45..15a27d5 100644
--- a/postpic/datareader/smileih5.py
+++ b/postpic/datareader/smileih5.py
@@ -177,7 +177,7 @@ def _modeexpansion_naiv(rawdata, theta=0, M=None):
         #F_total = F_total.swapaxes(0, 1)
         return F_total
     
-    def getComplex(self, key,theta=0):
+    def getExpanded(self, key,theta=0):
         
         array_list=[]
         modes=self.getAMmodes()[-1]
@@ -194,6 +194,12 @@ def getComplex(self, key,theta=0):
         #mod_complex_data = np.transpose(mod_complex_data, (0, 2, 1))     #Modified array of shape (Nmodes, Nr, Nx)
         return SmileiReader._modeexpansion_naiv(rawdata=mod_complex_data,theta=theta)
     
+    def check_pattern(self,key):
+        field_pattern=re.compile(r'^{}.*_mode_'.format(key))
+        match = [True if re.match(field_pattern, s) else False for s in list(self._data)]
+
+        return match
+    
     def data(self, key, **kwargs):
         '''
         should work with any key, that contains data, thus on every hdf5.Dataset,
@@ -201,11 +207,9 @@ def data(self, key, **kwargs):
         as a numpy array. Constant records will be detected and converted to
         a numpy array containing a single value only.
         '''
-        field_pattern=re.compile(r'^{}.*_mode_'.format(key))
-        match = [True if re.match(field_pattern, s) else False for s in list(self._data)]
 
-        if True in match:
-            return self.getComplex(key=key, **kwargs)
+        if key in self.getAMmodes()[0]:
+            return self.getExpanded(key=key, **kwargs)
         else:
             record = self._data[key]
         
@@ -250,30 +254,39 @@ def __str__(self):
     # override inherited method to count points after mode expansion
     def gridoffset(self, key, axis):
         axid = helper.axesidentify[axis]
+        
         if axid == 91:  # theta
             return 0
+        elif key in self.getAMmodes()[0] and axid in [0,90]:
+            key="{}_mode_0".format(key)
+            axid = int(axid/90)
+            return super(SmileiReader, self).gridoffset(key=key, axis=axid)
         else:
-            # r, theta, z
-            axidremap = {90: 0, 2: 1}[axid]
-            return super(SmileiReader, self).gridoffset(key, axidremap)
+            return super(SmileiReader, self).gridoffset(key, axis)
 
     # override inherited method to count points after mode expansion
     def gridspacing(self, key, axis):
         axid = helper.axesidentify[axis]
-        if axid == 91:  # theta
-            return 2 * np.pi / self.gridpoints(key, axis)
+        
+        if key in self.getAMmodes()[0] and axid in [0,90]:
+            key="{}_mode_0".format(key)
+            axid = int(axid/90)
+            return super(SmileiReader, self).gridspacing(key=key, axis=axid)
         else:
-            # r, theta, z
-            axidremap = {90: 0, 2: 1}[axid]
-            return super(SmileiReader, self).gridspacing(key, axidremap)
+            return super(SmileiReader, self).gridspacing(key, axis)
     
     # override inherited method to count points after mode expansion
     def gridpoints(self, key, axis):
         axid = helper.axesidentify[axis]
-        axid =int(axid/90)
-        (Nx, Nr) = self._data[key].shape
-        Nr=Nr/2
-        return (Nx,Nr)[axid]
+        
+        if key in self.getAMmodes()[0]:
+            key="{}_mode_0".format(key)
+            (Nx, Nr) = self._data[key].shape
+            Nr=Nr/2
+            axid =int(axid/90)
+            return (Nx,Nr)[axid]
+        else:
+            return super(SmileiReader, self).gridpoints(key=key,axis=axis)
 
 
 class SmileiSeries(Simulationreader_ifc):

From 8b9f174ee4b544a9c4b7cde20e26a5db05507d73 Mon Sep 17 00:00:00 2001
From: VinithSamson11 <109979733+VinithSamson11@users.noreply.github.com>
Date: Mon, 22 Jan 2024 04:11:40 +0530
Subject: [PATCH 06/10] Update smileih5.py

The doc-stings were added and the code is tested using run-tests.py
---
 postpic/datareader/smileih5.py | 163 +++++++++++++++++++--------------
 1 file changed, 93 insertions(+), 70 deletions(-)

diff --git a/postpic/datareader/smileih5.py b/postpic/datareader/smileih5.py
index 15a27d5..87416b2 100644
--- a/postpic/datareader/smileih5.py
+++ b/postpic/datareader/smileih5.py
@@ -75,6 +75,7 @@ def _getindexfile(path):
         _generateh5indexfile(indexfile, h5files)
     return indexfile
 
+
 class SmileiReader(OpenPMDreader):
     '''
     The Smilei reader can read  a single h5 file or combine the information of
@@ -120,18 +121,22 @@ def __init__(self, h5file, iteration=None):
             raise IOError("Iteration {} is in valid".format(iteration))
         else:
             self._iteration = int(iteration)
-        
+
         self._data = self._h5['/data/{:010d}/'.format(self._iteration)]
         self.attrs = self._data.attrs
 
-    #To return the number of AM modes and Field names in the dump.
-    def getAMmodes(self):    
-        strings=np.array(self._data)
-        mask =  np.array(["_mode_" in s for s in strings])
-        arr=strings[mask]
+    def getAMmodes(self):
+        '''
+        This method is used to get the prefix of field names and
+        No.of AM modes available in the dump.
+        And it works only for AM mode technique.
+        '''
+        strings = np.array(self._data)
+        mask = np.array(["_mode_" in s for s in strings])
+        arr = strings[mask]
         max_suffix = float('-inf')
         max_suffix_string = None
-        prefix_list=[]
+        prefix_list = []
 
         for i in arr:
             prefix, suffix = i.split('_mode_')
@@ -141,46 +146,61 @@ def getAMmodes(self):
                 max_suffix_string = i
             prefix_list.append(prefix)
 
-        availModes=[i for i in range(0,int(max_suffix_string[-1])+1)]
-        return [prefix_list, int(max_suffix_string[-1])+1, availModes]
-    
+        availModes = [i for i in range(0, int(max_suffix_string[-1])+1)]
+
+        # [field names prefix, available AM modes]
+        return [prefix_list, availModes]
+
     @staticmethod
-    def _modeexpansion_naiv(rawdata, theta=0, M=None):
+    def _modeexpansion_naiv(rawdata, theta=0):
         '''
-        converts to radial data using `modeexpansion`, possibly for multiple
-        theta at once.
+        This method performes mode expansion of the raw data (an array consisting of complex
+        numbers) for both single and multiple theta vaues.
+
+        The output array has the shape (No.of theta, Nx, Nr)
+
+        Args:
+            rawdata : numpy array
+            The elements of this array are complex numbers, this got structured from the
+            raw data dumped in h5 file through getExpanded(key, theta) function.
+
+            theta : float/integer OR list of floats/integer
+
+        Output F_total is an array of real numbers which has shape (Np.of theta, Nx, Nr),
+        this F_total is the real value summation of the fourier series.
         '''
         if np.asarray(theta).shape is ():
             # single theta
             theta = [theta]
-        # multiple theta
-        #data = np.asarray([SmileiReader._modeexpansion_naiv_single(rawdata, theta=t)
-        #                   for t in theta])
-        # switch from (theta, r, z) to (r, theta, z)
-        
-        Nt=len(theta)
-        (Nm, Nr, Nx) = rawdata.shape
-        F_total = np.zeros((Nt,Nr,Nx))                       
-        mode =[m for m in range(0,Nm)]
 
-        if M is not None:
-            if Nm ==1 and M not in SmileiReader.getAMmodes[-1]:
-                raise IOError('Mode {} not available in dump'.format(M))
+        Nt = len(theta)
+        (Nm, Nr, Nx) = rawdata.shape
+        F_total = np.zeros((Nt, Nr, Nx))
+        mode = [m for m in range(0, Nm)]
 
         for t in theta:
-            index=theta.index(t)
-            if M is None:
-                for m in mode:
-                    F_total[index] += np.real(rawdata[m])*np.cos(m*t)+np.imag(rawdata[m])*np.sin(m*t)
-            else:
-                F_total[index] += np.real(rawdata[0])*np.cos(M*t)+np.imag(rawdata[0])*np.sin(M*t)
-        #F_total = F_total.swapaxes(0, 1)
+            index = theta.index(t)
+
+            for m in mode:
+                F_total[index] += np.real(rawdata[m])*np.cos(m*t)
+                +np.imag(rawdata[m])*np.sin(m*t)
+
         return F_total
-    
-    def getExpanded(self, key,theta=0):
-        
-        array_list=[]
-        modes=self.getAMmodes()[-1]
+
+    def getExpanded(self, key, theta=0):
+        '''
+        getExpanded() method converts the raw data real number array from h5file into
+        a complex number array (This convertion is important while performing mode expansion)
+        and finally returns the mode expanded array.
+
+        This method takes input from the h5files dump which has the following format,
+        field array = [[real_1,imag_1,real_2,image_2,.....],...]
+        The shape of this field array is (Nx, 2x Nr)
+        After real->complex conversion, the field array shape takes the form (Nx, Nr)
+        This final array is fed into _modeexpansion_naiv method.
+        '''
+        array_list = []
+        modes = self.getAMmodes()[-1]
         for mode in modes:
 
             field_name = key+"_mode_"+str(mode)
@@ -189,30 +209,33 @@ def getExpanded(self, key,theta=0):
             reshaped_array = field_array.reshape(field_array_shape[0], field_array_shape[1]//2, 2)
             complex_array = reshaped_array[:, :, 0] + 1j * reshaped_array[:, :, 1]
             array_list.append(complex_array)
-        
-        mod_complex_data= np.stack(array_list,axis=0)                    #Modified array of shape (Nmodes, Nx, Nr)
-        #mod_complex_data = np.transpose(mod_complex_data, (0, 2, 1))     #Modified array of shape (Nmodes, Nr, Nx)
-        return SmileiReader._modeexpansion_naiv(rawdata=mod_complex_data,theta=theta)
-    
-    def check_pattern(self,key):
-        field_pattern=re.compile(r'^{}.*_mode_'.format(key))
-        match = [True if re.match(field_pattern, s) else False for s in list(self._data)]
-
-        return match
-    
+        # Modified array of shape (Nmodes, Nx, Nr)
+        mod_complex_data = np.stack(array_list, axis=0)
+
+        return SmileiReader._modeexpansion_naiv(rawdata=mod_complex_data, theta=theta)
+
     def data(self, key, **kwargs):
         '''
         should work with any key, that contains data, thus on every hdf5.Dataset,
         but not on hdf5.Group. Will extract the data, convert it to SI and return it
         as a numpy array. Constant records will be detected and converted to
         a numpy array containing a single value only.
+
+        If the key is in AM mode dump, then it performs the mode expansion.
+        The theta values for which we need to perform
+        mode expansion can be given as keyword args.
+
+        Example:
+            Data = data(key='El', theta=[0,pi/2,pi])
+            Now the Data array has the shape (3, Nx, Nr)
         '''
 
+        # checking whether the key is in AM mode dump
         if key in self.getAMmodes()[0]:
             return self.getExpanded(key=key, **kwargs)
         else:
             record = self._data[key]
-        
+
         if "value" in record.attrs:
             # constant data (a single int or float)
             ret = np.float64(record.attrs['value']) * record.attrs['unitSI']
@@ -220,7 +243,7 @@ def data(self, key, **kwargs):
             # array data
             ret = np.float64(record[()]) * record.attrs['unitSI']
         return ret
-        
+
     def _keyE(self, component, **kwargs):
         # Smilei deviates from openPMD standard: Ex instead of E/x
         axsuffix = {0: 'x', 1: 'y', 2: 'z', 90: 'r', 91: 't'}[helper.axesidentify[component]]
@@ -233,9 +256,9 @@ def _keyB(self, component, **kwargs):
 
     def _simgridkeys(self):
         # Smilei deviates from openPMD standard: Ex instead of E/x
-        return ['Ex', 'Ey', 'Ez','Er','El','Et',
-                'Bx', 'By', 'Bz','Br','Bl','Bt',
-                'Jx','Jy','Jz','Jr','Jl','Jt','Rho']
+        return ['Ex', 'Ey', 'Ez', 'Er', 'El', 'Et',
+                'Bx', 'By', 'Bz', 'Br', 'Bl', 'Bt',
+                'Jx', 'Jy', 'Jz', 'Jr', 'Jl', 'Jt', 'Rho']
 
     def getderived(self):
         '''
@@ -250,43 +273,43 @@ def getderived(self):
     def __str__(self):
         return '<SmileiReader at "{}" at iteration {:d}>'.format(self.dumpidentifier,
                                                                  self.timestep())
-    
-    # override inherited method to count points after mode expansion
+
+    # To get the offsets of the grid.
     def gridoffset(self, key, axis):
         axid = helper.axesidentify[axis]
-        
+
         if axid == 91:  # theta
             return 0
-        elif key in self.getAMmodes()[0] and axid in [0,90]:
-            key="{}_mode_0".format(key)
+        elif key in self.getAMmodes()[0] and axid in [0, 90]:
+            key = "{}_mode_0".format(key)
             axid = int(axid/90)
             return super(SmileiReader, self).gridoffset(key=key, axis=axid)
         else:
             return super(SmileiReader, self).gridoffset(key, axis)
 
-    # override inherited method to count points after mode expansion
+    # To get the grid spacing.
     def gridspacing(self, key, axis):
         axid = helper.axesidentify[axis]
-        
-        if key in self.getAMmodes()[0] and axid in [0,90]:
-            key="{}_mode_0".format(key)
+
+        if key in self.getAMmodes()[0] and axid in [0, 90]:
+            key = "{}_mode_0".format(key)
             axid = int(axid/90)
             return super(SmileiReader, self).gridspacing(key=key, axis=axid)
         else:
             return super(SmileiReader, self).gridspacing(key, axis)
-    
-    # override inherited method to count points after mode expansion
+
+    # To get the grid points
     def gridpoints(self, key, axis):
         axid = helper.axesidentify[axis]
-        
+
         if key in self.getAMmodes()[0]:
-            key="{}_mode_0".format(key)
+            key = "{}_mode_0".format(key)
             (Nx, Nr) = self._data[key].shape
-            Nr=Nr/2
-            axid =int(axid/90)
-            return (Nx,Nr)[axid]
+            Nr = Nr/2
+            axid = int(axid/90)
+            return (Nx, Nr)[axid]
         else:
-            return super(SmileiReader, self).gridpoints(key=key,axis=axis)
+            return super(SmileiReader, self).gridpoints(key=key, axis=axis)
 
 
 class SmileiSeries(Simulationreader_ifc):

From 4db59d34a32bfdbb6dbcb69c03fe41d9c601ecb1 Mon Sep 17 00:00:00 2001
From: Samson J Kennedy <vinith11022003@gmail.com>
Date: Wed, 24 Jan 2024 03:18:04 +0530
Subject: [PATCH 07/10] Updated the author names and # level 0 1 2 comments

---
 postpic/datareader/smileih5.py | 138 +++++++++++++++++----------------
 1 file changed, 73 insertions(+), 65 deletions(-)

diff --git a/postpic/datareader/smileih5.py b/postpic/datareader/smileih5.py
index 87416b2..1d81cec 100644
--- a/postpic/datareader/smileih5.py
+++ b/postpic/datareader/smileih5.py
@@ -16,6 +16,7 @@
 #
 # Stephan Kuschel, 2023
 # Carolin Goll, 2023
+# Vinith Samson J, 2024
 
 from __future__ import absolute_import, division, print_function, unicode_literals
 
@@ -125,32 +126,6 @@ def __init__(self, h5file, iteration=None):
         self._data = self._h5['/data/{:010d}/'.format(self._iteration)]
         self.attrs = self._data.attrs
 
-    def getAMmodes(self):
-        '''
-        This method is used to get the prefix of field names and
-        No.of AM modes available in the dump.
-        And it works only for AM mode technique.
-        '''
-        strings = np.array(self._data)
-        mask = np.array(["_mode_" in s for s in strings])
-        arr = strings[mask]
-        max_suffix = float('-inf')
-        max_suffix_string = None
-        prefix_list = []
-
-        for i in arr:
-            prefix, suffix = i.split('_mode_')
-            suffix_int = int(suffix)
-            if suffix_int > max_suffix:
-                max_suffix = suffix_int
-                max_suffix_string = i
-            prefix_list.append(prefix)
-
-        availModes = [i for i in range(0, int(max_suffix_string[-1])+1)]
-
-        # [field names prefix, available AM modes]
-        return [prefix_list, availModes]
-
     @staticmethod
     def _modeexpansion_naiv(rawdata, theta=0):
         '''
@@ -187,9 +162,39 @@ def _modeexpansion_naiv(rawdata, theta=0):
 
         return F_total
 
-    def getExpanded(self, key, theta=0):
+# --- Level 0 methods ---
+
+    def _listAMmodes(self):
         '''
-        getExpanded() method converts the raw data real number array from h5file into
+        This method is used to get the list of
+        [prefix of field names, No.of AM modes] available in the dump.
+        And it works only for AM mode technique.
+        '''
+        strings = np.array(self._data)
+        mask = np.array(["_mode_" in s for s in strings])
+        arr = strings[mask]
+        max_suffix = float('-inf')
+        max_suffix_string = None
+        prefix_list = []
+
+        for i in arr:
+            prefix, suffix = i.split('_mode_')
+            suffix_int = int(suffix)
+            if suffix_int > max_suffix:
+                max_suffix = suffix_int
+                max_suffix_string = i
+            prefix_list.append(prefix)
+
+        availModes = [i for i in range(0, int(max_suffix_string[-1])+1)]
+
+        # [field names prefix, available AM modes]
+        return [prefix_list, availModes]
+
+# --- Level 1 methods ---
+
+    def _getExpanded(self, key, theta=0):
+        '''
+        _getExpanded() method converts the raw data real number array from h5file into
         a complex number array (This convertion is important while performing mode expansion)
         and finally returns the mode expanded array.
 
@@ -200,7 +205,7 @@ def getExpanded(self, key, theta=0):
         This final array is fed into _modeexpansion_naiv method.
         '''
         array_list = []
-        modes = self.getAMmodes()[-1]
+        modes = self._listAMmodes()[-1]
         for mode in modes:
 
             field_name = key+"_mode_"+str(mode)
@@ -209,6 +214,7 @@ def getExpanded(self, key, theta=0):
             reshaped_array = field_array.reshape(field_array_shape[0], field_array_shape[1]//2, 2)
             complex_array = reshaped_array[:, :, 0] + 1j * reshaped_array[:, :, 1]
             array_list.append(complex_array)
+
         # Modified array of shape (Nmodes, Nx, Nr)
         mod_complex_data = np.stack(array_list, axis=0)
 
@@ -227,12 +233,12 @@ def data(self, key, **kwargs):
 
         Example:
             Data = data(key='El', theta=[0,pi/2,pi])
-            Now the Data array has the shape (3, Nx, Nr)
+            Now the Data will array have the shape (3, Nx, Nr)
         '''
 
         # checking whether the key is in AM mode dump
-        if key in self.getAMmodes()[0]:
-            return self.getExpanded(key=key, **kwargs)
+        if key in self._listAMmodes()[0]:
+            return self._getExpanded(key=key, **kwargs)
         else:
             record = self._data[key]
 
@@ -244,43 +250,13 @@ def data(self, key, **kwargs):
             ret = np.float64(record[()]) * record.attrs['unitSI']
         return ret
 
-    def _keyE(self, component, **kwargs):
-        # Smilei deviates from openPMD standard: Ex instead of E/x
-        axsuffix = {0: 'x', 1: 'y', 2: 'z', 90: 'r', 91: 't'}[helper.axesidentify[component]]
-        return 'E{}'.format(axsuffix)
-
-    def _keyB(self, component, **kwargs):
-        # Smilei deviates from openPMD standard: Bx instead of B/x
-        axsuffix = {0: 'x', 1: 'y', 2: 'z', 90: 'r', 91: 't'}[helper.axesidentify[component]]
-        return 'B{}'.format(axsuffix)
-
-    def _simgridkeys(self):
-        # Smilei deviates from openPMD standard: Ex instead of E/x
-        return ['Ex', 'Ey', 'Ez', 'Er', 'El', 'Et',
-                'Bx', 'By', 'Bz', 'Br', 'Bl', 'Bt',
-                'Jx', 'Jy', 'Jz', 'Jr', 'Jl', 'Jt', 'Rho']
-
-    def getderived(self):
-        '''
-        return all other fields dumped, except E and B.
-        '''
-        ret = []
-        self['.'].visit(ret.append)
-        # TODO: remove E and B fields and particles from list
-        ret.sort()
-        return ret
-
-    def __str__(self):
-        return '<SmileiReader at "{}" at iteration {:d}>'.format(self.dumpidentifier,
-                                                                 self.timestep())
-
     # To get the offsets of the grid.
     def gridoffset(self, key, axis):
         axid = helper.axesidentify[axis]
 
         if axid == 91:  # theta
             return 0
-        elif key in self.getAMmodes()[0] and axid in [0, 90]:
+        elif key in self._listAMmodes()[0] and axid in [0, 90]:
             key = "{}_mode_0".format(key)
             axid = int(axid/90)
             return super(SmileiReader, self).gridoffset(key=key, axis=axid)
@@ -291,7 +267,7 @@ def gridoffset(self, key, axis):
     def gridspacing(self, key, axis):
         axid = helper.axesidentify[axis]
 
-        if key in self.getAMmodes()[0] and axid in [0, 90]:
+        if key in self._listAMmodes()[0] and axid in [0, 90]:
             key = "{}_mode_0".format(key)
             axid = int(axid/90)
             return super(SmileiReader, self).gridspacing(key=key, axis=axid)
@@ -302,7 +278,7 @@ def gridspacing(self, key, axis):
     def gridpoints(self, key, axis):
         axid = helper.axesidentify[axis]
 
-        if key in self.getAMmodes()[0]:
+        if key in self._listAMmodes()[0]:
             key = "{}_mode_0".format(key)
             (Nx, Nr) = self._data[key].shape
             Nr = Nr/2
@@ -311,6 +287,38 @@ def gridpoints(self, key, axis):
         else:
             return super(SmileiReader, self).gridpoints(key=key, axis=axis)
 
+# --- Level 2 methods ---
+
+    def _keyE(self, component, **kwargs):
+        # Smilei deviates from openPMD standard: Ex instead of E/x
+        axsuffix = {0: 'x', 1: 'y', 2: 'z', 90: 'r', 91: 't'}[helper.axesidentify[component]]
+        return 'E{}'.format(axsuffix)
+
+    def _keyB(self, component, **kwargs):
+        # Smilei deviates from openPMD standard: Bx instead of B/x
+        axsuffix = {0: 'x', 1: 'y', 2: 'z', 90: 'r', 91: 't'}[helper.axesidentify[component]]
+        return 'B{}'.format(axsuffix)
+
+    def _simgridkeys(self):
+        # Smilei deviates from openPMD standard: Ex instead of E/x
+        return ['Ex', 'Ey', 'Ez', 'Er', 'El', 'Et',
+                'Bx', 'By', 'Bz', 'Br', 'Bl', 'Bt',
+                'Jx', 'Jy', 'Jz', 'Jr', 'Jl', 'Jt', 'Rho']
+
+    def getderived(self):
+        '''
+        return all other fields dumped, except E and B.
+        '''
+        ret = []
+        self['.'].visit(ret.append)
+        # TODO: remove E and B fields and particles from list
+        ret.sort()
+        return ret
+
+    def __str__(self):
+        return '<SmileiReader at "{}" at iteration {:d}>'.format(self.dumpidentifier,
+                                                                 self.timestep())
+
 
 class SmileiSeries(Simulationreader_ifc):
     '''

From bb405ddebc5fabba10e298f96cb191e5b3f75826 Mon Sep 17 00:00:00 2001
From: Samson J Kennedy <vinith11022003@gmail.com>
Date: Wed, 24 Jan 2024 03:19:32 +0530
Subject: [PATCH 08/10] Updated the latest version details and features
 implimented

---
 CHANGELOG.md | 27 +++++++++++++++++++++++++++
 1 file changed, 27 insertions(+)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 4c86801..e341faf 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -4,6 +4,33 @@ Changelog of postpic
 current master
 --------------
 
+v0.6
+----
+
+2023-01-23
+
+This is the first version which is compatible with `Smilei` PIC (https://smileipic.github.io/Smilei/) data format. The addon datareader code is `smileih5.py` (https://github.com/skuschel/postpic/blob/dev/postpic/datareader/smileih5.py). Might be slower in runtime compared to `happi` module, a default data processing module for smilei (https://github.com/SmileiPIC/Smilei/tree/master/happi). But soon will be optimised.
+
+**Highlights**
+
+* Quick and efficient way of importing h5file using index file method, so that multiple .h5 files can be imported at the same time.
+* Added support for files written by the smilei `version v4.8` (https://github.com/SmileiPIC/Smilei/releases/tag/v4.8).
+* Azimuthal Mode decomposition technique is implimented.
+* Can extract field data from both Cartesian and Azimuthal mode decomposed fields.
+
+**Incompatible adjustments to previous version**
+
+* The `data()`, `gridoffset()`, `gridpoints()`, `gridspacing()` and `_modeexpansion_naiv()` methods were altered from `openPMDh5.py` standards format due to different data dumping format of SmileiPIC.
+* Since the Azimuthal mode field data dump format were different from FBpic dump format, a new method `_getExpanded()` is introduced to structure the data dumped by SmileiPIC into an array of complex numbers which is then fed into `_modeexpansion_naiv()` method.
+
+**Other improvements and new features**
+
+* The fields (`both cartesian and Azimuthal mode`) can be accessed by method `data(key, **kwargs)`, the keys are in the format `El`, `Bx`, `Rho` etc.,
+* The species data also can be accessed by `x`, `Px`, `q`, `w` etc.,
+* Azimuthal Mode decomposition technique is implimented, by mentioning theta (either single value or list of angles) as kwargs in `data` method.
+* Gridoffset, gridpoints and gridspacing data can be accessed using `gridoffset(key, axis)`, `gridpoints(key, axis)` and `gridspacing(key, axis)` methods respectively.
+* List of Azimuthal fields names and Modes available in the dump can be obtained using `_listAMmodes()` method.
+
 
 v0.5
 ----

From 4c777cd2c804d36c2d7c48a50e8386a07ac5cc71 Mon Sep 17 00:00:00 2001
From: Samson J Kennedy <vinith11022003@gmail.com>
Date: Wed, 24 Jan 2024 04:11:11 +0530
Subject: [PATCH 09/10] Update CHANGELOG.md

---
 CHANGELOG.md | 22 ++++------------------
 1 file changed, 4 insertions(+), 18 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index e341faf..3330c1c 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -4,32 +4,18 @@ Changelog of postpic
 current master
 --------------
 
-v0.6
-----
-
-2023-01-23
-
-This is the first version which is compatible with `Smilei` PIC (https://smileipic.github.io/Smilei/) data format. The addon datareader code is `smileih5.py` (https://github.com/skuschel/postpic/blob/dev/postpic/datareader/smileih5.py). Might be slower in runtime compared to `happi` module, a default data processing module for smilei (https://github.com/SmileiPIC/Smilei/tree/master/happi). But soon will be optimised.
 
 **Highlights**
 
-* Quick and efficient way of importing h5file using index file method, so that multiple .h5 files can be imported at the same time.
-* Added support for files written by the smilei `version v4.8` (https://github.com/SmileiPIC/Smilei/releases/tag/v4.8).
-* Azimuthal Mode decomposition technique is implimented.
-* Can extract field data from both Cartesian and Azimuthal mode decomposed fields.
+* Add support to read the `Smilei` PIC (https://smileipic.github.io/Smilei/) data format in both cartesian and azimuthal geometry. Postpic uses a build in azimuthal mode expansion very similar to the one used for fbpic.
+* To read smilei data, postpic only relies on the hdf5 package and not smilei's happi module for data access. Paricle ID's (ParticleTracking as described by smilei) can be read directly from the hdf5. Happi requires to sort the IDs and write a new hdf5, which can be twice as big as the original dumps. Using postpic's access this step will be skipped and thus access is much faster (but by default with unordered particle IDs as in any other code).
+
 
 **Incompatible adjustments to previous version**
 
-* The `data()`, `gridoffset()`, `gridpoints()`, `gridspacing()` and `_modeexpansion_naiv()` methods were altered from `openPMDh5.py` standards format due to different data dumping format of SmileiPIC.
-* Since the Azimuthal mode field data dump format were different from FBpic dump format, a new method `_getExpanded()` is introduced to structure the data dumped by SmileiPIC into an array of complex numbers which is then fed into `_modeexpansion_naiv()` method.
 
-**Other improvements and new features**
 
-* The fields (`both cartesian and Azimuthal mode`) can be accessed by method `data(key, **kwargs)`, the keys are in the format `El`, `Bx`, `Rho` etc.,
-* The species data also can be accessed by `x`, `Px`, `q`, `w` etc.,
-* Azimuthal Mode decomposition technique is implimented, by mentioning theta (either single value or list of angles) as kwargs in `data` method.
-* Gridoffset, gridpoints and gridspacing data can be accessed using `gridoffset(key, axis)`, `gridpoints(key, axis)` and `gridspacing(key, axis)` methods respectively.
-* List of Azimuthal fields names and Modes available in the dump can be obtained using `_listAMmodes()` method.
+**Other improvements and new features**
 
 
 v0.5

From a1c763f3048157883087ff3a69cd4871b06fdcf9 Mon Sep 17 00:00:00 2001
From: Samson J Kennedy <vinith11022003@gmail.com>
Date: Fri, 26 Jan 2024 23:08:45 +0530
Subject: [PATCH 10/10] Percentage Error fixed

---
 postpic/datareader/smileih5.py | 25 ++++++++++++-------------
 1 file changed, 12 insertions(+), 13 deletions(-)

diff --git a/postpic/datareader/smileih5.py b/postpic/datareader/smileih5.py
index 1d81cec..f7b407b 100644
--- a/postpic/datareader/smileih5.py
+++ b/postpic/datareader/smileih5.py
@@ -144,23 +144,22 @@ def _modeexpansion_naiv(rawdata, theta=0):
         Output F_total is an array of real numbers which has shape (Np.of theta, Nx, Nr),
         this F_total is the real value summation of the fourier series.
         '''
-        if np.asarray(theta).shape is ():
-            # single theta
+        if np.array(theta).shape is ():
             theta = [theta]
 
-        Nt = len(theta)
-        (Nm, Nr, Nx) = rawdata.shape
-        F_total = np.zeros((Nt, Nr, Nx))
-        mode = [m for m in range(0, Nm)]
+        array_list = []
 
+        (Nm, Nx, Nr) = rawdata.shape
+        F_total = np.zeros((Nx, Nr))
+        mode = [m for m in range(0, Nm)]
         for t in theta:
-            index = theta.index(t)
 
             for m in mode:
-                F_total[index] += np.real(rawdata[m])*np.cos(m*t)
-                +np.imag(rawdata[m])*np.sin(m*t)
+                F_total += np.real(rawdata[m])*np.cos(m*t)+np.imag(rawdata[m])*np.sin(m*t)
+            array_list.append(F_total)
 
-        return F_total
+        mod_F_total = np.stack(array_list, axis=0)
+        return mod_F_total
 
 # --- Level 0 methods ---
 
@@ -209,7 +208,7 @@ def _getExpanded(self, key, theta=0):
         for mode in modes:
 
             field_name = key+"_mode_"+str(mode)
-            field_array = np.array(self._data[field_name])*(self._data[field_name].attrs['unitSI'])
+            field_array = np.array(self._data[field_name])
             field_array_shape = field_array.shape
             reshaped_array = field_array.reshape(field_array_shape[0], field_array_shape[1]//2, 2)
             complex_array = reshaped_array[:, :, 0] + 1j * reshaped_array[:, :, 1]
@@ -217,8 +216,8 @@ def _getExpanded(self, key, theta=0):
 
         # Modified array of shape (Nmodes, Nx, Nr)
         mod_complex_data = np.stack(array_list, axis=0)
-
-        return SmileiReader._modeexpansion_naiv(rawdata=mod_complex_data, theta=theta)
+        factor = self._data["{}_mode_0".format(key)].attrs['unitSI']
+        return SmileiReader._modeexpansion_naiv(rawdata=mod_complex_data, theta=theta)*factor
 
     def data(self, key, **kwargs):
         '''