Skip to content

Scattering models

This is the API reference for the scattering model part of the echoSMs package.

Each type of model is contained in a separate Python class (with name ending in Model), but with common calling signatures across all model classes, as defined in ScatterModelBase. There are also classes to provide ready access to the benchmark models and reference model definitions.

echosms.ScatterModelBase()

Bases: ABC

Base class for a class that provides a scattering model.

All scattering models should inherit from this class, have a name that ends with 'Model', and provide initialisation and calculate_ts_single() functions.

Initialise.

ATTRIBUTE DESCRIPTION
long_name

The long name of the model.

TYPE: str

short_name

A short version of the model's long name, typically an acronym.

TYPE: str

analytical_type

Whether the model implements an exact or an approximate model.

TYPE: str

boundary_types

The types of boundary conditions that the model provides, e.g., 'fixed rigid', 'pressure release', 'fluid filled'

TYPE: list[boundary_type]

shapes

The target shapes that the model can represent.

TYPE: list[str]

max_ka

An approximate maximum ka value that will result in accurate target strength results. Note that ka is often not the only parameter that determines the accuracy of the model (e.g., aspect ratio and incident angle can also affect the accuracy).

TYPE: float

no_expand_parameters

The model parameters that are not expanded into Pandas DataFrame columns or Xarray DataArray coordinates. They will instead end up as a dict in the DataFrame or DataArray attrs attribute.

TYPE: list[str]

Source code in src/echosms/scattermodelbase.py 
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
@abc.abstractmethod
def __init__(self):
    """Initialise.

    Attributes
    ----------
    long_name : str
        The long name of the model.
    short_name : str
        A short version of the model's long name, typically an acronym.
    analytical_type : str
        Whether the model implements an ``exact`` or an ``approximate`` model.
    boundary_types : list[boundary_type]
        The types of boundary conditions that the model provides, e.g., 'fixed rigid',
        'pressure release', 'fluid filled'
    shapes : list[str]
        The target shapes that the model can represent.
    max_ka : float
        An approximate maximum ka value that will result in accurate target strength results.
        Note that ka is often not the only parameter that determines the accuracy of the
        model (e.g., aspect ratio and incident angle can also affect the accuracy).
    no_expand_parameters : list[str]
        The model parameters that are not expanded into Pandas DataFrame columns or
        Xarray DataArray coordinates. They will instead end up as a dict in the DataFrame or
        DataArray `attrs` attribute.
    """
    self.long_name = ''
    self.short_name = ''
    self.analytical_type = ''
    self.boundary_types = []
    self.shapes = []
    self.max_ka = np.nan
    self.no_expand_parameters = []

calculate_ts(data, expand=False, inplace=False, multiprocess=False, progress=False)

Calculate the target strength (TS) for many parameters.

PARAMETER DESCRIPTION
data

Requirements for the different input data types are:

  • DataFrame: column names must match the function parameter names in calculate_ts_single(). One TS value will be calculated for each row in the DataFrame.
  • DataArray: dimension names must match the function parameter names in calculate_ts_single(). TS values will be calculated for all combinations of the coordinate variables.
  • dict: keys must match the function parameters in calculate_ts_single(). TS values will be calculated for all combinations of the dict values.

TYPE: dict | DataFrame | DataArray

multiprocess

Split the ts calculation across CPU cores. Multiprocessing is currently provided by mapply. For more sophisticated uses it may be preferred to use a multiprocessing package of your choice directly on the calculate_ts_single() method. See the code in this method (calculate_ts()) for an example.

TYPE: bool DEFAULT: False

expand

Only applicable if data is a dict. If True, will use as_dataframe() to expand the dict into a DataFrame with one column per dict key and return that, adding a column named ts for the results.

TYPE: bool DEFAULT: False

inplace

Only applicable if data is a DataFrame. If True, the results will be added to the input DataFrame in a column named ts. If a ts column already exists, it is overwritten.

TYPE: bool DEFAULT: False

progress

If True, will produce a progress bar while running models

TYPE: bool DEFAULT: False

RETURNS DESCRIPTION
None | list[float] | Series | DataFrame

The return type and value are determined by the type of the input variable (data) and the expand and inplace parameters:

  • dict input and expand=False returns a list of floats.
  • dict input and expand=True returns a DataFrame.
  • DataFrame input and inplace=False returns a Series.
  • DataFrame input and inplace=True modifies data and returns None.
  • DataArray input always modifies data and returns None.
Source code in src/echosms/scattermodelbase.py 
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
def calculate_ts(self, data: dict | pd.DataFrame | xr.DataArray, expand: bool=False,
                 inplace: bool=False, multiprocess: bool=False,
                 progress: bool=False) -> None | list[float] | pd.Series | pd.DataFrame:
    """Calculate the target strength (TS) for many parameters.

    Parameters
    ----------
    data :
        Requirements for the different input data types are:

        - **DataFrame**: column names must match the function parameter names in
          calculate_ts_single(). One TS value will be calculated for each row in the DataFrame.
        - **DataArray**: dimension names must match the function parameter names in
          calculate_ts_single(). TS values will be calculated for all combinations of the
          coordinate variables.
        - **dict**: keys must match the function parameters in calculate_ts_single().
          TS values will be calculated for all combinations of the dict values.

    multiprocess :
        Split the ts calculation across CPU cores. Multiprocessing is currently provided by
        [mapply](https://github.com/ddelange/mapply). For more
        sophisticated uses it may be preferred to use a multiprocessing package of your choice
        directly on the `calculate_ts_single()` method. See the code in this method
        (`calculate_ts()`) for an example.

    expand :
        Only applicable if `data` is a dict. If `True`, will use
        [`as_dataframe()`][echosms.utils.as_dataframe]
        to expand the dict into a DataFrame with one column per dict key
        and return that, adding a column named `ts` for the results.

    inplace :
        Only applicable if `data` is a DataFrame. If `True`, the results
        will be added to the input DataFrame in a column named `ts`. If a `ts` column
        already exists, it is overwritten.

    progress :
        If `True`, will produce a progress bar while running models

    Returns
    -------
    :
        The return type and value are determined by the type of the input variable (`data`) and
        the `expand` and `inplace` parameters:

        - dict input and `expand=False` returns a list of floats.
        - dict input and `expand=True` returns a DataFrame.
        - DataFrame input and `inplace=False` returns a Series.
        - DataFrame input and `inplace=True` modifies `data` and returns `None`.
        - DataArray input always modifies `data` and returns `None`.

    """
    match data:
        case dict():
            data_df = as_dataframe(data, self.no_expand_parameters)
        case pd.DataFrame():
            data_df = data
        case xr.DataArray():
            data_df = data.to_dataframe().reset_index()
            data_df.attrs = data.attrs
        case _:
            raise ValueError(f'Data type of {type(data)} is not supported'
                             ' (only dictionaries, Pandas DataFrames and'
                             ' Xarray DataArrays are).')

    self.validate_parameters(data_df)

    # Get the non-expandable model parameters
    p = data_df.attrs['parameters'] if 'parameters' in data_df.attrs else {}

    # Note: the args argument in the apply call below requires a tuple. data_df.attrs is a
    # dict and the default behaviour is to make a tuple using the dict keys. The trailing comma
    # and parenthesis instead causes the tuple to have one entry of the dict.

    if multiprocess:
        from mapply.mapply import mapply
        ts = mapply(data_df, self.__ts_helper, args=(p,), axis=1, progressbar=progress)
    else:  # this uses just one CPU
        if progress:
            tqdm.pandas(desc=self.short_name, unit=' models',
                        bar_format='{l_bar}{bar} [{n_fmt}/{total_fmt}; {rate_noinv_fmt}]')
            ts = data_df.progress_apply(self.__ts_helper, args=(p,), axis=1)
        else:
            ts = data_df.apply(self.__ts_helper, args=(p,), axis=1)

    match data:
        case dict() if expand:
            data_df['ts'] = ts
            return data_df
        case dict():
            return ts.to_list()
        case pd.DataFrame() if inplace:
            data_df['ts'] = ts
            return None
        case pd.DataFrame():
            return ts.rename('ts', inplace=True)
        case xr.DataArray():
            data.values = ts.to_numpy().reshape(data.shape)
            return None
        case _:
            raise AssertionError('This code should never be reached - unsupported input data '
                                 f'type of {type(data)}.')

calculate_ts_single() abstractmethod

Calculate the TS for one parameter set.

Source code in src/echosms/scattermodelbase.py 
250
251
252
@abc.abstractmethod
def calculate_ts_single(self):
    """Calculate the TS for one parameter set."""

validate_parameters(p) abstractmethod

Validate the model parameters.

PARAMETER DESCRIPTION
p

The model parameters.

TYPE: dict | DataFrame | DataArray

RAISES DESCRIPTION
ValueError

If any of the model parameters are invalid.
KeyError

If any required model parameters are not present.
Source code in src/echosms/scattermodelbase.py 
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
@abc.abstractmethod
def validate_parameters(self, p: dict | pd.DataFrame | xr.DataArray):
    """Validate the model parameters.

    Parameters
    ----------
    p :
        The model parameters.

    Raises
    ------
    ValueError
        If any of the model parameters are invalid.
    KeyError
        If any required model parameters are not present.
    """

echosms.DCMModel()

Bases: ScatterModelBase

Modal series deformed cylinder model (DCM).

This class contains methods to calculate acoustic scatter from finite straight cylinders with various boundary conditions.

Source code in src/echosms/dcmmodel.py 
17
18
19
20
21
22
23
24
def __init__(self):
    super().__init__()
    self.long_name = 'deformed cylinder model'
    self.short_name = 'dcm'
    self.analytical_type = 'approximate analytical'
    self.boundary_types = [bt.fixed_rigid, bt.pressure_release, bt.fluid_filled]
    self.shapes = ['finite cylinder']
    self.max_ka = 20  # [1]

calculate_ts(data, expand=False, inplace=False, multiprocess=False, progress=False)

Calculate the target strength (TS) for many parameters.

PARAMETER DESCRIPTION
data

Requirements for the different input data types are:

  • DataFrame: column names must match the function parameter names in calculate_ts_single(). One TS value will be calculated for each row in the DataFrame.
  • DataArray: dimension names must match the function parameter names in calculate_ts_single(). TS values will be calculated for all combinations of the coordinate variables.
  • dict: keys must match the function parameters in calculate_ts_single(). TS values will be calculated for all combinations of the dict values.

TYPE: dict | DataFrame | DataArray

multiprocess

Split the ts calculation across CPU cores. Multiprocessing is currently provided by mapply. For more sophisticated uses it may be preferred to use a multiprocessing package of your choice directly on the calculate_ts_single() method. See the code in this method (calculate_ts()) for an example.

TYPE: bool DEFAULT: False

expand

Only applicable if data is a dict. If True, will use as_dataframe() to expand the dict into a DataFrame with one column per dict key and return that, adding a column named ts for the results.

TYPE: bool DEFAULT: False

inplace

Only applicable if data is a DataFrame. If True, the results will be added to the input DataFrame in a column named ts. If a ts column already exists, it is overwritten.

TYPE: bool DEFAULT: False

progress

If True, will produce a progress bar while running models

TYPE: bool DEFAULT: False

RETURNS DESCRIPTION
None | list[float] | Series | DataFrame

The return type and value are determined by the type of the input variable (data) and the expand and inplace parameters:

  • dict input and expand=False returns a list of floats.
  • dict input and expand=True returns a DataFrame.
  • DataFrame input and inplace=False returns a Series.
  • DataFrame input and inplace=True modifies data and returns None.
  • DataArray input always modifies data and returns None.
Source code in src/echosms/scattermodelbase.py 
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
def calculate_ts(self, data: dict | pd.DataFrame | xr.DataArray, expand: bool=False,
                 inplace: bool=False, multiprocess: bool=False,
                 progress: bool=False) -> None | list[float] | pd.Series | pd.DataFrame:
    """Calculate the target strength (TS) for many parameters.

    Parameters
    ----------
    data :
        Requirements for the different input data types are:

        - **DataFrame**: column names must match the function parameter names in
          calculate_ts_single(). One TS value will be calculated for each row in the DataFrame.
        - **DataArray**: dimension names must match the function parameter names in
          calculate_ts_single(). TS values will be calculated for all combinations of the
          coordinate variables.
        - **dict**: keys must match the function parameters in calculate_ts_single().
          TS values will be calculated for all combinations of the dict values.

    multiprocess :
        Split the ts calculation across CPU cores. Multiprocessing is currently provided by
        [mapply](https://github.com/ddelange/mapply). For more
        sophisticated uses it may be preferred to use a multiprocessing package of your choice
        directly on the `calculate_ts_single()` method. See the code in this method
        (`calculate_ts()`) for an example.

    expand :
        Only applicable if `data` is a dict. If `True`, will use
        [`as_dataframe()`][echosms.utils.as_dataframe]
        to expand the dict into a DataFrame with one column per dict key
        and return that, adding a column named `ts` for the results.

    inplace :
        Only applicable if `data` is a DataFrame. If `True`, the results
        will be added to the input DataFrame in a column named `ts`. If a `ts` column
        already exists, it is overwritten.

    progress :
        If `True`, will produce a progress bar while running models

    Returns
    -------
    :
        The return type and value are determined by the type of the input variable (`data`) and
        the `expand` and `inplace` parameters:

        - dict input and `expand=False` returns a list of floats.
        - dict input and `expand=True` returns a DataFrame.
        - DataFrame input and `inplace=False` returns a Series.
        - DataFrame input and `inplace=True` modifies `data` and returns `None`.
        - DataArray input always modifies `data` and returns `None`.

    """
    match data:
        case dict():
            data_df = as_dataframe(data, self.no_expand_parameters)
        case pd.DataFrame():
            data_df = data
        case xr.DataArray():
            data_df = data.to_dataframe().reset_index()
            data_df.attrs = data.attrs
        case _:
            raise ValueError(f'Data type of {type(data)} is not supported'
                             ' (only dictionaries, Pandas DataFrames and'
                             ' Xarray DataArrays are).')

    self.validate_parameters(data_df)

    # Get the non-expandable model parameters
    p = data_df.attrs['parameters'] if 'parameters' in data_df.attrs else {}

    # Note: the args argument in the apply call below requires a tuple. data_df.attrs is a
    # dict and the default behaviour is to make a tuple using the dict keys. The trailing comma
    # and parenthesis instead causes the tuple to have one entry of the dict.

    if multiprocess:
        from mapply.mapply import mapply
        ts = mapply(data_df, self.__ts_helper, args=(p,), axis=1, progressbar=progress)
    else:  # this uses just one CPU
        if progress:
            tqdm.pandas(desc=self.short_name, unit=' models',
                        bar_format='{l_bar}{bar} [{n_fmt}/{total_fmt}; {rate_noinv_fmt}]')
            ts = data_df.progress_apply(self.__ts_helper, args=(p,), axis=1)
        else:
            ts = data_df.apply(self.__ts_helper, args=(p,), axis=1)

    match data:
        case dict() if expand:
            data_df['ts'] = ts
            return data_df
        case dict():
            return ts.to_list()
        case pd.DataFrame() if inplace:
            data_df['ts'] = ts
            return None
        case pd.DataFrame():
            return ts.rename('ts', inplace=True)
        case xr.DataArray():
            data.values = ts.to_numpy().reshape(data.shape)
            return None
        case _:
            raise AssertionError('This code should never be reached - unsupported input data '
                                 f'type of {type(data)}.')

calculate_ts_single(medium_c, medium_rho, a, b, theta, f, boundary_type, target_c=None, target_rho=None, validate_parameters=True, **kwargs)

Calculate the scatter from a finite cylinder using the modal series deformed cylinder model.

PARAMETER DESCRIPTION
medium_c

Sound speed in the fluid medium surrounding the target [m/s].

TYPE: float

medium_rho

Density of the fluid medium surrounding the target [kg/m³].

TYPE: float

a

Radius of the cylinderical target [m].

TYPE: float

b

Length of the cylinderical target [m].

TYPE: float

theta

Pitch angle to calculate the scattering as per the echoSMs coordinate system [°].

TYPE: float

f

Frequency to calculate the scattering at [Hz].

TYPE: float

boundary_type

The model type. Supported model types are given in the boundary_types class attribute.

TYPE: boundary_type

target_c

Sound speed in the fluid inside the sphere [m/s]. Only required for boundary_type of fluid_filled.

TYPE: None | float DEFAULT: None

target_rho

Density of the fluid inside the sphere [kg/m³]. Only required for boundary_type of fluid_filled.

TYPE: None | float DEFAULT: None

validate_parameters

Whether to validate the model parameters.

TYPE: bool DEFAULT: True

RETURNS DESCRIPTION
float

The target strength (re 1 m²) of the target [dB].
Notes

The class implements the code in Section B.1 of Jech et al. (2015).

References

Jech, J.M., Horne, J.K., Chu, D., Demer, D.A., Francis, D.T.I., Gorska, N., Jones, B., Lavery, A.C., Stanton, T.K., Macaulay, G.J., Reeder, D.B., Sawada, K., 2015. Comparisons among ten models of acoustic backscattering used in aquatic ecosystem research. Journal of the Acoustical Society of America 138, 3742–3764. https://doi.org/10.1121/1.4937607

Source code in src/echosms/dcmmodel.py 
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
def calculate_ts_single(self, medium_c: float, medium_rho: float, a: float, b: float,
                        theta: float, f: float, boundary_type: bt,
                        target_c: None | float=None, target_rho: None | float=None,
                        validate_parameters: bool=True,
                        **kwargs) -> float:
    """
    Calculate the scatter from a finite cylinder using the modal series deformed cylinder model.

    Parameters
    ----------
    medium_c :
        Sound speed in the fluid medium surrounding the target [m/s].
    medium_rho :
        Density of the fluid medium surrounding the target [kg/m³].
    a :
        Radius of the cylinderical target [m].
    b :
        Length of the cylinderical target [m].
    theta :
        Pitch angle to calculate the scattering as per the echoSMs
        [coordinate system](conventions.md#coordinate-systems) [°].
    f :
        Frequency to calculate the scattering at [Hz].
    boundary_type :
        The model type. Supported model types are given in the `boundary_types` class attribute.
    target_c :
        Sound speed in the fluid inside the sphere [m/s].
        Only required for `boundary_type` of ``fluid_filled``.
    target_rho :
        Density of the fluid inside the sphere [kg/m³].
        Only required for `boundary_type` of ``fluid_filled``.
    validate_parameters :
        Whether to validate the model parameters.

    Returns
    -------
    :
        The target strength (re 1 m²) of the target [dB].

    Notes
    -----
    The class implements the code in Section B.1 of Jech et al. (2015).

    References
    ----------
    Jech, J.M., Horne, J.K., Chu, D., Demer, D.A., Francis, D.T.I., Gorska, N., Jones, B.,
    Lavery, A.C., Stanton, T.K., Macaulay, G.J., Reeder, D.B., Sawada, K., 2015.
    Comparisons among ten models of acoustic backscattering used in aquatic ecosystem
    research. Journal of the Acoustical Society of America 138, 3742–3764.
    <https://doi.org/10.1121/1.4937607>
    """
    if validate_parameters:
        self.validate_parameters(locals())

    if theta == 0.0:
        return nan

    theta_rad = theta*pi/180.
    kL = wavenumber(medium_c, f)*b
    K = wavenumber(medium_c, f) * sin(theta_rad)
    Ka = K*a

    m = range(30)  # TODO this needs to vary with f

    match boundary_type:
        case bt.fixed_rigid:
            series = map(lambda m: (-1)**m * Neumann(m)*(jvp(m, Ka) / h1vp(m, Ka)), m)
        case bt.pressure_release:
            series = map(lambda m: (-1)**m * Neumann(m)*(jv(m, Ka) / hankel1(m, Ka)), m)
        case bt.fluid_filled:
            g = target_rho/medium_rho
            h = target_c/medium_c
            gh = g*h
            Kda = K/h*a

            def Cm(m):
                numer = (jvp(m, Kda)*yv(m, Ka)) / (jv(m, Kda)*jvp(m, Ka))\
                    - gh*(yvp(m, Ka)/jvp(m, Ka))
                denom = (jvp(m, Kda)*jv(m, Ka)) / (jv(m, Kda)*jvp(m, Ka)) - gh
                return numer/denom

            series = map(lambda m: 1j**(2*m) * Neumann(m) / (1 + 1j*Cm(m)), m)
        case _:
            raise ValueError(f'The {self.long_name} model does not support '
                             f'a model type of "{boundary_type}".')

    fbs = 1j*b/pi * (sin(kL*cos(theta_rad)) / (kL*cos(theta_rad))) * sum(series)
    return 20*log10(abs(fbs))  # ts

validate_parameters(params)

Validate the model parameters.

See here for calling details.

Source code in src/echosms/dcmmodel.py 
26
27
28
29
30
31
32
33
34
35
36
37
38
39
def validate_parameters(self, params):
    """Validate the model parameters.

    See [here][echosms.ScatterModelBase.validate_parameters] for calling details.
    """
    p = as_dict(params)
    super()._present_and_in(p, ['boundary_type'], self.boundary_types)
    super()._present_and_positive(p, ['medium_rho', 'medium_c', 'a', 'b', 'f'])

    types = np.unique(np.atleast_1d(p['boundary_type']))
    for t in types:
        if t == bt.fluid_filled:
            super()._present_and_positive(p, ['target_c', 'target_rho'],
                                          mask=p['boundary_type'] == t)

DWBA models

There are several models that use the distorted-wave Born approximation, documented below. There are also some functions to make cylinder and spheroid shapes for use in the DWBA models.

echosms.DWBAModel()

Bases: ScatterModelBase

Distorted-wave Born approximation (DWBA) scattering models.

This class calculates acoustic scatter from piecewise cylindical shapes with weakly scattering material contrasts. It can also add a stochastic component to provide the stochastic DWBA (SDWBA) model.

Source code in src/echosms/dwbamodel.py 
22
23
24
25
26
27
28
29
30
31
32
33
34
def __init__(self):
    super().__init__()
    self.long_name = 'distorted-wave Born approximation'
    self.short_name = 'dwba'
    self.analytical_type = 'approximate'
    self.boundary_types = [bt.fluid_filled]
    self.shapes = ['piecewise cylindical']
    self.max_ka = 20
    # The distorted wave Born approximation is increasingly inaccurate outside these limits:
    self.g_range = [0.95, 1.05]
    self.h_range = [0.95, 1.05]
    self.no_expand_parameters = ['a', 'rv_pos', 'rv_tan']
    self.rng = np.random.default_rng()  # for SDWBA

calculate_ts(data, expand=False, inplace=False, multiprocess=False, progress=False)

Calculate the target strength (TS) for many parameters.

PARAMETER DESCRIPTION
data

Requirements for the different input data types are:

  • DataFrame: column names must match the function parameter names in calculate_ts_single(). One TS value will be calculated for each row in the DataFrame.
  • DataArray: dimension names must match the function parameter names in calculate_ts_single(). TS values will be calculated for all combinations of the coordinate variables.
  • dict: keys must match the function parameters in calculate_ts_single(). TS values will be calculated for all combinations of the dict values.

TYPE: dict | DataFrame | DataArray

multiprocess

Split the ts calculation across CPU cores. Multiprocessing is currently provided by mapply. For more sophisticated uses it may be preferred to use a multiprocessing package of your choice directly on the calculate_ts_single() method. See the code in this method (calculate_ts()) for an example.

TYPE: bool DEFAULT: False

expand

Only applicable if data is a dict. If True, will use as_dataframe() to expand the dict into a DataFrame with one column per dict key and return that, adding a column named ts for the results.

TYPE: bool DEFAULT: False

inplace

Only applicable if data is a DataFrame. If True, the results will be added to the input DataFrame in a column named ts. If a ts column already exists, it is overwritten.

TYPE: bool DEFAULT: False

progress

If True, will produce a progress bar while running models

TYPE: bool DEFAULT: False

RETURNS DESCRIPTION
None | list[float] | Series | DataFrame

The return type and value are determined by the type of the input variable (data) and the expand and inplace parameters:

  • dict input and expand=False returns a list of floats.
  • dict input and expand=True returns a DataFrame.
  • DataFrame input and inplace=False returns a Series.
  • DataFrame input and inplace=True modifies data and returns None.
  • DataArray input always modifies data and returns None.
Source code in src/echosms/scattermodelbase.py 
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
def calculate_ts(self, data: dict | pd.DataFrame | xr.DataArray, expand: bool=False,
                 inplace: bool=False, multiprocess: bool=False,
                 progress: bool=False) -> None | list[float] | pd.Series | pd.DataFrame:
    """Calculate the target strength (TS) for many parameters.

    Parameters
    ----------
    data :
        Requirements for the different input data types are:

        - **DataFrame**: column names must match the function parameter names in
          calculate_ts_single(). One TS value will be calculated for each row in the DataFrame.
        - **DataArray**: dimension names must match the function parameter names in
          calculate_ts_single(). TS values will be calculated for all combinations of the
          coordinate variables.
        - **dict**: keys must match the function parameters in calculate_ts_single().
          TS values will be calculated for all combinations of the dict values.

    multiprocess :
        Split the ts calculation across CPU cores. Multiprocessing is currently provided by
        [mapply](https://github.com/ddelange/mapply). For more
        sophisticated uses it may be preferred to use a multiprocessing package of your choice
        directly on the `calculate_ts_single()` method. See the code in this method
        (`calculate_ts()`) for an example.

    expand :
        Only applicable if `data` is a dict. If `True`, will use
        [`as_dataframe()`][echosms.utils.as_dataframe]
        to expand the dict into a DataFrame with one column per dict key
        and return that, adding a column named `ts` for the results.

    inplace :
        Only applicable if `data` is a DataFrame. If `True`, the results
        will be added to the input DataFrame in a column named `ts`. If a `ts` column
        already exists, it is overwritten.

    progress :
        If `True`, will produce a progress bar while running models

    Returns
    -------
    :
        The return type and value are determined by the type of the input variable (`data`) and
        the `expand` and `inplace` parameters:

        - dict input and `expand=False` returns a list of floats.
        - dict input and `expand=True` returns a DataFrame.
        - DataFrame input and `inplace=False` returns a Series.
        - DataFrame input and `inplace=True` modifies `data` and returns `None`.
        - DataArray input always modifies `data` and returns `None`.

    """
    match data:
        case dict():
            data_df = as_dataframe(data, self.no_expand_parameters)
        case pd.DataFrame():
            data_df = data
        case xr.DataArray():
            data_df = data.to_dataframe().reset_index()
            data_df.attrs = data.attrs
        case _:
            raise ValueError(f'Data type of {type(data)} is not supported'
                             ' (only dictionaries, Pandas DataFrames and'
                             ' Xarray DataArrays are).')

    self.validate_parameters(data_df)

    # Get the non-expandable model parameters
    p = data_df.attrs['parameters'] if 'parameters' in data_df.attrs else {}

    # Note: the args argument in the apply call below requires a tuple. data_df.attrs is a
    # dict and the default behaviour is to make a tuple using the dict keys. The trailing comma
    # and parenthesis instead causes the tuple to have one entry of the dict.

    if multiprocess:
        from mapply.mapply import mapply
        ts = mapply(data_df, self.__ts_helper, args=(p,), axis=1, progressbar=progress)
    else:  # this uses just one CPU
        if progress:
            tqdm.pandas(desc=self.short_name, unit=' models',
                        bar_format='{l_bar}{bar} [{n_fmt}/{total_fmt}; {rate_noinv_fmt}]')
            ts = data_df.progress_apply(self.__ts_helper, args=(p,), axis=1)
        else:
            ts = data_df.apply(self.__ts_helper, args=(p,), axis=1)

    match data:
        case dict() if expand:
            data_df['ts'] = ts
            return data_df
        case dict():
            return ts.to_list()
        case pd.DataFrame() if inplace:
            data_df['ts'] = ts
            return None
        case pd.DataFrame():
            return ts.rename('ts', inplace=True)
        case xr.DataArray():
            data.values = ts.to_numpy().reshape(data.shape)
            return None
        case _:
            raise AssertionError('This code should never be reached - unsupported input data '
                                 f'type of {type(data)}.')

calculate_ts_single(medium_c, medium_rho, theta, phi, f, target_c, target_rho, a, rv_pos, rv_tan, phase_sd=0, num_runs=1, validate_parameters=True, **kwargs)

Distorted-wave Born approximation scattering model.

Implements the distorted-wave Born approximation and stochastic DWBA model for calculating the acoustic backscatter from weakly scattering bodies.

PARAMETER DESCRIPTION
medium_c

Sound speed in the fluid medium surrounding the target [m/s].

TYPE: float

medium_rho

Density of the fluid medium surrounding the target [kg/m³].

TYPE: float

theta

Pitch angle to calculate the scattering as per the echoSMs coordinate system [°].

TYPE: float

phi

Roll angle to calculate the scattering as per the echoSMs coordinate system [°].

TYPE: float

f

Frequency to run the model at [Hz]

TYPE: float

target_c

Sound speed in the fluid inside the target [m/s].

TYPE: float

target_rho

Density of the fluid inside the target [kg/m³].

TYPE: float

a

The radii of the discs that define the target shape [m].

TYPE: Iterable[float]

rv_pos

An interable of vectors of the 3D positions of the centre of each disc that defines the target shape. Each vector should have three values corresponding to the x, y, and z coordinates [m] of the disc centre.

TYPE: Iterable[ndarray]

rv_tan

An interable of unit vectors of the tangent to the target body axis at the points given in rv_pos. Each vector should have three values corresponding to the x, y, and z components of the tangent vector.

TYPE: Iterable[ndarray]

phase_sd

If non-zero, this model becomes the SDWBA (stochastic DWBA). A random phase is applied to each term in the DWBA integral, obtained from a Gaussian distribution centered on 0 with standard deviation of phase_sd [°].

TYPE: float DEFAULT: 0

num_runs

The number of times to run the SDWBA model. The mean TS (calculated in the linear domain) is returned. Intended to be used in conjunction with phase_sd.

TYPE: int DEFAULT: 1

validate_parameters

Whether to validate the model parameters.

TYPE: bool DEFAULT: True

RETURNS DESCRIPTION
float

The target strength (re 1 m²) [dB] of the target.
Notes

This class implements the method presented in Eqn (5) of Stanton et al. (1998). This caters to objects that are discretised into a set of adjacent discs with potentially offset centres.

The DWBA model density and sound speed values are often specified as ratios of the medium and target values (g & h). However, to maintain compatibility with other echoSMs models this implementation requires actual densities and sound speeds.

The SDWBA component is as per Demer & Conti (2003).

References

Demer, D. A., & Conti, S. G. (2003). Reconciling theoretical versus empirical target strengths of krill: Effects of phase variability on the distorted-wave Born approximation. ICES Journal of Marine Science, 60, 429-434. https://doi.org/10.1016/S1054-3139(03)00002-X

Stanton, T. K., Chu, D., & Wiebe, P. H. (1998). Sound scattering by several zooplankton groups. II. Scattering models. The Journal of the Acoustical Society of America, 103(1), 236-253. https://doi.org/10.1121/1.421110

Source code in src/echosms/dwbamodel.py 
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
def calculate_ts_single(self, medium_c: float, medium_rho: float, theta: float, phi: float,
                        f: float, target_c: float, target_rho: float,
                        a: Iterable[float], rv_pos: Iterable[np.ndarray],
                        rv_tan: Iterable[np.ndarray],
                        phase_sd: float=0, 
                        num_runs: int=1, validate_parameters: bool=True,
                        **kwargs) -> float:
    """Distorted-wave Born approximation scattering model.

    Implements the distorted-wave Born approximation and stochastic DWBA
    model for calculating the acoustic backscatter from weakly scattering bodies.

    Parameters
    ----------
    medium_c :
        Sound speed in the fluid medium surrounding the target [m/s].
    medium_rho :
        Density of the fluid medium surrounding the target [kg/m³].
    theta :
        Pitch angle to calculate the scattering as per the echoSMs
        [coordinate system](conventions.md#coordinate-systems) [°].
    phi :
        Roll angle to calculate the scattering as per the echoSMs
        [coordinate system](conventions.md#coordinate-systems) [°].
    f :
        Frequency to run the model at [Hz]
    target_c :
        Sound speed in the fluid inside the target [m/s].
    target_rho :
        Density of the fluid inside the target [kg/m³].
    a :
        The radii of the discs that define the target shape [m].
    rv_pos :
        An interable of vectors of the 3D positions of the centre of each disc that
        defines the target shape. Each vector should have three values corresponding to
        the _x_, _y_, and _z_ coordinates [m] of the disc centre.
    rv_tan :
        An interable of unit vectors of the tangent to the target body axis at
        the points given in `rv_pos`. Each vector should have three values corresponding to
        the _x_, _y_, and _z_ components of the tangent vector.
    phase_sd :
        If non-zero, this model becomes the SDWBA (stochastic DWBA). A random phase is
        applied to each term in the DWBA integral, obtained from a Gaussian distribution
        centered on 0 with standard deviation of `phase_sd` [°].
    num_runs :
        The number of times to run the SDWBA model. The mean TS (calculated in the
        linear domain) is returned. Intended to be used in conjunction with `phase_sd`.
    validate_parameters :
        Whether to validate the model parameters.

    Returns
    -------
    :
        The target strength (re 1 m²) [dB] of the target.

    Notes
    -----
    This class implements the method presented in Eqn (5) of Stanton et al. (1998). This
    caters to objects that are discretised into a set of adjacent discs with potentially
    offset centres.

    The DWBA model density and sound speed values are often specified as ratios of the medium
    and target values (g & h). However, to maintain compatibility with other echoSMs models
    this implementation requires actual densities and sound speeds.

    The SDWBA component is as per Demer & Conti (2003).

    References
    ----------
    Demer, D. A., & Conti, S. G. (2003). Reconciling theoretical versus empirical target
    strengths of krill: Effects of phase variability on the distorted-wave Born approximation.
    ICES Journal of Marine Science, 60, 429-434.
    <https://doi.org/10.1016/S1054-3139(03)00002-X>

    Stanton, T. K., Chu, D., & Wiebe, P. H. (1998). Sound scattering by several zooplankton
    groups. II. Scattering models. The Journal of the Acoustical Society of America, 103(1),
    236-253. <https://doi.org/10.1121/1.421110>
    """
    if validate_parameters:
        self.validate_parameters(locals())

    do_sdwba = False if phase_sd == 0.0 else True

    # The structure of this code follows closely the formulae in Stanton et al (1998). Where
    # relevant, the equation numbers from that paper are given.

    k1 = wavenumber(medium_c, f)
    k2 = wavenumber(target_c, f)
    kappa_1 = 1.0 / (medium_rho * medium_c**2)  # Eqn (4) for medium
    kappa_2 = 1.0 / (target_rho * target_c**2)  # Eqn (4) for target
    gamma_k = (kappa_2 - kappa_1) / kappa_1  # Eqn (2)
    gamma_rho = (target_rho - medium_rho) / target_rho  # Eqn (3)
    # Note: the (gamma_k - gamma_rho) term in Eqn (5) can be simplified using g & h to:
    # 1/gh^2 + 1/g - 2.

    # Calculate the distance between each disc using the disc position vectors. The Euclidean
    # distance is the L2 norm so use np.norm(). Could also use
    # scipy.spatial.distance.euclidean(), but that is apparently a lot slower.
    dist = np.array([np.linalg.norm(r1-r0) for r0, r1 in zip(rv_pos[0:-1], rv_pos[1:])])  # [m]

    # Thickness of each disc based on the distance between discs. The first and last
    # discs are treated differently.
    d_rv_pos = np.hstack((dist[0], dist[0:-1]/2 + dist[1:]/2, dist[-1]))  # [m]

    # Calculate direction of incident wave given theta and phi. The echoSMs convention
    # has the target rotating and the incident vector always being (0,0,1), but for the DWBA
    # we keep the body stationary and change the incident vector.
    rot = R.from_euler('ZYX', (0, theta-90, -phi), degrees=True)
    k_hat_i = rot.apply([0, 0, 1])  # needs to be a unit vector
    k_i2 = k_hat_i * k2  # incident vector with magnitude equal to wavenumber inside the target

    # This code is a little complex because it does both the DWBA and SDWBA
    phase_factors = np.ones(len(a))  # for DWBA
    runs = []

    for run in range(int(num_runs)):  # is only ever 1 run for the DWBA
        if do_sdwba:
            phase_factors = np.exp(1j*self.rng.normal(scale=radians(phase_sd), size=len(a)))

        integral = 0.0
        for a_, r_pos, d_r_pos, r_tan, p in zip(a, rv_pos, d_rv_pos, rv_tan, phase_factors):
            # The round() is here because sometimes the dot product gets values slightly outside
            # the [-1, 1] range (due to floating point inaccuracies) and cos will complain.
            cbeta_tilt = cos(pi/2 - acos(round(k_hat_i @ r_tan, 8)))

            # This is the integral part of Eqn (5) with addition of
            # Eqn (4) from Demer & Conti (2003) for the SDWBA part
            integral += (gamma_k-gamma_rho) * exp(2j*(k_i2@r_pos))\
                * a_*j1(2*k2*a_*cbeta_tilt) / cbeta_tilt * abs(d_r_pos) * p
        runs.append(integral)

    return 20*log10(np.mean(abs(k1/4.0*np.array(runs))))

validate_parameters(params)

Validate the model parameters.

See here for calling details.

Source code in src/echosms/dwbamodel.py 
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
def validate_parameters(self, params):
    """Validate the model parameters.

    See [here][echosms.ScatterModelBase.validate_parameters] for calling details.
    """
    p = as_dict(params)
    super()._present_and_positive(p, ['medium_rho', 'medium_c', 'target_rho', 'target_c', 'f'])

    g = p['target_rho'] / p['medium_rho']
    h = p['target_c'] / p['medium_c']

    if (np.any(g < self.g_range[0])) or np.any((g > self.g_range[1])):
        warnings.warn('Ratio of target and medium densities (g) is outside the DWBA limits.')

    if (np.any(h < self.h_range[0])) or np.any((h > self.h_range[1])):
        warnings.warn('Ratio of target and medium sound speeds (h) are '
                      'outside the DWBA limits.')

    if not np.all([isclose(1.0, np.linalg.norm(v)) for v in p['rv_tan']]):
        raise ValueError('All vectors in rv_tan must be of unit length.')

echosms.PTDWBAModel()

Bases: ScatterModelBase

Phase-tracking distorted-wave Born approximation scattering model.

Source code in src/echosms/ptdwbamodel.py 
14
15
16
17
18
19
20
21
22
def __init__(self):
    super().__init__()
    self.long_name = 'phase-tracking distorted-wave Born approximation'
    self.short_name = 'pt-dwba'
    self.analytical_type = 'approximate'
    self.boundary_types = [bt.fluid_filled]
    self.shapes = ['unrestricted voxel-based']
    self.max_ka = 20
    self.no_expand_parameters = ['volume', 'voxel_size', 'rho', 'c']

calculate_ts(data, expand=False, inplace=False, multiprocess=False, progress=False)

Calculate the target strength (TS) for many parameters.

PARAMETER DESCRIPTION
data

Requirements for the different input data types are:

  • DataFrame: column names must match the function parameter names in calculate_ts_single(). One TS value will be calculated for each row in the DataFrame.
  • DataArray: dimension names must match the function parameter names in calculate_ts_single(). TS values will be calculated for all combinations of the coordinate variables.
  • dict: keys must match the function parameters in calculate_ts_single(). TS values will be calculated for all combinations of the dict values.

TYPE: dict | DataFrame | DataArray

multiprocess

Split the ts calculation across CPU cores. Multiprocessing is currently provided by mapply. For more sophisticated uses it may be preferred to use a multiprocessing package of your choice directly on the calculate_ts_single() method. See the code in this method (calculate_ts()) for an example.

TYPE: bool DEFAULT: False

expand

Only applicable if data is a dict. If True, will use as_dataframe() to expand the dict into a DataFrame with one column per dict key and return that, adding a column named ts for the results.

TYPE: bool DEFAULT: False

inplace

Only applicable if data is a DataFrame. If True, the results will be added to the input DataFrame in a column named ts. If a ts column already exists, it is overwritten.

TYPE: bool DEFAULT: False

progress

If True, will produce a progress bar while running models

TYPE: bool DEFAULT: False

RETURNS DESCRIPTION
None | list[float] | Series | DataFrame

The return type and value are determined by the type of the input variable (data) and the expand and inplace parameters:

  • dict input and expand=False returns a list of floats.
  • dict input and expand=True returns a DataFrame.
  • DataFrame input and inplace=False returns a Series.
  • DataFrame input and inplace=True modifies data and returns None.
  • DataArray input always modifies data and returns None.
Source code in src/echosms/scattermodelbase.py 
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
def calculate_ts(self, data: dict | pd.DataFrame | xr.DataArray, expand: bool=False,
                 inplace: bool=False, multiprocess: bool=False,
                 progress: bool=False) -> None | list[float] | pd.Series | pd.DataFrame:
    """Calculate the target strength (TS) for many parameters.

    Parameters
    ----------
    data :
        Requirements for the different input data types are:

        - **DataFrame**: column names must match the function parameter names in
          calculate_ts_single(). One TS value will be calculated for each row in the DataFrame.
        - **DataArray**: dimension names must match the function parameter names in
          calculate_ts_single(). TS values will be calculated for all combinations of the
          coordinate variables.
        - **dict**: keys must match the function parameters in calculate_ts_single().
          TS values will be calculated for all combinations of the dict values.

    multiprocess :
        Split the ts calculation across CPU cores. Multiprocessing is currently provided by
        [mapply](https://github.com/ddelange/mapply). For more
        sophisticated uses it may be preferred to use a multiprocessing package of your choice
        directly on the `calculate_ts_single()` method. See the code in this method
        (`calculate_ts()`) for an example.

    expand :
        Only applicable if `data` is a dict. If `True`, will use
        [`as_dataframe()`][echosms.utils.as_dataframe]
        to expand the dict into a DataFrame with one column per dict key
        and return that, adding a column named `ts` for the results.

    inplace :
        Only applicable if `data` is a DataFrame. If `True`, the results
        will be added to the input DataFrame in a column named `ts`. If a `ts` column
        already exists, it is overwritten.

    progress :
        If `True`, will produce a progress bar while running models

    Returns
    -------
    :
        The return type and value are determined by the type of the input variable (`data`) and
        the `expand` and `inplace` parameters:

        - dict input and `expand=False` returns a list of floats.
        - dict input and `expand=True` returns a DataFrame.
        - DataFrame input and `inplace=False` returns a Series.
        - DataFrame input and `inplace=True` modifies `data` and returns `None`.
        - DataArray input always modifies `data` and returns `None`.

    """
    match data:
        case dict():
            data_df = as_dataframe(data, self.no_expand_parameters)
        case pd.DataFrame():
            data_df = data
        case xr.DataArray():
            data_df = data.to_dataframe().reset_index()
            data_df.attrs = data.attrs
        case _:
            raise ValueError(f'Data type of {type(data)} is not supported'
                             ' (only dictionaries, Pandas DataFrames and'
                             ' Xarray DataArrays are).')

    self.validate_parameters(data_df)

    # Get the non-expandable model parameters
    p = data_df.attrs['parameters'] if 'parameters' in data_df.attrs else {}

    # Note: the args argument in the apply call below requires a tuple. data_df.attrs is a
    # dict and the default behaviour is to make a tuple using the dict keys. The trailing comma
    # and parenthesis instead causes the tuple to have one entry of the dict.

    if multiprocess:
        from mapply.mapply import mapply
        ts = mapply(data_df, self.__ts_helper, args=(p,), axis=1, progressbar=progress)
    else:  # this uses just one CPU
        if progress:
            tqdm.pandas(desc=self.short_name, unit=' models',
                        bar_format='{l_bar}{bar} [{n_fmt}/{total_fmt}; {rate_noinv_fmt}]')
            ts = data_df.progress_apply(self.__ts_helper, args=(p,), axis=1)
        else:
            ts = data_df.apply(self.__ts_helper, args=(p,), axis=1)

    match data:
        case dict() if expand:
            data_df['ts'] = ts
            return data_df
        case dict():
            return ts.to_list()
        case pd.DataFrame() if inplace:
            data_df['ts'] = ts
            return None
        case pd.DataFrame():
            return ts.rename('ts', inplace=True)
        case xr.DataArray():
            data.values = ts.to_numpy().reshape(data.shape)
            return None
        case _:
            raise AssertionError('This code should never be reached - unsupported input data '
                                 f'type of {type(data)}.')

calculate_ts_single(volume, theta, phi, f, voxel_size, rho, c, validate_parameters=True, **kwargs)

Phase-tracking distorted-wave Born approximation scattering model.

Implements the phase-tracking distorted-wave Born approximation model for calculating the acoustic backscatter from weakly scattering bodies.

PARAMETER DESCRIPTION
volume

The object to be modelled as a 3D volume of voxels. Array contents should be 0 for the surrounding medium, then increasing by 1 for each additional material type (i.e., 1, 2, 3, etc). volume should be arranged as per the echoSMs coordinate system, where

  • axis 0 (rows) is the x-axis
  • axis 1 (columns) is the y-axis
  • axis 2: (slices) is the z-axis

Increasing axes indices correspond to increasing x, y, and z values.

TYPE: ndarray[int]

theta

Pitch angle to calculate the scattering as per the echoSMs coordinate system [°].

TYPE: float

phi

Roll angle to calculate the scattering as per the echoSMs coordinate system [°].

TYPE: float

f

Frequency to run the model at [Hz]

TYPE: float

voxel_size

The size of the voxels in volume [m], ordered (x, y, z). This code assumes that the voxels are cubes so y and z are currently irrelevant.

TYPE: Iterable[float]

rho

Densities of each material. Must have at least the same number of entries as unique integers in volume [kg/m³].

TYPE: Iterable[float]

c

Sound speed of each material. Must have at least the same number of entries as unique integers in volume [m/s].

TYPE: Iterable[float]

validate_parameters

Whether to validate the model parameters.

TYPE: bool DEFAULT: True

RETURNS DESCRIPTION
float

The target strength (re 1 m²) [dB] of the target.
Notes

This class implements the method presented in Jones et. al. (2009). The code is based closely on the Matlab code in Jones (2006).

References

Jones, B. A. (2006). Acoustic scattering of broadband echolocation signals from prey of Blainville's beaked whales: Modeling and analysis. Master of Science, Massachusetts Institute of Technology. https://doi.org/10.1575/1912/1283

Jones, B. A., Lavery, A. C., & Stanton, T. K. (2009). Use of the distorted wave Born approximation to predict scattering by inhomogeneous objects: Application to squid. The Journal of the Acoustical Society of America, 125(1), 73-88. https://doi.org/10.1121/1.3021298

Source code in src/echosms/ptdwbamodel.py 
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
def calculate_ts_single(self, volume: np.ndarray[int], theta: float, phi: float,
                        f: float, voxel_size: Iterable[float],
                        rho: Iterable[float], c: Iterable[float],
                        validate_parameters: bool=True, **kwargs) -> float:
    """Phase-tracking distorted-wave Born approximation scattering model.

    Implements the phase-tracking distorted-wave Born approximation
    model for calculating the acoustic backscatter from weakly scattering bodies.

    Parameters
    ----------
    volume :
        The object to be modelled as a 3D volume of voxels. Array contents should be 0
        for the surrounding medium, then increasing by 1 for each additional material
        type (i.e., 1, 2, 3, etc). `volume` should be arranged as per the echoSMs
        [coordinate system](conventions.md#coordinate-systems), where

        - axis 0 (rows) is the _x_-axis
        - axis 1 (columns) is the _y_-axis
        - axis 2: (slices) is the _z_-axis

        Increasing axes indices correspond to increasing _x_, _y_, and _z_ values.

    theta :
        Pitch angle to calculate the scattering as per the echoSMs
        [coordinate system](conventions.md#coordinate-systems) [°].

    phi :
        Roll angle to calculate the scattering as per the echoSMs
        [coordinate system](conventions.md#coordinate-systems) [°].

    f :
        Frequency to run the model at [Hz]

    voxel_size :
        The size of the voxels in `volume` [m], ordered (_x_, _y_, _z_).
        This code assumes that the voxels are cubes so _y_ and _z_ are currently irrelevant.

    rho :
        Densities of each material. Must have at least the same number of entries as unique
        integers in `volume` [kg/m³].

    c :
        Sound speed of each material. Must have at least the same number of entries as unique
        integers in `volume` [m/s].
    validate_parameters : 
        Whether to validate the model parameters.

    Returns
    -------
    : float
        The target strength (re 1 m²) [dB] of the target.

    Notes
    -----
    This class implements the method presented in Jones et. al. (2009). The code is
    based closely on the Matlab code in Jones (2006).

    References
    ----------
    Jones, B. A. (2006). Acoustic scattering of broadband echolocation signals
    from prey of Blainville's beaked whales: Modeling and analysis. Master of Science,
    Massachusetts Institute of Technology. <https://doi.org/10.1575/1912/1283>

    Jones, B. A., Lavery, A. C., & Stanton, T. K. (2009). Use of the distorted
    wave Born approximation to predict scattering by inhomogeneous objects:
    Application to squid. The Journal of the Acoustical Society of America,
    125(1), 73-88. <https://doi.org/10.1121/1.3021298>
    """
    if validate_parameters:
        self.validate_parameters(locals())

    # Make sure things are numpy arrays
    rho = np.atleast_1d(rho)
    c = np.atleast_1d(c)
    voxel_size = np.array(voxel_size)

    # volume of the voxels [m^3]
    dv = voxel_size.prod()

    # input parameter checks
    if not len(volume.shape) == 3:
        raise TypeError('The volume input variable must be 3-dimensional.')

    if not voxel_size.shape[0] == 3:
        raise TypeError('The voxel_size input variable must contain 3 items.')

    if not np.any(voxel_size > 0):
        raise ValueError('All voxel_size values must be positive.')

    if f < 0.0:
        raise ValueError('The f input variable must contain only positive values.')

    if (theta < -0.0) or (theta > 180.0):
        raise ValueError('The theta (pitch) angle must be between -180.0 and +180.0')

    if (phi < -180.0) or (phi > 180.0):
        raise ValueError('The phi (roll) angle must be between -180.0 and +180.0')

    if volume.min() != 0:
        raise ValueError('The volume input variable must contain zeros.')

    categories = np.unique(volume)
    if not len(categories == (volume.max() + 1)):
        raise ValueError('The integers in volume must include all values in the series '
                         '(0, 1, 2, ..., n), where n is the largest integer in volume.')

    if not len(rho) >= len(categories):
        raise ValueError('The target_rho variable must contain at least as many values as '
                         'unique integers in the volume variable.')

    if not len(c) >= len(categories):
        raise ValueError('The target_c variable must contain at least as many values '
                         'as unique integers in the volume variable.')

    # density and sound speed ratios for all object materials
    g = rho[1:] / rho[0]
    h = c[1:] / c[0]

    # Do the pitch and roll rotations

    # Convert echoSMs rotation angles (which are intrinsic) into extrinsic as
    # that is what ndimage.rotate() below uses.
    if phi == 0.0:  # short circuit the coordinate transformation if we can
        pitch = theta-90
        roll = 0.0
    else:
        rot = R.from_euler('ZYX', (0, theta-90, -phi), degrees=True)
        # for backscatter we don't care about yaw
        _, pitch, roll = rot.as_euler('zyz', degrees=True)

    v = ndimage.rotate(volume, pitch, axes=(0, 2), order=0)
    v = ndimage.rotate(v, roll, axes=(1, 2), order=0)

    categories = np.unique(v)  # or just take the max?

    # wavenumbers in the various media
    k = 2.0*np.pi * f / c

    # DWBA coefficients
    # amplitudes in media 1,2,...,n
    Cb = 1.0/(g * h**2) + 1.0/g - 2.0  # gamma_kappa - gamma_rho
    Ca = k[0]**2 * Cb / (4.0*np.pi)  # summation coefficient

    # Differential phase for each voxel.
    dph = np.zeros(v.shape)
    masks = []
    for i, category in enumerate(categories):
        masks.append(np.isin(v, category))
        dph[masks[i]] = k[i] * voxel_size[0]
    masks.pop(0)  # don't need to keep the category[0] mask

    # cumulative summation of phase along the z-direction
    phase = dph.cumsum(axis=2) - dph/2.0
    dA = np.zeros(phase.shape, dtype=np.complex128)

    # differential phases for each voxel
    for i, m in enumerate(masks):
        dA[m] = Ca[i] * np.exp(2.0*1j*phase[m]) * dv

    # Convert to TS
    return 20.0 * np.log10(np.abs(dA.sum()))

validate_parameters(params)

Validate the model parameters.

See here for calling details.

Source code in src/echosms/ptdwbamodel.py 
24
25
26
27
28
29
def validate_parameters(self, params):
    """Validate the model parameters.

    See [here][echosms.ScatterModelBase.validate_parameters] for calling details.
    """
    p = as_dict(params)

echosms.dwbautils

Miscellaneous functions for the DWBA models.

DWBAdata()

Example datasets for the SDWBA and DWBA models.

Source code in src/echosms/dwbautils.py 
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
def __init__(self):
    # Load in the shapes data
    self.file = Path(__file__).parent/Path('resources')/Path('DWBA_shapes.toml')
    with open(self.file, 'rb') as f:
        try:
            shapes = tomllib.load(f)
        except tomllib.TOMLDecodeError as e:
            raise SyntaxError(f'Error while parsing file "{self.defs_filename.name}"') from e

    # Put the shapes into a dict of DWBAorganism().
    self.dwba_models = {}
    for s in shapes['shape']:
        mx = max(s['x'])
        s['x'] = np.array(s['x'])-mx  # put head at x=0
        # Estimate rv_tan from a spline through (x,y,z).
        tck, u = splprep([s['x'], s['y'], s['z']])
        rv_tan = np.vstack(splev(u, tck, der=1))
        # Make sure rv_tan holds only unit vectors
        n = np.linalg.norm(np.vstack(rv_tan), axis=0)
        rv_tan = (rv_tan / n).T

        # Convert the x, y, z lists into a 2D array with one row for each (x,y,z) point
        rv_pos = np.vstack((np.array(s['x']), np.array(s['y']), np.array(s['z']))).T

        organism = DWBAorganism(rv_pos, np.array(s['a']), np.array(s['g']), np.array(s['h']),
                                s['name'], s.get('source', ''), s.get('note', ''), rv_tan,
                                s['aphiaid'], s['length'], s['vernacular'],)
        self.dwba_models[s['name']] = organism
as_dict()

DWBA model shapes as a dict.

RETURNS DESCRIPTION
dict

All the DWBA model shapes. The dataset name is the dict key and the value is an instance of DWBAorganism.
Source code in src/echosms/dwbautils.py 
283
284
285
286
287
288
289
290
291
292
293
def as_dict(self) -> dict:
    """DWBA model shapes as a dict.

    Returns
    -------
    :
        All the DWBA model shapes. The dataset name is the dict key and the value is an
        instance of `DWBAorganism`.

    """
    return self.dwba_models
model(name)

DWBA model shape with requested name.

PARAMETER DESCRIPTION
name

The name of a DWBA model shape.

TYPE: str

RETURNS DESCRIPTION
DWBAorganism

An instance of DWBAorganism or None if there is no model with name.
Source code in src/echosms/dwbautils.py 
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
def model(self, name: str) -> DWBAorganism:
    """DWBA model shape with requested name.

    Parameters
    ----------
    name :
        The name of a DWBA model shape.

    Returns
    -------
    :
        An instance of `DWBAorganism` or None if there is no model with `name`.

    """
    try:
        return self.dwba_models[name]
    except KeyError:
        return None
names()

Available DWBA model names.

Source code in src/echosms/dwbautils.py 
279
280
281
def names(self):
    """Available DWBA model names."""
    return [*self.dwba_models]

DWBAorganism(rv_pos, a, g, h, name, source='', note='', rv_tan=None, aphiaid=1, length=0.0, vernacular_name='') dataclass

DWBA shape and property class to represent an organism.

ATTRIBUTE DESCRIPTION
rv_pos

An iterable of vectors of the 3D positions of the centre of each disc that defines the target shape. Each vector should have three values corresponding to the x, y, and z coordinates [m] of the disc centre, as per the echoSMs coordinate system.

TYPE: iterable[ndarray]

a

The radii [m] of each disc.

TYPE: ndarray

g

The density contrast between medium and organism (organism divied by medium).

TYPE: ndarray

h

The sound speed contrast betwwen medium and organism (organism divied by medium).

TYPE: ndarray

source

A link to the source of the data.

TYPE: str

note

Information about the data.

TYPE: str

rv_tan

An iterable of unit vectors of the tangent to the body axis at the points given by (x, y, z). Each vector has three values corresponding to the x, y, and z components of the tangent vector [m]. If not given, unit vectors pointing along the positive x-axis are used.

TYPE: ndarray

plot()

Do a simple plot of the DWBA model data.

Source code in src/echosms/dwbautils.py 
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
def plot(self):
    """Do a simple plot of the DWBA model data."""
    import matplotlib.pyplot as plt
    fig, axs = plt.subplots(2, 1, layout='compressed')
    x = self.rv_pos[:, 0]*1e3
    outline1 = (-self.rv_pos[:, 1] + self.a)*1e3
    outline2 = (-self.rv_pos[:, 1] - self.a)*1e3
    axs[0].plot(x, -self.rv_pos[:, 1]*1e3, '.-', c='C0')
    axs[0].plot(x, outline1, c='C1')
    axs[0].plot(x, outline2, c='C1')
    axs[0].plot([x[0], x[0]], [outline1[0], outline2[0]], c='C1')
    axs[0].plot([x[-1], x[-1]], [outline1[-1], outline2[-1]], c='C1')
    axs[0].set_title('Dorsal', loc='left', fontsize=8)
    axs[0].axis('scaled')
    axs[0].xaxis.set_inverted(True)

    outline1 = (-self.rv_pos[:, 2] + self.a)*1e3
    outline2 = (-self.rv_pos[:, 2] - self.a)*1e3
    axs[1].plot(x, -self.rv_pos[:, 2]*1e3, '.-', c='C0')
    axs[1].plot(x, outline1, c='C1')
    axs[1].plot(x, outline2, c='C1')
    axs[1].plot([x[0], x[0]], [outline1[0], outline2[0]], c='C1')
    axs[1].plot([x[-1], x[-1]], [outline1[-1], outline2[-1]], c='C1')
    axs[1].set_title('Lateral', loc='left', fontsize=8)
    axs[1].axis('scaled')
    axs[1].xaxis.set_inverted(True)

    plt.suptitle(self.name)
    plt.show()

create_dwba_cylinder(radius, length, spacing=0.0001)

Create shape description variables for the DWBA model for cylinders.

The shape descriptions are essentially a set of discs and their orientation.

PARAMETER DESCRIPTION
radius

The radius [m] of the cylinder.

TYPE: float

length

The length [m] of the cylinder.

TYPE: float

spacing

The spacing [m] between successive discs.

TYPE: float DEFAULT: 0.0001

RETURNS DESCRIPTION
rv_pos

An iterable of vectors of the 3D positions of the centre of each disc that defines the cylinder. Each vector has three values corresponding to the x, y, and z coordinates [m] of the disc centre.

TYPE: iterable[ndarray]

rv_tan

An iterable of unit vectors of the tangent to the cylinder body axis at the points given in rv_pos. Each vector has three values corresponding to the x, y, and z components of the tangent vector.

TYPE: iterable[ndarray]

a

The radii [m] of the discs that define the spheroid.

TYPE: iterable

Source code in src/echosms/dwbautils.py 
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
def create_dwba_cylinder(radius: float, length: float, spacing: float = 0.0001):
    """Create shape description variables for the DWBA model for cylinders.

    The shape descriptions are essentially a set of discs and their orientation.

    Parameters
    ----------
    radius :
        The radius [m] of the cylinder.
    length :
        The length [m] of the cylinder.
    spacing :
        The spacing [m] between successive discs.

    Returns
    -------
    rv_pos : iterable[np.ndarray]
        An iterable of vectors of the 3D positions of the centre of each disc that
        defines the cylinder. Each vector has three values corresponding to
        the _x_, _y_, and _z_ coordinates [m] of the disc centre.
    rv_tan : iterable[np.ndarray]
        An iterable of unit vectors of the tangent to the cylinder body axis at
        the points given in `rv_pos`. Each vector has three values corresponding to
        the _x_, _y_, and _z_ components of the tangent vector.
    a : iterable
        The radii [m] of the discs that define the spheroid.
    """
    pos = np.linspace(0, length, int(round(length/spacing)))
    rv_pos = [np.array([x, 0, 0]) for x in pos]
    rv_tan = [np.array([1, 0, 0])] * len(pos)
    a = [radius] * len(pos)

    return rv_pos, rv_tan, a

create_dwba_from_xyza(x, y, z, a, name, g=1.0, h=1.0, source='', note='')

Create a DWBAorganism instance from shape data.

Converts a centreline and radius definiton of the DWBA shape into that required by the echoSMs implementation of the DWBA (centreline, tangential, and radii vectors).

PARAMETER DESCRIPTION
x

x-coordinates [m] of the centreline of the DWBA shape as per the echoSMs coordinate system.

TYPE: Iterable[float]

y

y-coordinates [m] of the centreline of the DWBA shape as per the echoSMs coordinate system.

TYPE: Iterable[float]

z

z-coordinates [m] of the centreline of the DWBA shape as per the echoSMs coordinate system.

TYPE: Iterable[float]

a

radius [m] of the DWBA shape at each centreline (x,y,z) position.

TYPE: Iterable[float]

name

A name for the organism.

TYPE: str

source

A link/URL/DOI or description for the source of the data.

TYPE: str DEFAULT: ''

note

Notes about the organism or data.

TYPE: str DEFAULT: ''

g

A single value of g, the ratio of organism density divided by the medium density. This is applied to all parts of the shape.

TYPE: float DEFAULT: 1.0

h

A single value of h, the ratio of organism sound speed divided by the medium sound speed. This is applied to all parts of the shape.

TYPE: float DEFAULT: 1.0

RETURNS DESCRIPTION
An instance of DWBAorganism.
Notes

Here is an example of how to use this function to read in .sat files from the ZooScatR package and convert them into the format required by the echoSMs DWBA model.

import pandas as pd
from echosms import create_dwba_from_xyza

filepath = 'shape.sat'

s = pd.read_csv(filepath, delimiter=' ', names=['y', 'x', 'a'])

# Adjust the data to match the echoSMs units
s /= 1000  # convert mm to m

# .sat files don't have a z coordinate but echoSMs requires it
s['z'] = 0.0

# An example of flipping left to right to match the echoSMS coordinate system
s['x'] = max(s['x']) - s['x']

shape = create_dwba_from_xyza(s['x'], s['y'], s['z'], s['a'], name=filepath, g=1.05, h=1.05)
Source code in src/echosms/dwbautils.py 
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
def create_dwba_from_xyza(x, y, z, a, name: str, g: float = 1.0, h: float = 1.0,
                          source: str = '', note: str = ''):
    """Create a DWBAorganism instance from shape data.

    Converts a centreline and radius definiton of the DWBA shape into
    that required by the echoSMs implementation of the DWBA (centreline, tangential, and
    radii vectors).

    Parameters
    ----------
    x : Iterable[float]
        x-coordinates [m] of the centreline of the DWBA shape as per the echoSMs
        [coordinate system](conventions.md#coordinate-systems).
    y : Iterable[float]
        y-coordinates [m] of the centreline of the DWBA shape as per the echoSMs
        [coordinate system](conventions.md#coordinate-systems).
    z : Iterable[float]
        z-coordinates [m] of the centreline of the DWBA shape as per the echoSMs
        [coordinate system](conventions.md#coordinate-systems).
    a : Iterable[float]
        radius [m] of the DWBA shape at each centreline (x,y,z) position.
    name :
        A name for the organism.
    source :
        A link/URL/DOI or description for the source of the data.
    note :
        Notes about the organism or data.
    g :
        A single value of g, the ratio of organism density divided by the
        medium density. This is applied to all parts of the shape.
    h :
        A single value of h, the ratio of organism sound speed divided by
        the medium sound speed. This is applied to all parts of the shape.

    Returns
    -------
        An instance of DWBAorganism.

    Notes
    -----
    Here is an example of how to use this function to read in .sat files
    from the [ZooScatR](https://github.com/AustralianAntarcticDivision/ZooScatR)
    package and convert them into the format required
    by the echoSMs DWBA model.

    ```py
    import pandas as pd
    from echosms import create_dwba_from_xyza

    filepath = 'shape.sat'

    s = pd.read_csv(filepath, delimiter=' ', names=['y', 'x', 'a'])

    # Adjust the data to match the echoSMs units
    s /= 1000  # convert mm to m

    # .sat files don't have a z coordinate but echoSMs requires it
    s['z'] = 0.0

    # An example of flipping left to right to match the echoSMS coordinate system
    s['x'] = max(s['x']) - s['x']

    shape = create_dwba_from_xyza(s['x'], s['y'], s['z'], s['a'], name=filepath, g=1.05, h=1.05)

    ```
    """
    # Estimate rv_tan from a spline through (x,y).
    tck, u = splprep([x, y, z])
    rv_tan = np.vstack(splev(u, tck, der=1))
    # Make sure rv_tan holds only unit vectors
    n = np.linalg.norm(np.vstack(rv_tan), axis=0)
    rv_tan = (rv_tan / n).T

    # Convert the x, y, and z into a 2D array and get one row for each (x,y,z) point.
    rv_pos = np.vstack((x, y, z)).T

    gg = np.full((1, rv_pos.shape[0]), g)
    hh = np.full((1, rv_pos.shape[0]), h)

    return DWBAorganism(rv_pos, a, gg, hh, name, source, note, rv_tan)

create_dwba_spheroid(major_radius, minor_radius, spacing=0.0001)

Create shape description variables for the DWBA model for spheres and prolate spheroids.

The shape descriptions are essentially a set of discs and their orientation.

Notes

Currently only supports prolate spheroids and spheres (set major_radius and minor_radius to the same value to get a sphere).

PARAMETER DESCRIPTION
major_radius

The major radius [m] of the spheroid.

TYPE: float

minor_radius

The minor radius [m] of the spheroid.

TYPE: float

spacing

The spacing [m] between successive discs.

TYPE: float DEFAULT: 0.0001

RETURNS DESCRIPTION
rv_pos

An iterable of vectors of the 3D positions of the centre of each disc that defines the spheroid. Each vector has three values corresponding to the x, y, and z coordinates [m] of the disc centre.

TYPE: iterable[ndarray]

rv_tan

An iterable of unit vectors of the tangent to the target body axis at the points given in rv_pos. Each vector has three values corresponding to the x, y, and z components of the tangent vector.

TYPE: iterable[ndarray]

a

The radii [m] of the discs that define the spheroid.

TYPE: iterable

Source code in src/echosms/dwbautils.py 
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
def create_dwba_spheroid(major_radius: float, minor_radius: float, spacing: float = 0.0001):
    """Create shape description variables for the DWBA model for spheres and prolate spheroids.

    The shape descriptions are essentially a set of discs and their orientation.

    Notes
    -----
    Currently only supports prolate spheroids and spheres (set `major_radius` and `minor_radius`
    to the same value to get a sphere).

    Parameters
    ----------
    major_radius :
        The major radius [m] of the spheroid.
    minor_radius :
        The minor radius [m] of the spheroid.
    spacing :
        The spacing [m] between successive discs.

    Returns
    -------
    rv_pos : iterable[np.ndarray]
        An iterable of vectors of the 3D positions of the centre of each disc that
        defines the spheroid. Each vector has three values corresponding to
        the _x_, _y_, and _z_ coordinates [m] of the disc centre.
    rv_tan : iterable[np.ndarray]
        An iterable of unit vectors of the tangent to the target body axis at
        the points given in `rv_pos`. Each vector has three values corresponding to
        the _x_, _y_, and _z_ components of the tangent vector.
    a : iterable
        The radii [m] of the discs that define the spheroid.
    """
    v = np.linspace(0, np.pi, int(round(2*major_radius/spacing)))
    a = minor_radius*np.sin(v)  # radius at points along the spheroid
    x = major_radius-major_radius*np.cos(v)  # shift so that origin is at one end

    # List of disc position vectors
    rv_pos = [np.array([i, 0, 0]) for i in x]

    # List of tangent vectors to sphere axis (all the same for spheroids)
    rv_tan = [np.array([1, 0, 0])] * len(x)

    return rv_pos, rv_tan, a

echosms.ESModel()

Bases: ScatterModelBase

Elastic sphere (ES) scattering model.

This class calculates acoustic backscatter from elastic spheres.

Source code in src/echosms/esmodel.py 
17
18
19
20
21
22
23
24
def __init__(self):
    super().__init__()
    self.long_name = 'elastic sphere'
    self.short_name = 'es'
    self.analytical_type = 'exact'
    self.boundary_types = [bt.elastic]
    self.shapes = ['sphere']
    self.max_ka = 20  # [1]

calculate_ts(data, expand=False, inplace=False, multiprocess=False, progress=False)

Calculate the target strength (TS) for many parameters.

PARAMETER DESCRIPTION
data

Requirements for the different input data types are:

  • DataFrame: column names must match the function parameter names in calculate_ts_single(). One TS value will be calculated for each row in the DataFrame.
  • DataArray: dimension names must match the function parameter names in calculate_ts_single(). TS values will be calculated for all combinations of the coordinate variables.
  • dict: keys must match the function parameters in calculate_ts_single(). TS values will be calculated for all combinations of the dict values.

TYPE: dict | DataFrame | DataArray

multiprocess

Split the ts calculation across CPU cores. Multiprocessing is currently provided by mapply. For more sophisticated uses it may be preferred to use a multiprocessing package of your choice directly on the calculate_ts_single() method. See the code in this method (calculate_ts()) for an example.

TYPE: bool DEFAULT: False

expand

Only applicable if data is a dict. If True, will use as_dataframe() to expand the dict into a DataFrame with one column per dict key and return that, adding a column named ts for the results.

TYPE: bool DEFAULT: False

inplace

Only applicable if data is a DataFrame. If True, the results will be added to the input DataFrame in a column named ts. If a ts column already exists, it is overwritten.

TYPE: bool DEFAULT: False

progress

If True, will produce a progress bar while running models

TYPE: bool DEFAULT: False

RETURNS DESCRIPTION
None | list[float] | Series | DataFrame

The return type and value are determined by the type of the input variable (data) and the expand and inplace parameters:

  • dict input and expand=False returns a list of floats.
  • dict input and expand=True returns a DataFrame.
  • DataFrame input and inplace=False returns a Series.
  • DataFrame input and inplace=True modifies data and returns None.
  • DataArray input always modifies data and returns None.
Source code in src/echosms/scattermodelbase.py 
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
def calculate_ts(self, data: dict | pd.DataFrame | xr.DataArray, expand: bool=False,
                 inplace: bool=False, multiprocess: bool=False,
                 progress: bool=False) -> None | list[float] | pd.Series | pd.DataFrame:
    """Calculate the target strength (TS) for many parameters.

    Parameters
    ----------
    data :
        Requirements for the different input data types are:

        - **DataFrame**: column names must match the function parameter names in
          calculate_ts_single(). One TS value will be calculated for each row in the DataFrame.
        - **DataArray**: dimension names must match the function parameter names in
          calculate_ts_single(). TS values will be calculated for all combinations of the
          coordinate variables.
        - **dict**: keys must match the function parameters in calculate_ts_single().
          TS values will be calculated for all combinations of the dict values.

    multiprocess :
        Split the ts calculation across CPU cores. Multiprocessing is currently provided by
        [mapply](https://github.com/ddelange/mapply). For more
        sophisticated uses it may be preferred to use a multiprocessing package of your choice
        directly on the `calculate_ts_single()` method. See the code in this method
        (`calculate_ts()`) for an example.

    expand :
        Only applicable if `data` is a dict. If `True`, will use
        [`as_dataframe()`][echosms.utils.as_dataframe]
        to expand the dict into a DataFrame with one column per dict key
        and return that, adding a column named `ts` for the results.

    inplace :
        Only applicable if `data` is a DataFrame. If `True`, the results
        will be added to the input DataFrame in a column named `ts`. If a `ts` column
        already exists, it is overwritten.

    progress :
        If `True`, will produce a progress bar while running models

    Returns
    -------
    :
        The return type and value are determined by the type of the input variable (`data`) and
        the `expand` and `inplace` parameters:

        - dict input and `expand=False` returns a list of floats.
        - dict input and `expand=True` returns a DataFrame.
        - DataFrame input and `inplace=False` returns a Series.
        - DataFrame input and `inplace=True` modifies `data` and returns `None`.
        - DataArray input always modifies `data` and returns `None`.

    """
    match data:
        case dict():
            data_df = as_dataframe(data, self.no_expand_parameters)
        case pd.DataFrame():
            data_df = data
        case xr.DataArray():
            data_df = data.to_dataframe().reset_index()
            data_df.attrs = data.attrs
        case _:
            raise ValueError(f'Data type of {type(data)} is not supported'
                             ' (only dictionaries, Pandas DataFrames and'
                             ' Xarray DataArrays are).')

    self.validate_parameters(data_df)

    # Get the non-expandable model parameters
    p = data_df.attrs['parameters'] if 'parameters' in data_df.attrs else {}

    # Note: the args argument in the apply call below requires a tuple. data_df.attrs is a
    # dict and the default behaviour is to make a tuple using the dict keys. The trailing comma
    # and parenthesis instead causes the tuple to have one entry of the dict.

    if multiprocess:
        from mapply.mapply import mapply
        ts = mapply(data_df, self.__ts_helper, args=(p,), axis=1, progressbar=progress)
    else:  # this uses just one CPU
        if progress:
            tqdm.pandas(desc=self.short_name, unit=' models',
                        bar_format='{l_bar}{bar} [{n_fmt}/{total_fmt}; {rate_noinv_fmt}]')
            ts = data_df.progress_apply(self.__ts_helper, args=(p,), axis=1)
        else:
            ts = data_df.apply(self.__ts_helper, args=(p,), axis=1)

    match data:
        case dict() if expand:
            data_df['ts'] = ts
            return data_df
        case dict():
            return ts.to_list()
        case pd.DataFrame() if inplace:
            data_df['ts'] = ts
            return None
        case pd.DataFrame():
            return ts.rename('ts', inplace=True)
        case xr.DataArray():
            data.values = ts.to_numpy().reshape(data.shape)
            return None
        case _:
            raise AssertionError('This code should never be reached - unsupported input data '
                                 f'type of {type(data)}.')

calculate_ts_single(medium_c, medium_rho, a, f, target_longitudinal_c, target_transverse_c, target_rho, validate_parameters=True, **kwargs)

Calculate the backscatter from an elastic sphere for one set of parameters.

PARAMETER DESCRIPTION
medium_c

Sound speed in the fluid medium surrounding the sphere [m/s].

TYPE: float

medium_rho

Density of the fluid medium surrounding the sphere [kg/m³].

TYPE: float

a

Radius of the sphere [m].

TYPE: float

f

Frequency to calculate the scattering at [Hz].

TYPE: float

target_longitudinal_c

Longitudinal sound speed in the material inside the sphere [m/s].

TYPE: float

target_transverse_c

Transverse sound speed in the material inside the sphere [m/s].

TYPE: float

target_rho

Density of the material inside the sphere [kg/m³].

TYPE: float

validate_parameters

Whether to validate the model parameters.

TYPE: bool DEFAULT: True

RETURNS DESCRIPTION
float

The target strength (re 1 m²) of the sphere [dB].
Notes

The class implements the code in MacLennan (1981).

References

MacLennan, D. N. (1981). The Theory of Solid Spheres as Sonar Calibration Targets. Scottish Fisheries Research Report Number 22. Department of Agriculture and Fisheries for Scotland. https://archive.org/details/scottish-fisheries-research-reports_1981_22

Source code in src/echosms/esmodel.py 
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
def calculate_ts_single(self, medium_c: float, medium_rho: float, a: float, f: float,
                        target_longitudinal_c: float, target_transverse_c: float,
                        target_rho: float, validate_parameters: bool=True,
                        **kwargs) -> float:
    """
    Calculate the backscatter from an elastic sphere for one set of parameters.

    Parameters
    ----------
    medium_c :
        Sound speed in the fluid medium surrounding the sphere [m/s].
    medium_rho :
        Density of the fluid medium surrounding the sphere [kg/m³].
    a :
        Radius of the sphere [m].
    f :
        Frequency to calculate the scattering at [Hz].
    target_longitudinal_c :
        Longitudinal sound speed in the material inside the sphere [m/s].
    target_transverse_c :
        Transverse sound speed in the material inside the sphere [m/s].
    target_rho :
        Density of the material inside the sphere [kg/m³].
    validate_parameters :
        Whether to validate the model parameters.

    Returns
    -------
    :
        The target strength (re 1 m²) of the sphere [dB].

    Notes
    -----
    The class implements the code in MacLennan (1981).

    References
    ----------
    MacLennan, D. N. (1981). The Theory of Solid Spheres as Sonar Calibration Targets.
    Scottish Fisheries Research Report Number 22. Department of Agriculture and Fisheries
    for Scotland. <https://archive.org/details/scottish-fisheries-research-reports_1981_22>
    """
    if validate_parameters:
        self.validate_parameters(locals())

    q = wavenumber(medium_c, f)*a
    q1 = q*medium_c/target_longitudinal_c
    q2 = q*medium_c/target_transverse_c
    alpha = 2. * (target_rho/medium_rho) * (target_transverse_c/medium_c)**2
    beta = (target_rho/medium_rho) * (target_longitudinal_c/medium_c)**2 - alpha

    # Use n instead of l (ell) because l looks like 1.
    def S(n):
        A2 = (n**2 + n-2) * spherical_jn(n, q2) + q2**2 * spherical_jnpp(n, q2)
        A1 = 2*n*(n+1) * (q1*spherical_jn(n, q1, True) - spherical_jn(n, q1))
        B2 = A2*q1**2 * (beta*spherical_jn(n, q1) - alpha*spherical_jnpp(n, q1))\
            - A1*alpha * (spherical_jn(n, q2) - q2*spherical_jn(n, q2, True))
        B1 = q * (A2*q1*spherical_jn(n, q1, True) - A1*spherical_jn(n, q2))
        eta_n = atan(-(B2*spherical_jn(n, q, True) - B1*spherical_jn(n, q))
                     / (B2*spherical_yn(n, q, True) - B1*spherical_yn(n, q)))

        return (-1)**n * (2*n+1) * sin(eta_n) * exp(1j*eta_n)

    # Estimate the number of terms to use in the summation
    n_max = round(q+10)
    tol = 1e-10  # somewhat arbitrary
    while abs(S(n_max)) > tol:
        n_max += 10

    if n_max > 200:
        warn('TS results may be inaccurate because the modal series required a large '
             f'number ({n_max}) of terms to converge.')

    n = range(n_max)

    f_inf = -2.0/q * sum(map(S, n))

    return 10*log10(a**2 * abs(f_inf)**2 / 4.0)

validate_parameters(params)

Validate the model parameters.

See here for calling details.

Source code in src/echosms/esmodel.py 
26
27
28
29
30
31
32
33
34
35
def validate_parameters(self, params):
    """Validate the model parameters.

    See [here][echosms.ScatterModelBase.validate_parameters] for calling details.
    """
    p = as_dict(params)
    super()._present_and_in(p, ['boundary_type'], self.boundary_types)
    super()._present_and_positive(p, ['medium_rho', 'medium_c', 'a', 'f',
                                      'target_longitudinal_c',
                                      'target_transverse_c', 'target_rho'])

echosms.HPModel()

Bases: ScatterModelBase

High-pass (HP) scattering model.

Source code in src/echosms/hpmodel.py 
11
12
13
14
15
16
17
18
def __init__(self):
    super().__init__()
    self.long_name = 'high pass'
    self.short_name = 'hp'
    self.analytical_type = 'approximate'
    self.boundary_types = [bt.fluid_filled, bt.elastic, bt.fixed_rigid]
    self.shapes = ['sphere', 'prolate spheroid', 'cylinder', 'bent cylinder']
    self.max_ka = 20  # [1]

calculate_ts(data, expand=False, inplace=False, multiprocess=False, progress=False)

Calculate the target strength (TS) for many parameters.

PARAMETER DESCRIPTION
data

Requirements for the different input data types are:

  • DataFrame: column names must match the function parameter names in calculate_ts_single(). One TS value will be calculated for each row in the DataFrame.
  • DataArray: dimension names must match the function parameter names in calculate_ts_single(). TS values will be calculated for all combinations of the coordinate variables.
  • dict: keys must match the function parameters in calculate_ts_single(). TS values will be calculated for all combinations of the dict values.

TYPE: dict | DataFrame | DataArray

multiprocess

Split the ts calculation across CPU cores. Multiprocessing is currently provided by mapply. For more sophisticated uses it may be preferred to use a multiprocessing package of your choice directly on the calculate_ts_single() method. See the code in this method (calculate_ts()) for an example.

TYPE: bool DEFAULT: False

expand

Only applicable if data is a dict. If True, will use as_dataframe() to expand the dict into a DataFrame with one column per dict key and return that, adding a column named ts for the results.

TYPE: bool DEFAULT: False

inplace

Only applicable if data is a DataFrame. If True, the results will be added to the input DataFrame in a column named ts. If a ts column already exists, it is overwritten.

TYPE: bool DEFAULT: False

progress

If True, will produce a progress bar while running models

TYPE: bool DEFAULT: False

RETURNS DESCRIPTION
None | list[float] | Series | DataFrame

The return type and value are determined by the type of the input variable (data) and the expand and inplace parameters:

  • dict input and expand=False returns a list of floats.
  • dict input and expand=True returns a DataFrame.
  • DataFrame input and inplace=False returns a Series.
  • DataFrame input and inplace=True modifies data and returns None.
  • DataArray input always modifies data and returns None.
Source code in src/echosms/scattermodelbase.py 
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
def calculate_ts(self, data: dict | pd.DataFrame | xr.DataArray, expand: bool=False,
                 inplace: bool=False, multiprocess: bool=False,
                 progress: bool=False) -> None | list[float] | pd.Series | pd.DataFrame:
    """Calculate the target strength (TS) for many parameters.

    Parameters
    ----------
    data :
        Requirements for the different input data types are:

        - **DataFrame**: column names must match the function parameter names in
          calculate_ts_single(). One TS value will be calculated for each row in the DataFrame.
        - **DataArray**: dimension names must match the function parameter names in
          calculate_ts_single(). TS values will be calculated for all combinations of the
          coordinate variables.
        - **dict**: keys must match the function parameters in calculate_ts_single().
          TS values will be calculated for all combinations of the dict values.

    multiprocess :
        Split the ts calculation across CPU cores. Multiprocessing is currently provided by
        [mapply](https://github.com/ddelange/mapply). For more
        sophisticated uses it may be preferred to use a multiprocessing package of your choice
        directly on the `calculate_ts_single()` method. See the code in this method
        (`calculate_ts()`) for an example.

    expand :
        Only applicable if `data` is a dict. If `True`, will use
        [`as_dataframe()`][echosms.utils.as_dataframe]
        to expand the dict into a DataFrame with one column per dict key
        and return that, adding a column named `ts` for the results.

    inplace :
        Only applicable if `data` is a DataFrame. If `True`, the results
        will be added to the input DataFrame in a column named `ts`. If a `ts` column
        already exists, it is overwritten.

    progress :
        If `True`, will produce a progress bar while running models

    Returns
    -------
    :
        The return type and value are determined by the type of the input variable (`data`) and
        the `expand` and `inplace` parameters:

        - dict input and `expand=False` returns a list of floats.
        - dict input and `expand=True` returns a DataFrame.
        - DataFrame input and `inplace=False` returns a Series.
        - DataFrame input and `inplace=True` modifies `data` and returns `None`.
        - DataArray input always modifies `data` and returns `None`.

    """
    match data:
        case dict():
            data_df = as_dataframe(data, self.no_expand_parameters)
        case pd.DataFrame():
            data_df = data
        case xr.DataArray():
            data_df = data.to_dataframe().reset_index()
            data_df.attrs = data.attrs
        case _:
            raise ValueError(f'Data type of {type(data)} is not supported'
                             ' (only dictionaries, Pandas DataFrames and'
                             ' Xarray DataArrays are).')

    self.validate_parameters(data_df)

    # Get the non-expandable model parameters
    p = data_df.attrs['parameters'] if 'parameters' in data_df.attrs else {}

    # Note: the args argument in the apply call below requires a tuple. data_df.attrs is a
    # dict and the default behaviour is to make a tuple using the dict keys. The trailing comma
    # and parenthesis instead causes the tuple to have one entry of the dict.

    if multiprocess:
        from mapply.mapply import mapply
        ts = mapply(data_df, self.__ts_helper, args=(p,), axis=1, progressbar=progress)
    else:  # this uses just one CPU
        if progress:
            tqdm.pandas(desc=self.short_name, unit=' models',
                        bar_format='{l_bar}{bar} [{n_fmt}/{total_fmt}; {rate_noinv_fmt}]')
            ts = data_df.progress_apply(self.__ts_helper, args=(p,), axis=1)
        else:
            ts = data_df.apply(self.__ts_helper, args=(p,), axis=1)

    match data:
        case dict() if expand:
            data_df['ts'] = ts
            return data_df
        case dict():
            return ts.to_list()
        case pd.DataFrame() if inplace:
            data_df['ts'] = ts
            return None
        case pd.DataFrame():
            return ts.rename('ts', inplace=True)
        case xr.DataArray():
            data.values = ts.to_numpy().reshape(data.shape)
            return None
        case _:
            raise AssertionError('This code should never be reached - unsupported input data '
                                 f'type of {type(data)}.')

calculate_ts_single(shape, medium_c, a, f, boundary_type, medium_rho=None, target_c=None, target_rho=None, theta=None, L=None, rho_c=None, irregular=False, validate_parameters=True, **kwargs)

Calculate the backscatter using the high pass model for one set of parameters.

PARAMETER DESCRIPTION
shape

The shape to model. Must be one of the shapes given in the shapes variable.

TYPE: str

medium_c

Sound speed in the fluid medium surrounding the target [m/s].

TYPE: float

medium_rho

Density of the fluid medium surrounding the target [kg/m³]. Not required when boundary_type is fixed_rigid.

TYPE: None | float DEFAULT: None

target_c

Longitudinal sound speed in the material inside the target [m/s]. Not required when boundary_type is fixed_rigid.

TYPE: None | float DEFAULT: None

target_rho

Density of the material inside the target [kg/m³]. Not required when boundary_type is fixed_rigid.

TYPE: None | float DEFAULT: None

a

Radius of the sphere, length of semi-minor axis of the prolate spheriod, or cylindrical radius of the straight or bent cylinder [m].

TYPE: float

f

Frequency to calculate the scattering at [Hz].

TYPE: float

boundary_type

The boundary type for the model.

TYPE: boundary_type

theta

Pitch angle to calculate the scattering as per the echoSMs coordinate system [°]. Only required for the straight cylinder shape.

TYPE: None | float DEFAULT: None

L

Total length of the prolate spheroid and straight cylinder, or arc length of the bent cylinder [m]. Only required for prolate spheroid, cylinder, and bent cylinder shapes.

TYPE: None | float DEFAULT: None

rho_c

Radius of curvature of the axis of the bent cylinder [m]. Only required for the bent cylinder shape.

TYPE: None | float DEFAULT: None

irregular

Set to True if the modelled object is not exactly a sphere, prolate spheroid, straight or uniformly beny cylinder.

TYPE: bool DEFAULT: False

validate_parameters

Whether to validate the model parameters.

TYPE: bool DEFAULT: True

RETURNS DESCRIPTION
float

The target strength (re 1 m²) of the sphere [dB].
Notes

The class implements the high-pass model in Stanton (1989) for spheres, prolate spheroids, cylinders, and bent cylinders with fluid filled, elastic, and rigid fixed boundary conditions. There are several restrictions on valid input parameters, so a careful reading of Stanton (1989) is recommended.

The theta angle convention used in Stanton (1989) is the same as the echoSMs coordinate system convention.

Stanton (1989) also provides parameters for gas-filled shapes, but more prior knowledge is required about the gas for useful results (e.g., damping characteristics of the gas and medium) and these have not been implemented here. There are many other models available that accurately model gas-filled shapes such that the lack of the high-pass gas-filled model should not be missed.

References

Stanton, T. K. (1989). Simple approximate formulas for backscattering of sound by spherical and elongated objects. The Journal of the Acoustical Society of America, 86(4), 1499-1510. https://doi.org/10.1121/1.398711

Source code in src/echosms/hpmodel.py 
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
def calculate_ts_single(self, shape: str, medium_c: float, a: float, f: float,
                        boundary_type: bt, medium_rho: None | float=None,
                        target_c: None | float=None, target_rho: None | float=None,
                        theta: None | float=None,
                        L: None | float=None, rho_c: None | float=None,
                        irregular: bool=False,
                        validate_parameters: bool=True, **kwargs) -> float:
    """
    Calculate the backscatter using the high pass model for one set of parameters.

    Parameters
    ----------
    shape :
        The shape to model. Must be one of the shapes given in the `shapes` variable.
    medium_c :
        Sound speed in the fluid medium surrounding the target [m/s].
    medium_rho :
        Density of the fluid medium surrounding the target [kg/m³]. Not required when
        `boundary_type` is `fixed_rigid`.
    target_c :
        Longitudinal sound speed in the material inside the target [m/s]. Not required when
        `boundary_type` is `fixed_rigid`.
    target_rho :
        Density of the material inside the target [kg/m³]. Not required when
        `boundary_type` is `fixed_rigid`.
    a :
        Radius of the sphere, length of semi-minor axis of the prolate spheriod, or cylindrical
        radius of the straight or bent cylinder [m].
    f :
        Frequency to calculate the scattering at [Hz].
    boundary_type :
        The boundary type for the model.
    theta :
        Pitch angle to calculate the scattering as per the echoSMs
        [coordinate system](conventions.md#coordinate-systems) [°]. Only required for the straight cylinder shape.
    L :
        Total length of the prolate spheroid and straight cylinder, or arc length of
        the bent cylinder [m]. Only required for prolate spheroid, cylinder, and bent cylinder
        shapes.
    rho_c :
        Radius of curvature of the axis of the bent cylinder [m]. Only required for the
        bent cylinder shape.
    irregular :
        Set to `True` if the modelled object is not exactly a sphere, prolate spheroid,
        straight or uniformly beny cylinder.
    validate_parameters :
        Whether to validate the model parameters.

    Returns
    -------
    :
        The target strength (re 1 m²) of the sphere [dB].

    Notes
    -----
    The class implements the high-pass model in Stanton (1989) for spheres, prolate spheroids,
    cylinders, and bent cylinders with fluid filled, elastic, and rigid fixed boundary
    conditions. There are several restrictions on valid input parameters, so a careful
    reading of Stanton (1989) is recommended.

    The theta angle convention used in Stanton (1989) is the same as the echoSMs
    [coordinate system convention](conventions.md#coordinate-systems).

    Stanton (1989) also provides parameters for gas-filled shapes, but more
    prior knowledge is required about the gas for useful results (e.g., damping
    characteristics of the gas and medium) and these have not been implemented here. There
    are many other models available that accurately model gas-filled shapes such that the
    lack of the high-pass gas-filled model should not be missed.

    References
    ----------
    Stanton, T. K. (1989). Simple approximate formulas for backscattering of sound
    by spherical and elongated objects. The Journal of the Acoustical Society of
    America, 86(4), 1499-1510.
    <https://doi.org/10.1121/1.398711>
    """
    if validate_parameters:
        self.validate_parameters(locals())

    if boundary_type == bt.fixed_rigid:
        # just need something large
        g = 1e20
        h = 1e20
    else:
        g = target_rho/medium_rho
        h = target_c/medium_c

    k = wavenumber(medium_c, f)

    G = 1.0
    F = 1.0

    def alpha_pic(g, h):
        return (1-g*h*h)/(2*g*h*h) + (1-g)/(1+g)

    match shape:
        case 'sphere':
            alpha_pis = (1-g*h*h)/(3*g*h*h) + (1-g)/(1+2*g)
            R = (g*h-1)/(g*h+1)
            if irregular:
                match boundary_type:
                    case bt.fluid_filled:
                        F = 40 * (k*a)**(-0.4)
                        G = 1-0.8*exp(-2.5*(k*a-2.25)**2)
                    case bt.elastic:
                        F = 15 * (k*a)**(-1.9)
                    case bt.fixed_rigid:
                        F = 15 * (k*a)**(-1.9)

            sigma_bs = a*a * (k*a)**4 * alpha_pis**2 * G\
                / (1 + 4*(k*a)**4 * alpha_pis**2/(R**2 * F))
        case 'prolate spheroid':
            a_pic = alpha_pic(g, h)
            if irregular:
                match boundary_type:
                    case bt.fluid_filled:
                        F = 2.5 * (k*a)**(1.65)
                        G = 1-0.8*exp(-2.5*(k*a-2.3)**2)
                    case bt.elastic:
                        F = 1.8 * (k*a)**(-0.4)
                    case bt.fixed_rigid:
                        F = 1.8 * (k*a)**(-0.4)

            sigma_bs = 1/9 * L*L * (k*a)**4 * a_pic**2 * G\
                / (1 + 16/9*(k*a)**4 * a_pic**2/(R**2 * F))
        case 'cylinder':
            theta = radians(theta)
            a_pic = alpha_pic(g, h)
            s = sin(k*L*cos(theta)) / (k*L*cos(theta))
            Ka = k*sin(theta)*a
            if irregular:
                match boundary_type:
                    case bt.fluid_filled:
                        F = 3 * (k*a)**(0.65)
                        G = 1-0.8*exp(-2.5*(k*a-2.0)**2)
                    case bt.elastic:
                        F = 3.5 * (k*a)**(-1.0)
                    case bt.fixed_rigid:
                        F = 3.5 * (k*a)**(-1.0)

            sigma_bs = 0.25 * L*L * (Ka)**4 * a_pic**2 * s*s * G\
                / (1 + pi*(Ka)**4 * a_pic**2/(R**2 * F))
        case 'bent cylinder':
            a_pic = alpha_pic(g, h)
            H = 1.
            if irregular:
                match boundary_type:
                    case bt.fluid_filled:
                        F = 3.0 * (k*a)**(0.65)
                        G = 1-0.8*exp(-2.5*(k*a-2.0)**2)
                    case bt.elastic:
                        F = 2.5 * (k*a)**(-1.0)
                    case bt.fixed_rigid:
                        F = 2.5 * (k*a)**(-1.0)

            sigma_bs = 0.25 * L*L * (k*a)**4 * a_pic**2 * H*H*G\
                / (1 + L*L*(k*a)**4 * a_pic**2 * H*H/(rho_c*a*R**2 * F))

    return 10*log10(sigma_bs)

validate_parameters(params)

Validate the model parameters.

See here for calling details.

Source code in src/echosms/hpmodel.py 
20
21
22
23
24
25
26
27
28
29
30
31
32
33
def validate_parameters(self, params):
    """Validate the model parameters.

    See [here][echosms.ScatterModelBase.validate_parameters] for calling details.
    """
    p = as_dict(params)
    super()._present_and_positive(p, ['medium_c', 'a', 'f'])

    if not p['shape'].isin(self.shapes).all():
        raise ValueError('The shape parameter must be one of: ' + ', '.join(self.shapes))

    if not p['boundary_type'].isin(self.boundary_types).all():
        raise ValueError('The boundary_type parameter must be one of: ' +
                         ', '.join(self.boundary_types))

echosms.KAModel()

Bases: ScatterModelBase

Kirchhoff approximation (KA) scattering model.

This class calculates acoustic scatter from arbitrary surfaces.

Source code in src/echosms/kamodel.py 
17
18
19
20
21
22
23
24
25
def __init__(self):
    super().__init__()
    self.long_name = 'Kirchhoff approximation'
    self.short_name = 'ka'
    self.analytical_type = 'approximate'
    self.boundary_types = [bt.pressure_release]
    self.shapes = ['closed surfaces']
    self.max_ka = 20  # [1]
    self.no_expand_parameters = ['mesh']

calculate_ts(data, expand=False, inplace=False, multiprocess=False, progress=False)

Calculate the target strength (TS) for many parameters.

PARAMETER DESCRIPTION
data

Requirements for the different input data types are:

  • DataFrame: column names must match the function parameter names in calculate_ts_single(). One TS value will be calculated for each row in the DataFrame.
  • DataArray: dimension names must match the function parameter names in calculate_ts_single(). TS values will be calculated for all combinations of the coordinate variables.
  • dict: keys must match the function parameters in calculate_ts_single(). TS values will be calculated for all combinations of the dict values.

TYPE: dict | DataFrame | DataArray

multiprocess

Split the ts calculation across CPU cores. Multiprocessing is currently provided by mapply. For more sophisticated uses it may be preferred to use a multiprocessing package of your choice directly on the calculate_ts_single() method. See the code in this method (calculate_ts()) for an example.

TYPE: bool DEFAULT: False

expand

Only applicable if data is a dict. If True, will use as_dataframe() to expand the dict into a DataFrame with one column per dict key and return that, adding a column named ts for the results.

TYPE: bool DEFAULT: False

inplace

Only applicable if data is a DataFrame. If True, the results will be added to the input DataFrame in a column named ts. If a ts column already exists, it is overwritten.

TYPE: bool DEFAULT: False

progress

If True, will produce a progress bar while running models

TYPE: bool DEFAULT: False

RETURNS DESCRIPTION
None | list[float] | Series | DataFrame

The return type and value are determined by the type of the input variable (data) and the expand and inplace parameters:

  • dict input and expand=False returns a list of floats.
  • dict input and expand=True returns a DataFrame.
  • DataFrame input and inplace=False returns a Series.
  • DataFrame input and inplace=True modifies data and returns None.
  • DataArray input always modifies data and returns None.
Source code in src/echosms/scattermodelbase.py 
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
def calculate_ts(self, data: dict | pd.DataFrame | xr.DataArray, expand: bool=False,
                 inplace: bool=False, multiprocess: bool=False,
                 progress: bool=False) -> None | list[float] | pd.Series | pd.DataFrame:
    """Calculate the target strength (TS) for many parameters.

    Parameters
    ----------
    data :
        Requirements for the different input data types are:

        - **DataFrame**: column names must match the function parameter names in
          calculate_ts_single(). One TS value will be calculated for each row in the DataFrame.
        - **DataArray**: dimension names must match the function parameter names in
          calculate_ts_single(). TS values will be calculated for all combinations of the
          coordinate variables.
        - **dict**: keys must match the function parameters in calculate_ts_single().
          TS values will be calculated for all combinations of the dict values.

    multiprocess :
        Split the ts calculation across CPU cores. Multiprocessing is currently provided by
        [mapply](https://github.com/ddelange/mapply). For more
        sophisticated uses it may be preferred to use a multiprocessing package of your choice
        directly on the `calculate_ts_single()` method. See the code in this method
        (`calculate_ts()`) for an example.

    expand :
        Only applicable if `data` is a dict. If `True`, will use
        [`as_dataframe()`][echosms.utils.as_dataframe]
        to expand the dict into a DataFrame with one column per dict key
        and return that, adding a column named `ts` for the results.

    inplace :
        Only applicable if `data` is a DataFrame. If `True`, the results
        will be added to the input DataFrame in a column named `ts`. If a `ts` column
        already exists, it is overwritten.

    progress :
        If `True`, will produce a progress bar while running models

    Returns
    -------
    :
        The return type and value are determined by the type of the input variable (`data`) and
        the `expand` and `inplace` parameters:

        - dict input and `expand=False` returns a list of floats.
        - dict input and `expand=True` returns a DataFrame.
        - DataFrame input and `inplace=False` returns a Series.
        - DataFrame input and `inplace=True` modifies `data` and returns `None`.
        - DataArray input always modifies `data` and returns `None`.

    """
    match data:
        case dict():
            data_df = as_dataframe(data, self.no_expand_parameters)
        case pd.DataFrame():
            data_df = data
        case xr.DataArray():
            data_df = data.to_dataframe().reset_index()
            data_df.attrs = data.attrs
        case _:
            raise ValueError(f'Data type of {type(data)} is not supported'
                             ' (only dictionaries, Pandas DataFrames and'
                             ' Xarray DataArrays are).')

    self.validate_parameters(data_df)

    # Get the non-expandable model parameters
    p = data_df.attrs['parameters'] if 'parameters' in data_df.attrs else {}

    # Note: the args argument in the apply call below requires a tuple. data_df.attrs is a
    # dict and the default behaviour is to make a tuple using the dict keys. The trailing comma
    # and parenthesis instead causes the tuple to have one entry of the dict.

    if multiprocess:
        from mapply.mapply import mapply
        ts = mapply(data_df, self.__ts_helper, args=(p,), axis=1, progressbar=progress)
    else:  # this uses just one CPU
        if progress:
            tqdm.pandas(desc=self.short_name, unit=' models',
                        bar_format='{l_bar}{bar} [{n_fmt}/{total_fmt}; {rate_noinv_fmt}]')
            ts = data_df.progress_apply(self.__ts_helper, args=(p,), axis=1)
        else:
            ts = data_df.apply(self.__ts_helper, args=(p,), axis=1)

    match data:
        case dict() if expand:
            data_df['ts'] = ts
            return data_df
        case dict():
            return ts.to_list()
        case pd.DataFrame() if inplace:
            data_df['ts'] = ts
            return None
        case pd.DataFrame():
            return ts.rename('ts', inplace=True)
        case xr.DataArray():
            data.values = ts.to_numpy().reshape(data.shape)
            return None
        case _:
            raise AssertionError('This code should never be reached - unsupported input data '
                                 f'type of {type(data)}.')

calculate_ts_single(medium_c, theta, phi, f, mesh, boundary_type, validate_parameters=True, **kwargs)

Calculate the scatter using the ka model for one set of parameters.

PARAMETER DESCRIPTION
medium_c

Sound speed in the fluid medium surrounding the target [m/s].

TYPE: float

theta

Pitch angle to calculate the scattering as per the echoSMs coordinate system [°].

TYPE: float

phi

Roll angle to calculate the scattering as per the echoSMs coordinate system [°].

TYPE: float

f

Frequency to calculate the scattering at [Hz].

TYPE: float

mesh

The triangular mesh that defines the scattering surface. This parameter must provide attributes with names of:

  • triangles_center (the position of the centre of each triangular face [m]),
  • face_normals (the outward-pointing unit normals for each triangular face),
  • area_faces (the area of each triangular face [m²]).

A suitable library for creating and manipulating triangular meshes is trimesh. Trimesh will accept the usual nodes/facets definition of a mesh and calculate the above attributes automatically.

TYPE: Any

boundary_type

The boundary type. Supported types are given in the boundary_types class variable.

TYPE: boundary_type

validate_parameters

Whether to validate the model parameters.

TYPE: bool DEFAULT: True

RETURNS DESCRIPTION
float

The target strength (re 1 m²) of the target [dB].
Notes

The class implements the code in Foote (1985).

References

Foote, K. G. (1985). Rather-high-frequency sound scattering of swimbladdered fish. The Journal of the Acoustical Society of America, 78(2), 688–700. https://doi.org/10.1121/1.392438

Source code in src/echosms/kamodel.py 
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
def calculate_ts_single(self, medium_c: float, theta: float, phi: float, f: float,
                        mesh: Any, boundary_type: bt,
                        validate_parameters: bool=True, **kwargs) -> float:
    """
    Calculate the scatter using the ka model for one set of parameters.

    Parameters
    ----------
    medium_c :
        Sound speed in the fluid medium surrounding the target [m/s].
    theta :
        Pitch angle to calculate the scattering as per the echoSMs
        [coordinate system](conventions.md#coordinate-systems) [°].
    phi :
        Roll angle to calculate the scattering as per the echoSMs
        [coordinate system](conventions.md#coordinate-systems) [°].
    f :
        Frequency to calculate the scattering at [Hz].
    mesh :
        The triangular mesh that defines the scattering surface. This parameter must provide
        attributes with names of:

        - `triangles_center` (the position of the centre of each triangular face [m]),
        - `face_normals` (the outward-pointing unit normals for each triangular face),
        - `area_faces` (the area of each triangular face [m²]).

        A suitable library for creating and manipulating triangular meshes
        is [trimesh](https://trimesh.org). Trimesh will accept the usual nodes/facets
        definition of a mesh and calculate the above attributes automatically.
    boundary_type :
        The boundary type. Supported types are given in the `boundary_types` class variable.
    validate_parameters :
        Whether to validate the model parameters.

    Returns
    -------
    :
        The target strength (re 1 m²) of the target [dB].

    Notes
    -----
    The class implements the code in Foote (1985).

    References
    ----------
    Foote, K. G. (1985). Rather-high-frequency sound scattering of swimbladdered fish.
    The Journal of the Acoustical Society of America, 78(2), 688–700.
    <https://doi.org/10.1121/1.392438>

    """
    if validate_parameters:
        self.validate_parameters(locals())

    if boundary_type not in self.boundary_types:
        raise ValueError(f'The {self.long_name} model does not support '
                         f'a model type of "{boundary_type}".')

    # This model keeps the organism fixed and varies the incident wave vector. So need
    # to convert the theta and phi echoSMs coordinate sytem Tait-Bryan angles
    # into an (x,y,z) vector.

    # Acoustic wave incident vector and its' norm
    rot = R.from_euler('ZYX', (0, theta-90, -phi), degrees=True)
    k_norm = rot.as_matrix() @ np.array([[0, 0, 1]]).T
    k = k_norm * wavenumber(medium_c, f)

    r = mesh.triangles_center  # position vector of each surface element
    dS = mesh.area_faces.reshape((-1, 1))  # [m^2]

    kn_nn = mesh.face_normals @ k_norm

    fbs = 1./wavelength(medium_c, f)\
        * np.sum(np.exp(2j*r @ k) * np.heaviside(kn_nn, 0.5) * kn_nn * dS)

    return 10*log10(abs(fbs)**2)  # ts

validate_parameters(params)

Validate the model parameters.

See here for calling details.

Source code in src/echosms/kamodel.py 
27
28
29
30
31
32
33
34
def validate_parameters(self, params):
    """Validate the model parameters.

    See [here][echosms.ScatterModelBase.validate_parameters] for calling details.
    """
    p = as_dict(params)
    super()._present_and_in(p, ['boundary_type'], self.boundary_types)
    super()._present_and_positive(p, ['medium_c', 'f'])

KRM model & utilities

echosms.KRMModel()

Bases: ScatterModelBase

Kirchhoff ray mode (KRM) scattering model.

Source code in src/echosms/krmmodel.py 
30
31
32
33
34
35
36
37
38
def __init__(self):
    super().__init__()
    self.long_name = 'Kirchhoff ray mode'
    self.short_name = 'krm'
    self.analytical_type = 'approximate'
    self.boundary_types = [bt.fluid_filled]
    self.shapes = ['closed surfaces']
    self.max_ka = 20  # [1]
    self.no_expand_parameters = ['bodies']

calculate_ts(data, expand=False, inplace=False, multiprocess=False, progress=False)

Calculate the target strength (TS) for many parameters.

PARAMETER DESCRIPTION
data

Requirements for the different input data types are:

  • DataFrame: column names must match the function parameter names in calculate_ts_single(). One TS value will be calculated for each row in the DataFrame.
  • DataArray: dimension names must match the function parameter names in calculate_ts_single(). TS values will be calculated for all combinations of the coordinate variables.
  • dict: keys must match the function parameters in calculate_ts_single(). TS values will be calculated for all combinations of the dict values.

TYPE: dict | DataFrame | DataArray

multiprocess

Split the ts calculation across CPU cores. Multiprocessing is currently provided by mapply. For more sophisticated uses it may be preferred to use a multiprocessing package of your choice directly on the calculate_ts_single() method. See the code in this method (calculate_ts()) for an example.

TYPE: bool DEFAULT: False

expand

Only applicable if data is a dict. If True, will use as_dataframe() to expand the dict into a DataFrame with one column per dict key and return that, adding a column named ts for the results.

TYPE: bool DEFAULT: False

inplace

Only applicable if data is a DataFrame. If True, the results will be added to the input DataFrame in a column named ts. If a ts column already exists, it is overwritten.

TYPE: bool DEFAULT: False

progress

If True, will produce a progress bar while running models

TYPE: bool DEFAULT: False

RETURNS DESCRIPTION
None | list[float] | Series | DataFrame

The return type and value are determined by the type of the input variable (data) and the expand and inplace parameters:

  • dict input and expand=False returns a list of floats.
  • dict input and expand=True returns a DataFrame.
  • DataFrame input and inplace=False returns a Series.
  • DataFrame input and inplace=True modifies data and returns None.
  • DataArray input always modifies data and returns None.
Source code in src/echosms/scattermodelbase.py 
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
def calculate_ts(self, data: dict | pd.DataFrame | xr.DataArray, expand: bool=False,
                 inplace: bool=False, multiprocess: bool=False,
                 progress: bool=False) -> None | list[float] | pd.Series | pd.DataFrame:
    """Calculate the target strength (TS) for many parameters.

    Parameters
    ----------
    data :
        Requirements for the different input data types are:

        - **DataFrame**: column names must match the function parameter names in
          calculate_ts_single(). One TS value will be calculated for each row in the DataFrame.
        - **DataArray**: dimension names must match the function parameter names in
          calculate_ts_single(). TS values will be calculated for all combinations of the
          coordinate variables.
        - **dict**: keys must match the function parameters in calculate_ts_single().
          TS values will be calculated for all combinations of the dict values.

    multiprocess :
        Split the ts calculation across CPU cores. Multiprocessing is currently provided by
        [mapply](https://github.com/ddelange/mapply). For more
        sophisticated uses it may be preferred to use a multiprocessing package of your choice
        directly on the `calculate_ts_single()` method. See the code in this method
        (`calculate_ts()`) for an example.

    expand :
        Only applicable if `data` is a dict. If `True`, will use
        [`as_dataframe()`][echosms.utils.as_dataframe]
        to expand the dict into a DataFrame with one column per dict key
        and return that, adding a column named `ts` for the results.

    inplace :
        Only applicable if `data` is a DataFrame. If `True`, the results
        will be added to the input DataFrame in a column named `ts`. If a `ts` column
        already exists, it is overwritten.

    progress :
        If `True`, will produce a progress bar while running models

    Returns
    -------
    :
        The return type and value are determined by the type of the input variable (`data`) and
        the `expand` and `inplace` parameters:

        - dict input and `expand=False` returns a list of floats.
        - dict input and `expand=True` returns a DataFrame.
        - DataFrame input and `inplace=False` returns a Series.
        - DataFrame input and `inplace=True` modifies `data` and returns `None`.
        - DataArray input always modifies `data` and returns `None`.

    """
    match data:
        case dict():
            data_df = as_dataframe(data, self.no_expand_parameters)
        case pd.DataFrame():
            data_df = data
        case xr.DataArray():
            data_df = data.to_dataframe().reset_index()
            data_df.attrs = data.attrs
        case _:
            raise ValueError(f'Data type of {type(data)} is not supported'
                             ' (only dictionaries, Pandas DataFrames and'
                             ' Xarray DataArrays are).')

    self.validate_parameters(data_df)

    # Get the non-expandable model parameters
    p = data_df.attrs['parameters'] if 'parameters' in data_df.attrs else {}

    # Note: the args argument in the apply call below requires a tuple. data_df.attrs is a
    # dict and the default behaviour is to make a tuple using the dict keys. The trailing comma
    # and parenthesis instead causes the tuple to have one entry of the dict.

    if multiprocess:
        from mapply.mapply import mapply
        ts = mapply(data_df, self.__ts_helper, args=(p,), axis=1, progressbar=progress)
    else:  # this uses just one CPU
        if progress:
            tqdm.pandas(desc=self.short_name, unit=' models',
                        bar_format='{l_bar}{bar} [{n_fmt}/{total_fmt}; {rate_noinv_fmt}]')
            ts = data_df.progress_apply(self.__ts_helper, args=(p,), axis=1)
        else:
            ts = data_df.apply(self.__ts_helper, args=(p,), axis=1)

    match data:
        case dict() if expand:
            data_df['ts'] = ts
            return data_df
        case dict():
            return ts.to_list()
        case pd.DataFrame() if inplace:
            data_df['ts'] = ts
            return None
        case pd.DataFrame():
            return ts.rename('ts', inplace=True)
        case xr.DataArray():
            data.values = ts.to_numpy().reshape(data.shape)
            return None
        case _:
            raise AssertionError('This code should never be reached - unsupported input data '
                                 f'type of {type(data)}.')

calculate_ts_single(medium_c, medium_rho, theta, f, organism, high_ka_medium='body', low_ka_medium='body', validate_parameters=True, **kwargs)

Calculate the scatter using the Kirchhoff ray mode model for one set of parameters.

Warning

The mode solution (low ka) part of this model has not yet been verified to give correct results.

PARAMETER DESCRIPTION
medium_c

Sound speed in the fluid medium surrounding the target [m/s].

TYPE: float

medium_rho

Density in the fluid medium surrounding the target [kg/m³]

TYPE: float

theta

Pitch angle to calculate the scattering at, as per the echoSMs coordinate system [°].

TYPE: float

f

Frequency to calculate the scattering at [Hz].

TYPE: float

organism

The shapes that make up the model. This is typically a shape for the body and zero or more enclosed shapes that repesent internal parts of the organism.

TYPE: KRMorganism

high_ka_medium

If set to body the sound speed and density of the organism body is used for the fluid surrounding any inclusions. If set to anything else (e.g., water) the sound speed and density given by medium_c and medium_rho are used. This parameter applies to the Kirchhoff approximation part of the model (i.e., high ka) and corresponds to the use (or not) of the approximation given in Clay & Horne (1994) on the line immediately below Eqn (13): k_b ≈ k at low contrast.

TYPE: str DEFAULT: 'body'

low_ka_medium

If set to body the sound speed and density of the organism body is used for the fluid surrounding any inclusions. If set to anything else (e.g., water) the sound speed and density given by medium_c and medium_rho are used. This parameter applies to the mode solution part of the model (i.e., low ka) and corresponds to the use (or not) of the approximation given in Clay & Horne (1994) on the line immediately below Eqn (13): k_b ≈ k at low contrast.

TYPE: str DEFAULT: 'body'

validate_parameters

Whether to validate the model parameters.

TYPE: bool DEFAULT: True

RETURNS DESCRIPTION
float

The target strength (re 1 m²) of the target [dB].
Notes

The class implements the code in Clay & Horne (1994) and when ka < 0.15 uses Clay (1992).

The high_ka_medium and low_ka_medium parameters allow the user to select which medium surrounds the inclusions (e.g., the swimbladder) - the fish body or the water surrounding the fish body. The equations in Clay & Horne (1994) used the body but included a sentence saying that the water could be used at low contrast (between the water and the body). A later paper (Horne & Jech, 1999) used water for low ka (the mode solution) and the body for higher ka (the Kirchhoff approximation). Some open-source KRM model codes always use the water.

References

Clay, C. S. (1992). Composite ray-mode approximations for backscattered sound from gas-filled cylinders and swimbladders. The Journal of the Acoustical Society of America, 92(4), 2173–2180. https://doi.org/10.1121/1.405211

Clay, C. S., & Horne, J. K. (1994). Acoustic models of fish: The Atlantic cod (Gadus morhua). The Journal of the Acoustical Society of America, 96(3), 1661–1668. https://doi.org/10.1121/1.410245

Horne, J. K., & J. M. Jech. (1999). Multi–frequency estimates of fish abundance: constraints of rather high frequencies. ICES Journal of Marine Science, 56 (2), 184–199. https://doi.org/10.1006/jmsc.1998.0432

Source code in src/echosms/krmmodel.py 
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
def calculate_ts_single(self, medium_c: float, medium_rho: float, theta: float,
                        f: float, organism: KRMorganism, high_ka_medium: str='body',
                        low_ka_medium: str='body',
                        validate_parameters: bool=True, **kwargs) -> float:
    """
    Calculate the scatter using the Kirchhoff ray mode model for one set of parameters.

    Warning
    --------
    The mode solution (low _ka_) part of this model has not yet been verified to give
    correct results.

    Parameters
    ----------
    medium_c :
        Sound speed in the fluid medium surrounding the target [m/s].
    medium_rho :
        Density in the fluid medium surrounding the target [kg/m³]
    theta :
        Pitch angle to calculate the scattering at, as per the echoSMs
        [coordinate system](conventions.md#coordinate-systems) [°].
    f :
        Frequency to calculate the scattering at [Hz].
    organism :
        The shapes that make up the model. This is typically a shape for the body and zero or
        more enclosed shapes that repesent internal parts of the organism.
    high_ka_medium :
        If set to `body` the sound speed and density of the organism body is used for
        the fluid surrounding any inclusions. If set to anything else (e.g., `water`)
        the sound speed and density given by `medium_c` and `medium_rho` are used.
        This parameter applies to the Kirchhoff approximation part of
        the model (i.e., high _ka_) and corresponds to the use (or not) of the
        approximation given in Clay & Horne (1994) on the line immediately below Eqn (13):
        _k_b ≈ k at low contrast_.
    low_ka_medium :
        If set to `body` the sound speed and density of the organism body is used for
        the fluid surrounding any inclusions. If set to anything else (e.g., `water`)
        the sound speed and density given by `medium_c` and `medium_rho` are used.
        This parameter applies to the mode solution part of the model (i.e., low _ka_)
        and corresponds to the use (or not) of the approximation given in Clay & Horne (1994)
        on the line immediately below Eqn (13): _k_b ≈ k at low contrast_.
    validate_parameters :
        Whether to validate the model parameters.

    Returns
    -------
    :
        The target strength (re 1 m²) of the target [dB].

    Notes
    -----
    The class implements the code in Clay & Horne (1994) and when _ka_ < 0.15 uses Clay (1992).

    The `high_ka_medium` and `low_ka_medium` parameters allow the user to select which
    medium surrounds the inclusions (e.g., the swimbladder) - the fish body or
    the water surrounding the fish body. The equations in Clay & Horne (1994) used the body
    but included a sentence saying that the water could be used at low contrast
    (between the water and the body). A later paper
    (Horne & Jech, 1999) used water for low _ka_ (the mode solution) and the body for
    higher _ka_ (the Kirchhoff approximation). Some open-source KRM model codes always
    use the water.

    References
    ----------
    Clay, C. S. (1992). Composite ray-mode approximations for backscattered sound from
    gas-filled cylinders and swimbladders. The Journal of the Acoustical Society of
    America, 92(4), 2173–2180.
    <https://doi.org/10.1121/1.405211>

    Clay, C. S., & Horne, J. K. (1994). Acoustic models of fish: The Atlantic cod
    (_Gadus morhua_). The Journal of the Acoustical Society of America, 96(3), 1661–1668.
    <https://doi.org/10.1121/1.410245>

    Horne, J. K., & J. M. Jech. (1999). Multi–frequency estimates of fish abundance:
    constraints of rather high frequencies. ICES Journal of Marine Science, 56 (2), 184–199.
    <https://doi.org/10.1006/jmsc.1998.0432>
    """
    if validate_parameters:
        self.validate_parameters(locals())

    theta = radians(theta)

    body = organism.body

    k = wavenumber(medium_c, f)
    k_b = wavenumber(body.c, f)

    # Reflection coefficient between water and body
    R_wb = (body.rho*body.c - medium_rho*medium_c)\
        / (body.rho*body.c + medium_rho*medium_c)
    TwbTbw = 1-R_wb**2  # Eqn (15)

    sl = []  # scattering lengths for inclusions
    for incl in organism.inclusions:
        # Reflection coefficient between body and inclusion
        # The paper gives R_bc in terms of g & h, but it can also be done in the
        # same manner as R_wb above.
        gp = incl.rho / body.rho  # p is 'prime' to fit with paper notation
        hp = incl.c / body.c

        R_bc = (gp*hp-1) / (gp*hp+1)  # Eqn (9)

        # Equivalent radius of inclusion (as per Part A of paper)
        a_e = sqrt(incl.volume() / (pi * incl.length()))

        # Choose which modelling approach to use
        if k*a_e < 0.15:  # Do the mode solution for the inclusion
            if low_ka_medium != 'body':
                gp = incl.rho / medium_rho
                hp = incl.c / medium_c
            sl.append(self._mode_solution(1/gp, 1/hp, k, a_e, incl.length(), theta))
        elif incl.boundary == bt.pressure_release:
            kk = k_b if high_ka_medium == 'body' else k
            sl.append(self._soft_KA(incl, k, kk, R_bc, TwbTbw, theta))
        elif incl.boundary == bt.fluid_filled:
            kk = k_b if high_ka_medium == 'body' else k
            sl.append(self._fluid_KA(incl, k, kk, R_bc, TwbTbw, theta))
        else:
            raise ValueError(f'Unsupported boundary of "{incl.boundary}" for KRM inclusion')

    # Do the Kirchhoff-ray approximation for the body. This is always done as a fluid.
    body_sl = self._fluid_KA(body, k, k_b, R_wb, TwbTbw, theta)

    return 20*log10(abs(body_sl + sum(sl)))

validate_parameters(params)

Validate the model parameters.

See here for calling details.

Source code in src/echosms/krmmodel.py 
40
41
42
43
44
45
46
47
48
49
def validate_parameters(self, params):
    """Validate the model parameters.

    See [here][echosms.ScatterModelBase.validate_parameters] for calling details.
    """
    p = as_dict(params)
    super()._present_and_positive(p, ['medium_c', 'f'])

    if np.any(np.atleast_1d(p['theta']) < 65) or np.any(np.atleast_1d(p['theta']) > 115):
        raise KeyError('Incidence angle(s) (theta) are outside 65 to 115°')

echosms.KRMdata()

Example datasets for the KRM model.

Source code in src/echosms/krmdata.py 
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
def __init__(self):
    # Load in the NOAA KRM shapes data
    self.file = Path(__file__).parent/Path('resources')/Path('KRM_shapes.toml')
    with open(self.file, 'rb') as f:
        try:
            shapes = tomllib.load(f)
        except tomllib.TOMLDecodeError as e:
            raise SyntaxError(f'Error while parsing file "{self.defs_filename.name}"') from e

    # Put the shapes into a dict of KRMorganism(). Use some default values for sound speed and
    # density
    self.krm_models = {}
    for s in shapes['shape']:
        # These KRM data have the head pointing in the -ve x direction,
        # opposite to the echoSMs coordinate convetion, so fix the
        # x-coordinates here when ingesting the data. And set the posterior end
        # of the organism to have x=0
        m = max(s['x_b'])
        body = KRMshape(bt.fluid_filled, -np.array(s['x_b']), np.array(s['w_b']),
                        np.array(s['z_bU']), np.array(s['z_bL']),
                        s['body_c'], s['body_rho'])
        swimbladder = KRMshape(bt.pressure_release, -np.array(s['x_sb']), np.array(s['w_sb']),
                               np.array(s['z_sbU']), np.array(s['z_sbL']),
                               s['swimbladder_c'], s['swimbladder_rho'])
        self.krm_models[s['name']] = KRMorganism(s['name'], s['source'],
                                                 body, [swimbladder],
                                                 s['aphiaid'], s['length'],
                                                 s['vernacular'])

as_dict()

KRM model shapes as a dict.

RETURNS DESCRIPTION
dict

All the KRM model shapes. The dataset name is the dict key and the value is an instance of KRMorganism.
Source code in src/echosms/krmdata.py 
156
157
158
159
160
161
162
163
164
165
166
def as_dict(self) -> dict:
    """KRM model shapes as a dict.

    Returns
    -------
    :
        All the KRM model shapes. The dataset name is the dict key and the value is an instance
        of `KRMorganism`.

    """
    return self.krm_models

model(name)

KRM model shape with requested name.

PARAMETER DESCRIPTION
name

The name of a KRM model shape.

TYPE: str

RETURNS DESCRIPTION
KRMorganism

An instance of KRMorganism or None if there is no model with name.
Source code in src/echosms/krmdata.py 
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
def model(self, name: str) -> KRMorganism:
    """KRM model shape with requested name.

    Parameters
    ----------
    name :
        The name of a KRM model shape.

    Returns
    -------
    :
        An instance of `KRMorganism` or None if there is no model with `name`.

    """
    try:
        return self.krm_models[name]
    except KeyError:
        return None

names()

Available KRM model names.

Source code in src/echosms/krmdata.py 
152
153
154
def names(self):
    """Available KRM model names."""
    return [*self.krm_models]

ts(name) staticmethod

KRM model TS from model name.

PARAMETER DESCRIPTION
name

The name of a KRM model shape.

TYPE: str

RETURNS DESCRIPTION
ndarray

The TS (re 1 m²) for some default model parameters [dB] or None if no TS data are available.
Source code in src/echosms/krmdata.py 
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
@staticmethod
def ts(name: str) -> np.ndarray:
    """KRM model TS from model `name`.

    Parameters
    ----------
    name :
        The name of a KRM model shape.

    Returns
    -------
    :
        The TS (re 1 m²) for some default model parameters [dB] or None if no TS data
        are available.

    """
    # Sometimes there will be TS results for the model (available for testing of the
    # model), so load them in if present.
    tsfile = Path(__file__).parent/Path('resources')/Path('NOAA_KRM_ts_' + name + '.csv')

    if tsfile.exists():
        return pd.read_csv(tsfile)

    return None

echosms.KRMorganism(name, source, body, inclusions, aphiaid=1, length=0.0, vernacular_name='') dataclass

KRM body and inclusion shape(s).

ATTRIBUTE DESCRIPTION
name

A name for the organism.

TYPE: str

source

A link to or description of the source of the organism data.

TYPE: str

body

The shape that represents the organism's body.

TYPE: KRMshape

inclusions

The shapes that are internal to the organism (e.g., swimbladder, backbone, etc)

TYPE: List[KRMshape]

aphiaid

The aphiaID of the organism

TYPE: int

length

The length of the organism (m)

TYPE: float

vernacular_name

A vernacular name of the organism

TYPE: str

plot()

Plot of organism shape.

Source code in src/echosms/krmdata.py 
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
def plot(self):
    """Plot of organism shape."""
    import matplotlib.pyplot as plt

    plt.plot(self.body.x*1e3, self.body.z_U*1e3, self.body.x*1e3, self.body.z_L*1e3, c='black')
    for i in [0, -1]:  # close the ends of the shape
        plt.plot([self.body.x[i]*1e3]*2, [self.body.z_U[i]*1e3, self.body.z_L[i]*1e3],
                 c='black')

    for s in self.inclusions:
        c = 'C0' if s.boundary == bt.fluid_filled else 'C1'
        plt.plot(s.x*1e3, s.z_U*1e3, s.x*1e3, s.z_L*1e3, c=c)
        for i in [0, -1]:  # close the ends of the shape
            plt.plot([s.x[i]*1e3]*2, [s.z_U[i]*1e3, s.z_L[i]*1e3], c=c)

    plt.gca().set_aspect('equal')
    plt.gca().xaxis.set_inverted(True)
    plt.title(self.name)
    plt.show()

echosms.KRMshape(boundary, x, w, z_U, z_L, c, rho) dataclass

KRM shape and property class.

ATTRIBUTE DESCRIPTION
boundary

The shape bounday condition - either pressure_release or fluid_filled.

TYPE: boundary_type

x

The x-axis coordinates [m].

TYPE: ndarray

w

Width of the shape [m].

TYPE: ndarray

z_U

Distance from the axis to the upper surface of the shape [m].

TYPE: ndarray

z_L

Distance from the axis to the lower surface of the shape [m].

TYPE: ndarray

c

Sound speed in the shape [m/s].

TYPE: float

rho

Density of the shape material [kg/m³].

TYPE: float

length()

Length of the shape.

RETURNS DESCRIPTION
float

The length of the shape [m].
Source code in src/echosms/krmdata.py 
58
59
60
61
62
63
64
65
66
def length(self) -> float:
    """Length of the shape.

    Returns
    -------
    :
        The length of the shape [m].
    """
    return self.x[-1] - self.x[0]

volume()

Volume of the shape.

RETURNS DESCRIPTION
float

The volume of the shape [m³].
Source code in src/echosms/krmdata.py 
46
47
48
49
50
51
52
53
54
55
56
def volume(self) -> float:
    """Volume of the shape.

    Returns
    -------
    :
        The volume of the shape [m³].
    """
    thickness = np.diff(self.x)
    thickness = np.append(thickness, thickness[1])
    return np.sum(np.pi * (self.z_U - self.z_L) * self.w * thickness)

echosms.MSSModel()

Bases: ScatterModelBase

Modal series solution (MSS) scattering model.

This class calculates acoustic scatter from spheres and shells with various boundary conditions, as listed in the boundary_types class attribute.

Source code in src/echosms/mssmodel.py 
17
18
19
20
21
22
23
24
25
26
def __init__(self):
    super().__init__()
    self.long_name = 'modal series solution'
    self.short_name = 'mss'
    self.analytical_type = 'exact'
    self.boundary_types = [bt.fixed_rigid, bt.pressure_release, bt.fluid_filled,
                           bt.fluid_shell_fluid_interior,
                           bt.fluid_shell_pressure_release_interior]
    self.shapes = ['sphere']
    self.max_ka = 20  # [1]

calculate_ts(data, expand=False, inplace=False, multiprocess=False, progress=False)

Calculate the target strength (TS) for many parameters.

PARAMETER DESCRIPTION
data

Requirements for the different input data types are:

  • DataFrame: column names must match the function parameter names in calculate_ts_single(). One TS value will be calculated for each row in the DataFrame.
  • DataArray: dimension names must match the function parameter names in calculate_ts_single(). TS values will be calculated for all combinations of the coordinate variables.
  • dict: keys must match the function parameters in calculate_ts_single(). TS values will be calculated for all combinations of the dict values.

TYPE: dict | DataFrame | DataArray

multiprocess

Split the ts calculation across CPU cores. Multiprocessing is currently provided by mapply. For more sophisticated uses it may be preferred to use a multiprocessing package of your choice directly on the calculate_ts_single() method. See the code in this method (calculate_ts()) for an example.

TYPE: bool DEFAULT: False

expand

Only applicable if data is a dict. If True, will use as_dataframe() to expand the dict into a DataFrame with one column per dict key and return that, adding a column named ts for the results.

TYPE: bool DEFAULT: False

inplace

Only applicable if data is a DataFrame. If True, the results will be added to the input DataFrame in a column named ts. If a ts column already exists, it is overwritten.

TYPE: bool DEFAULT: False

progress

If True, will produce a progress bar while running models

TYPE: bool DEFAULT: False

RETURNS DESCRIPTION
None | list[float] | Series | DataFrame

The return type and value are determined by the type of the input variable (data) and the expand and inplace parameters:

  • dict input and expand=False returns a list of floats.
  • dict input and expand=True returns a DataFrame.
  • DataFrame input and inplace=False returns a Series.
  • DataFrame input and inplace=True modifies data and returns None.
  • DataArray input always modifies data and returns None.
Source code in src/echosms/scattermodelbase.py 
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
def calculate_ts(self, data: dict | pd.DataFrame | xr.DataArray, expand: bool=False,
                 inplace: bool=False, multiprocess: bool=False,
                 progress: bool=False) -> None | list[float] | pd.Series | pd.DataFrame:
    """Calculate the target strength (TS) for many parameters.

    Parameters
    ----------
    data :
        Requirements for the different input data types are:

        - **DataFrame**: column names must match the function parameter names in
          calculate_ts_single(). One TS value will be calculated for each row in the DataFrame.
        - **DataArray**: dimension names must match the function parameter names in
          calculate_ts_single(). TS values will be calculated for all combinations of the
          coordinate variables.
        - **dict**: keys must match the function parameters in calculate_ts_single().
          TS values will be calculated for all combinations of the dict values.

    multiprocess :
        Split the ts calculation across CPU cores. Multiprocessing is currently provided by
        [mapply](https://github.com/ddelange/mapply). For more
        sophisticated uses it may be preferred to use a multiprocessing package of your choice
        directly on the `calculate_ts_single()` method. See the code in this method
        (`calculate_ts()`) for an example.

    expand :
        Only applicable if `data` is a dict. If `True`, will use
        [`as_dataframe()`][echosms.utils.as_dataframe]
        to expand the dict into a DataFrame with one column per dict key
        and return that, adding a column named `ts` for the results.

    inplace :
        Only applicable if `data` is a DataFrame. If `True`, the results
        will be added to the input DataFrame in a column named `ts`. If a `ts` column
        already exists, it is overwritten.

    progress :
        If `True`, will produce a progress bar while running models

    Returns
    -------
    :
        The return type and value are determined by the type of the input variable (`data`) and
        the `expand` and `inplace` parameters:

        - dict input and `expand=False` returns a list of floats.
        - dict input and `expand=True` returns a DataFrame.
        - DataFrame input and `inplace=False` returns a Series.
        - DataFrame input and `inplace=True` modifies `data` and returns `None`.
        - DataArray input always modifies `data` and returns `None`.

    """
    match data:
        case dict():
            data_df = as_dataframe(data, self.no_expand_parameters)
        case pd.DataFrame():
            data_df = data
        case xr.DataArray():
            data_df = data.to_dataframe().reset_index()
            data_df.attrs = data.attrs
        case _:
            raise ValueError(f'Data type of {type(data)} is not supported'
                             ' (only dictionaries, Pandas DataFrames and'
                             ' Xarray DataArrays are).')

    self.validate_parameters(data_df)

    # Get the non-expandable model parameters
    p = data_df.attrs['parameters'] if 'parameters' in data_df.attrs else {}

    # Note: the args argument in the apply call below requires a tuple. data_df.attrs is a
    # dict and the default behaviour is to make a tuple using the dict keys. The trailing comma
    # and parenthesis instead causes the tuple to have one entry of the dict.

    if multiprocess:
        from mapply.mapply import mapply
        ts = mapply(data_df, self.__ts_helper, args=(p,), axis=1, progressbar=progress)
    else:  # this uses just one CPU
        if progress:
            tqdm.pandas(desc=self.short_name, unit=' models',
                        bar_format='{l_bar}{bar} [{n_fmt}/{total_fmt}; {rate_noinv_fmt}]')
            ts = data_df.progress_apply(self.__ts_helper, args=(p,), axis=1)
        else:
            ts = data_df.apply(self.__ts_helper, args=(p,), axis=1)

    match data:
        case dict() if expand:
            data_df['ts'] = ts
            return data_df
        case dict():
            return ts.to_list()
        case pd.DataFrame() if inplace:
            data_df['ts'] = ts
            return None
        case pd.DataFrame():
            return ts.rename('ts', inplace=True)
        case xr.DataArray():
            data.values = ts.to_numpy().reshape(data.shape)
            return None
        case _:
            raise AssertionError('This code should never be reached - unsupported input data '
                                 f'type of {type(data)}.')

calculate_ts_single(medium_c, medium_rho, a, f, boundary_type, target_c=None, target_rho=None, shell_c=None, shell_rho=None, shell_thickness=None, validate_parameters=True, **kwargs)

Calculate the scatter using the mss model for one set of parameters.

PARAMETER DESCRIPTION
medium_c

Sound speed in the fluid medium surrounding the target [m/s].

TYPE: float

medium_rho

Density of the fluid medium surrounding the target [kg/m³].

TYPE: float

a

Radius of the spherical target [m].

TYPE: float

f

Frequency to calculate the scattering at [Hz].

TYPE: float

boundary_type

The boundary type. Supported types are given in the boundary_types class variable.

TYPE: boundary_type

target_c

Sound speed in the fluid inside the sphere [m/s]. Only required for boundary_type of fluid_filled.

TYPE: None | float DEFAULT: None

target_rho

Density of the fluid inside the sphere [kg/m³]. Only required for boundary_type of fluid_filled.

TYPE: None | float DEFAULT: None

shell_c

Sound speed in the spherical shell [m/s]. Only required for boundary_types that include a fluid shell.

TYPE: None | float DEFAULT: None

shell_rho

Density in the spherical shell [kg/m³]. Only required for boundary_types that include a fluid shell.

TYPE: None | float DEFAULT: None

shell_thickness

Thickness of the spherical shell [m]. This value is subtracted from a to give the radius of the interior sphere. Only required for boundary_types that include a fluid shell.

TYPE: None | float DEFAULT: None

validate_parameters

Whether to validate the model parameters.

TYPE: bool DEFAULT: True

RETURNS DESCRIPTION
float

The target strength (re 1 m²) of the target [dB].
Notes

The class implements the code in Section A.1 of Jech et al. (2015).

References

Jech, J.M., Horne, J.K., Chu, D., Demer, D.A., Francis, D.T.I., Gorska, N., Jones, B., Lavery, A.C., Stanton, T.K., Macaulay, G.J., Reeder, D.B., Sawada, K., 2015. Comparisons among ten models of acoustic backscattering used in aquatic ecosystem research. Journal of the Acoustical Society of America 138, 3742–3764. https://doi.org/10.1121/1.4937607

Source code in src/echosms/mssmodel.py 
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
def calculate_ts_single(self, medium_c: float, medium_rho: float, a: float, f: float,
                        boundary_type: bt,
                        target_c: None | float=None, target_rho: None | float=None,
                        shell_c: None | float=None, shell_rho: None | float=None,
                        shell_thickness: None | float=None,
                        validate_parameters: bool=True,
                        **kwargs) -> float:
    """
    Calculate the scatter using the mss model for one set of parameters.

    Parameters
    ----------
    medium_c :
        Sound speed in the fluid medium surrounding the target [m/s].
    medium_rho :
        Density of the fluid medium surrounding the target [kg/m³].
    a :
        Radius of the spherical target [m].
    f :
        Frequency to calculate the scattering at [Hz].
    boundary_type :
        The boundary type. Supported types are given in the `boundary_types` class variable.
    target_c :
        Sound speed in the fluid inside the sphere [m/s].
        Only required for `boundary_type` of ``fluid_filled``.
    target_rho :
        Density of the fluid inside the sphere [kg/m³].
        Only required for `boundary_type` of ``fluid_filled``.
    shell_c :
        Sound speed in the spherical shell [m/s].
        Only required for `boundary_type`s that include a fluid shell.
    shell_rho :
        Density in the spherical shell [kg/m³].
        Only required for `boundary_type`s that include a fluid shell.
    shell_thickness :
        Thickness of the spherical shell [m]. This value is subtracted from ``a`` to give
        the radius of the interior sphere.
        Only required for `boundary_type`s that include a fluid shell.
    validate_parameters :
        Whether to validate the model parameters.

    Returns
    -------
    :
        The target strength (re 1 m²) of the target [dB].

    Notes
    -----
    The class implements the code in Section A.1 of Jech et al. (2015).

    References
    ----------
    Jech, J.M., Horne, J.K., Chu, D., Demer, D.A., Francis, D.T.I., Gorska, N.,
    Jones, B., Lavery, A.C., Stanton, T.K., Macaulay, G.J., Reeder, D.B., Sawada, K., 2015.
    Comparisons among ten models of acoustic backscattering used in aquatic ecosystem
    research. Journal of the Acoustical Society of America 138, 3742–3764.
    <https://doi.org/10.1121/1.4937607>
    """
    if validate_parameters:
        self.validate_parameters(locals())

    k0 = wavenumber(medium_c, f)
    ka = k0*a
    n = np.arange(0, round(ka+20))

    match boundary_type:
        case bt.fixed_rigid:
            A = list(map(lambda x: -spherical_jn(x, ka, True) / h1(x, ka, True), n))
        case bt.pressure_release:
            A = list(map(lambda x: -spherical_jn(x, ka) / h1(x, ka), n))
        case bt.fluid_filled:
            k1a = wavenumber(target_c, f)*a
            gh = target_rho/medium_rho * target_c/medium_c

            def Cn_fr(n):
                return\
                    ((spherical_jn(n, k1a, True)*spherical_yn(n, ka))
                        / (spherical_jn(n, k1a)*spherical_jn(n, ka, True))
                        - gh*(spherical_yn(n, ka, True)/spherical_jn(n, ka, True)))\
                    / ((spherical_jn(n, k1a, True)*spherical_jn(n, ka))
                       / (spherical_jn(n, k1a)*spherical_jn(n, ka, True))-gh)

            A = -1/(1 + 1j*np.asarray(list(map(Cn_fr, n)), dtype=complex))
        case bt.fluid_shell_fluid_interior:
            b = a - shell_thickness

            g21 = shell_rho / medium_rho
            h21 = shell_c / medium_c
            g32 = target_rho / shell_rho
            h32 = target_c / shell_c

            k1a = wavenumber(medium_c, f) * a
            k2 = wavenumber(shell_c, f)
            k3b = wavenumber(target_c, f) * b

            def Cn_fsfi(n):
                (b1, b2, a11, a21, a12, a22, a32, a13, a23, a33) =\
                    MSSModel.__eqn9(n, k1a, g21, h21, k2*a, k2*b, k3b, h32, g32)
                return (b1*a22*a33 + a13*b2*a32 - a12*b2*a33 - b1*a23*a32)\
                    / (a11*a22*a33 + a13*a21*a32 - a12*a21*a33 - a11*a23*a32)

            A = list(map(Cn_fsfi, n))
        case bt.fluid_shell_pressure_release_interior:
            b = a - shell_thickness

            g21 = shell_rho / medium_rho
            h21 = shell_c / medium_c

            k1a = wavenumber(medium_c, f) * a
            k2 = wavenumber(shell_c, f)
            ksa = k2 * a  # ksa is used in the paper, but isn't that the same as k2a?

            def Cn_fspri(n):
                (b1, b2, d1, d2, a11, a21) = MSSModel.__eqn10(n, k1a, g21, h21, ksa, k2*a, k2*b)
                return (b1*d2-d1*b2) / (a11*d2-d1*a21)

            A = list(map(Cn_fspri, n))
        case _:
            raise ValueError(f'The {self.long_name} model does not support '
                             f'a model type of "{boundary_type}".')

    fbs = -1j/k0 * np.sum((-1)**n * (2*n+1) * A)
    return 20*log10(abs(fbs))  # ts

validate_parameters(params)

Validate the model parameters.

See here for calling details.

Source code in src/echosms/mssmodel.py 
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
def validate_parameters(self, params):
    """Validate the model parameters.

    See [here][echosms.ScatterModelBase.validate_parameters] for calling details.

    """
    p = as_dict(params)
    super()._present_and_in(p, ['boundary_type'], self.boundary_types)
    super()._present_and_positive(p, ['medium_rho', 'a', 'f'])

    types = np.unique(np.atleast_1d(p['boundary_type']))
    for t in types:
        mask = p['boundary_type'] == t
        match t:
            case bt.fluid_filled:
                super()._present_and_positive(p, ['target_c', 'target_rho'], mask=mask)
            case bt.fluid_shell_fluid_interior:
                super()._present_and_positive(p, ['target_c', 'target_rho', 'shell_c',
                                                  'shell_rho', 'shell_thickness'], mask=mask)
            case bt.fluid_shell_pressure_release_interior:
                super()._present_and_positive(p, ['shell_c', 'shell_rho', 'shell_thickness'],
                                              mask=mask)

echosms.PSMSModel()

Bases: ScatterModelBase

Prolate spheroidal modal series (PSMS) scattering model.

Note

The fluid filled boundary type implementation is under development and is of limited use at the moment.

Source code in src/echosms/psmsmodel.py 
18
19
20
21
22
23
24
25
def __init__(self):
    super().__init__()
    self.long_name = 'prolate spheroidal modal series'
    self.short_name = 'psms'
    self.analytical_type = 'exact'
    self.boundary_types = [bt.fixed_rigid, bt.pressure_release, bt.fluid_filled]
    self.shapes = ['prolate spheroid']
    self.max_ka = 10  # [1]

calculate_ts(data, expand=False, inplace=False, multiprocess=False, progress=False)

Calculate the target strength (TS) for many parameters.

PARAMETER DESCRIPTION
data

Requirements for the different input data types are:

  • DataFrame: column names must match the function parameter names in calculate_ts_single(). One TS value will be calculated for each row in the DataFrame.
  • DataArray: dimension names must match the function parameter names in calculate_ts_single(). TS values will be calculated for all combinations of the coordinate variables.
  • dict: keys must match the function parameters in calculate_ts_single(). TS values will be calculated for all combinations of the dict values.

TYPE: dict | DataFrame | DataArray

multiprocess

Split the ts calculation across CPU cores. Multiprocessing is currently provided by mapply. For more sophisticated uses it may be preferred to use a multiprocessing package of your choice directly on the calculate_ts_single() method. See the code in this method (calculate_ts()) for an example.

TYPE: bool DEFAULT: False

expand

Only applicable if data is a dict. If True, will use as_dataframe() to expand the dict into a DataFrame with one column per dict key and return that, adding a column named ts for the results.

TYPE: bool DEFAULT: False

inplace

Only applicable if data is a DataFrame. If True, the results will be added to the input DataFrame in a column named ts. If a ts column already exists, it is overwritten.

TYPE: bool DEFAULT: False

progress

If True, will produce a progress bar while running models

TYPE: bool DEFAULT: False

RETURNS DESCRIPTION
None | list[float] | Series | DataFrame

The return type and value are determined by the type of the input variable (data) and the expand and inplace parameters:

  • dict input and expand=False returns a list of floats.
  • dict input and expand=True returns a DataFrame.
  • DataFrame input and inplace=False returns a Series.
  • DataFrame input and inplace=True modifies data and returns None.
  • DataArray input always modifies data and returns None.
Source code in src/echosms/scattermodelbase.py 
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
def calculate_ts(self, data: dict | pd.DataFrame | xr.DataArray, expand: bool=False,
                 inplace: bool=False, multiprocess: bool=False,
                 progress: bool=False) -> None | list[float] | pd.Series | pd.DataFrame:
    """Calculate the target strength (TS) for many parameters.

    Parameters
    ----------
    data :
        Requirements for the different input data types are:

        - **DataFrame**: column names must match the function parameter names in
          calculate_ts_single(). One TS value will be calculated for each row in the DataFrame.
        - **DataArray**: dimension names must match the function parameter names in
          calculate_ts_single(). TS values will be calculated for all combinations of the
          coordinate variables.
        - **dict**: keys must match the function parameters in calculate_ts_single().
          TS values will be calculated for all combinations of the dict values.

    multiprocess :
        Split the ts calculation across CPU cores. Multiprocessing is currently provided by
        [mapply](https://github.com/ddelange/mapply). For more
        sophisticated uses it may be preferred to use a multiprocessing package of your choice
        directly on the `calculate_ts_single()` method. See the code in this method
        (`calculate_ts()`) for an example.

    expand :
        Only applicable if `data` is a dict. If `True`, will use
        [`as_dataframe()`][echosms.utils.as_dataframe]
        to expand the dict into a DataFrame with one column per dict key
        and return that, adding a column named `ts` for the results.

    inplace :
        Only applicable if `data` is a DataFrame. If `True`, the results
        will be added to the input DataFrame in a column named `ts`. If a `ts` column
        already exists, it is overwritten.

    progress :
        If `True`, will produce a progress bar while running models

    Returns
    -------
    :
        The return type and value are determined by the type of the input variable (`data`) and
        the `expand` and `inplace` parameters:

        - dict input and `expand=False` returns a list of floats.
        - dict input and `expand=True` returns a DataFrame.
        - DataFrame input and `inplace=False` returns a Series.
        - DataFrame input and `inplace=True` modifies `data` and returns `None`.
        - DataArray input always modifies `data` and returns `None`.

    """
    match data:
        case dict():
            data_df = as_dataframe(data, self.no_expand_parameters)
        case pd.DataFrame():
            data_df = data
        case xr.DataArray():
            data_df = data.to_dataframe().reset_index()
            data_df.attrs = data.attrs
        case _:
            raise ValueError(f'Data type of {type(data)} is not supported'
                             ' (only dictionaries, Pandas DataFrames and'
                             ' Xarray DataArrays are).')

    self.validate_parameters(data_df)

    # Get the non-expandable model parameters
    p = data_df.attrs['parameters'] if 'parameters' in data_df.attrs else {}

    # Note: the args argument in the apply call below requires a tuple. data_df.attrs is a
    # dict and the default behaviour is to make a tuple using the dict keys. The trailing comma
    # and parenthesis instead causes the tuple to have one entry of the dict.

    if multiprocess:
        from mapply.mapply import mapply
        ts = mapply(data_df, self.__ts_helper, args=(p,), axis=1, progressbar=progress)
    else:  # this uses just one CPU
        if progress:
            tqdm.pandas(desc=self.short_name, unit=' models',
                        bar_format='{l_bar}{bar} [{n_fmt}/{total_fmt}; {rate_noinv_fmt}]')
            ts = data_df.progress_apply(self.__ts_helper, args=(p,), axis=1)
        else:
            ts = data_df.apply(self.__ts_helper, args=(p,), axis=1)

    match data:
        case dict() if expand:
            data_df['ts'] = ts
            return data_df
        case dict():
            return ts.to_list()
        case pd.DataFrame() if inplace:
            data_df['ts'] = ts
            return None
        case pd.DataFrame():
            return ts.rename('ts', inplace=True)
        case xr.DataArray():
            data.values = ts.to_numpy().reshape(data.shape)
            return None
        case _:
            raise AssertionError('This code should never be reached - unsupported input data '
                                 f'type of {type(data)}.')

calculate_ts_single(medium_c, medium_rho, a, b, theta, f, boundary_type, target_c=None, target_rho=None, validate_parameters=True)

Prolate spheroid modal series (PSMS) solution model.

PARAMETER DESCRIPTION
medium_c

Sound speed in the fluid medium surrounding the target [m/s].

TYPE: float

medium_rho

Density of the fluid medium surrounding the target [kg/m³].

TYPE: float

a

Prolate spheroid major axis radius [m].

TYPE: float

b

Prolate spheroid minor axis radius [m].

TYPE: float

theta

Pitch angle to calculate the scattering as per the echoSMs coordinate system [°].

TYPE: float

f

Frequency to calculate the scattering at [Hz].

TYPE: float

boundary_type

The model type. Supported model types are given in the boundary_types class variable.

TYPE: boundary_type

target_c

Sound speed in the fluid inside the target [m/s]. Only required for boundary_type of fluid_filled.

TYPE: None | float DEFAULT: None

target_rho

Density of the fluid inside the target [kg/m³]. Only required for boundary_type of fluid_filled.

TYPE: None | float DEFAULT: None

validate_parameters

Whether to validate the input parameters.

TYPE: bool DEFAULT: True

RETURNS DESCRIPTION
float

The target strength (re 1 m²) of the target [dB].
Notes

The backscattered target strength of a pressure release or fluid-filled prolate spheroid is calculated using the PSMS method of Furusawa (1988) and corrections in Furusawa et al. (1994).

References

Furusawa, M. (1988). "Prolate spheroidal models for predicting general trends of fish target strength," J. Acoust. Soc. Jpn. 9, 13-24. https://doi.org/10.1250/ast.9.13

Furusawa, M., Miyanohana, Y., Ariji, M., and Sawada, Y. (1994). “Prediction of krill target strength by liquid prolate spheroid model,” Fish. Sci., 60, 261-265. https://doi.org/10.2331/FISHSCI.60.261

Source code in src/echosms/psmsmodel.py 
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
def calculate_ts_single(self, medium_c: float, medium_rho: float, a: float, b: float,
                        theta: float, f: float, boundary_type: bt,
                        target_c: None | float=None, target_rho: None | float=None,
                        validate_parameters: bool=True) -> float:
    """Prolate spheroid modal series (PSMS) solution model.

    Parameters
    ----------
    medium_c :
        Sound speed in the fluid medium surrounding the target [m/s].
    medium_rho :
        Density of the fluid medium surrounding the target [kg/m³].
    a :
        Prolate spheroid major axis radius [m].
    b :
        Prolate spheroid minor axis radius [m].
    theta :
        Pitch angle to calculate the scattering as per the echoSMs
        [coordinate system](conventions.md#coordinate-systems) [°].
    f :
        Frequency to calculate the scattering at [Hz].
    boundary_type :
        The model type. Supported model types are given in the `boundary_types` class variable.
    target_c :
        Sound speed in the fluid inside the target [m/s].
        Only required for `boundary_type` of ``fluid_filled``.
    target_rho :
        Density of the fluid inside the target [kg/m³].
        Only required for `boundary_type` of ``fluid_filled``.
    validate_parameters :
        Whether to validate the input parameters.

    Returns
    -------
    :
        The target strength (re 1 m²) of the target [dB].

    Notes
    -----
    The backscattered target strength of a pressure release or fluid-filled prolate spheroid
    is calculated using the PSMS method of Furusawa (1988) and corrections in
    Furusawa et al. (1994).

    References
    ----------
    Furusawa, M. (1988). "Prolate spheroidal models for predicting general
        trends of fish target strength," J. Acoust. Soc. Jpn. 9, 13-24.
        <https://doi.org/10.1250/ast.9.13>

    Furusawa, M., Miyanohana, Y., Ariji, M., and Sawada, Y. (1994).
        “Prediction of krill target strength by liquid prolate spheroid
        model,” Fish. Sci., 60, 261-265. <https://doi.org/10.2331/FISHSCI.60.261>
    """
    if validate_parameters:
        self.validate_parameters(locals())

    if boundary_type not in self.boundary_types:
        raise ValueError(f'The {self.long_name} model does not support '
                         f'a model type of "{boundary_type}".')

    xim = (1.0 - (b/a)**2)**(-0.5)
    q = a/xim  # semi-focal length

    km = wavenumber(medium_c, f)
    hm = km*q

    if boundary_type == bt.fluid_filled:
        g = target_rho / medium_rho
        ht = wavenumber(target_c, f)*q
        simplified = False
        # Note: we can implement simpler equations if sound speeds are
        # similar between the medium and the target. The simplified
        # equations are quicker, so it is worth to do.
        # But, it appears that target_c ≈ medium_c is not the only requirement for
        # the simplification to work well - see section 4.1.1 in:
        #
        # Gonzalez, J. D., Lavia, E. F., Blanc, S., & Prario, I. (2016).
        # Acoustic scattering by prolate and oblate liquid spheroids.
        # Proceedings of the 22nd International Congress on Acoustics.
        # Acoustics for the 21st Century, Buenos Aires.
        #
        # So, for the moment, the simplification is turned off.

        # if abs(1.0-target_c/medium_c) <= 0.01:
        #    simplified = True

    # Phi, the roll angle is fixed for this implementation
    phi_inc = np.pi  # incident direction
    phi_sca = np.pi + phi_inc  # scattered direction

    theta_inc = np.deg2rad(theta)  # incident direction
    theta_sca = np.pi - theta_inc  # scattered direction

    # Approximate limits on the summations. These come from Furusawa (1988), but other
    # implementations of this code use more complex functions to calculate the maximum orders
    # of spheroidal wave functions to calculate. It is also feasible to work to a defined
    # convergence level. This is a potenital future improvement.
    m_max = int(np.ceil(2*km*b))  # +1
    n_max = int(m_max + np.ceil(hm/2))  # +3

    f_sca = 0.0
    for m in range(m_max+1):
        epsilon_m = Neumann(m)
        cos_term = np.cos(m*(phi_sca - phi_inc))

        if boundary_type == bt.fluid_filled and not simplified:
            Am = PSMSModel._fluidfilled(m, n_max, hm, ht, xim, g, theta_inc)

        for n_i, n in enumerate(range(m, n_max+1)):
            # Use the normalisation offered by spheroidalwavefunctions.pro_ang1() because
            # it removes the need to calculate a normalisation factor (N_mn in Furusawa, 1998).
            # This is because the pro_ang1(norm=True) norm is unity.
            Smn_inc, _ = pro_ang1(m, n, hm, np.cos(theta_inc), norm=True)
            Smn_sca, _ = pro_ang1(m, n, hm, np.cos(theta_sca), norm=True)

            # FYI, the code used to use the Meixner-Schäfke normalisation scheme for the
            # angular function of the first kind (eqn 21.7.11 in Abramowitz, M., and Stegun,
            # I. A. (1964). Handbook of Mathematical Functions with Formulas, Graphs, and
            # Mathematical Tables # (Dover, New York), 10th ed.
            # N_mn = 2/(2*n+1) * factorial(n+m) / factorial(n-m)

            R1m, dR1m = pro_rad1(m, n, hm, xim)
            R2m, dR2m = pro_rad2(m, n, hm, xim)

            match boundary_type:
                case bt.fluid_filled:
                    if simplified:
                        E1, E3 = PSMSModel._fluidfilled_Emn(m, n, n, hm, ht, xim, g)
                        Amn = -E1/E3
                    else:
                        Amn = Am[n_i][0]
                case bt.pressure_release:
                    Amn = -R1m/(R1m + 1j*R2m)
                case bt.fixed_rigid:
                    Amn = -dR1m/(dR1m + 1j*dR2m)

            f_sca += epsilon_m * Smn_inc * Amn * Smn_sca * cos_term

    return 20*np.log10(np.abs(-2j / km * f_sca))

validate_parameters(params)

Validate the model parameters.

See here for calling details.

Source code in src/echosms/psmsmodel.py 
27
28
29
30
31
32
33
34
35
36
37
38
39
40
def validate_parameters(self, params):
    """Validate the model parameters.

    See [here][echosms.ScatterModelBase.validate_parameters] for calling details.
    """
    p = as_dict(params)
    super()._present_and_in(p, ['boundary_type'], self.boundary_types)
    super()._present_and_positive(p, ['medium_c', 'medium_rho', 'a', 'b', 'f'])

    types = np.unique(np.atleast_1d(p['boundary_type']))
    for t in types:
        if t == bt.fluid_filled:
            super()._present_and_positive(p, ['target_c', 'target_rho'],
                                          mask=p['boundary_type'] == t)