#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Scattering calculation (:mod:`hydrac.model.scattering.scattering`)
==================================================================
.. autoclass:: Scattering
:members:
:private-members:
"""
import numpy as np
from hydrac.util.parameters import Parameters
from .thq import Thq
from .dwba import DWBA
from .emp import Emp
from .hp import Hp
from hydrac.util.display import Display
class InstrAmbiguity(Exception):
pass
[docs]class Scattering(object):
"""
Scattering class.
This class computes the scattering model given the input properties of
particles specified by the class :class:`hydrac.model.particle.Particle`,
and several other inputs such as the carrier frequency of the instruments,
the water parameters (ex. sound velocity in water)...
Two major cases are considered regarding the inversion procedure:
**The case of singlefrequency inversion**: In this case, the mean radius of
the particles is entered as input. The model is computed at high
resolution for display purpose and potential averaging over size
distribnution if one is indicated, but only one single value will be
used.
**The case of multifrequency inversion**: One of the multifrequency
inversion methods (NNLS, see :mod:`hydrac.inversion.multifrequence.nnls`)
considers an undetermined inversion system, where the number a size classes
is defined, the concentrations of which correspond to the system unknown.
In this case, in addition to the high-resolution model, it is convenient to
compute the model for each ka values corresponding to the instrument
frequencies and the target sizes of the particles. As such, the model is
computed at two distinct resolutions.
The "model" returned by the Scattering class is two-fold:
- The form function :math:`f`
- The total scattering cross-section :math:`\\chi`, from which the
attenuation due to scattering can be evaluated.
This class can be used as a standalone if one is interested in computing
a specific model, and export it elsewhere for other purposes. The Typical
usage for a standalone package to compute models (here shown for a model
parameters stored in the database, copepDWBA)::
# Define the particle properties with a Particle object
>>>P = hydrac.model.particle.Particle('copepDWBA',1)
# Define the water properties with a Water object
>>>W = hydrac.model.water.Water(preset='sea')
# Compute the target strength and attenuation
>>>S=hydrac.model.scattering.scattering.Scattering(P,
W,
freq=[1e6, 4e6],
tag_M='multi',
mode_depl='Mooring',
am=0.06e-3,
aM=2.8e-3)
# Access the form function
>>>S.param.f_inf
Within hydrac processing chain, the Scattering class takes the data from
the hydroacoustic object as inputs. It is solely called when performing
the inversion in the package :mod:`hdrac.inversion.inversion.Inversion`.
All the results are stored into a ``scattering`` modified dictionnary
(see :mod:`hydrac.util.parameters`).
Parameters
----------
particle : object
Particle object from :class:`hydrac.model.particle.Particle`
water : object
Water object from :class:`hydrac.model.water.Water`
ori : float
Angle of observation => angle between incident and oscerved scattered
wave (:math:`\\pi for backscattering`)
freq : float (numpy.array or list)
Target frequencies of the hydroacoustic instrument
tag_M : str, {'mono','multi'}
Tag specifying whether the instrument is multi or single frequency
gen_model : bool
Save the model in netCDF as a standalone
mode_depl : str, {'Mooring', 'Cast', ...}
Deployment strategy (see :mod:`hydrac.instrument.mod_deployment`)
am, aM : float, float
Minimum and maximum radius for model computation
dsigma : float
Distribution width if the user wants the model to be function of mean
radius. Not used for multifrequency inversions.
spacing : str, {'lin', 'log'}
Spacing type of the radii vector (log to better sample fine particles)
UND : int
Underdetermination factor for the multifrequency inversion system.
"""
def __init__(self, particle, water, ori=-np.pi, freq=None, tag_M=None,
gen_model=None, mode_depl=None, am=30e-6, aM=600e-6,
dsigma=None, spacing='log', UND=4):
self.particle = particle
self.water = water
self.ori = ori
self.scattering = Parameters({})
self.mode_depl = mode_depl
self.cw = self.water.mean_TSCR(self.water.paramw.cw, self.mode_depl)
self.rho = self.water.mean_TSCR(self.water.paramw.rho, self.mode_depl)
self.scattering.g = self.particle.param.rhos / self.rho
self.scattering.h = self.particle.param.c / self.cw
self.HR_tag = 0
if tag_M is None:
pass
else:
self.scattering.freq = freq
self.scattering.tag_M = tag_M
if tag_M == 'mono':
self.scattering.freq = np.unique(freq)
if len(self.scattering.freq) > 1:
raise InstrAmbiguity('More than one frequency ' +
'for one instrument supposed ' +
'to use a sole frequency ')
elif tag_M == 'multi':
self.scattering.a_UND, self.scattering.ka_UND,\
self.scattering.f_infUND, self.scattering.c_infUND,\
self.scattering.aij = self.RES_model(am, aM, spacing, UND)
self.scattering.a, self.scattering.ka, self.scattering.f_inf,\
self.scattering.c_inf = self.RES_model(
am, aM, 'log', na=100)[0: -1]
if 'DWBA' in self.particle.name:
self.dsigma = None
# self.scattering.a_HR = self.scattering.a
# self.scattering.ka_HR = self.scattering.ka
# self.scattering.f_infHR = self.scattering.f_inf
# self.scattering.c_infHR = self.scattering.c_inf
self.scattering.a_HR,\
self.scattering.ka_HR,\
self.scattering.f_infHR,\
self.scattering.c_infHR = self.RES_model(am,
aM,
'log',
na=500)[0: -1]
else:
self.scattering.a_HR,\
self.scattering.ka_HR,\
self.scattering.f_infHR,\
self.scattering.c_infHR = self.RES_model(am, aM, 'log')[0: -1]
if dsigma is not None:
self.scattering.dsigma = dsigma
self.scattering.ka0, self.scattering.f_0,\
self.scattering.c_0 = self.ensemble_average_model(
self.scattering.ka, dsigma=self.scattering.dsigma, na=1e4)
self.scattering.ka0_HR, self.scattering.f_0HR,\
self.scattering.c_0HR = self.ensemble_average_model(
self.scattering.ka_HR,
dsigma=self.scattering.dsigma, na=1e4)
a0_HR = []
a0 = []
k_tmp = 2 * np.pi * np.array(self.scattering.freq) / self.cw
[a0_HR.append(a / b)
for a, b in zip(self.scattering.ka0_HR, k_tmp)]
[a0.append(a/b)
for a, b in zip(self.scattering.ka0, k_tmp)]
self.scattering.a0_HR = a0_HR[0]
self.scattering.a0 = a0[0]
self.tmp_del()
[docs] def tmp_del(self):
"""
Deletes temporary model information
"""
try:
del self.scattering.atmp, self.scattering.katmp, \
self.scattering.f_inftmp, self.scattering.c_inftmp
except AttributeError:
pass
[docs] def ensemble_average_model(self, xx, dsigma, freq=[1e6], na=1e4,
am=1e-6, aM=0.01, space='lin'):
"""
Computes the ensemble averaged model over the size distribution if the
dsigma input of the class is not None.
Parameters
----------
xx : float
Taget :math:`ka_0` vector, with :math:`a_0` the mean radius, similar
to the ``ka`` variable
dsigma : float
Distribution width
freq : float
Reference frequency (single frequency case)
na : float (numpy.array or list)
Number of integration points for the size distribution
am, aM : object
Min and max radius to be considered for the size distribution
space : str
Spacing of the new ka vector
Returns
-------
axlng : float
Taget :math:`ka_0` vector, with :math:`a_0` the mean radius
aflng : float
Ensemble averaged form function, function of :math:`ka_0`
achilng : float
Ensemble averaged total scattering cross-section , function of
:math:`ka_0`
"""
a, ka, fe, chi = self.RES_model(am, aM, space=space,
na=na, freq=freq)[0: -1]
self.tmp_del()
nf, nz = np.shape(xx)
axlng = []
aflng = []
achilng = []
nz = np.shape(xx)[1]
x = ka.T
dxT = np.tile(np.diff(ka).T, (1, nf))
chiT = np.tile(chi[:, 0:-1].T, (1, nf))
feT = np.tile(fe[:, 0:-1].T, (1, nf))
xT = np.tile(ka.T[0: -1], (1, nf))
# lognormal distribution
sigma = [dsigma * xo.T for xo in list(xx.T)]
mu = [np.log((xo ** 2) / np.sqrt(xo ** 2 + sigma ** 2))
for xo, sigma in zip(xx.T, sigma)]
sigman = [np.sqrt(np.log((sigma / xo) ** 2 + 1))
for xo, sigma in zip(xx.T, sigma)]
plognorm = [(1 / (x[0:-1] * s * np.sqrt(2 * np.pi))) *
np.exp(-((np.log(x[0: -1]) - m) * (np.log(x[0: -1]) - m)) /
(2 * s ** 2)) for s, m in zip(sigman, mu)]
# new ka0 vector, equal to ka in this case
[axlng.append(np.sum(dxT * xT * p, axis=0)) for p in plognorm]
# integration over the form function and tot. scatt. cross section
ax3 = [np.sum(dxT * xT ** 3 * p, axis=0) for p in plognorm]
axf = [np.sum(dxT * xT ** 2 * feT ** 2 * p, axis=0) for p in plognorm]
axchi = [np.sum(xT ** 2 * chiT * dxT * p, axis=0) for p in plognorm]
# final result building
[aflng.append(np.sqrt(a * b * (1 / c)))
for a, b, c in zip(axlng, axf, ax3)]
[achilng.append(a * b * (1 / c)) for a, b, c in zip(axlng, axchi, ax3)]
return np.array(axlng).T, np.array(aflng).T, np.array(achilng).T
[docs] def RES_model(self, am, aM, space='log', UND=None, na=None, freq=None):
"""
Calculates the scattering model by calling the right method
Parameters
----------
am, aM : float
Taget :math:`ka_0` vector, with :math:`a_0` the mean radius, similar
to the math:`ka` variable
space : float
Distribution width
UND : float
Reference frequency (single frequency case)
na : float (numpy.array or list)
Number of integration points for the size distribution
freq : object
Returns
-------
atmp : float
Radius vector
katmp : float
:math:`ka` vector
f_inftmp : float
Intrinsic form function, function of :math:`ka`
c_inftmp : float
Intrinsic total scattering cross-section, function of :math:`ka`
aij : float
Backscattering cross-section, function of :math:`ka`
"""
if UND is None:
if na is None:
na = 1000
else:
na = np.int(na)
else:
try:
na = len(self.scattering.freq) * UND
except TypeError:
na = 1
if freq is None:
freq = self.scattering.freq
if space == 'log':
am = np.log10(am)
aM = np.log10(aM)
print('uiui1111')
if self.particle.mtype == 'THQ':
if 'DWBA' in self.particle.name:
DWBA(self, freq, np.logspace(am, aM, na), self.ori)
else:
Thq(self, freq, np.logspace(am, aM, na), self.ori)
elif self.particle.mtype == 'EMP':
Emp(self, freq, np.logspace(am, aM, na), self.ori)
elif self.particle.mtype == 'HP':
Hp(self, freq, np.logspace(am, aM, na))
print('uiui')
elif space == 'lin':
if self.particle.mtype == 'THQ':
Thq(self, freq, np.linspace(am, aM, na), self.ori)
elif self.particle.mtype == 'EMP':
Emp(self, freq, np.linspace(am, aM, na), self.ori)
elif self.particle.mtype == 'HP':
Hp(self, freq, np.linspace(am, aM, na))
print('lala {}'.format(np.shape(self.scattering.f_inftmp)))
if 'DWBA' in self.particle.name:
aij = np.multiply((self.particle.param.shape.L_a *
self.scattering.atmp) ** 2,
np.abs(self.scattering.f_inftmp)**2) / 16
elif self.particle.shape in ['CYL', 'bCYL', 'PRO']:
aij = np.multiply((self.particle.param.L_a *
self.scattering.atmp) ** 2,
np.abs(self.scattering.f_inftmp)**2) / 16
else:
aij = np.multiply(self.scattering.atmp ** 2,
np.abs(self.scattering.f_inftmp)**2) / 4
return self.scattering.atmp, self.scattering.katmp,\
self.scattering.f_inftmp, self.scattering.c_inftmp, aij
[docs] def scatt_func_display(self, preset='basic', *args, **kwargs):
"""
Display the scattering model. Typical usage::
# default view as a (1row, 2columns) subplot
>>>S.scatt_func_display()
# view as a (2row, 1column) subplot
>>>S.scatt_func_display(subplot=[2,1])
"""
if 'subplot' not in kwargs.keys():
kwargs.update({'subplot': (1, 2)})
if preset == 'basic':
if self.scattering.tag_M == 'mono':
Display(self.scattering.ka_HR,
np.abs(self.scattering.f_infHR),
self.scattering.ka_HR,
self.scattering.c_infHR,
subplot=kwargs['subplot'],
plt_type=['loglog', 'loglog'],
xlabel=['ka', 'ka'],
color=['r', 'b', 'r', 'b', 'r', 'b'])
elif self.scattering.tag_M == 'multi':
Display(self.scattering.ka_HR[0, :],
np.abs(self.scattering.f_infHR[0, :]),
self.scattering.ka_HR[0, :],
self.scattering.c_infHR[0, :],
subplot=kwargs['subplot'],
plt_type=['loglog', 'loglog'],
xlabel=['ka', 'ka'],
color=['r', 'b', 'r', 'b', 'r', 'b'])
# add a functionality in Display class to close all figures...
def name_model(self):
return self.particle._descrip + '_' +\
self.particle.mtype + '_' + self.particle.shape
[docs] def calc_f(self, func, ka):
"""
This method is called by each class called by the present Scattering
class to compute the model for each ka vectors in case of
multifrequency instruments.
Parameters
----------
func : class method
method used to compute the desired scattering model
ka : float
:math:`ka` vector
"""
self.scattering.ind = 1
toto = []
[toto.append(func(self.particle.param, ka[u, :]))
for u in range(0, np.size(ka, 0))]
a, b = zip(*toto)
return list(a), list(b)