flake8 code cleanups.

gprMax/geometry_outputs.py

@@ -58,7 +58,7 @@ class GeometryView(object):
         self.dz = dz
         self.basefilename = filename
         self.fileext = fileext
-        
+
         if self.fileext == '.vti':
             # Calculate number of cells according to requested sampling for geometry view
             self.vtk_xscells = round_value(self.xs / self.dx)
@@ -71,7 +71,7 @@ class GeometryView(object):
             self.vtk_nycells = self.vtk_yfcells - self.vtk_yscells
             self.vtk_nzcells = self.vtk_zfcells - self.vtk_zscells
             self.datawritesize = int(np.dtype(np.uint32).itemsize * self.vtk_nxcells * self.vtk_nycells * self.vtk_nzcells) + 2 * (int(np.dtype(np.int8).itemsize * self.vtk_nxcells * self.vtk_nycells * self.vtk_nzcells))
-        
+
         elif self.fileext == '.vtp':
             self.vtk_numpoints = (self.nx + 1) * (self.ny + 1) * (self.nz + 1)
             self.vtk_numpoint_components = 3

gprMax/gprMax.py

@@ -391,7 +391,7 @@ def run_model(args, modelrun, numbermodelruns, inputfile, usernamespace):
         if all(value == 0 for value in G.pmlthickness.values()):
             if G.messages:
                 print('PML boundaries: switched off')
-            pass # If all the PMLs are switched off don't need to build anything
+            pass  # If all the PMLs are switched off don't need to build anything
         else:
             if G.messages:
                 if all(value == G.pmlthickness['xminus'] for value in G.pmlthickness.values()):
@@ -482,7 +482,7 @@ def run_model(args, modelrun, numbermodelruns, inputfile, usernamespace):
             geometryview.write_vtk(modelrun, numbermodelruns, G, pbar)
             pbar.close()
     if G.geometryobjectswrite:
-        
+
         for i, geometryobject in enumerate(G.geometryobjectswrite):
             pbar = tqdm(total=geometryobject.datawritesize, unit='byte', unit_scale=True, desc='Writing geometry object file {} of {}, {}'.format(i + 1, len(G.geometryobjectswrite), os.path.split(geometryobject.filename)[1]), ncols=get_terminal_width() - 1, file=sys.stdout, disable=G.tqdmdisable)
             geometryobject.write_hdf5(G, pbar)

gprMax/grid.py

@@ -78,12 +78,12 @@ class FDTDGrid(Grid):
         self.title = ''
         self.messages = True
         self.tqdmdisable = False
-        
+
-        # Threshold (dB) down from maximum power (0dB) of main frequency used to calculate highest frequency for disperion analysis
+        # Threshold (dB) down from maximum power (0dB) of main frequency used to calculate highest frequency for disperion analysis
         self.highestfreqthres = 60
         # Maximum allowable percentage physical phase-velocity phase error
         self.maxnumericaldisp = 2
-        
+
         self.nx = 0
         self.ny = 0
         self.nz = 0
@@ -157,7 +157,7 @@ def dispersion_analysis(G):
 
     # Physical phase velocity error (percentage); grid sampling density; material with maximum permittivity; maximum frequency of interest
     results = {'deltavp': False, 'N': False, 'material': False, 'maxfreq': False}
-    
+
     # Find maximum frequency
     maxfreqs = []
     for waveform in G.waveforms:
@@ -221,19 +221,19 @@ def dispersion_analysis(G):
 
         # Minimum wavelength
         minwavelength = minvelocity / results['maxfreq']
-        
+
         # Maximum spatial step
         delta = max(G.dx, G.dy, G.dz)
-        
+
         # Courant stability factor
         S = (c * G.dt) / delta
-        
+
         # Grid sampling density
         results['N'] = minwavelength / delta
-        
+
         # Numerical phase velocity
         vp = np.pi / (results['N'] * np.arcsin((1 / S) * np.sin((np.pi * S) / results['N'])))
-        
+
         # Physical phase velocity error (percentage)
         results['deltavp'] = (((vp * c) - c) / c) * 100
 

gprMax/input_cmd_funcs.py

@@ -38,7 +38,7 @@ Coordinate_tuple = namedtuple('Coordinate', ['x', 'y', 'z'])
 
 class Coordinate(Coordinate_tuple):
     """Subclass of a namedtuple where __str__ outputs 'x y z'"""
-    
+
     def __str__(self):
         return '{:g} {:g} {:g}'.format(self.x, self.y, self.z)
 
@@ -53,7 +53,7 @@ def command(cmd, *parameters):
     Returns:
         s (str): the printed string
     """
-    
+
     # remove Nones
     filtered = filter(lambda x: x is not None, parameters)
     # convert to str
@@ -76,11 +76,11 @@ def command(cmd, *parameters):
 
 def rotate90_point(x, y, rotate90origin=()):
     """Rotates a point 90 degrees CCW in the x-y plane.
-        
+
     Args:
         x, y (float): Coordinates.
         rotate90origin (tuple): x, y origin for 90 degree CCW rotation in x-y plane.
-        
+
     Returns:
         xrot, yrot (float): Rotated coordinates.
     """
@@ -158,7 +158,7 @@ def domain(x, y, z):
     Returns:
         domain (Coordinate): Namedtuple of the extent of the domain.
     """
-    
+
     domain = Coordinate(x, y, z)
     command('domain', domain)
 
@@ -174,7 +174,7 @@ def dx_dy_dz(x, y, z):
     Returns:
         dx_dy_dz (float): Tuple of the spatial resolutions.
     """
-    
+
     dx_dy_dz = Coordinate(x, y, z)
     command('dx_dy_dz', dx_dy_dz)
 
@@ -190,7 +190,7 @@ def time_window(time_window):
     Returns:
         time_window (float): Duration of simulation.
     """
-    
+
     command('time_window', time_window)
 
     return time_window
@@ -206,7 +206,7 @@ def material(permittivity, conductivity, permeability, magconductivity, name):
         magconductivity (float): Magnetic loss of the material.
         name (str): Material identifier.
     """
-    
+
     command('material', permittivity, conductivity, permeability, magconductivity, name)
 
 
@@ -222,7 +222,7 @@ def geometry_view(xs, ys, zs, xf, yf, zf, dx, dy, dz, filename, type='n'):
     Returns:
         s, f, d (tuple): 3 namedtuple Coordinate for the start, finish coordinates and spatial discretisation
     """
-    
+
     s = Coordinate(xs, ys, zs)
     f = Coordinate(xf, yf, zf)
     d = Coordinate(dx, dy, dz)
@@ -243,7 +243,7 @@ def snapshot(xs, ys, zs, xf, yf, zf, dx, dy, dz, time, filename):
     Returns:
         s, f, d (tuple): 3 namedtuple Coordinate for the start, finish coordinates and spatial discretisation
     """
-    
+
     s = Coordinate(xs, ys, zs)
     f = Coordinate(xf, yf, zf)
     d = Coordinate(dx, dy, dz)
@@ -269,14 +269,13 @@ def edge(xs, ys, zs, xf, yf, zf, material, rotate90origin=()):
     Returns:
         s, f (tuple): 2 namedtuple Coordinate for the start and finish coordinates
     """
-    
+
     if rotate90origin:
         if xs == xf:
             polarisation = 'y'
         else:
             polarisation = 'x   '
         xs, ys, xf, yf = rotate90_edge(xs, ys, xf, yf, polarisation, rotate90origin)
-            
 
     s = Coordinate(xs, ys, zs)
     f = Coordinate(xf, yf, zf)
@@ -296,10 +295,10 @@ def plate(xs, ys, zs, xf, yf, zf, material, rotate90origin=()):
     Returns:
         s, f (tuple): 2 namedtuple Coordinate for the start and finish coordinates
     """
-    
+
     if rotate90origin:
         xs, ys, xf, yf = rotate90_plate(xs, ys, xf, yf, rotate90origin)
-    
+
     s = Coordinate(xs, ys, zs)
     f = Coordinate(xf, yf, zf)
     command('plate', s, f, material)
@@ -320,12 +319,12 @@ def triangle(x1, y1, z1, x2, y2, z2, x3, y3, z3, thickness, material, averaging=
     Returns:
         v1, v2, v3 (tuple): 3 namedtuple Coordinate for the vertices
     """
-    
+
     if rotate90origin:
         x1, y1 = rotate90_point(x1, y1, rotate90origin)
         x2, y2 = rotate90_point(x2, y2, rotate90origin)
         x3, y3 = rotate90_point(x3, y3, rotate90origin)
-    
+
     v1 = Coordinate(x1, y1, z1)
     v2 = Coordinate(x2, y2, z2)
     v3 = Coordinate(x3, y3, z3)
@@ -346,10 +345,10 @@ def box(xs, ys, zs, xf, yf, zf, material, averaging='', rotate90origin=()):
     Returns:
         s, f (tuple): 2 namedtuple Coordinate for the start and finish coordinates
     """
-    
+
     if rotate90origin:
         xs, ys, xf, yf = rotate90_plate(xs, ys, xf, yf, rotate90origin)
-    
+
     s = Coordinate(xs, ys, zs)
     f = Coordinate(xf, yf, zf)
     command('box', s, f, material, averaging)
@@ -369,7 +368,7 @@ def sphere(x, y, z, radius, material, averaging=''):
     Returns:
         c (tuple): namedtuple Coordinate for the center of the sphere
     """
-    
+
     c = Coordinate(x, y, z)
     command('sphere', c, radius, material, averaging)
 
@@ -389,11 +388,11 @@ def cylinder(x1, y1, z1, x2, y2, z2, radius, material, averaging='', rotate90ori
     Returns:
         c1, c2 (tuple): 2 namedtuple Coordinate for the centres of the two faces of the cylinder.
     """
-    
+
     if rotate90origin:
         x1, y1 = rotate90_point(x1, y1, rotate90origin)
         x2, y2 = rotate90_point(x2, y2, rotate90origin)
-    
+
     c1 = Coordinate(x1, y1, z1)
     c2 = Coordinate(x2, y2, z2)
     command('cylinder', c1, c2, radius, material, averaging)
@@ -414,7 +413,7 @@ def cylindrical_sector(axis, ctr1, ctr2, t1, t2, radius, startingangle, sweptang
         material (str): Material identifier(s).
         averaging (str): Turn averaging on or off.
     """
-    
+
     command('cylindrical_sector', axis, ctr1, ctr2, t1, t2, radius, startingangle, sweptangle, material, averaging)
 
 
@@ -427,7 +426,7 @@ def excitation_file(file1):
     Returns:
         file1 (str): filename
     """
-    
+
     command('excitation_file', file1)
 
     return file1
@@ -445,7 +444,7 @@ def waveform(shape, amplitude, frequency, identifier):
     Returns:
         identifier (str): is an identifier for the waveform used to assign it to a source.
     """
-    
+
     command('waveform', shape, amplitude, frequency, identifier)
 
     return identifier
@@ -466,7 +465,7 @@ def hertzian_dipole(polarisation, f1, f2, f3, identifier, t0=None, t_remove=None
     Returns:
         coordinates (tuple): namedtuple Coordinate of the source location
     """
-    
+
     if rotate90origin:
         if polarisation == 'x':
             xf = f1 + dxdy[0]
@@ -476,7 +475,7 @@ def hertzian_dipole(polarisation, f1, f2, f3, identifier, t0=None, t_remove=None
             xf = f1
             yf = f2 + dxdy[1]
             newpolarisation = 'x'
-        
+
         f1, f2, xf, yf = rotate90_edge(f1, f2, xf, yf, polarisation, rotate90origin)
         polarisation = newpolarisation
 
@@ -502,7 +501,7 @@ def magnetic_dipole(polarisation, f1, f2, f3, identifier, t0=None, t_remove=None
     Returns:
         coordinates (tuple): namedtuple Coordinate of the source location
     """
-    
+
     if rotate90origin:
         if polarisation == 'x':
             xf = f1 + dxdy[0]
@@ -512,7 +511,7 @@ def magnetic_dipole(polarisation, f1, f2, f3, identifier, t0=None, t_remove=None
             xf = f1
             yf = f2 + dxdy[1]
             newpolarisation = 'x'
-        
+
         f1, f2, xf, yf = rotate90_edge(f1, f2, xf, yf, polarisation, rotate90origin)
         polarisation = newpolarisation
 
@@ -539,7 +538,7 @@ def voltage_source(polarisation, f1, f2, f3, resistance, identifier, t0=None, t_
     Returns:
         coordinates (tuple): namedtuple Coordinate of the source location
     """
-    
+
     if rotate90origin:
         if polarisation == 'x':
             xf = f1 + dxdy[0]
@@ -549,7 +548,7 @@ def voltage_source(polarisation, f1, f2, f3, resistance, identifier, t0=None, t_
             xf = f1
             yf = f2 + dxdy[1]
             newpolarisation = 'x'
-        
+
         f1, f2, xf, yf = rotate90_edge(f1, f2, xf, yf, polarisation, rotate90origin)
         polarisation = newpolarisation
 
@@ -562,7 +561,7 @@ def voltage_source(polarisation, f1, f2, f3, resistance, identifier, t0=None, t_
 
 def transmission_line(polarisation, f1, f2, f3, resistance, identifier, t0=None, t_remove=None, dxdy=None, rotate90origin=()):
     """Prints the #transmission_line: polarisation, f1, f2, f3, resistance, identifier, [t0, t_remove]
-        
+
     Args:
         polarisation (str):  is the polarisation of the source and can be 'x', 'y', or 'z'.
         f1 f2 f3 (float): are the coordinates (x,y,z) of the source in the model.
@@ -572,11 +571,11 @@ def transmission_line(polarisation, f1, f2, f3, resistance, identifier, t0=None,
         t_remove (float): is a time to remove the source.
         dxdy (float): Tuple of x-y spatial resolutions. For rotation purposes only.
         rotate90origin (tuple): x, y origin for 90 degree CCW rotation in x-y plane.
-        
+
     Returns:
         coordinates (tuple): namedtuple Coordinate of the source location
     """
-    
+
     if rotate90origin:
         if polarisation == 'x':
             xf = f1 + dxdy[0]
@@ -586,14 +585,14 @@ def transmission_line(polarisation, f1, f2, f3, resistance, identifier, t0=None,
             xf = f1
             yf = f2 + dxdy[1]
             newpolarisation = 'x'
-        
+
         f1, f2, xf, yf = rotate90_edge(f1, f2, xf, yf, polarisation, rotate90origin)
         polarisation = newpolarisation
-    
+
     c = Coordinate(f1, f2, f3)
     # since command ignores None, this is safe:
     command('transmission_line', polarisation, str(c), resistance, identifier, t0, t_remove)
-    
+
     return c
 
 
@@ -611,7 +610,7 @@ def rx(x, y, z, identifier=None, to_save=None, polarisation=None, dxdy=None, rot
     Returns:
         coordinates (tuple): namedtuple Coordinate of the receiver location
     """
-    
+
     if rotate90origin:
         if polarisation == 'x':
             xf = x + dxdy[0]
@@ -620,7 +619,7 @@ def rx(x, y, z, identifier=None, to_save=None, polarisation=None, dxdy=None, rot
             xf = x
             yf = y + dxdy[1]
         x, y, xf, yf = rotate90_edge(x, y, xf, yf, polarisation, rotate90origin)
-    
+
     c = Coordinate(x, y, z)
     command('rx', str(c), identifier, ' '.join(to_save))
 

gprMax/input_cmds_file.py

@@ -173,7 +173,7 @@ def check_cmd_names(processedlines, checkessential=True):
         cmd = processedlines[lindex].split(':')
         cmdname = cmd[0]
         cmdparams = cmd[1]
-        
+
         # Check if there is space between command name and parameters, i.e. check first character of parameter string
         if ' ' not in cmdparams[0]:
             raise CmdInputError('There must be a space between the command name and parameters in ' + processedlines[lindex])

gprMax/input_cmds_geometry.py

@@ -65,7 +65,7 @@ def process_geometrycmds(geometry, G):
 
             matstr = os.path.splitext(os.path.split(matfile)[1])[0]
             numexistmaterials = len(G.materials)
-            
+
             # Read materials from file
             with open(matfile, 'r') as f:
                 # Strip out any newline characters and comments that must begin with double hashes
@@ -73,7 +73,7 @@ def process_geometrycmds(geometry, G):
 
             # Check validity of command names
             singlecmdsimport, multicmdsimport, geometryimport = check_cmd_names(materials, checkessential=False)
-            
+
             # Process parameters for commands that can occur multiple times in the model
             process_multicmds(multicmdsimport, G)
 
@@ -873,7 +873,7 @@ def process_geometrycmds(geometry, G):
             # Search and process any modifiers for the fractal box
             for object in geometry:
                 tmp = object.split()
-                
+
                 if tmp[0] == '#add_surface_roughness:':
                     if len(tmp) < 13:
                         raise CmdInputError("'" + ' '.join(tmp) + "'" + ' requires at least twelve parameters')

gprMax/input_cmds_multiuse.py

@@ -502,11 +502,11 @@ def process_multicmds(multicmds, G):
             m.se = se
             m.mr = float(tmp[2])
             m.sm = float(tmp[3])
-            
+
             # Set material averaging to False if infinite conductivity, i.e. pec
             if m.se == float('inf'):
                 m.averagable = False
-            
+
             if G.messages:
                 tqdm.write('Material {} with epsr={:g}, sig={:g} S/m; mur={:g}, sig*={:g} S/m created.'.format(m.ID, m.er, m.se, m.mr, m.sm))
 

gprMax/materials.py

@@ -156,7 +156,7 @@ def process_materials(G):
         else:
             materialsdata = [['\nID', '\nName', '\nType', '\neps_r', 'sigma\n[S/m]', '\nDelta eps_r', 'tau\n[s]', 'omega\n[Hz]', 'delta\n[Hz]', 'gamma\n[Hz]', '\nmu_r', 'sigma*\n[S/m]', 'Dielectric\nsmoothable']]
 
-    for material in G.materials:       
+    for material in G.materials:
         # Calculate update coefficients for material
         material.calculate_update_coeffsE(G)
         material.calculate_update_coeffsH(G)
@@ -181,23 +181,23 @@ def process_materials(G):
             materialtext.append('{:g}'.format(material.se))
             if Material.maxpoles > 0:
                 if 'debye' in material.type:
-                    materialtext.append(', '.join('{:g}'.format(deltaer) for deltaer in material.deltaer))
+                    materialtext.append(', '.join('{:g}'.format(deltaer for deltaer in material.deltaer)))
-                    materialtext.append(', '.join('{:g}'.format(tau) for tau in material.tau))
+                    materialtext.append(', '.join('{:g}'.format(tau for tau in material.tau)))
                     materialtext.append('')
                     materialtext.append('')
                     materialtext.append('')
                 elif 'lorentz' in material.type:
-                    materialtext.append(', '.join('{:g}'.format(deltaer) for deltaer in material.deltaer))
+                    materialtext.append(', '.join('{:g}'.format(deltaer for deltaer in material.deltaer)))
                     materialtext.append('')
-                    materialtext.append(', '.join('{:g}'.format(tau) for tau in material.tau))
+                    materialtext.append(', '.join('{:g}'.format(tau for tau in material.tau)))
-                    materialtext.append(', '.join('{:g}'.format(alpha) for alpha in material.alpha))
+                    materialtext.append(', '.join('{:g}'.format(alpha for alpha in material.alpha)))
                     materialtext.append('')
                 elif 'drude' in material.type:
                     materialtext.append('')
                     materialtext.append('')
-                    materialtext.append(', '.join('{:g}'.format(tau) for tau in material.tau))
+                    materialtext.append(', '.join('{:g}'.format(tau for tau in material.tau)))
                     materialtext.append('')
-                    materialtext.append(', '.join('{:g}'.format(alpha) for tau in material.alpha))
+                    materialtext.append(', '.join('{:g}'.format(alpha for alpha in material.alpha)))
                 else:
                     materialtext.append('')
                     materialtext.append('')

gprMax/pml.py

@@ -262,7 +262,7 @@ class PML(object):
 def build_pmls(G, pbar):
     """This function builds instances of the PML and calculates the initial parameters and coefficients including setting profile
         (based on underlying material er and mr from solid array).
-        
+
     Args:
         G (class): Grid class instance - holds essential parameters describing the model.
         pbar (class): Progress bar class instance.
@@ -317,6 +317,6 @@ def build_pmls(G, pbar):
                         summr += material.mr
                 averageer = sumer / (G.nx * G.ny)
                 averagemr = summr / (G.nx * G.ny)
-            
+
             pml.calculate_update_coeffs(averageer, averagemr, G)
             pbar.update()

gprMax/sources.py

@@ -364,7 +364,7 @@ class PlaneWave(Source):
         """
 
         super(Source, self).__init__()
-        
+
         # Coordinates defining Huygen's surface
         self.xs = 0
         self.xf = 0
@@ -372,54 +372,54 @@ class PlaneWave(Source):
         self.yf = 0
         self.zs = 0
         self.zf = 0
-        
+
         # Spherical coordinates defining incident unit wavevector (k)
-        self.theta = 0 # 0 <= theta <= 180
+        self.theta = 0  # 0 <= theta <= 180
-        self.phi = 0 # 0 <= phi <= 360
+        self.phi = 0  # 0 <= phi <= 360
-        
+
         # Angle that incident electric field makes with k cross z
-        self.psi = 0 # 0 <= psi <= 360
+        self.psi = 0  # 0 <= psi <= 360
-    
+
     def calculate_origin(self, G):
         """Calculate origin of TF/SF interface with incident wavefront."""
-        
+
         if self.theta >= 0 and self.theta <= 90:
             if self.phi >= 0 and self.phi <= 90:
                 self.xcoordorigin = 0
                 self.ycoordorigin = 0
                 self.zcoordorigin = 0
-            
+
             elif self.phi > 90 and self.phi <= 180:
                 self.xcoordorigin = G.nx
                 self.ycoordorigin = 0
                 self.zcoordorigin = 0
-            
+
             elif self.phi > 180 and self.phi <= 270:
                 self.xcoordorigin = G.nx
                 self.ycoordorigin = G.ny
                 self.zcoordorigin = 0
-            
+
             elif self.phi > 270 and self.phi <= 360:
                 self.xcoordorigin = 0
                 self.ycoordorigin = G.ny
                 self.zcoordorigin = 0
-    
+
         elif self.theta > 90 and self.theta <= 180:
             if self.phi >= 0 and self.phi <= 90:
                 self.xcoordorigin = 0
                 self.ycoordorigin = 0
                 self.zcoordorigin = G.nz
-            
+
             elif self.phi > 90 and self.phi <= 180:
                 self.xcoordorigin = G.nx
                 self.ycoordorigin = 0
                 self.zcoordorigin = G.nz
-            
+
             elif self.phi > 180 and self.phi <= 270:
                 self.xcoordorigin = G.nx
                 self.ycoordorigin = G.ny
                 self.zcoordorigin = G.nz
-            
+
             elif self.phi > 270 and self.phi <= 360:
                 self.xcoordorigin = 0
                 self.ycoordorigin = G.ny
@@ -427,31 +427,20 @@ class PlaneWave(Source):
 
     def calculate_vector_components(self):
         """Calculate components of incident fields."""
-        
+
         self.theta = np.deg2rad(self.theta)
         self.phi = np.deg2rad(self.phi)
         self.psi = np.deg2rad(self.psi)
 
         # Components of incident unit wavevector
-        self.kx = np.sin(theta) * np.cos(phi)
+        self.kx = np.sin(self.theta) * np.cos(self.phi)
-        self.ky = np.sin(theta) * np.sin(phi)
+        self.ky = np.sin(self.theta) * np.sin(self.phi)
-        self.kz = np.cos(theta)
+        self.kz = np.cos(self.theta)
 
         # Components of incident field vectors
-        self.Exinc = np.cos(psi) * np.sin(phi) - np.sin(psi) * np.cos(theta) * np.cos(phi)
+        self.Exinc = np.cos(self.psi) * np.sin(self.phi) - np.sin(self.psi) * np.cos(self.theta) * np.cos(self.phi)
-        self.Eyinc = -np.cos(psi) * np.cos(phi) - np.sin(psi) * np.cos(theta) * np.sin(phi)
+        self.Eyinc = -np.cos(self.psi) * np.cos(self.phi) - np.sin(self.psi) * np.cos(self.theta) * np.sin(self.phi)
-        self.Ezinc = np.sin(psi) * np.sin(theta)
+        self.Ezinc = np.sin(self.psi) * np.sin(self.theta)
-        self.Hxinc = np.sin(psi) * np.sin(phi) + np.cos(psi) * np.cos(theta) * np.cos(phi)
+        self.Hxinc = np.sin(self.psi) * np.sin(self.phi) + np.cos(self.psi) * np.cos(self.theta) * np.cos(self.phi)
-        self.Hyinc = -np.sin(psi) * np.cos(phi) + np.cos(psi) * np.cos(theta) * np.sin(phi)
+        self.Hyinc = -np.sin(self.psi) * np.cos(self.phi) + np.cos(self.psi) * np.cos(self.theta) * np.sin(self.phi)
-        self.Hzinc = -np.cos(psi) * np.sin(theta)
+        self.Hzinc = -np.cos(self.psi) * np.sin(self.theta)
-
-
-
-
-
-
-
-
-
-
-