# -*- coding: utf-8 -*-
"""
@author: Rochlitz.R
"""
import numpy as np
[docs]
class PHC:
"""
PyhedCalculations (PHC) class for reference solution calculated with
**pyhed** sub-module of **COMET** library.
Methods
-------
- calc_reference()
calucalte semi-analytic reference solution for given 1D CSEM setup-
Usage
-----
Load a config file (JSON format) from a custEM modeling to obtain the Tx
information and calculate results for E and/or H on specified coordinates.
>>> from custEM.misc import pyhed_calculations as ph
>>> calculator = ph.PHC("CONFIG_FILE")
>>> calculator.calc_reference(np.array([0., 100., 0.]), 'EH', max_procs=4)
"""
def __init__(self, config_file):
"""
Initialize PHC class utilizing exisiting configuration file.
Required arguments
------------------
- config_file, type str
specified model configuration file of custEM results in JSON format
"""
if type(config_file) is str:
import json
with open(config_file, 'r') as cfile:
self.__dict__.update(json.load(cfile))
[docs]
def calc_reference(self, coords, EH, freq=None, max_procs=1):
"""
Calculate semi-analytic reference solution for 1D CSEm setups.
Required arguments
------------------
- coords, type array of shape (:, 3)
array of target coordinates
- EH, type str
specify either 'E' for electric or 'H' for magnetic fields
Keyword arguments
-----------------
- max_procs = 1, type int
allowed number of processes for **pyhed** internal multiprocessing
"""
print('... calculating ' + self.s_type +
' source reference solution ...')
from comet import pyhed as ph
if self.s_type == 'hed':
source = ph.loop.buildDipole(self.origin, length=self.length,
angle=np.deg2rad(self.azimuth))
elif self.s_type == 'line':
source = ph.loop.buildLine(self.start, self.stop,
num_segs=self.n_segs)
elif self.s_type == 'loop_c':
source = ph.loop.buildCircle(self.R, P=self.origin,
num_segs=self.n_segs)
elif self.s_type == 'loop_r':
S_1 = [self.start[0], self.start[1], 0.0]
S_2 = [self.stop[0], self.start[1], 0.0]
S_3 = [self.stop[0], self.stop[1], 0.0]
S_4 = [self.start[0], self.stop[1], 0.0]
source = ph.loop.buildRectangle(np.array((S_1, S_2, S_3, S_4)),
num_segs=self.n_segs)
elif self.s_type == 'path':
print('Notice! Only first Tx is considered for calculating '
'pyhed solutions for specified path source. Continuing ...')
source = ph.loop.buildLoop(np.array(self.tx[0]),
num_segs=self.n_segs,
grounded=self.grounding[0])
source.setLoopMesh(coords)
source.config.rho = 1. / np.array(
self.sigma[1:1+self.n_layers]).ravel()
source.config.d = self.layer_thicknesses
if freq is None:
source.config.f = self.omega / (2.0 * np.pi)
else:
source.config.f = freq
source.config.current = self.current
source.config.ftype = EH
data = source.calculate(interpolate=False, maxCPUCount=max_procs)
return(data.T)