#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Distorted Wave Born Approximation model (:mod:`hydrac.model.scattering.dwba`)
=============================================================================
.. autoclass:: DWBA
:members:
:private-members:
"""
from scipy import special
from scipy.integrate import quad
import numpy as np
from transonic import boost, jit
class Memoize:
def __init__(self, f):
self.f = f
self.memo = {}
def __call__(self, *args):
if args not in self.memo:
self.memo[args] = self.f(*args)
# Warning: You may wish to do a deepcopy here if returning objects
return self.memo[args]
cos = Memoize(np.cos)
jv = Memoize(special.jv)
A1 = "float64[:]"
A2 = "float64[:,:]"
@boost
def DWBAintegrand(s: float, k1: A1,rpos1: A2,rpos2: A2,a1: A1, a2: A1, betatilt: A1, g1: float, h1: float, g2: float, h2: float):
#def DWBAintegrand(self, s, k1, rpos1, rpos2, a1, a2, betatilt, g1, h1, g2, h2):
"""
Integration function. Needs to be simplifyed in the future to make it
low level function.
"""
i = np.complex(0, 1)
rposx = s * (rpos2[0, 0] - rpos1[0, 0]) + rpos1[0, 0]
rposy = s * (rpos2[1, 0] - rpos1[1, 0]) + rpos1[1, 0]
rposz = s * (rpos2[2, 0] - rpos1[2, 0]) + rpos1[2, 0]
rpos = np.array([rposx, rposy, rposz]).T
a = s * (a2 - a1) + a1
g = s * (g2 - g1) + g1
h = s * (h2 - h1) + h1
gamgam = 1 / (g * h * h) + 1 / g - 2
if np.abs(np.abs(betatilt) - np.pi / 2) < 1e-10:
# limiting case for end-on incidence
bessy = np.linalg.norm(k1) * a / h
else:
bessy = special.jv(1, 2 * np.linalg.norm(k1) * a / h *
np.cos(betatilt)) / np.cos(betatilt)
return np.linalg.norm(k1) / 4 * gamgam * np.exp(2 * i * np.dot(k1.T, rpos) / h) * a * bessy * np.linalg.norm(rpos2 - rpos1)
class ModelNonExistent(Exception):
pass
#@boost
[docs]class DWBA(object):
"""DWBA class.
This class is meant to describe the scattering properties of fluid-like
elongated objects with circular cross secton whatever their shape and
physical properties.
Parameters
----------
obj : object
Scattering object from :class:
`hydrac.model.scattering.scattering.Scattering`
f : float (numpy.array)
Instrument frequencies
a : float (numpy.array)
Size vector
theta : float
Object orientation in the indicdent field
"""
def __init__(self, obj, f, a, theta):
self.obj = obj
self.f = f
self.obj.scattering.atmp = a
self.obj.scattering.theta = theta
self.obj.scattering.katmp = 2 * np.pi * np.outer(f, a) / self.obj.cw
model_index = self.obj.particle.param.simu.model_index
if self.obj.particle.nature == 'ORG':
if self.obj.particle.param.simu.model == 1:
if model_index == 1:
self.obj.scattering.f_inftmp,\
self.obj.scattering.c_inftmp = np.array(
self.obj.calc_f(self.dwba_bentcylinder,
self.obj.scattering.katmp))
elif model_index == 2:
pass
elif model_index == 3:
self.obj.scattering.f_inftmp,\
self.obj.scattering.c_inftmp = np.array(
self.obj.calc_f(self.dwba_ellispoid,
self.obj.scattering.katmp))
else:
raise ModelNonExistent('Model parameters cannot be' +
'interpreted : model_index ' +
'should be 1 or 2 or 3 but is' +
'here equal to {}' +
''.format(model_index))
else:
pass
[docs] def dwba_ellispoid(self, param, ka):
"""
This method computes the DWBA model for ellisposidal fluid-like
particles. This model is valid regardless of the particle body
properties yet hasn't been implemented here for simplicity.
Note that here, a is the semi-minor axis of the animal.
Parameters
----------
param : object
Modified dictionnary from :mod:`hydrac.util.parameters` containing
the particles parameters.
ka : float (numpy.ndarray)
vector of products of wave number and particle size; where the model
is computed.
Returns
-------
f_inf : float (numpy.array)
Model form function :math:`f_{\infty}`
c_inf : float (numpy.array)
Model total scattering cross-secton :math:`\\chi_{\infty}`
"""
# shape parameters
e_ca = param.shape.L_a / 2
e_cb = param.shape.L_b / 2
e_ba = e_cb / e_ca
param.shape.L_a_ave = 2 * np.sqrt(e_ca * e_cb)
# simulation parameter
# nb of ka point
n = len(ka)
# different incident angle
ang = np.array(param.simu.ang)
try:
m = len(ang)
except TypeError:
m = 1
th = ang * np.pi / 180
# azimuthal angle in radians
phi = param.shape.phi * np.pi / 180
ka0 = np.atleast_2d(ka).T
coef = np.atleast_2d((ka0 * ka0) * e_ba * e_ca / param.shape.L_a)
g = param.phy.g1
h = param.phy.hL
Cb = (1-g * h * h) / (g * h * h) - (g-1) / g
mux = np.atleast_2d(2 * np.sin(th) * np.cos(phi))
muy = np.atleast_2d(2 * np.sin(th) * np.sin(phi))
muz = np.atleast_2d(2 * np.cos(th))
# prepare matrices to speed up computations
Ka = np.tile(ka0, (1, m))
Coef = np.tile(coef, (1, m))
Mux = np.tile(mux, (n, 1))
Muy = np.tile(muy, (n, 1))
Muz = np.tile(muz, (n, 1))
Mu = np.sqrt(Mux**2 + e_ba * e_ba * Muy**2 + e_ca * e_ca * Muz**2)
Arg = (Ka / h) * Mu + np.spacing(1)
if np.size(Ka[:, 0]) == 1000:
print('rui{}'.format(Arg[500, :]))
tmp = []
[tmp.append(special.spherical_jn(1, Arg[:, j]) / Arg[:, j])
for j in range(m)]
tmp = np.array(tmp).T
f = Cb * Coef * tmp
if param.simu.aveA_flag == 1:
orient_ave_para = [param.shape.ang, param.shape.dang]
f1 = np.atleast_2d(self.orient_ave(ang,
f,
param.simu.aveA_PDF,
orient_ave_para))
else:
f1 = np.sqrt(f * np.conj(f))
# L_simu = param.shape.L_a * self.obj.scattering.atmp
f2 = np.squeeze(f1) * 4
q = np.zeros(np.shape(f2), dtype=np.float)
return f2, q
[docs] def dwba_bentcylinder(self, param, ka):
"""
This method computes the DWBA model for elongated fluid-like
particles. This model is valid regardless of the particle body
properties yet if it is possible to make the body properties vary,
this requires additionnal method to input the body properties along
the lengthwise axis of the animal. Similarly, the lengthwise radius
can vary, yet this requires additional method to input these values.
Maybe think of inputing this information using a txt file in the
future.
Note that here, a is the radius of the circular cross section, that can
vary along the length of the animal.
This method performs an integration of the scattering form function
over the length of the animal. This requires a lot of computing time
and the integration scheme (:func:`DWBAintegrand`) needs to be fastend
in the future.
Parameters
----------
param : object
Modified dictionnary from :mod:`hydrac.util.parameters` containing
the particles parameters.
ka : float (numpy.ndarray)
vector of products of wave number and particle size; where the model
is computed.
Returns
-------
f_inf : float (numpy.array)
Model form function :math:`f_{\infty}`
c_inf : float (numpy.array)
Model total scattering cross-secton :math:`\\chi_{\infty}`
"""
L_a = param.shape.L_a
param.simu.ni = param.simu.min_ni
# of integration points
n_int = np.int(param.simu.ni)
# different incident angle
ang = np.array(param.simu.ang)
ka0 = np.atleast_2d(ka).T
a0 = np.atleast_2d(self.obj.scattering.atmp).T
k = np.mean(ka0 / a0)
g = param.phy.g1
h = param.phy.hL
rho_L = param.shape.rho_L
order = param.shape.order
gamma = 0.5 / rho_L
ratio = 2*rho_L
z = np.linspace(-1, 1, n_int)
taper = np.atleast_2d(np.sqrt(1 - (z**order)))
z = np.sin(gamma) * z
x = (1-np.sqrt(1 - z**2))
x = ratio * x
z = ratio * z
a0_t = np.dot(taper.T, a0.T)
gvec = np.zeros((n_int, 1)) + g
hvec = np.zeros((n_int, 1)) + h
rposnorm = np.array([z, x, x*0]).T
L_simu = L_a * a0
f = np.zeros((len(ang), len(a0)))
for i in range(0, len(a0)):
avec = np.atleast_2d(a0_t[:, i]).T
L = L_a * a0[i, 0]
rposvec = rposnorm * L / 2
fbsvector = []
for phi in ang:
phirad = phi * np.pi / 180
k1 = k * np.array([np.cos(phirad), np.sin(phirad), 0]).T
fbs = self.findfbs(avec, rposvec, gvec, hvec, k1)
fbsvector.append(fbs)
f[:, i] = np.squeeze(np.array(np.abs(fbsvector))**2)
f = f.T / L_simu / L_simu
# f = f.T
if param.simu.aveA_flag == 1:
orient_ave_para = [param.shape.ang, param.shape.dang]
f1 = np.atleast_2d(self.orient_ave(ang,
f,
param.simu.aveA_PDF,
orient_ave_para))
else:
f1 = np.sqrt(f * np.conj(f))
f2 = np.squeeze(np.sqrt(f1)) * 4
q = np.zeros(np.shape(f2), dtype=np.float)
return np.abs(f2), q
[docs] def findfbs(self, avec, rposvec, gvec, hvec, k1):
"""
Calculates the form function over the length of the animal for one
particular ka value and one particular angle of orientation since the
form function may be averaged over angle distribution.
"""
(n, m) = np.shape(avec)
fbs = 0
for ii in range(0, n-1):
a1 = avec[ii]
a2 = avec[ii+1]
rpos1 = np.atleast_2d(rposvec[ii, :]).T.copy()
rpos2 = np.atleast_2d(rposvec[ii+1, :]).T.copy()
g1 = gvec[ii, 0]
g2 = gvec[ii+1, 0]
h1 = hvec[ii, 0]
h2 = hvec[ii+1, 0]
alphatilt = np.arccos(np.dot(k1.T, (rpos2 - rpos1)) /
(np.linalg.norm(k1) *
np.linalg.norm(rpos2 - rpos1)))
betatilt = np.abs(alphatilt - np.pi / 2)
fbs = fbs + self.complex_quadrature(lambda s:
DWBAintegrand(s,
k1,
rpos1,
rpos2,
a1,
a2,
betatilt,
g1,
h1,
g2,
h2), 0, 1)[0]
#
# fbs = fbs + self.complex_quadrature(lambda s:
# self.DWBAintegrand(s,
# k1,
# rpos1,
# rpos2,
# a1,
# a2,
# betatilt,
# g1,
# h1,
# g2,
# h2), 0, 1)[0]
return fbs
## @boost
## def DWBAintegrand(self, s: float, k1: A1,rpos1: A2,rpos2: A2,a1: A1, a2: A1, betatilt: A1, g1: float, h1: float, g2: float, h2: float):
# @jit
# def DWBAintegrand(self, s, k1, rpos1, rpos2, a1, a2, betatilt, g1, h1, g2, h2):
# """
# Integration function. Needs to be simplifyed in the future to make it
# low level function.
# """
# i = np.complex(0, 1)
# rposx = s * (rpos2[0, 0] - rpos1[0, 0]) + rpos1[0, 0]
# rposy = s * (rpos2[1, 0] - rpos1[1, 0]) + rpos1[1, 0]
# rposz = s * (rpos2[2, 0] - rpos1[2, 0]) + rpos1[2, 0]
# rpos = np.array([rposx, rposy, rposz]).T
# a = s * (a2 - a1) + a1
# g = s * (g2 - g1) + g1
# h = s * (h2 - h1) + h1
#
# gamgam = 1 / (g * h * h) + 1 / g - 2
# if np.abs(np.abs(betatilt) - np.pi / 2) < 1e-10:
# # limiting case for end-on incidence
# bessy = np.linalg.norm(k1) * a / h
# else:
# bessy = special.jv(1, 2 * np.linalg.norm(k1) * a / h *
# np.cos(betatilt)) / np.cos(betatilt)
#
# return np.linalg.norm(k1) / 4 * gamgam * np.exp(2 * i * np.dot(k1.T, rpos) / h) * a * bessy * np.linalg.norm(rpos2 - rpos1)
[docs] def complex_quadrature(self, func, a, b, **kwargs):
"""
Calculates the integral of a complex function, since scipy.quad
doesn't allow it.
Parameters
----------
func : python function
Integrand function
a, b : float (numpy.ndarray)
Boundaries for intergation
Returns
-------
f_int : complex128
Integrated function
"""
def real_func(x):
return np.real(func(x))
def imag_func(x):
return np.imag(func(x))
real_integral = quad(real_func, a, b, **kwargs)
imag_integral = quad(imag_func, a, b, **kwargs)
return (real_integral[0] + 1j*imag_integral[0],
real_integral[1:], imag_integral[1:])
[docs] def buildpos(self, param):
"""
Method not used here but will be to design the animage shape, hence the
radius variation along the lengthwise axis of the animal.
"""
n_int = np.int(param.simu.ni)
order = np.int(param.shape.order)
rho_L = param.shape.rho_L
L_a = param.shape.L_a
# uniformly bent cylinder and regularly tapered
gamma = 0.5 / rho_L
ratio = 2 * rho_L
z = np.linspace(-1, 1, n_int)
taper = np.atleast_2d(np.sqrt(1 - np.linspace(-1, 1, n_int)**order))
# normalized by rho - radius of curvature
z = np.sin(gamma) * z
x = (1 - np.sqrt(1 - z * z))
# normalized by L/2
x = np.atleast_2d(ratio * x).T
z = np.atleast_2d(ratio * z).T
taper2 = np.reshape(taper, (1, n_int))
th_tilt = np.zeros((n_int, 1))
r_pos = np.sqrt(x * x + z * z)
gamma_t = np.arctan2(z, x)
dx = np.diff(x[:, 0]) + np.spacing(1)
dz = np.diff(z[:, 0])
alpha_t = np.zeros(np.shape(gamma_t))
alpha_t[0:-1, 0] = np.arctan(dz / dx)
alpha_t[-1, 0] = np.arctan(dz[n_int-2] / dx[n_int-2])
indx1 = np.where(alpha_t < 0)[0]
if len(indx1) > 0:
th_tilt[indx1] = alpha_t[indx1] + np.pi / 2
indx2 = np.where(alpha_t >= 0)[0]
if len(indx2) > 0:
th_tilt[indx2] = alpha_t[indx2] - np.pi / 2
dr1 = np.sqrt(dx * dx + dz * dz)
dr = np.zeros(np.shape(alpha_t))
dr[0, 0] = dr1[0]
dr[1:, 0] = dr1
return r_pos, th_tilt, dr, gamma_t, taper, x, z
[docs] def orient_ave(self, ang, f, pdf_type, ori_para):
"""
Performs average over the angle distribution.
Parameters
----------
ang : float
Angles of orientation for which the model has been computed
f : float (numpy.ndarray)
Form function
pdf_type : str, {uniform, gaussian}
Type of pdf to consider, input by the user during model parameters
definition in :mod:`hydrac.model.particle`
ori_para : float
Standard deviation and mean values
Returns
-------
outy : float (numpy.array)
Averaged model
"""
# n = points in freq., m = points in orientation
n, m = np.shape(f)
if m == 1:
outy = np.sqrt(f * np.conj(f))
return outy
if pdf_type == 1:
pdf = np.ones(m) / m
else:
dang = ang[1] - ang[0]
angm = ori_para[0]
angstd = ori_para[1]
pdf = dang * np.exp(-0.5 * (ang - angm) * (ang-angm) /
(angstd * angstd)
) / (np.sqrt(2 * np.pi) * angstd)
outy = np.sqrt(np.dot((f * np.conj(f)), pdf))
return outy