#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Thoery based model (:mod:`hydrac.model.scattering.thq`)
=======================================================
.. autoclass:: Thq
:members:
:private-members:
"""
from scipy import special
import numpy as np
[docs]class Thq(object):
"""Theoretical models class.
This class deals with analytical exact models from modal series
solutions. It handles spherical and straight cylinders shapes, of
organic or mineral nature.
The class calls the method
:func:`hydrac.model.scattering.scattering.Scattering.calc_f` to compute
the model for each :math:`ka` vectors generated by each instrument
frequency.
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
Angle between the particle principal axis (ex. cylinder axis) and the
incident pressure field
"""
def __init__(self, obj, f, a, theta):
self.obj = obj
self.obj.scattering.atmp = a
self.obj.scattering.theta = theta
self.obj.scattering.katmp = 2 * np.pi * np.outer(f, a) / self.obj.cw
if self.obj.particle.nature == 'ORG':
self.obj.scattering.f_inftmp,\
self.obj.scattering.c_inftmp = np.array(
self.obj.calc_f(self.thq_fluid_func,
self.obj.scattering.katmp))
else:
if self.obj.particle.shape == 'SPH':
self.obj.scattering.f_inftmp,\
self.obj.scattering.c_inftmp = np.array(
self.obj.calc_f(self.thq_elast_func,
self.obj.scattering.katmp))
elif self.obj.particle.shape == 'CYL':
self.obj.scattering.f_inftmp,\
self.obj.scattering.c_inftmp = np.array(
self.obj.calc_f(self.thq_cyl_elast_func,
self.obj.scattering.katmp))
[docs] def thq_fluid_func(self, param, ka):
""" Method to compute the modal series solution for a perfect fluid
sphere (See Stanton 1988).
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}`
"""
p = np.zeros((len(ka), ), dtype=np.float)
q = np.zeros((len(ka), ), dtype=np.float)
g = self.obj.scattering.g
h = self.obj.scattering.h
for tutu in range(0, len(ka)):
cc = []
j = [special.spherical_jn(merci, ka[tutu])
for merci in range(0, int(param['mod_num']) + 1)]
jh = [special.spherical_jn(merci, ka[tutu] / h)
for merci in range(0, int(param['mod_num']) + 1)]
y = [special.spherical_yn(merci, ka[tutu])
for merci in range(0, int(param['mod_num']) + 1)]
cc.append(jh[1] * y[0] / (j[1] * jh[0]) - g * h * y[1] / j[1])
cc[0] = cc[0] / (jh[1] * j[0] / (j[1] * jh[0]) - g * h)
alpha = j[0] - 2 * j[2]
alphah = jh[0] - 2 * (jh[2])
beta = y[0] - 2 * y[2]
cc.append(alphah * y[1] / (alpha * jh[1]) - g * h * beta / alpha)
cc[1] = cc[1] / (alphah * j[1] / (alpha * jh[1]) - g * h)
for n in range(3, np.int(param['mod_num']) + 1):
m = n-1
alpha = m * j[m - 1] - (m + 1) * j[m - 1 + 2]
alphah = m * jh[m - 1] - (m + 1) * jh[m - 1 + 2]
beta = m * y[m - 1] - (m + 1) * y[m - 1 + 2]
cc.append(alphah * y[m - 1 + 1] / (alpha * jh[m - 1 + 1]) -
g * h * beta / alpha)
cc[n - 1] = cc[n - 1] / (alphah * j[m - 1 + 1] /
(alpha * jh[m - 1 + 1]) - g * h)
sumr = 0
sumi = 0
p0 = 1
for n in range(1, np.int(param['mod_num']) + 1):
m = n-1
con = ((-1) ** m) * (2 * m + 1) / (1 + cc[m] * cc[m])
sumr = sumr + con
sumi = sumi + con * cc[m]
p0 = np.sqrt(sumr * sumr + sumi * sumi)
p[tutu] = p0
p[tutu] = 2 * p[tutu] / ka[tutu]
del cc
print(m)
return p, q
[docs] def thq_elast_func(self, param, ka):
""" Method to compute the modal series solution for a perfect elastic
sphere (See Faran 1951).
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}`
"""
cp = param.cp
cs = param.cs
rhos = param.rhos
# Ligne d'apres permet computation du modele
# pour chaque données de temp salinité...
# ce qui peut etre très intéressant pour la suite :
# rho=self.obj.water.paramw.rho;c=self.obj.cw
rho = self.obj.rho
c = self.obj.cw
x = ka
xt = x * c / cs
xl = x * c / cp
L = np.zeros(len(x), dtype=np.complex)
CHI = np.zeros(len(x), dtype=np.complex)
for zz in range(0, len(x)):
jx = [special.spherical_jn(merci, x[zz])
for merci in range(0, int(param['mod_num']) + 1)]
jxD = [special.spherical_jn(merci, x[zz], derivative=True)
for merci in range(0, int(param['mod_num']) + 1)]
nx = [special.spherical_yn(merci, x[zz])
for merci in range(0, int(param['mod_num']) + 1)]
nxD = [special.spherical_yn(merci, x[zz], derivative=True)
for merci in range(0, int(param['mod_num']) + 1)]
jxl = [special.spherical_jn(merci, xl[zz])
for merci in range(0, int(param['mod_num']) + 1)]
jxDl = [special.spherical_jn(merci, xl[zz], derivative=True)
for merci in range(0, int(param['mod_num']) + 1)]
jxt = [special.spherical_jn(merci, xt[zz])
for merci in range(0, int(param['mod_num']) + 1)]
jxDt = [special.spherical_jn(merci, xt[zz], derivative=True)
for merci in range(0, int(param['mod_num']) + 1)]
for loop in range(0, param['mod_num']):
beta1 = -xt[zz] ** 2 * (rho / rhos) * jx[loop]
beta2 = x[zz] * jxD[loop]
alpha11 = xt[zz] ** 2 * (rho / rhos) * (
jx[loop] + np.complex(0, 1) * nx[loop])
alpha21 = -x[zz] * (jxD[loop] + np.complex(0, 1) * nxD[loop])
alpha12 = (2 * loop * (loop + 1) - xt[zz] ** 2
) * jxl[loop] - 4 * xl[zz] * jxDl[loop]
alpha22 = xl[zz] * jxDl[loop]
alpha32 = 2 * (jxl[loop] - xl[zz] * jxDl[loop])
alpha13 = (2 * loop * (loop + 1)) * (
xt[zz] * jxDt[loop] - jxt[loop])
alpha23 = loop * (loop + 1) * jxt[loop]
alpha33 = 2 * xt[zz] * jxDt[loop] + (
xt[zz] ** 2 - 2 * (loop + 1) * loop + 2) * jxt[loop]
bn = np.linalg.det([[beta1, alpha12, alpha13],
[beta2, alpha22, alpha23],
[0, alpha32, alpha33]]
) / np.linalg.det(
[[alpha11, alpha12, alpha13],
[alpha21, alpha22, alpha23],
[0, alpha32, alpha33]])
L[zz] = L[zz] + (2 / (np.complex(0, 1) * x[zz])) * (
((-1) ** loop) * (2 * loop + 1) * bn)
CHI[zz] = CHI[zz] + (2 / (x[zz] ** 2)) * (
2 * loop + 1) * np.abs(bn) ** 2
print('fhu {}'.format(np.shape(L)))
return np.abs(L), CHI
[docs] def thq_cyl_elast_func(self, param, ka):
""" Method to compute the modal series solution for a perfect elastic
cylinder of finite length (See Stanton 1988b).
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}`
"""
cp = param.cp
cs = param.cs
rhos = param.rhos
# Ligne d'apres permet computation du modele
# pour chaque données de temp salinité...
# ce qui peut etre très intéressant pour la suite :
# rho=self.obj.water.paramw.rho;c=self.obj.cw
rho = self.obj.rho
c = self.obj.cw
x = ka
xt = x * c / cs
xl = x * c / cp
L = np.zeros(len(x), dtype=np.complex)
CHI = np.zeros(len(x), dtype=np.complex)
theta = self.obj.scattering.theta * np.pi / 180
theta = 0.5 * np.pi
L_a = param.L_a
Long = L_a * self.obj.scattering.atmp
delta = x * Long * np.cos(theta)
for zz in range(0, len(x)):
jx = [special.jn(merci, x[zz]*np.sin(theta))
for merci in range(0, int(param['mod_num']) + 1)]
jxD = [special.jvp(merci, x[zz]*np.sin(theta))
for merci in range(0, int(param['mod_num']) + 1)]
nx = [special.yn(merci, x[zz]*np.sin(theta))
for merci in range(0, int(param['mod_num']) + 1)]
nxD = [special.yvp(merci, x[zz]*np.sin(theta))
for merci in range(0, int(param['mod_num']) + 1)]
jxl = [special.jn(merci, xl[zz]*np.sin(theta))
for merci in range(0, int(param['mod_num']) + 1)]
jxDl = [special.jvp(merci, xl[zz]*np.sin(theta))
for merci in range(0, int(param['mod_num']) + 1)]
jxt = [special.jn(merci, xt[zz]*np.sin(theta))
for merci in range(0, int(param['mod_num']) + 1)]
jxDt = [special.jvp(merci, xt[zz]*np.sin(theta))
for merci in range(0, int(param['mod_num']) + 1)]
for loop in range(0, param['mod_num']):
dm = -jx[loop] / nx[loop]
am = -x[zz] * jxD[loop] / jx[loop]
betam = -x[zz] * nxD[loop] / nx[loop]
am1 = -xl[zz] * jxDl[loop] / jxl[loop]
am2 = -xt[zz] * jxDt[loop] / jxt[loop]
dz1 = ((-xt[zz]**2) / 2) * (((am1 / (am1+1)) - (loop**2) /
(am2+loop**2-0.5*xt[zz]**2)) /
(((am1+loop**2-0.5*xt[zz]**2) /
(am1+1))
- ((loop**2) * (1 + am2) /
(am2 + loop**2 -
0.5 * xt[zz]**2))))
phim = (-rho / rhos) * dz1
tanetam = dm * (phim + am) / (phim + betam)
etam = np.arctan(tanetam)
if loop == 0:
bm = np.sin(etam) * np.exp(-np.complex(0, 1) * etam)
else:
bm = 2 * np.sin(etam) * np.exp(-np.complex(0, 1) * etam)
L[zz] = L[zz] -\
(Long[zz] / np.pi) *\
np.sin(delta[zz] / delta[zz]) *\
bm * np.cos((loop)*np.pi)
return 4 * np.abs(L) / Long, CHI