diff --git a/acoustics/__init__.py b/acoustics/__init__.py index 0de5434b..b0306d91 100644 --- a/acoustics/__init__.py +++ b/acoustics/__init__.py @@ -24,7 +24,6 @@ import acoustics.reflection import acoustics.room import acoustics.signal -import acoustics.turbulence #import acoustics.utils import acoustics.weighting diff --git a/acoustics/_turbulence.py b/acoustics/_turbulence.py deleted file mode 100644 index b4bdada5..00000000 --- a/acoustics/_turbulence.py +++ /dev/null @@ -1,558 +0,0 @@ -""" -This file contains all abstract base classes related to atmospheric turbulence. -""" - -import numpy as np -import matplotlib.pyplot as plt -import abc -import six - -from scipy.special import gamma -from scipy.special import iv as bessel # Modified Bessel function of the first kind. - - -class abstractstaticmethod(staticmethod): - """ - Abstract static method. - """ - __slots__ = () - def __init__(self, function): - super(abstractstaticmethod, self).__init__(function) - function.__isabstractmethod__ = True - __isabstractmethod__ = True - -@six.add_metaclass(abc.ABCMeta) -class Spectrum(object): - """ - Abstract turbulence spectrum. - """ - - wavenumber_resolution = None - """ - Wavenumber resolution - """ - - max_mode_order = None - """ - Maximum amount of modes to consider. - """ - - - #_required_attributes = ['x', 'y', 'z', 'wavenumber_resolution', 'spatial_resolution', 'max_mode_order'] - _required_attributes = ['wavenumber_resolution', 'max_mode_order'] - - def _construct(self, attributes, *args, **kwargs): - for attr in attributes: - try: - setattr(self, attr, kwargs[attr]) - except KeyError: - raise ValueError(attr) - #raise ValueError("Requires " + attr + ".") - - def __init__(self, *args, **kwargs): - - self._sizes = np.zeros(3) - """ - Length of the field in each direction. - """ - - missing = list() - - for cls in reversed(self.__class__.__mro__): - try: - self._construct(getattr(cls, '_required_attributes'), args, **kwargs) - except AttributeError as error: - pass - except ValueError as error: - missing.append(str(error)) - if missing: - raise ValueError("Missing arguments: " + str(set(missing))) - - @property - def modes(self): - """ - Vector of modes. - """ - return np.arange(0, self.max_mode_order) - - @property - def wavenumber(self): - """ - Vector of wavenumbers. - """ - return self.wavenumber_resolution * self.modes - - @abc.abstractmethod - def mode_amplitude(): - pass - - @abc.abstractmethod - def spectral_density(self): - """ - Spectral density of this object. - """ - pass - - @abstractstaticmethod - def spectral_density_function(): - """ - The spectral density function that is used in this model. - """ - pass - - - ##def correlation(self): - ##""" - ##Correlation for this object. - ##""" - ##return self.correlation_function(self.mu_0, self.a, self.r) - - ##def structure(self): - ##""" - ##Structure for this object. - ##""" - ##return self.structure_function(self.mu_0, self.a, self.r) - - ##def spectral_density(self): - ##""" - ##Spectral density for this object. - ##""" - ##return self.spectral_density_function(self.mu_0, self.a, self.wavenumber) - #return getattr(self, 'spectral_density_function_'+self.ndim+'d')(self.mu_0, self.a, self.wavenumber) - - - -@six.add_metaclass(abc.ABCMeta) -class Spectrum1D(Spectrum): - """ - Abstract class for one-dimensional turbulence spectra. - """ - - NDIM = 1 - """ - Amount of dimensions. - """ - - -@six.add_metaclass(abc.ABCMeta) -class Spectrum2D(Spectrum): - """ - Abstract class for two-dimensional turbulence spectra. - """ - - NDIM = 2 - """ - Amount of dimensions. - """ - - #def __init__(self, *args, **kwargs): - - #super(Spectrum2D, self).__init__(*args, **kwargs) - #attributes = ['r', 'z'] - #self._construct(attributes, args, kwargs) - - - #def __init__(self, wavenumber, max_mode_order, r, z, *args, **kwargs): - #Spectrum.__init__(self, wavenumber, max_mode_order) - - _max_mode_order = None - - _required_attributes = ['plane'] - - - @property - def wavenumber_resolution(self): - return self._wavenumber_resolution - - @wavenumber_resolution.setter - def wavenumber_resolution(self, x): - self._wavenumber_resolution = x - self.randomize() - - @property - def max_mode_order(self): - return self._max_mode_order - - @max_mode_order.setter - def max_mode_order(self, x): - self._max_mode_order = x - self.randomize() - - - #@property - #def plane(self): - #""" - #Tuple indicating the plane that is modelled. - #""" - #return self._sizes.astype(bool) - - - def mode_amplitude(self): - """ - Mode amplitudes :math:`G(\\mathbf{k})`. - - :rtype: A `n`-dimensional array where `n` is equal to the amount of dimensions of `k_n`. - - The mode amplitudes are calculating using - - .. math:: G (\\mathbf{k}_n ) = \\sqrt{4 \\pi \\Delta k F(\\mathbf{k}_n) \\mathbf{k}_n} - - where :math:`\\mathbf{k}_n` are the wavenumber, :math:`\\Delta k` the wavenumber resolution, - and :math:`F` the spectral density. - - See Salomons, below equation J.24. - - """ - n = np.arange(0, self.max_mode_order) - return np.sqrt(4.0 * np.pi * self.wavenumber_resolution * self.spectral_density() * self.wavenumber) - - #def _field(self): - """Numba version??""" - - #r_vector = np.arange(0.0, self.r, self.spatial_resolution) - #z_vector = np.arange(0.0, self.z, self.spatial_resolution) - - #shape = (len(self.r_vector), len(self.z_vector), len(self.wavenumber)) - #mu = np.zeros(shape) - - #for n in range(self.max_mode_order): - #for r in r_vector: - #for z in z_vector: - #for k in self.wavenumber: - #mu += - - #return mu - - - def randomize(self): - """ - Create new random values for :math:`\\theta_n` and :math:`\\alpha_n`. - - :rtype: self - - .. warning:: This function should always be called before :meth:`field` when a new random field should be generated. - - .. note:: This function is called whenever :attr:`max_mode_order` or :attr:`wavenumber_resolution` is changed. - - """ - self.alpha = np.random.random_sample(self.max_mode_order) * np.pi # Create random alpha_n - self.theta = np.random.random_sample(self.max_mode_order) * np.pi # Create random alpha_n - return self - - - def plot_mode_amplitudes(self, filename=None): - """ - Calculate and plot mode amplitudes. - """ - fig = plt.figure() - ax = fig.add_subplot(111) - #ax.set_title("Mode {}".format(n)) - - ax.semilogx(self.wavenumber, self.mode_amplitude()) - ax.set_xlabel(r'$k$ in $\mathrm{m}^{-1}$') - ax.set_ylabel(r'$G$') - ax.grid() - ax.set_title('Mode amplitude as function of wavenumber') - - if filename: - fig.savefig(filename) - else: - fig.show() - - - def plot_structure(self): - """ - Plot the structure function. - """ - pass - - def plot_spectral_density(self, filename=None): - """ - Plot the spectral density. - """ - - fig = plt.figure() - ax = fig.add_subplot(111) - - ax.loglog(self.wavenumber, self.spectral_density()) - ax.set_xlabel(r'$k$ in $\mathrm{m}^{-1}$') - ax.set_ylabel(r'$F$') - ax.grid() - ax.set_title('Spectral density as function of wavenumber') - - - if filename: - fig.savefig(filename) - else: - fig.show() - - - -@six.add_metaclass(abc.ABCMeta) -class Spectrum3D(Spectrum): - """ - Abstract class for one-dimensional turbulence spectra. - """ - - NDIM = 3 - """ - Amount of dimensions. - """ - - -@six.add_metaclass(abc.ABCMeta) -class GaussianTemp(Spectrum): - """ - Abstract class for Gaussian spectrum when only temperature fluctuations are considered. - """ - - _required_attributes = ['a', 'mu_0'] - - - a = None - """ - Characteristic length :math:`a`. - """ - - mu_0 = None - """ - The standard deviation of the refractive-index fluctuation :math:`\\mu` is :math:`\\mu_0`. - """ - - def spectral_density(self): - return self.spectral_density_function(self.wavenumber, self.a, self.mu_0) - - - @staticmethod - def correlation_function(r, a, mu_0): - """ - Correlation function :math:`B(r)`. - - :param r: :math:`r` - :param a: :math:`a` - :param mu_0: :math:`\\mu_0` - - The correlation function is given by - - .. math:: B(r) = \\mu_0^2 \\exp{(-r^2/a^2)} - - See Salomons, equation I.28. - """ - return mu_0**2.0 * np.exp(-r**2.0/a**2.0) - - @staticmethod - def structure_function(r, a, mu_0): - """ - Structure function :math:`D(r)`. - - :param r: :math:`r` - :param a: :math:`a` - :param mu_0: :math:`\\mu_0` - - The structure function is given by - - .. math:: D(r) = 2 \\mu_0^2 \\left[ 1 - \\exp{(-r^2/a^2)} \\right] - - See Salomons, equation I.29. - """ - return 2.0 * mu_0**2.0 * (1.0 - np.exp(-r**2/a**2)) - - -@six.add_metaclass(abc.ABCMeta) -class KolmogorovTemp(object): - """ - Abstract class for Kolmogorov spectrum when only temperature fluctuations are considered. - """ - - def spectral_density(self): - return self.spectral_density_function(self.wavenumber, self.C) - - - @staticmethod - def correlation_function(): - """ - Correlation function is not defined for Kolmogorov spectrum. - """ - raise AttributeError("Correlation function is not defined for Kolmogorov spectrum.") - - @staticmethod - def structure_function(r, C, p=2.0/3.0): - """ - Structure function :math:`D(r)`. - - :param r: :math:`r` - :param C: :math:`C` - - The structure function is given by - - .. math:: D(r) = C^2 r^p - - where :math:`p = 2/3`. - - See Salomons, equation I.34. - """ - return C**2.0 * r**p - - -@six.add_metaclass(abc.ABCMeta) -class VonKarmanTemp(Spectrum): - """ - Abstract class for Von Karman spectrum when only temperature fluctuations are considered. - """ - - _required_attributes = ['a', 'mu_0'] - - - a = None - """ - Characteristic length :math:`a`. - """ - - mu_0 = None - """ - The standard deviation of the refractive-index fluctuation :math:`\\mu` is :math:`\\mu_0`. - """ - - def spectral_density(self): - return self.spectral_density_function(self.wavenumber, self.a, self.mu_0) - - - @staticmethod - def correlation_function(r, a, mu_0): - """ - Correlation function :math:`B(r)`. - - :param r: :math:`r` - :param a: :math:`a` - :param mu_0: :math:`\\mu_0` - - The correlation function is given by - - .. math:: B(r) = \\mu_0^2 \\frac{2^{2/3}}{\\Gamma(1/3)} \\left(\\frac{r}{a}\\right)^{1/3} K_{1/3} \\left(\\frac{r}{a}\\right) - - See Salomons, equation I.39. - """ - return mu_0**2.0 * 2.0**(2.0/3.0) / gamma(1.0/3.0) * (r/a)**(1.0/3.0) * bessel(1.0/3.0, r/a) - - @staticmethod - def structure_function(r, a, mu_0, smaller_than_factor=0.1): - """ - Structure function :math:`D(r)`. - - :param r: :math:`r` - :param a: :math:`a` - :param mu_0: :math:`\\mu_0` - :param smaller_than_factor: Factor - - - .. math:: D(r) = 2 \\mu_0^2 \\left[ 1 - \\frac{2^{2/3}}{\\Gamma(1/3)} \\left(\\frac{r}{a}\\right)^{1/3} K_{1/3} \\left(\\frac{r}{a}\\right) \\right] - - When :math:`a \\ll r`, or 'r < smaller_than_factor * a' - - .. math:: D(r) \\approx \\mu_0^2 \\frac{\\sqrt{\\pi}}{\\Gamma(7/6)} \\left( \\frac{r}{a} \\right)^{2/3} - - See Salomons, equation I.40. - """ - return (r < smaller_than_factor * a) * \ - ( mu_0**2.0 * np.sqrt(np.pi)/gamma(7.0/6.0) * (r/a)**(2.0/3.0) ) + \ - (r >= smaller_than_factor * a) * \ - ( mu_0**2.0 * (1.0 - 2.0**(2.0/3.0) / gamma(1.0/3.0) * (r/a)**(1.0/3.0) * bessel(1.0/3.0, r/a) ) ) - - #if r < smaller_than_factor * a: - #return mu_0**2.0 * np.sqrt(np.pi)/gamma(7.0/6.0) * (r/a)**(2.0/3.0) - #else: - #return mu_0**2.0 * (1.0 - 2.0**(2.0/3.0) / gamma(1.0/3.0) * (r/a)**(1.0/3.0) * bessel(1.0/3.0, r/a) ) - - - - -@six.add_metaclass(abc.ABCMeta) -class GaussianTempWind(object): - """ - Abstract class for Gaussian spectrum when both temperature and wind fluctuations are considered. - """ - - _required_attributes = ['a', 'sigma_T', 'T_0', 'sigma_nu', 'c_0'] - - - a = None - """ - Characteristic length :math:`a`. - """ - - - @staticmethod - def r(x, y, z): - """ - Distance :math:`r`. - - :param x: x - :param y: y - :param z: z - - .. math:: r = \\sqrt{x^2 + y^2 + z^2} - - """ - return (x**2.0 + y**2.0 + z**2.0) - - - @staticmethod - def rho(y, z): - """ - Distance :math:`\\rho`. - - :param y: y - :param z: z - - .. math:: \\rho = \\sqrt{y^2 + z^2} - - """ - return (z**2.0 + y**2.0)**0.5 - - - def spectral_density(self): - return self.spectral_density_function(self.wavenumber, self.theta, tuple(self.plane), self.a, self.sigma_T, self.T_0, self.sigma_nu, self.c_0) - - @staticmethod - def correlation_function(r, a, sigma_T, T_0, sigma_nu, c_0, rho): - """ - Correlation function :math:`B(r)`. - - - .. math:: B(x,y,z) = \\left[ \\frac{\\sigma_T^2}{4 T_0)^2} + \\frac{\\sigma_{\\nu}^2}{c_0^2} \\left( 1 - \\frac{\\rho^2}{a^2} \\right) \\right] \\exp{\\left( -r^2/a^2 \\right)} - - See Salomons, equation I.48. - """ - return (sigma_T/(2.0*T_0))**2.0 + (sigma_nu/c_0)**2.0 * (1.0 - (rho/a)**2.0) * np.exp(-(r/a)**2.0) - - - - -@six.add_metaclass(abc.ABCMeta) -class KolmogorovTempWind(object): - """ - Abstract class for Kolmogorov spectrum when both temperature and wind fluctuations are considered. - """ - pass - -@six.add_metaclass(abc.ABCMeta) -class VonKarmanTempWind(object): - """ - Abstract class for Von Karman spectrum when both temperature and wind fluctuations are considered. - """ - - _required_attributes = ['c_0', 'T_0', 'C_v', 'C_T', 'L'] - - CONSTANT_A = 5.0 / (18.0 * np.pi * gamma(1.0/3.0)) - """ - Constant `A`. - - .. math:: A = 5 / [ 18 \\pi \\Gamma(1/3) ] \\approx 0.0330 - - """ - - def spectral_density(self): - return self.spectral_density_function(self.wavenumber, self.theta, tuple(self.plane), self.c_0, self.T_0, self.C_v, self.C_T, self.L, self.CONSTANT_A) - - - diff --git a/acoustics/turbulence.py b/acoustics/turbulence.py deleted file mode 100644 index ae5f9179..00000000 --- a/acoustics/turbulence.py +++ /dev/null @@ -1,553 +0,0 @@ -""" -Turbulence -========== - -Turbulence in the atmosphere affects wave propagation. -This module contains implementations of models that can be used to create random turbulence fields. - -References are made to the book 'Computational Atmospheric Acoustics' by 'Erik M. Salomons', published in 2001. - -================ -Abstract classes -================ -.. inheritance-diagram:: acoustics._turbulence -.. automodule:: acoustics._turbulence - -========== -Turbulence -========== - -.. inheritance-diagram:: acoustics.turbulence -""" - - -import matplotlib.pyplot as plt -from scipy.special import gamma -from scipy.special import iv as bessel # Modified Bessel function of the first kind. - -from ._turbulence import * - -import numexpr as ne - -#try: - #import numba -#except ImportError: - #pass - -class Gaussian1DTemp(GaussianTemp, Spectrum1D): - """ - One-dimensional Gaussian turbulence spectrum supporting temperature fluctuations. - """ - - @staticmethod - def spectral_density_function(k, a, mu_0): - """ - One-dimensional spectral density :math:`V(k)`. - - :param k: :math:`k` - :param a: :math:`a` - :param mu_0: :math:`\\mu_0` - - The spectral density function is given by - - .. math:: V(k) = \\mu_0^2 \\frac{a}{2\\sqrt{\\pi}} \\exp{(-k^2a^2/4)} - - See Salomons, equation I.30. - """ - return mu_0**2.0 * a / (2.0 * np.sqrt(np.pi)) * np.exp(-k**2.0 * a**2.0 / 4) - - -class Kolmogorov1DTemp(KolmogorovTemp, Spectrum1D): - """ - One-dimensional Kolmogorov turbulence spectrum supporting temperature fluctuations. - """ - - @staticmethod - def spectral_density_function(k, C, p=2.0/3.0): - """ - One-dimensional spectral density density :math:`V(k)`. - - :param k: :math:`k` - :param C: :math:`C` - :param p: :math:`p = 2/3` - - The spectral density function is given by - - .. math:: V(k) = C^2 \\frac{\\Gamma(p+1)}{2\\pi} \\sin{\\left(\\frac{1}{2}\\pi p\\right)} \\left| k \\right|^{-p-1} - - See Salomons, equation I.35. - """ - return C**2.0 * gamma(p+1.0)/(2.0*np.pi) * np.sin(0.5*np.pi*p)*np.abs(k)**(-p-1.0) - - -class VonKarman1DTemp(VonKarmanTemp, Spectrum1D): - """ - One-dimensional Von Karman turbulence spectrum supporting temperature fluctuations. - """ - - @staticmethod - def spectral_density_function(k, a, mu_0): - """ - One-dimensional spectral density function :math:`V(k)`. - - :param r: :math:`r` - :param a: :math:`a` - :param mu_0: :math:`\\mu_0` - - The spectral density function is given by - - .. math:: V(k) = \\mu_0 \\frac{\\Gamma(5/6)}{\\Gamma(1/3) \\pi^{1/2}} \\frac{a}{\\left( 1 + k^2 a^2 \\right)^{5/6}} - - See Salomons, equation I.41. - """ - return mu_0 * gamma(5.0/6.0) / (gamma(1.0/3.0)*np.sqrt(np.pi)) * a / (1.0 + k**2.0 * a**2.0)**(5.0/6.0) - - -class Gaussian2DTemp(GaussianTemp, Spectrum2D): - """ - Two-dimensional Gaussian turbulence spectrum supporting temperature fluctuations. - """ - - @staticmethod - def spectral_density_function(k, a, mu_0): - """ - Two-dimensional spectral density :math:`F(k)`. - - .. math:: F(k) = \\mu_0^2 \\frac{ a^2 }{4 \\pi} \\exp{(-k^2 a^2 / 4)} - - See Salomons, equation I.31. - """ - return mu_0**2.0 * a**2.0 / (4.0 *np.pi) * np.exp(-k**2.0 * a**2.0 / 4) - - -class Kolmogorov2DTemp(KolmogorovTemp, Spectrum2D): - """ - Two-dimensional Kolmogorov turbulence spectrum supporting temperature fluctuations. - """ - - @staticmethod - def spectral_density_function(k, C, p=2.0/3.0): - """ - Two-dimensional spectraldensity density :math:`F(k)`. - - .. math:: F(k) = C^2 \\frac{\\Gamma^2(0.5 p + 1) 2^p}{2 \\pi^2} \\sin{\\left(\\frac{1}{2}\\pi p\\right)} \\left| k \\right|^{-p-2} - - See Salomons, equation I.36. - """ - return C**2.0 * gamma(0.5*p+1.0)*2.0**p/(2.0*np.pi**2.0) * np.sin(0.5*np.pi*p)*np.abs(k)**(-p-2.0) - - -class VonKarman2DTemp(VonKarmanTemp, Spectrum2D): - """ - Two-dimensional Von Karman turbulence spectrum supporting temperature fluctuations. - """ - - @staticmethod - def spectral_density_function(k, a, mu_0): - """ - Two-dimensional spectral density function :math:`F(k)`. - - .. math:: F(k) = \\mu_0^2 \\frac{\\Gamma(8/6)}{\\Gamma(1/3) \\pi} \\frac{a^2}{\\left( 1 + k^2 a^2 \\right)^{8/6}} - - See Salomons, equation I.42. - """ - return mu_0**2.0 * gamma(8.0/6.0) / (gamma(1.0/3.0)*np.pi) * a**2 / (1.0 + k**2.0 * a**2.0)**(8.0/6.0) - - -class Gaussian3DTemp(GaussianTemp, Spectrum3D): - """ - Three-dimensional Gaussian turbulence spectrum supporting temperature fluctuations. - """ - - @staticmethod - def spectral_density_function(k, a, mu_0): - """ - Three-dimensional spectral density :math:`\\Phi(k)`. - - .. math:: \\Phi(k) = \\mu_0^2 \\frac{a^3}{8\\pi^{3/2}} \\exp{(-k^2 a^2 /4)} - - See Salomons, equation I.32. - """ - return mu_0**2.0 * a**3.0 * np.exp(-k**2.0*a**2.0 / 4) - - -class Kolmogorov3DTemp(KolmogorovTemp, Spectrum3D): - """ - Three-dimensional Kolmogorov turbulence spectrum supporting temperature fluctuations. - """ - - @staticmethod - def spectral_density_function(k, C, p=2.0/3.0): - """ - Three-dimensional spectral density density :math:`\\Phi(k)`. - - .. math:: \\Phi(k) = C^2 \\frac{\\Gamma(p+2)}{4\\pi^2} \\sin{\\left(\\frac{1}{2}\\pi p\\right)} \\left| k \\right|^{-p-3} - - See Salomons, equation I.37. - """ - return C**2.0 * gamma(p+2.0)/(4.0*np.pi**2.0) * np.sin(0.5*np.pi*p)*np.abs(k)**(-p-3.0) - -class VonKarman3DTemp(VonKarmanTemp, Spectrum3D): - """ - Three-dimensional Von Karman turbulence spectrum supporting temperature fluctuations. - """ - - @staticmethod - def spectral_density_function(k, a, mu_0): - """ - Three-dimensional spectral density function :math:`\\Phi(k)`. - - .. math:: \\Phi(k) = \\mu_0 \\frac{\\Gamma(11/6)}{\\Gamma(1/3) \\pi^{3/2}} \\frac{a^3}{\\left( 1 + k^2 a^2 \\right)^{11/6}} - - See Salomons, equation I.43. - """ - return mu_0 * gamma(11.0/6.0) / (gamma(1.0/3.0)*np.pi**(1.5)) * a**3 / (1.0 + k**2.0 * a**2.0)**(11.0/6.0) - - - -class Gaussian2DTempWind(GaussianTempWind, Spectrum2D): - """ - Two-dimensional Gaussian turbulence spectrum supporting temperature and wind fluctuations. - """ - - - - - @staticmethod - def spectral_density_function(k, theta, plane, a, sigma_T, T_0, sigma_mu, c_0): - """ - Two-dimensional spectral density function :math:`F(k)`. - - - :param k: :math:`k` - :param plane: Tuple indicating which planes to consider. - - - The spectral density is calculated according to - - .. math:: F(k_x, k_y) = F(k_x, k_z) = \\frac{a^2}{4 \\pi} \\left( \\frac{\\sigma_T^2}{4 T_0^2} + \\frac{\\sigma_{\\mu}^2 [k_z^2 a^2 + 2] }{4 c_0^2} \\right) \\exp{(-k^2 a^2 / 4)} - - or - - .. math:: F(k_y, k_z) = \\frac{a^2}{4 \\pi} \\left( \\frac{\\sigma_T^2}{4 T_0^2} + \\frac{\\sigma_{\\mu}^2 [k^2 a^2 + 2] }{4 c_0^2} \\right) \\exp{(-k^2 a^2 / 4)} - - depending on the chosen plane. - - - See Salomons, page 215, and equation I.49 and I.50. - """ - - if plane == (1,0,1): # xz-plane - k_x = k * np.cos(theta) - k_z = k * np.sin(theta) - k = (k_x**2.0 + k_z**2.0)**(0.5) - return a**2.0/(4.0*np.pi) * ( (sigma_T/(2.0*T_0))**2.0 + sigma_mu**2.0/(4.0*c_0**2.0)*(k_z**2.0*a**2.0+1)) * np.exp(-k**2.0*a**2.0/4) - - elif plane == (1,1,0): # xy-plane - k_x = k * np.cos(theta) - k_y = k * np.sin(theta) - k = (k_x**2.0 + k_y**2.0)**(0.5) - return a**2.0/(4.0*np.pi) * ( (sigma_T/(2.0*T_0))**2.0 + sigma_mu**2.0/(4.0*c_0**2.0)*(k_y**2.0*a**2.0+1)) * np.exp(-k**2.0*a**2.0/4) - - elif plane == (0,1,1): # yz-plane - k_y = k * np.cos(theta) - k_z = k * np.sin(theta) - k = (k_y**2.0 + k_z**2.0)**(0.5) - return a**2.0/(4.0*np.pi) * ( (sigma_T/(2.0*T_0))**2.0 + sigma_mu**2.0/(4.0*c_0**2.0)*(k**2.0*a**2.0+1)) * np.exp(-k**2.0*a**2.0/4) - - else: - raise ValueError("Incorrect wavenumbers given.") - - - -class Kolmogorov2DTempWind(KolmogorovTempWind, Spectrum2D): - """ - Two-dimensional Kolmogorov turbulence spectrum support temperature and wind fluctuations. - """ - pass - -class VonKarman2DTempWind(VonKarmanTempWind, Spectrum2D): - """ - Two-dimensional Von Karman turbulence spectrum supporting temperature and wind fluctuations. - """ - - @staticmethod - def spectral_density_function(k, theta, plane, c_0, T_0, C_v, C_T, L, A): - """ - Two-dimensional spectral density function :math:`F(k)`. - - :param k: Wavenumber :math:`k` - :param c_0: :math:`c_0` - :param T_0: :math:`T_0` - :param C_v: :math:`C_v` - :param C_T: :math:`C_T` - :param L: :math:`L` - - See Salomons, equation I.53. - """ - K_0 = 2.0 * np.pi / L - - if plane == (1,0,1): # xz-plane - k_var = k*np.sin(theta) - - elif plane == (1,1,0): # xy-plane - k_var = k*np.sin(theta) - - elif plane == (0,1,1): # yz-plane - k_var = k - - f1 = A / (k**2.0 + K_0**2.0 )**(8.0/6.0) - f2 = gamma(1.0/2.0)*gamma(8.0/6.0) / gamma(11.0/6.0) * C_T**2.0/(4.0*T_0**2.0) - f3 = gamma(3.0/2.0)*gamma(8.0/6.0)/gamma(17.0/6.0) + k_var**2.0/(k**2.0+K_0**2.0) * gamma(1.0/2.0)*gamma(14.0/6.0)/gamma(17.0/6.0) - f4 = 22.0*C_v**2.0/(12.0*c_0**2.0) - - return f1 * (f2 + f3 * f4) - -class Gaussian3DTempWind(GaussianTempWind, Spectrum3D): - """ - Three-dimensional Von Karman turbulence spectrum supporting temperature and wind fluctiations. - """ - - #@staticmethod - #def spectral_density_function(): - #""" - #Three-dimensional spectral density function :math:`\\Phi(k_x, k_y, k_z)`. - - #See Salomons, I.51. - #""" - #raise NotImplementedError - -class Comparison(object): - """ - Compare turbulence spectra. - """ - - def __init__(self, items): - - self.items = items - """ - Turbulence spectra. - """ - - def plot_mode_amplitudes(self, filename=None): - """ - Create a plot of the mode amplitudes for all turbulence spectra. - """ - - fig = plt.figure() - - ax = fig.add_subplot(111) - - for item in self.items: - ax.plot(item.wavenumber, item.mode_amplitude(), label=item.__class__.__name__) - - ax.set_xlabel(r'$k$ in $\mathrm{m}^{-1}$') - ax.set_ylabel(r'$G$') - ax.grid() - - if filename: - fig.savefig(filename) - else: - fig.show() - - - - def plot_spectral_density(self, filename=None): - """ - Plot the spectral density. - """ - - fig = plt.figure() - - ax = fig.add_subplot(111) - - for item in self.items: - ax.loglog(item.wavenumber, item.spectral_density(), label=item.__class__.__name__) - - ax.set_xlabel(r'$k$ in $\mathrm{m}^{-1}$') - ax.set_ylabel(r'$F$') - ax.grid() - ax.legend() - - if filename: - fig.savefig(filename) - else: - fig.show() - - -def _mu(G, r_mesh, k_nr, z_mesh, k_nz, alpha_n): - return ne.evaluate("G * cos(r_mesh * k_nr + z_mesh * k_nz + alpha_n)") - -def _generate(r, z, delta_k, mode_amplitudes, modes, theta, alpha): - - mu = np.zeros((len(r), len(z)), dtype='float64') - - r_mesh, z_mesh = np.meshgrid(r, z) - r_mesh = r_mesh.T - z_mesh = z_mesh.T - #r_mesh, z_mesh = np.meshgrid(r, z, indexing='ij') - - for n, G, theta_n, alpha_n in zip(modes, mode_amplitudes, theta, alpha): - #for i in range(len(modes)): #zip(modes, mode_amplitudes, theta, alpha): - #n = modes[i] - #G = mode_amplitudes[i] - #theta_n = theta[i] - #alpha_n = alpha[i] - - k_n = n * delta_k - - k_nr = k_n * np.cos(theta_n) # Wavenumber component - k_nz = k_n * np.sin(theta_n) # Wavenumber component - - #mu_n = G * np.cos(r_mesh * k_nr + z_mesh * k_nz + alpha_n) - - mu_n = _mu(G, r_mesh, k_nr, z_mesh, k_nz, alpha_n) - mu += mu_n - - return mu - - -class Field2D(object): - """ - Refractive index field. - """ - - mu = None - """ - Refractive index. - """ - - def __init__(self, x, y, z, spatial_resolution, spectrum): - - self.x = x - """ - Size of field in x-direction. - """ - self.y = y - """ - Size of field in y-direction. - """ - self.z = z - """ - Size of field in z-direction. - """ - self.spatial_resolution = spatial_resolution - """ - Spatial resolution. - """ - self.spectrum = spectrum - """ - Spectrum. - """ - - #try: - #self._generate = numba.autojit(_generate) - #except NameError: - self._generate = _generate - - - def randomize(self): - """ - Create new random values. This is a shortcut to :meth:`Spectrum2D.randomize()`. - """ - self.spectrum.randomize() - return self - - def generate(self): - r = self.x - z = self.z - r = np.arange(np.ceil(r/self.spatial_resolution)) * self.spatial_resolution - z = np.arange(np.ceil(z/self.spatial_resolution)) * self.spatial_resolution - delta_k = self.spectrum.wavenumber_resolution - - self.mu = self._generate(r, z, delta_k, self.spectrum.mode_amplitude(), self.spectrum.modes, self.spectrum.theta, self.spectrum.alpha) - - return self - - - #def generate(self): - #""" - #Create a random realization of the refractive-index fluctuations. To actually create a random field, call :meth:`randomize` first. - - #.. math:: \\mu(r) = \\sqrt{4\\pi \\Delta k} \\sum_n \\cos{\\left( \\mathbf{k}_n \cdot \\mathbf{r} + \\alpha_n \\right)} \\sqrt{F(\\mathbf{k_n} k_n} - - #""" - - #r = self.x - #z = self.z - - #r = np.arange(np.ceil(r/self.spatial_resolution)) * self.spatial_resolution - #z = np.arange(np.ceil(z/self.spatial_resolution)) * self.spatial_resolution - - ##r = np.arange(0.0, r, self.spatial_resolution) - ##z = np.arange(0.0, z, self.spatial_resolution) - - #delta_k = self.spectrum.wavenumber_resolution - - #mu = list() - - #mode_amplitudes = self.spectrum.mode_amplitude() - - #for n, G, theta_n, alpha_n in zip(self.spectrum.modes, mode_amplitudes, self.spectrum.theta, self.spectrum.alpha): - - #k_n = n * delta_k - - #k_nr = k_n * np.cos(theta_n) # Wavenumber component - #k_nz = k_n * np.sin(theta_n) # Wavenumber component - - ##r_mesh, z_mesh = np.meshgrid(r, z, indexing='ij') - #r_mesh, z_mesh = np.meshgrid(r, z) - #r_mesh = r_mesh.T - #z_mesh = z_mesh.T - - ##k_n_v = np.vstack( , k_n * np.sin(theta)) - - #mu_n = G * np.cos(r_mesh * k_nr + z_mesh * k_nz + alpha_n) - #mu.append(mu_n) - - #self.mu = sum(mu) - #return self - - #def plot_correlation(self): - #""" - #Plot the correlation. - #""" - - #fig = plt.figure(figsize=(16,12) - - - def plot(self, filename=None): - """ - Plot the field. - """ - - if self.mu is None: - raise ValueError("Need to calculate the refractive index first.") - - r = self.x - z = self.z - r = np.arange(np.ceil(r/self.spatial_resolution)) * self.spatial_resolution - z = np.arange(np.ceil(z/self.spatial_resolution)) * self.spatial_resolution - - fig = plt.figure(figsize=(16, 12), dpi=80) - ax = fig.add_subplot(111, aspect='equal') - ax.set_title("Refractive-index field") - plot = ax.pcolormesh(r, z, self.mu.T) - ax.set_xlabel(r'$r$ in m') - ax.set_ylabel(r'$z$ in m') - ax.set_xlim(r[0], r[-1]) - ax.set_ylim(z[0], z[-1]) - - orientation = 'horizontal' if self.x > self.z else 'vertical' - - c = fig.colorbar(plot, orientation=orientation, pad=0.06) - c.set_label(r'Refractive-index fluctuation $\mu$') - plt.tight_layout() - if filename: - fig.savefig(filename, bbox_inches='tight') - else: - fig.show() - - - - diff --git a/examples/salomons_figure_I3.py b/examples/salomons_figure_I3.py deleted file mode 100644 index d0771f26..00000000 --- a/examples/salomons_figure_I3.py +++ /dev/null @@ -1,36 +0,0 @@ -""" -This script reproduces figure I.3 from Salomons. -""" - -from acoustics.turbulence import Gaussian2DTempWind, VonKarman2DTempWind - - -def main(): - - - a = 1.0 # Correlation length - - - sigma_T = np.sqrt(1e-5) / 2.0 - T_0 = 1.0 - mu_0 = np.sqrt(1e-5) - sigma_nu = 0.0 - K_0 = 10.0 - L = 2.0 * np.pi / K_0 - C_T = np.sqrt(6e-7) - c_0 = 1.0 - C_v = np.sqrt(2e-6) - l_0 = 0.001 - k_max = 5.48 / l_0 - k_x = 0.0 - - g = Gaussian3DTempWind() - - - - - - - -if __name__ == '__main__': - main() \ No newline at end of file diff --git a/examples/salomons_figure_J1.py b/examples/salomons_figure_J1.py deleted file mode 100644 index fcbbfa19..00000000 --- a/examples/salomons_figure_J1.py +++ /dev/null @@ -1,75 +0,0 @@ -""" -This script reproduces figure J.1 from Salomons. -""" -import numpy as np -from acoustics.turbulence import Gaussian2DTempWind, VonKarman2DTempWind, Comparison, Field2D - -def main(): - - - """Given parameters.""" - N = 200 # Amount of modes - k_max = 10.0 # Maxmimum wavenumber - - """Parameters for Gaussian spectrum.""" - a = 1.0 # Correlation length - sigma_T = np.sqrt(1e-5) # Standard deviation in temperature - T_0 = 1.0 # Temperature - sigma_nu = 0.0 # Standard deviation in wind - - - """Parameters for Von Karman spectrum.""" - K_0 = 1.0/10.0 - L = 2.0 * np.pi / K_0 # Size of largest eddy. - C_T = np.sqrt(1e-7) # - T_0 = 1.0 # - C_v = 0.001 - c_0 = 1.0 - - - """Other settings.""" - wavenumber_resolution = k_max / N - spatial_resolution = 0.05 # We don't need it for the calculations but we do need it to create an instance. - - x = 20.0 - y = 0.0 - z = 40.0 - - - g = Gaussian2DTempWind(max_mode_order=N, - a=a, - sigma_T=sigma_T, - T_0=T_0, - sigma_nu=sigma_nu, - c_0=c_0, - wavenumber_resolution=wavenumber_resolution, - plane=(1,0,1) - ) - - vk = VonKarman2DTempWind(max_mode_order=N, - L=L, - C_T=C_T, - T_0=T_0, - C_v=C_v, - c_0=c_0, - wavenumber_resolution=wavenumber_resolution, - plane=(1,0,1) - ) - - - c = Comparison([g, vk]) - - """The following figure is a reproduction of Salomons J.1 figure.""" - c.plot_mode_amplitudes('salomons_figure_J1.png') - - """We can additionally calculate turbulent fields according to these two spectra.""" - field_g = Field2D(x=x, y=y, z=z, spatial_resolution=spatial_resolution, spectrum=g) - field_g.generate().plot('Gaussian2DTempWind_field.png') - - field_vk = Field2D(x=x, y=y, z=z, spatial_resolution=spatial_resolution, spectrum=vk) - field_vk.generate().plot('VonKarman2DTempWind_field.png') - - c.plot_spectral_density('Gaussian2DTempWind_and_VonKarman2DTempWind_spectral_density.png') - -if __name__ == '__main__': - main() \ No newline at end of file diff --git a/examples/turbulence.py b/examples/turbulence.py deleted file mode 100644 index 7ad5b06a..00000000 --- a/examples/turbulence.py +++ /dev/null @@ -1,59 +0,0 @@ -import numpy as np -from acoustics.turbulence import Gaussian2DTemp, VonKarman2DTemp, Comparison, Field2D - -def main(): - - mu_0 = np.sqrt(10.0**(-6)) - correlation_length = 1.0 # Typical correlation length for Gaussian spectrum. - - x = 20.0 - y = 0.0 - z = 40.0 - - plane = (1,0,1) - - #f_resolution = wavenumber_resolution / (2.0*np.pi) - - spatial_resolution = 0.05 - - - N = 100 - - min_wavenumber = 0.01 - max_wavenumber = 10.0 - - - wavenumber_resolution = (max_wavenumber - min_wavenumber) / N - - - """Create an object to describe an Gaussian turbulence spectrum.""" - g = Gaussian2DTemp(plane=plane, a=correlation_length, mu_0=mu_0, wavenumber_resolution=wavenumber_resolution, max_mode_order=N) - - """Create an object to describe a VonKarman turbulence spectrum.""" - s = VonKarman2DTemp(plane=plane, a=correlation_length, mu_0=mu_0, wavenumber_resolution=wavenumber_resolution, max_mode_order=N) - - g.plot_mode_amplitudes('Gaussian2DTemp_mode_amplitudes.png') - s.plot_mode_amplitudes('VonKarman2DTemp_mode_amplitudes.png') - - c = Comparison([g, s]) - - c.plot_mode_amplitudes('Gaussian2DTemp_and_VonKarman2DTemp_mode_amplitudes.png') - - - field_g = Field2D(x=x, y=y, z=z, spatial_resolution=spatial_resolution, spectrum=g) - field_s = Field2D(x=x, y=y, z=z, spatial_resolution=spatial_resolution, spectrum=s) - - field_g.generate().plot('Gaussian2DTemp_field.png') - field_s.generate().plot('VonKarman2DTemp_field.png') - - - - - - - - - - -if __name__ == '__main__': - main() \ No newline at end of file diff --git a/tests/test_turbulence.py b/tests/test_turbulence.py deleted file mode 100644 index 43301065..00000000 --- a/tests/test_turbulence.py +++ /dev/null @@ -1,77 +0,0 @@ -""" -Just testing whether an example executes without generating exceptions... -""" -import unittest -import numpy as np - -from acoustics.turbulence import Gaussian2DTempWind, VonKarman2DTempWind, Comparison, Field2D - - -class TurbulenceCase(unittest.TestCase): - - def test_some(self): - """Given parameters.""" - N = 20 # Amount of modes - k_max = 10.0 # Maxmimum wavenumber - - """Parameters for Gaussian spectrum.""" - a = 1.0 # Correlation length - sigma_T = np.sqrt(1e-5) # Standard deviation in temperature - T_0 = 1.0 # Temperature - sigma_nu = 0.0 # Standard deviation in wind - - - """Parameters for Von Karman spectrum.""" - K_0 = 1.0/10.0 - L = 2.0 * np.pi / K_0 # Size of largest eddy. - C_T = np.sqrt(1e-7) # - T_0 = 1.0 # - C_v = 0.001 - c_0 = 1.0 - - - """Other settings.""" - wavenumber_resolution = k_max / N - spatial_resolution = 0.05 # We don't need it for the calculations but we do need it to create an instance. - - x = 20.0 - y = 0.0 - z = 40.0 - - plane=(1,0,1) - - - g = Gaussian2DTempWind(max_mode_order=N, - a=a, - sigma_T=sigma_T, - T_0=T_0, - sigma_nu=sigma_nu, - c_0=c_0, - wavenumber_resolution=wavenumber_resolution, - plane=plane - ) - - vk = VonKarman2DTempWind(max_mode_order=N, - L=L, - C_T=C_T, - T_0=T_0, - C_v=C_v, - c_0=c_0, - wavenumber_resolution=wavenumber_resolution, - plane=plane - ) - - field_g = Field2D(x=x, y=y, z=z, spectrum=g, spatial_resolution=spatial_resolution) - field_vk = Field2D(x=x, y=y, z=z, spectrum=vk, spatial_resolution=spatial_resolution) - c = Comparison([g, vk]) - - field_g.generate() - field_vk.generate() - - - -if __name__ == '__main__': - unittest.main() - - - \ No newline at end of file