qpms/argproc.py

   1 '''
   2 Common snippets for argument processing in command line scripts.
   3 '''
   4
   5 import argparse
   6 import sys
   7 import warnings
   8
   9 def flatten(S):
  10     if S == []:
  11         return S
  12     if isinstance(S[0], list):
  13         return flatten(S[0]) + flatten(S[1:])
  14     return S[:1] + flatten(S[1:])
  15
  16 def make_action_sharedlist(opname, listname):
  17     class opAction(argparse.Action):
  18         def __call__(self, parser, args, values, option_string=None):
  19             if (not hasattr(args, listname)) or getattr(args, listname) is None:
  20                 setattr(args, listname, list())
  21             getattr(args, listname).append((opname, values))
  22     return opAction
  23
  24 def make_dict_action(argtype=None, postaction='store', first_is_key=True):
  25     class DictAction(argparse.Action):
  26         #def __init__(self, option_strings, dest, nargs=None, **kwargs):
  27         #    if nargs is not None:
  28         #    raise ValueError("nargs not allowed")
  29         #    super(DictAction, self).__init__(option_strings, dest, **kwargs)
  30         def __call__(self, parser, namespace, values, option_string=None):
  31             if first_is_key: # For the labeled versions
  32                 key = values[0]
  33                 vals = values[1:]
  34             else: # For the default values
  35                 key = None
  36                 vals = values
  37             if argtype is not None:
  38                 if (first_is_key and self.nargs == 2) or (not first_is_key and self.nargs == 1):
  39                     vals = argtype(vals[0]) # avoid having lists in this case
  40                 else:
  41                     vals = [argtype(val) for val in vals]
  42             ledict = getattr(namespace, self.dest, {})
  43             if ledict is None:
  44                 ledict = {}
  45             if postaction=='store':
  46                 ledict[key] = vals
  47             elif postaction=='append':
  48                 lelist = ledict.get(key, [])
  49                 lelist.append(vals)
  50                 ledict[key] = lelist
  51             setattr(namespace, self.dest, ledict)
  52     return DictAction
  53
  54
  55 class ArgumentProcessingError(Exception):
  56     pass
  57
  58 class AppendTupleAction(argparse.Action):
  59     ''' A variation on the 'append' builtin action from argparse, but uses tuples for the internal groupings instead of lists '''
  60     def __call__(self, parser, args, values, option_string=None):
  61         if (not hasattr(args, self.dest)) or getattr(args, self.dest) is None:
  62             setattr(args, self.dest, list())
  63         getattr(args, self.dest).append(tuple(values))
  64
  65 def float_range(string):
  66     """Tries to parse a string either as one individual float value
  67     or one of the following patterns:
  68
  69     first:last:increment
  70     first:last|steps
  71     first:last
  72
  73     (The last one is equivalent to first:last|50.)
  74     Returns either float or numpy array.
  75     """
  76     try:
  77         res = float(string)
  78         return res
  79     except ValueError:
  80         import re
  81         steps = None
  82         match = re.match(r's?([^:]+):([^|]+)\|(.+)', string)
  83         if match:
  84             steps = int(match.group(3))
  85         else:
  86             match = re.match(r's?([^:]+):([^:]+):(.+)', string)
  87             if match:
  88                 increment = float(match.group(3))
  89             else:
  90                 match = re.match(r's?([^:]+):(.+)', string)
  91                 if match:
  92                     steps = 50
  93                 else:
  94                     argparse.ArgumentTypeError('Invalid float/sequence format: "%s"' % string)
  95         first = float(match.group(1))
  96         last = float(match.group(2))
  97         import numpy as np
  98         if steps is not None:
  99             return np.linspace(first, last, num=steps)
 100         else:
 101             return np.arange(first, last, increment)
 102
 103 def int_or_None(string):
 104     """Tries to parse a string either as an int or None (if it contains only whitespaces)"""
 105     try:
 106         return int(string)
 107     except ValueError as ve:
 108         if string.strip() == '':
 109             return None
 110         else:
 111             raise ve
 112
 113 def sslice(string):
 114     """Tries to parse a string either as one individual int value
 115     or one of the following patterns:
 116
 117     first:last:increment
 118     first:last
 119
 120     first, last and increment must be parseable as ints
 121     or be empty (then
 122
 123     In each case, 's' letter can be prepended to the whole string to avoid
 124     argparse interpreting this as a new option (if the argument contains
 125     '-' or '+').
 126
 127     Returns either int or slice containing ints or Nones.
 128     """
 129     if string[0] == 's':
 130         string = string[1:]
 131     try:
 132         res = int(string)
 133         return res
 134     except ValueError:
 135         import re
 136         match = re.match(r'([^:]*):([^:]*):(.*)', string)
 137         if match:
 138             step = int_or_None(match.group(3))
 139         else:
 140             match = re.match(r'([^:]*):(.*)', string)
 141             if match:
 142                 step = None
 143             else:
 144                 argparse.ArgumentTypeError('Invalid int/slice format: "%s"' % string)
 145         start = int_or_None(match.group(1))
 146         stop = int_or_None(match.group(2))
 147         return slice(start, stop, step)
 148
 149 def sfloat(string):
 150     '''Tries to match a float, or a float with prepended 's'
 151
 152     Used as a workaraound for argparse's negative number matcher, which does not recognize
 153     scientific notation.
 154     '''
 155     try:
 156         res = float(string)
 157     except ValueError as exc:
 158         if string[0] == 's':
 159             res = float(string[1:])
 160         else: raise exc
 161     return res
 162
 163 def sint(string):
 164     '''Tries to match an int, or an int with prepended 's'
 165
 166     Used as a workaraound for argparse's negative number matcher if '+' is used as a
 167     prefix
 168     '''
 169     try:
 170         res = int(string)
 171     except ValueError as exc:
 172         if string[0] == 's':
 173             res = int(string[1:])
 174         else: raise exc
 175     return res
 176
 177 def material_spec(string):
 178     """Tries to parse a string as a material specification, i.e. a
 179     real or complex number or one of the string in built-in Lorentz-Drude models.
 180
 181     Tries to interpret the string as 1) float, 2) complex, 3) Lorentz-Drude key.
 182     Raises argparse.ArgumentTypeError on failure.
 183     """
 184     from .cymaterials import lorentz_drude
 185     if string in lorentz_drude.keys():
 186         return string
 187     else:
 188         try: lemat = float(string)
 189         except ValueError:
 190             try: lemat = complex(string)
 191             except ValueError as ve:
 192                 raise argparse.ArgumentTypeError("Material specification must be a supported material name %s, or a number" % (str(lorentz_drude.keys()),)) from ve
 193         return lemat
 194
 195 class ArgParser:
 196     ''' Common argument parsing engine for QPMS python CLI scripts. '''
 197
 198     def __add_planewave_argparse_group(ap):
 199         pwgrp = ap.add_argument_group('Incident wave specification', """
 200         Incident wave direction is given in terms of ISO polar and azimuthal angles θ, φ,
 201         which translate into cartesian coordinates as r̂ = (x, y, z) = (sin(θ) cos(φ), sin(θ) sin(φ), cos(θ)).
 202
 203         Wave polarisation is given in terms of parameters ψ, χ, where ψ is the angle between a polarisation
 204         ellipse axis and meridian tangent θ̂, and tg χ determines axes ratio;
 205         the electric field in the origin is then
 206
 207         E⃗ = cos(χ) (cos(ψ) θ̂ + sin(ψ) φ̂) + i sin(χ) (sin(ψ) θ̂ + cos(ψ) φ̂).
 208
 209         All the angles are given as multiples of π/2.
 210         """ # TODO EXAMPLES
 211         )
 212         pwgrp.add_argument("-φ", "--phi", type=float, default=0,
 213             help='Incident wave asimuth in multiples of π/2.')
 214         pwgrp.add_argument("-θ", "--theta", type=float_range, default=0,
 215             help='Incident wave polar angle in multiples of π/2. This might be a sequence in format FIRST:LAST:INCREMENT.')
 216         pwgrp.add_argument("-ψ", "--psi", type=float, default=0,
 217             help='Angle between polarisation ellipse axis and meridian tangent θ̂ in multiples of π/2.')
 218         pwgrp.add_argument("-χ", "--chi", type=float, default=0,
 219             help='Polarisation parameter χ in multiples of π/2. 0 for linear, 0.5 for circular pol.')
 220
 221     def __add_manyparticle_argparse_group(ap):
 222         mpgrp = ap.add_argument_group('Many particle specification', "TODO DOC")
 223         mpgrp.add_argument("-p", "--position", nargs='+', action=make_dict_action(argtype=sfloat, postaction='append',
 224             first_is_key=False), help="Particle positions, cartesion coordinates (default particle properties)")
 225         mpgrp.add_argument("+p", "++position", nargs='+', action=make_dict_action(argtype=sfloat, postaction='append',
 226             first_is_key=True), help="Particle positions, cartesian coordinates (labeled)")
 227         mpgrp.add_argument("-L", "--lMax", nargs=1, default={},
 228                 action=make_dict_action(argtype=int, postaction='store', first_is_key=False,),
 229                 help="Cutoff multipole degree (default)")
 230         mpgrp.add_argument("+L", "++lMax", nargs=2,
 231                 action=make_dict_action(argtype=int, postaction='store', first_is_key=True,),
 232             help="Cutoff multipole degree (labeled)")
 233         mpgrp.add_argument("-m", "--material", nargs=1, default={},
 234                 action=make_dict_action(argtype=material_spec, postaction='store', first_is_key=False,),
 235                 help='particle material (Au, Ag, ... for Lorentz-Drude or number for constant refractive index) (default)')
 236         mpgrp.add_argument("+m", "++material", nargs=2,
 237                 action=make_dict_action(argtype=material_spec, postaction='store', first_is_key=True,),
 238                 help='particle material (Au, Ag, ... for Lorentz-Drude or number for constant refractive index) (labeled)')
 239         mpgrp.add_argument("-r", "--radius", nargs=1, default={},
 240                 action=make_dict_action(argtype=float, postaction='store', first_is_key=False,),
 241                 help='particle radius (sphere or cylinder; default)')
 242         mpgrp.add_argument("+r", "++radius", nargs=2,
 243                 action=make_dict_action(argtype=float, postaction='store', first_is_key=True,),
 244                 help='particle radius (sphere or cylinder; labeled)')
 245         mpgrp.add_argument("-H", "--height", nargs=1, default={},
 246                 action=make_dict_action(argtype=float, postaction='store', first_is_key=False,),
 247                 help='particle radius (cylinder; default)')
 248         mpgrp.add_argument("+H", "++height", nargs=2,
 249                 action=make_dict_action(argtype=float, postaction='store', first_is_key=True,),
 250                 help='particle radius (cylinder; labeled)')
 251
 252     atomic_arguments = {
 253             'rectlattice2d_periods': lambda ap: ap.add_argument("-p", "--period", type=float, nargs='+', required=True, help='square/rectangular lattice periods', metavar=('px','[py]')),
 254             'rectlattice2d_counts': lambda ap: ap.add_argument("--size", type=int, nargs=2, required=True, help='rectangular array size (particle column, row count)', metavar=('NCOLS', 'NROWS')),
 255             'single_frequency_eV': lambda ap: ap.add_argument("-f", "--eV", type=float, required=True, help='radiation angular frequency in eV'),
 256             'multiple_frequency_eV_optional': lambda ap: ap.add_argument("-f", "--eV", type=float, nargs='*', help='radiation angular frequency in eV (additional)'),
 257             'seq_frequency_eV': lambda ap: ap.add_argument("-F", "--eV-seq", type=float, nargs=3, required=True, help='uniform radiation angular frequency sequence in eV', metavar=('FIRST', 'INCREMENT', 'LAST')),
 258             'real_frequencies_eV_ng': lambda ap: ap.add_argument("-f", "--eV", type=float_range, nargs=1, action='append', required=True, help='Angular frequency (or angular frequency range) in eV'), # nargs='+', action='extend' would be better, but action='extend' requires python>=3.8
 259             'single_material': lambda ap: ap.add_argument("-m", "--material", help='particle material (Au, Ag, ... for Lorentz-Drude or number for constant refractive index)', type=material_spec, required=True),
 260             'single_radius': lambda ap: ap.add_argument("-r", "--radius", type=float, required=True, help='particle radius (sphere or cylinder)'),
 261             'single_height': lambda ap: ap.add_argument("-H", "--height", type=float, help='cylindrical particle height; if not provided, particle is assumed to be spherical'),
 262             'single_kvec2': lambda ap: ap.add_argument("-k", '--kx-lim', nargs=2, type=sfloat, required=True, help='k vector', metavar=('KX_MIN', 'KX_MAX')),
 263             'kpi': lambda ap: ap.add_argument("--kpi", action='store_true', help="Indicates that the k vector is given in natural units instead of SI, i.e. the arguments given by -k shall be automatically multiplied by pi / period (given by -p argument)"),
 264             'bg_real_refractive_index': lambda ap: ap.add_argument("-n", "--refractive-index", type=float, default=1., help='background medium strictly real refractive index'),
 265             'bg_analytical': lambda ap: ap.add_argument("-B", "--background", type=material_spec, default=1., help="Background medium specification (constant real or complex refractive index, or supported material label)"),
 266             'single_lMax': lambda ap: ap.add_argument("-L", "--lMax", type=int, required=True, default=3, help='multipole degree cutoff'),
 267             'single_lMax_extend': lambda ap: ap.add_argument("--lMax-extend", type=int, required=False, default=6, help='multipole degree cutoff for T-matrix calculation (cylindrical particles only'),
 268             'outfile': lambda ap: ap.add_argument("-o", "--output", type=str, required=False, help='output path (if not provided, will be generated automatically)'), # TODO consider type=argparse.FileType('w')
 269             'plot_out': lambda ap: ap.add_argument("-O", "--plot-out", type=str, required=False, help="path to plot output (optional)"),
 270             'plot_do': lambda ap: ap.add_argument("-P", "--plot", action='store_true', help="if -p not given, plot to a default path"),
 271             'lattice2d_basis': lambda ap: ap.add_argument("-b", "--basis-vector", nargs='+', action=AppendTupleAction, help="basis vector in xy-cartesian coordinates (two required)", required=True, type=sfloat, dest='basis_vectors', metavar=('X', 'Y')),
 272             'planewave_pol_angles': __add_planewave_argparse_group,
 273             'multi_particle': __add_manyparticle_argparse_group,
 274     }
 275
 276     feature_sets_available = { # name : (description, dependencies, atoms not in other dependencies, methods called after parsing, "virtual" features provided)
 277             'const_real_background': ("Background medium with constant real refractive index", (), ('bg_real_refractive_index',), ('_eval_const_background_epsmu',), ('background', 'background_analytical')),
 278             'background' : ("Most general background medium specification currently supported", ('background_analytical',), (), (), ()),
 279             'background_analytical' : ("Background medium model holomorphic for 'reasonably large' complex frequency areas", (), ('bg_analytical',), ('_eval_analytical_background_epsmugen',), ('background',)),
 280             'single_particle': ("Single particle definition (shape [currently spherical or cylindrical]) and materials, incl. background)", ('background',), ('single_material', 'single_radius', 'single_height', 'single_lMax_extend'), ('_eval_single_tmgen',), ()),
 281             'multi_particle': ("One or more particle definition (shape [curently spherical or cylindrical]), materials, and positions)", ('background',), ('multi_particle',), ('_process_multi_particle',), ()),
 282             'single_lMax': ("Single particle lMax definition", (), ('single_lMax',), (), ()),
 283             'single_omega': ("Single angular frequency", (), ('single_frequency_eV',), ('_eval_single_omega',), ()),
 284             'omega_seq': ("Equidistant real frequency range with possibility of adding individual frequencies", (), ('seq_frequency_eV', 'multiple_frequency_eV_optional',), ('_eval_omega_seq',), ()),
 285             'omega_seq_real_ng': ("Equidistant real frequency ranges or individual frequencies (new syntax)", (), ('real_frequencies_eV_ng',), ('_eval_omega_seq_real_ng',), ()),
 286             'lattice2d': ("Specification of a generic 2d lattice (spanned by the x,y axes)", (), ('lattice2d_basis',), ('_eval_lattice2d',), ()),
 287             'rectlattice2d': ("Specification of a rectangular 2d lattice; conflicts with lattice2d", (), ('rectlattice2d_periods',), ('_eval_rectlattice2d',), ()),
 288             'rectlattice2d_finite': ("Specification of a rectangular 2d lattice; conflicts with lattice2d", ('rectlattice2d',), ('rectlattice2d_counts',), (), ()),
 289             'planewave': ("Specification of a normalised plane wave (typically used for scattering) with a full polarisation state", (), ('planewave_pol_angles',), ("_process_planewave_angles",), ()),
 290     }
 291
 292
 293     def __init__(self, features=[]):
 294         prefix_chars = '+-' if 'multi_particle' in features else '-'
 295         self.ap = argparse.ArgumentParser(prefix_chars=prefix_chars)
 296         self.features_enabled = set()
 297         self.call_at_parse_list = []
 298         self.parsed = False
 299         for feat in features:
 300             self.add_feature(feat)
 301         self._emg_register = {} # EpsMuGenerator dictionary to avoid recreating equivalent instances; filled by _add_emg()
 302         self._tmg_register = {} # TMatrixGenerator dictionary to avoid recreating equivalent instances; filled by _add_tmg()
 303         self._bspec_register = {} # Dictionary of used BaseSpecs to keep the equivalent instances unique; filled by _add_bspec()
 304
 305     def _add_emg(self, emgspec):
 306         """Looks up whether if an EpsMuGenerator from given material_spec has been already registered, and if not, creates a new one"""
 307         from .cymaterials import EpsMu, EpsMuGenerator, lorentz_drude
 308         if emgspec in self._emg_register.keys():
 309             return self._emg_register[emgspec]
 310         else:
 311             if isinstance(emgspec, (float, complex)):
 312                 emg = EpsMuGenerator(EpsMu(emgspec**2))
 313             else:
 314                 emg = EpsMuGenerator(lorentz_drude[emgspec])
 315             self._emg_register[emgspec] = emg
 316             return emg
 317
 318     def _add_tmg(self, tmgspec):
 319         """Looks up whether if a T-matrix from given T-matrix specification tuple has been already registered, and if not, creates a new one
 320
 321         T-matrix specification shall be of the form
 322         (bg_material_spec, fg_material_spec, shape_spec) where shape_spec is
 323         (radius, height, lMax_extend)
 324         """
 325         if tmgspec in self._tmg_register.keys():
 326             return self._tmg_register[tmgspec]
 327         else:
 328             from .cytmatrices import TMatrixGenerator
 329             bgspec, fgspec, (radius, height, lMax_extend) = tmgspec
 330             bg = self._add_emg(bgspec)
 331             fg = self._add_emg(fgspec)
 332             if height is None:
 333                 tmgen = TMatrixGenerator.sphere(bg, fg, radius)
 334             else:
 335                 tmgen = TMatrixGenerator.cylinder(bg, fg, radius, height, lMax_extend=lMax_extend)
 336             self._tmg_register[tmgspec] = tmgen
 337             return tmgen
 338
 339     def _add_bspec(self, key):
 340         if key in self._bspec_register.keys():
 341             return self._bspec_register[key]
 342         else:
 343             from .cybspec import BaseSpec
 344             if isinstance(key, BaseSpec):
 345                 bspec = key
 346             elif isinstance(key, int):
 347                 bspec = self._add_bspec(BaseSpec(lMax=key))
 348             else: raise TypeError("Can't register this as a BaseSpec")
 349             self._bspec_register[key] = bspec
 350             return bspec
 351
 352     def add_feature(self, feat):
 353         if feat not in self.features_enabled:
 354             if feat not in ArgParser.feature_sets_available:
 355                 raise ValueError("Unknown ArgParser feature: %s" % feat)
 356             #resolve dependencies
 357             _, deps, atoms, atparse, provides_virtual = ArgParser.feature_sets_available[feat]
 358             for dep in deps:
 359                 self.add_feature(dep)
 360             for atom in atoms: # maybe check whether that atom has already been added sometimes in the future?
 361                 ArgParser.atomic_arguments[atom](self.ap)
 362             for methodname in atparse:
 363                 self.call_at_parse_list.append(methodname)
 364             self.features_enabled.add(feat)
 365             for feat_virt in provides_virtual:
 366                 self.features_enabled.add(feat_virt)
 367
 368     def add_argument(self, *args, **kwargs):
 369         '''Add a custom argument directly to the standard library ArgParser object'''
 370         return self.ap.add_argument(*args, **kwargs)
 371
 372     def add_argument_group(self, *args, **kwargs):
 373         '''Add a custom argument group directly to the standard library ArgParser object'''
 374         return self.ap.add_argument_group(*args, **kwargs)
 375
 376     def parse_args(self, process_data = True,  *args, **kwargs):
 377         self.args = self.ap.parse_args(*args, **kwargs)
 378         if process_data:
 379             for method in self.call_at_parse_list:
 380                 try:
 381                     getattr(self, method)()
 382                 except ArgumentProcessingError:
 383                     err = sys.exc_info()[1]
 384                     self.ap.error(str(err))
 385         return self.args
 386
 387     def __getattr__(self, name):
 388         return getattr(self.args, name)
 389
 390
 391     # Methods to initialise the related data structures:
 392
 393     def _eval_const_background_epsmu(self): # feature: const_real_background
 394         self.args.background = self.args.refractive_index
 395         self._eval_analytical_background_epsmugen()
 396
 397     def _eval_analytical_background_epsmugen(self): # feature: background_analytical
 398         a = self.args
 399         from .cymaterials import EpsMu
 400         if isinstance(a.background, (float, complex)):
 401             self.background_epsmu = EpsMu(a.background**2)
 402         self.background_emg = self._add_emg(a.background)
 403
 404     def _eval_single_tmgen(self): # feature: single_particle
 405         a = self.args
 406         from .cymaterials import EpsMuGenerator, lorentz_drude
 407         from .cytmatrices import TMatrixGenerator
 408         self.foreground_emg = self._add_emg(a.material)
 409         self.tmgen = self._add_tmg((a.background, a.material, (a.radius, a.height, a.lMax_extend)))
 410         self.bspec = self._add_bspec(a.lMax)
 411
 412     def _eval_single_omega(self): # feature: single_omega
 413         from .constants import eV, hbar
 414         self.omega = self.args.eV * eV / hbar
 415
 416     def _eval_omega_seq(self): # feature: omega_seq
 417         import numpy as np
 418         from .constants import eV, hbar
 419         start, step, stop = self.args.eV_seq
 420         self.omegas = np.arange(start, stop, step)
 421         if self.args.eV:
 422             self.omegas = np.concatenate((self.omegas, np.array(self.args.eV)))
 423             self.omegas.sort()
 424         self.omegas *= eV/hbar
 425
 426     def _eval_omega_seq_real_ng(self): # feature: omega_seq_real_ng
 427         import numpy as np
 428         from .constants import eV, hbar
 429         eh = eV / hbar
 430         self.omegas = [omega_eV * eh for  omega_eV in flatten(self.args.eV)]
 431         self.omega_max = max(om if isinstance(om, float) else max(om) for om in self.omegas)
 432         self.omega_min = min(om if isinstance(om, float) else min(om) for om in self.omegas)
 433         self.omega_singles = [om for om in self.omegas if isinstance(om, float)]
 434         self.omega_ranges = [om for om in self.omegas if not isinstance(om, float)]
 435         self.omega_descr = ("%geV" % (self.omega_max / eh)) if (self.omega_max == self.omega_min) else (
 436                 "%g–%geV" % (self.omega_min / eh, self.omega_max / eh))
 437         self.allomegas = []
 438         for om in self.omegas:
 439             if isinstance(om, float):
 440                 self.allomegas.append(om)
 441             else:
 442                 self.allomegas.extend(om)
 443         self.allomegas = np.unique(self.allomegas)
 444
 445
 446     def _eval_lattice2d(self): # feature: lattice2d
 447         l = len(self.args.basis_vectors)
 448         if l != 2: raise ValueError('Two basis vectors must be specified (have %d)' % l)
 449         from .qpms_c import lll_reduce
 450         self.direct_basis = lll_reduce(self.args.basis_vectors, delta=1.)
 451         import numpy as np
 452         self.reciprocal_basis1 = np.linalg.inv(self.direct_basis.T)
 453         self.reciprocal_basis2pi = 2 * np.pi * self.reciprocal_basis1
 454
 455     def _eval_rectlattice2d(self): # feature: rectlattice2d
 456         a = self.args
 457         l = len(a.period)
 458         if (l == 1): # square lattice
 459             a.period = (a.period[0], a.period[0])
 460         else:
 461             a.period = (a.period[0], a.period[1])
 462             if (l > 2):
 463                 raise ValueError("At most two lattice periods allowed for a rectangular lattice (got %d)" % l)
 464
 465         import numpy as np
 466         a.basis_vectors = [(a.period[0], 0.), (0., a.period[1])]
 467         self.direct_basis = np.array(a.basis_vectors)
 468         self.reciprocal_basis1 = np.linalg.inv(self.direct_basis.T)
 469         self.reciprocal_basis2pi = 2 * np.pi * self.reciprocal_basis1
 470
 471     def _process_planewave_angles(self): #feature: planewave
 472         import math
 473         pi2 = math.pi/2
 474         a = self.args
 475         a.chi = a.chi * pi2
 476         a.psi = a.psi * pi2
 477         a.theta = a.theta * pi2
 478         a.phi = a.phi * pi2
 479
 480     def _process_multi_particle(self): # feature: multi_particle
 481         a = self.args
 482         self.tmspecs = {}
 483         self.tmgens = {}
 484         self.bspecs = {}
 485         self.positions = {}
 486         pos13, pos23, pos33 = False, False, False # used to
 487         if len(a.position.keys()) == 0:
 488             warnings.warn("No particle position (-p or +p) specified, assuming single particle in the origin / single particle per unit cell!")
 489             a.position[None] = [(0.,0.,0.)]
 490         for poslabel in a.position.keys():
 491             try:
 492                 lMax = a.lMax.get(poslabel, False) or a.lMax[None]
 493                 radius = a.radius.get(poslabel, False) or a.radius[None]
 494                 # Height is "inherited" only together with radius
 495                 height = a.height.get(poslabel, None) if poslabel in a.radius.keys() else a.height.get(None, None)
 496                 if hasattr(a, 'lMax_extend'):
 497                     lMax_extend = a.lMax_extend.get(poslabel, False) or a.lMax_extend.get(None, False) or None
 498                 else:
 499                     lMax_extend = None
 500                 material = a.material.get(poslabel, False) or a.material[None]
 501             except (TypeError, KeyError) as exc:
 502                 if poslabel is None:
 503                     raise ArgumentProcessingError("Unlabeled particles' positions (-p) specified, but some default particle properties are missing (--lMax, --radius, and --material have to be specified)") from exc
 504                 else:
 505                     raise ArgumentProcessingError(("Incomplete specification of '%s'-labeled particles: you must"
 506                         "provide at least ++lMax, ++radius, ++material arguments with the label, or the fallback arguments"
 507                         "--lMax, --radius, --material.")%(str(poslabel),)) from exc
 508             tmspec = (a.background, material, (radius, height, lMax_extend))
 509             self.tmspecs[poslabel] = tmspec
 510             self.tmgens[poslabel] = self._add_tmg(tmspec)
 511             self.bspecs[poslabel] = self._add_bspec(lMax)
 512             poslist_cured = []
 513             for pos in a.position[poslabel]:
 514                 if len(pos) == 1:
 515                     pos_cured = (0., 0., pos[0])
 516                     pos13 = True
 517                 elif len(pos) == 2:
 518                     pos_cured = (pos[0], pos[1], 0.)
 519                     pos23 = True
 520                 elif len(pos) == 3:
 521                     pos_cured = pos
 522                     pos33 = True
 523                 else:
 524                     raise argparse.ArgumentTypeError("Each -p / +p argument requires 1 to 3 cartesian coordinates")
 525                 poslist_cured.append(pos_cured)
 526             self.positions[poslabel] = poslist_cured
 527         if pos13 and pos23:
 528             warnings.warn("Both 1D and 2D position specifications used. The former are interpreted as z coordinates while the latter as x, y coordinates")
 529
 530     def get_particles(self):
 531         """Creates a list of Particle instances that can be directly used in ScatteringSystem.create().
 532
 533         Assumes that self._process_multi_particle() has been already called.
 534         """
 535         from .qpms_c import Particle
 536         plist = []
 537         for poslabel, poss in self.positions.items():
 538             t = self.tmgens[poslabel]
 539             bspec = self.bspecs[poslabel]
 540             plist.extend([Particle(pos, t, bspec=bspec) for pos in poss])
 541         return plist
 542
 543 #TODO perhaps move into another module
 544 def annotate_pdf_metadata(pdfPages, scriptname=None, keywords=None, author=None, title=None, subject=None, **kwargs):
 545     """Adds QPMS version-related metadata to a matplotlib PdfPages object
 546
 547     Use before closing the PDF file.
 548     """
 549     from .qpms_c import qpms_library_version
 550     d = pdfPages.infodict()
 551     d['Creator'] = "QPMS%s (lib rev. %s), https://qpms.necada.org" % (
 552         "" if scriptname is None else (" "+scriptname), qpms_library_version())
 553     if author is not None:
 554         d['Author'] = author
 555     if title is not None:
 556         d['Title'] = title
 557     if subject is not None:
 558         d['Subject'] = subject
 559     if keywords is not None:
 560         d['Keywords'] = ' '.join(keywords)
 561     d.update(kwargs)