Source code for pygimli.physics.gravimetry.MagneticsModelling

"""Magnetics forward operator."""
import numpy as np
import pygimli as pg
from .kernel import SolveGravMagHolstein

[docs] class MagneticsModelling(pg.frameworks.MeshModelling): """Magnetics modelling operator using Holstein (2007)."""
[docs] def __init__(self, mesh=None, points=None, cmp=["TFA"], igrf=[50, 13], foot=None): """Setup forward operator. Parameters ---------- mesh : pygimli:mesh tetrahedral or hexahedral mesh points : list|array of (x, y, z) measuring points cmp : list of str component of: gx, gy, gz, TFA, Bx, By, Bz, Bxy, Bxz, Byy, Byz, Bzz igrf : list|array of size 3 or 7 international geomagnetic reference field, either [D, I, H, X, Y, Z, F] - declination, inclination, horizontal field, X/Y/Z components, total field OR [X, Y, Z] - X/Y/Z components [lat, lon] - latitude, longitude (automatic IGRF) """ # check if components do not contain g! super().__init__() self._refineH2 = False # self.createRefinedForwardMesh(refine=False, pRefine=False) self.mesh_ = mesh self.sensorPositions = points self.components = cmp self.igrf = None if hasattr(igrf, "__iter__"): if len(igrf) == 2: # lat lon import pyIGRF self.igrf = pyIGRF.igrf_value(*igrf) else: self.igrf = igrf self.footprint = foot self.kernel = None self.J = pg.matrix.BlockMatrix() if self.mesh_ is not None: self.setMesh(self.mesh_)
[docs] def computeKernel(self): """Compute the kernel.""" points = np.column_stack([self.sensorPositions[:, 1], self.sensorPositions[:, 0], -np.abs(self.sensorPositions[:, 2])]) self.kernel = SolveGravMagHolstein(self.mesh().NED(), pnts=points, igrf=self.igrf, cmp=self.components, foot=self.footprint) self.J = pg.matrix.BlockMatrix() self.Ki = [] self.Ji = [] for iC in range(self.kernel.shape[1]): self.Ki.append(np.squeeze(self.kernel[:, iC, :])) self.Ji.append(pg.matrix.NumpyMatrix(self.Ki[-1])) self.J.addMatrix(self.Ji[-1], iC*self.kernel.shape[0], 0) self.J.recalcMatrixSize() self.setJacobian(self.J)
# better move the latter to # self.createKernel # def setMesh(self, mesh): # self.createKernel(mesh)
[docs] def response(self, model): """Compute forward response.""" if self.kernel is None: self.computeKernel() return
[docs] def createJacobian(self, model): """Do nothing as this is a linear problem.""" pass