Merge pull request #10 from skuschel/vsimreader

Vsimreader

.gitignore

@@ -1,5 +1,6 @@
 *~
 *.pyc
+*.odt
 .project
 .pydevproject
 *.egg-info

postpic/_const.py

@@ -28,7 +28,9 @@
                        'PY': 4, 'Py': 4, 'py': 4, 4: 4,
                        'PZ': 5, 'Pz': 5, 'pz': 5, 9: 9,
                        'weight': 9, 'w': 9, 10: 10,
-                       'id': 10, 'ID': 10})
+                       'id': 10, 'ID': 10,
+                       'mass': 11, 'm': 11, 'Mass': 11,
+                       'charge': 12, 'c': 12, 'Charge': 12, 'q': 12})
 
 # Some static functions
 

postpic/analyzer/analyzer.py

@@ -89,8 +89,8 @@ class SpeciesIdentifier(PhysicalConstants):
                   'gold20': True}
 
     #  unit: amu
-    _masslistelement = {'H': 1, 'He': 4, 'C': 12, 'N': 14, 'O': 16, 'F': 19,
-                        'Ne': 20.2, 'Al': 27, 'Si': 28, 'Ar': 40, 'Au': 197}
+    _masslistelement = {'H': 1, 'He': 4, 'Li': 6.9, 'C': 12, 'N': 14, 'O': 16, 'F': 19,
+                        'Ne': 20.2, 'Al': 27, 'Si': 28, 'Ar': 40, 'Rb': 85.5, 'Au': 197}
 
     @staticmethod
     def isejected(species):
@@ -105,7 +105,8 @@ def ision(cls, species):
     @classmethod
     def identifyspecies(cls, species):
         """
-        Returns a dictionary contining particle informations.
+        Returns a dictionary containing particle informations deduced from
+        the species name.
         The following keys in the dictionary will always be present:
         name   species name string
         mass    kg (SI)
@@ -125,15 +126,20 @@ def identifyspecies(cls, species):
 
         # Regex for parsing ion species name.
         # See docsting for valid examples
-        regex = '(?P<prae>(.*_)*)(?P<name>(ionc(?P<c1>\d+)m(?P<m2>\d+)|' \
-                'ionm(?P<m1>\d+)c(?P<c2>\d+))|(?P<electron>[Ee]le[ck]tron)|' \
-                '(?P<elem>[A-Za-z]+)(?P<elem_c>\d*))(?P<suffix>[a-z]*)'
+        regex = '(?P<prae>(.*_)*)' \
+                '((?P<elem>((?!El)[A-Z][a-z]?))(?P<elem_c>\d*)|' \
+                '(?P<ionmc>(ionc(?P<c1>\d+)m(?P<m2>\d+)|ionm(?P<m1>\d+)c(?P<c2>\d+)))' \
+                ')?' \
+                '(?P<plus>(Plus)*)' \
+                '(?P<electron>[Ee]le[ck])?' \
+                '(?P<suffix>\w*?)$'
         r = re.match(regex, s)
         if r is None:
             raise Exception('Species ' + str(s) +
                             ' does not match regex name pattern: ' +
                             str(regex))
         regexdict = r.groupdict()
+        # print(regexdict)
 
         # recognize anz prae and add dictionary key
         if regexdict['prae']:
@@ -142,44 +148,52 @@ def identifyspecies(cls, species):
                 if not key == '':
                     ret[key] = True
 
-        # Name Element + charge state: C1, C6, F2, F9, Au20, Pb34a
+        # Excluding patterns start here
+        # 1) Name Element + charge state: C1, C6, F2, F9, Au20, Pb34a
         if regexdict['elem']:
-            try:
-                ret['mass'] = float(cls._masslistelement[regexdict['elem']]) * \
-                    1836.2 * cls.me
+            ret['mass'] = float(cls._masslistelement[regexdict['elem']]) * \
+                1836.2 * cls.me
+            if regexdict['elem_c'] == '':
+                ret['charge'] = 0
+            else:
                 ret['charge'] = float(regexdict['elem_c']) * cls.qe
+            ret['ision'] = True
+        # 2) ionmc like
+        elif regexdict['ionmc']:
+            if regexdict['c1']:
+                ret['mass'] = float(regexdict['m2']) * 1836.2 * cls.me
+                ret['charge'] = float(regexdict['c1']) * cls.qe
+                ret['ision'] = True
+
+            if regexdict['c2']:
+                ret['mass'] = float(regexdict['m1']) * 1836.2 * cls.me
+                ret['charge'] = float(regexdict['c2']) * cls.qe
                 ret['ision'] = True
-            except KeyError:
-                # this pattern will also match, if name is defined in masslist,
-                # so just ignore if key is not found.
-                pass
 
+        # charge may be given via plus
+        if regexdict['plus']:
+            charge = len(regexdict['plus'])/4 * cls.qe
+            if ret['charge'] == 0:
+                ret['charge'] = charge
+
+        # Elektron can be appended to any ion name, so overwrite.
         if regexdict['electron']:
             ret['mass'] = cls.me
             ret['charge'] = -1 * cls.qe
             ret['ision'] = False
 
-        if regexdict['c1']:
-            ret['mass'] = float(regexdict['m2']) * 1836.2 * cls.me
-            ret['charge'] = float(regexdict['c1']) * cls.qe
-            ret['ision'] = True
-
-        if regexdict['c2']:
-            ret['mass'] = float(regexdict['m1']) * 1836.2 * cls.me
-            ret['charge'] = float(regexdict['c2']) * cls.qe
-            ret['ision'] = True
-
         # simply added to _masslist and _chargelist
         # this should not be used anymore
-        if regexdict['name'] in cls._masslist:
-            ret['mass'] = float(cls._masslist[regexdict['name']]) * cls.me
-        if regexdict['name'] in cls._chargelist:
-            ret['charge'] = float(cls._chargelist[regexdict['name']] * cls.qe)
-        if regexdict['name'] in cls._isionlist:
-            ret['ision'] = cls._isionlist[regexdict['name']]
+        # set only if property is not already set
+        if 'mass' not in ret and regexdict['suffix'] in cls._masslist:
+            ret['mass'] = float(cls._masslist[regexdict['suffix']]) * cls.me
+        if 'charge' not in ret and regexdict['suffix'] in cls._chargelist:
+            ret['charge'] = float(cls._chargelist[regexdict['suffix']] * cls.qe)
+        if 'ision' not in ret and regexdict['suffix'] in cls._isionlist:
+            ret['ision'] = cls._isionlist[regexdict['suffix']]
 
         if not (('mass' in ret) and ('charge' in ret) and ('ision' in ret)):
-            raise Error('species ' + species + ' not recognized.')
+            raise Exception('species ' + species + ' not recognized.')
 
         return ret
 

postpic/analyzer/particles.py

@@ -37,9 +37,12 @@ class _SingleSpeciesAnalyzer(object):
     def __init__(self, dumpreader, species):
         self.species = species
         self._dumpreader = dumpreader
-        self._idfy = identifyspecies(species)
-        self._mass = self._idfy['mass']  # SI
-        self._charge = self._idfy['charge']  # SI
+        self._mass = dumpreader.getSpecies(species, 'mass')
+        self._charge = dumpreader.getSpecies(species, 'charge')
+        if self._mass is None or self._charge is None:
+            self._idfy = identifyspecies(species)
+            self._mass = self._idfy['mass']  # SI
+            self._charge = self._idfy['charge']  # SI
         self.compresslog = []
         self._compressboollist = None
         self._cache = {}
@@ -55,8 +58,8 @@ def __getitem__(self, key):
         else:
             ret = self._dumpreader.getSpecies(self.species, key)
             if ret is not None:
-                ret = np.asfarray(ret, dtype='float64')
-                if self._compressboollist is not None:
+                ret = np.float64(ret)
+                if not isinstance(ret, float) and self._compressboollist is not None:
                     ret = ret[self._compressboollist]  # avoid executing this line too often.
                     self._cache[key] = ret
                     # if memomry is low, caching could be skipped entirely.
@@ -119,17 +122,30 @@ def uncompress(self):
     # --- Only very basic functions
 
     def __len__(self):  # = number of particles
-        w = self['weight']
-        if w is None:
-            ret = 0
-        else:
-            ret = len(w)
+        # find a valid dataset to count number of paricles
+        if self._compressboollist is not None:
+            return len(self._compressboollist)
+        for key in ['weight', 'x', 'px', 'y', 'py', 'z', 'pz']:
+            data = self[key]
+            try:
+                # len(3) will yield a TypeError, len([3]) returns 1
+                ret = len(data)
+                break
+            except(TypeError):
+                pass
         return ret
 
 # --- These functions are for practical use. Return None if not dumped.
 
     def weight(self):  # np.float64(np.array([4.3])) == 4.3 may cause error
-        return self['weight']
+        w = self['weight']
+        # on constant weight, self['weight'] may return a single scalar,
+        # which is converted to float by __getitem__
+        if isinstance(w, float):
+            ret = np.repeat(w, len(self))
+        else:
+            ret = w
+        return ret
 
     def mass(self):  # SI
         return np.repeat(self._mass, len(self))

postpic/datareader/__init__.py

@@ -314,13 +314,18 @@ def chooseCode(code):
     Args:
       code : string
         Possible options are:
-          - "EPOCH": .sdf files written by EPOCH1D, EPOCH2D or EPOCH3D.
           - "DUMMY": dummy class creating fake data.
+          - "EPOCH": .sdf files written by EPOCH1D, EPOCH2D or EPOCH3D.
+          - "VSIM": .hdf5 files written by VSim.
     '''
     if code in ['EPOCH', 'epoch', 'EPOCH1D', 'EPOCH2D', 'EPOCH3D']:
         from epochsdf import Sdfreader, Visitreader
         setdumpreadercls(Sdfreader)
         setsimreadercls(Visitreader)
+    elif code in ['VSim', 'VSIM', 'vsim']:
+        from vsimhdf5 import Hdf5reader, VSimReader
+        setdumpreadercls(Hdf5reader)
+        setsimreadercls(VSimReader)
     elif code in ['DUMMY', 'dummy']:
         from dummy import Dummyreader, Dummysim
         setdumpreadercls(Dummyreader)