Skip to content

API reference

This is the API reference for 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. There are also utility functions.

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.

Attributes:

Name Type Description
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[str]

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.
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[str]
        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.

Parameters:

Name Type Description Default
data Pandas DataFrame, Xarray DataArray or dict

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.
 required
multiprocess bool

Split the ts calculation across CPU cores. Multiprocessing is currently provided by mapply with little customisation. 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.

 False
expand bool

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.

 False
inplace bool

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.

 False
progress bool

If True, will produce a progress bar while running models

 False

Returns:

Type Description
None, list[float], Series, or 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
def calculate_ts(self, data, expand=False, inplace=False, multiprocess=False, progress=False):
    """Calculate the target strength (TS) for many parameters.

    Parameters
    ----------
    data : Pandas DataFrame, Xarray DataArray or dict
        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 : bool
        Split the ts calculation across CPU cores. Multiprocessing is currently provided by
        [mapply](https://github.com/ddelange/mapply) with little customisation. 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 : bool
        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 : bool
        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 : bool
        If `True`, will produce a progress bar while running models

    Returns
    -------
    : None, list[float], Series, or 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`.

    """
    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)
    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 
236
237
238
@abc.abstractmethod
def calculate_ts_single(self):
    """Calculate the TS for one parameter set."""

validate_parameters(p) abstractmethod

Validate the model parameters.

Parameters:

Name Type Description Default
p dict | DataFrame | DataArray

The model parameters.

 required

Raises:

Type 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 
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
@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 
19
20
21
22
23
24
25
26
def __init__(self):
    super().__init__()
    self.long_name = 'deformed cylinder model'
    self.short_name = 'dcm'
    self.analytical_type = 'approximate analytical'
    self.boundary_types = ['fixed rigid', 'pressure release', '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.

Parameters:

Name Type Description Default
data Pandas DataFrame, Xarray DataArray or dict

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.
 required
multiprocess bool

Split the ts calculation across CPU cores. Multiprocessing is currently provided by mapply with little customisation. 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.

 False
expand bool

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.

 False
inplace bool

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.

 False
progress bool

If True, will produce a progress bar while running models

 False

Returns:

Type Description
None, list[float], Series, or 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
def calculate_ts(self, data, expand=False, inplace=False, multiprocess=False, progress=False):
    """Calculate the target strength (TS) for many parameters.

    Parameters
    ----------
    data : Pandas DataFrame, Xarray DataArray or dict
        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 : bool
        Split the ts calculation across CPU cores. Multiprocessing is currently provided by
        [mapply](https://github.com/ddelange/mapply) with little customisation. 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 : bool
        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 : bool
        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 : bool
        If `True`, will produce a progress bar while running models

    Returns
    -------
    : None, list[float], Series, or 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`.

    """
    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)
    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.

Parameters:

Name Type Description Default
medium_c float

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

 required
medium_rho float

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

 required
a float

Radius of the cylinderical target [m].

 required
b float

Length of the cylinderical target [m].

 required
theta float

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

 required
f float

Frequency to calculate the scattering at [Hz].

 required
boundary_type str

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

 required
target_c float

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

 None
target_rho float

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

 None
validate_parameters bool

Whether to validate the model parameters.

 True

Returns:

Type 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
def calculate_ts_single(self, 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.

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

    Returns
    -------
    : 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>
    """
    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 'fixed rigid':
            series = map(lambda m: (-1)**m * Neumann(m)*(jvp(m, Ka) / h1vp(m, Ka)), m)
        case 'pressure release':
            series = map(lambda m: (-1)**m * Neumann(m)*(jv(m, Ka) / hankel1(m, Ka)), m)
        case '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 
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'])

    for bt in np.atleast_1d(p['boundary_type']):
        if bt == 'fluid filled':
            super()._present_and_positive(p, ['target_c', 'target_rho'])

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 
21
22
23
24
25
26
27
28
29
30
31
32
33
def __init__(self):
    super().__init__()
    self.long_name = 'distorted-wave Born approximation'
    self.short_name = 'dwba'
    self.analytical_type = 'approximate'
    self.boundary_types = ['weakly scattering']
    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.

Parameters:

Name Type Description Default
data Pandas DataFrame, Xarray DataArray or dict

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.
 required
multiprocess bool

Split the ts calculation across CPU cores. Multiprocessing is currently provided by mapply with little customisation. 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.

 False
expand bool

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.

 False
inplace bool

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.

 False
progress bool

If True, will produce a progress bar while running models

 False

Returns:

Type Description
None, list[float], Series, or 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
def calculate_ts(self, data, expand=False, inplace=False, multiprocess=False, progress=False):
    """Calculate the target strength (TS) for many parameters.

    Parameters
    ----------
    data : Pandas DataFrame, Xarray DataArray or dict
        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 : bool
        Split the ts calculation across CPU cores. Multiprocessing is currently provided by
        [mapply](https://github.com/ddelange/mapply) with little customisation. 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 : bool
        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 : bool
        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 : bool
        If `True`, will produce a progress bar while running models

    Returns
    -------
    : None, list[float], Series, or 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`.

    """
    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)
    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.

Parameters:

Name Type Description Default
medium_c float

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

 required
medium_rho float

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

 required
theta float

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

 required
phi float

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

 required
f float

Frequency to run the model at [Hz]

 required
target_c float

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

 required
target_rho float

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

 required
a iterable

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

 required
rv_pos iterable[ndarray]

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.

 required
rv_tan iterable[ndarray]

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.

 required
phase_sd float

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 [°].

 0
num_runs int

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.

 1
validate_parameters bool

Whether to validate the model parameters.

 True

Returns:

Type 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 
 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
def calculate_ts_single(self, 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) -> 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 : float
        Sound speed in the fluid medium surrounding the target [m/s].
    medium_rho : float
        Density of the fluid medium surrounding the target [kg/m³].
    theta : float
        Pitch angle to calculate the scattering as per the echoSMs
        [coordinate system](https://ices-tools-dev.github.io/echoSMs/
        conventions/#coordinate-systems) [°].
    phi : float
        Roll angle to calculate the scattering as per the echoSMs
        [coordinate system](https://ices-tools-dev.github.io/echoSMs/
        conventions/#coordinate-systems) [°].
    f : float
        Frequency to run the model at [Hz]
    target_c : float
        Sound speed in the fluid inside the target [m/s].
    target_rho : float
        Density of the fluid inside the target [kg/m³].
    a : iterable
        The radii of the discs that define the target shape [m].
    rv_pos : iterable[np.ndarray]
        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 : iterable[np.ndarray]
        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 : float
        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 : int
        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 : bool
        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 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 
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
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 
13
14
15
16
17
18
19
20
21
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 = ['weakly scattering']
    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.

Parameters:

Name Type Description Default
data Pandas DataFrame, Xarray DataArray or dict

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.
 required
multiprocess bool

Split the ts calculation across CPU cores. Multiprocessing is currently provided by mapply with little customisation. 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.

 False
expand bool

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.

 False
inplace bool

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.

 False
progress bool

If True, will produce a progress bar while running models

 False

Returns:

Type Description
None, list[float], Series, or 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
def calculate_ts(self, data, expand=False, inplace=False, multiprocess=False, progress=False):
    """Calculate the target strength (TS) for many parameters.

    Parameters
    ----------
    data : Pandas DataFrame, Xarray DataArray or dict
        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 : bool
        Split the ts calculation across CPU cores. Multiprocessing is currently provided by
        [mapply](https://github.com/ddelange/mapply) with little customisation. 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 : bool
        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 : bool
        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 : bool
        If `True`, will produce a progress bar while running models

    Returns
    -------
    : None, list[float], Series, or 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`.

    """
    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)
    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.

Parameters:

Name Type Description Default
volume Numpy ndarray[int]

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.

 required
theta float

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

 required
phi float

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

 required
f float

Frequency to run the model at [Hz]

 required
voxel_size iterable[float]

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.

 required
rho iterable[float]

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

 required
c iterable[float]

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

 required
validate_parameters bool

Whether to validate the model parameters.

 True

Returns:

Type 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 
 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
 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, 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.

    Parameters
    ----------
    volume : Numpy ndarray[int]
        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](https://ices-tools-dev.github.io/echoSMs/
        conventions/#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 : float
        Pitch angle to calculate the scattering as per the echoSMs
        [coordinate system](https://ices-tools-dev.github.io/echoSMs/
        conventions/#coordinate-systems) [°].

    phi : float
        Roll angle to calculate the scattering as per the echoSMs
        [coordinate system](https://ices-tools-dev.github.io/echoSMs/
        conventions/#coordinate-systems) [°].

    f : float
        Frequency to run the model at [Hz]

    voxel_size : iterable[float]
        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 : iterable[float]
        Densities of each material. Must have at least the same number of entries as unique
        integers in `volume` [kg/m³].

    c : iterable[float]
        Sound speed of each material. Must have at least the same number of entries as unique
        integers in `volume` [m/s].
    validate_parameters : bool
        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 
23
24
25
26
27
28
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 
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 __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 SDWBAorganism().
    self.dwba_models = {}
    for s in shapes['shape']:
        # 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)
        self.dwba_models[s['name']] = organism
as_dict()

DWBA model shapes as a dict.

Returns:

Type 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 
185
186
187
188
189
190
191
192
193
194
195
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.

Parameters:

Name Type Description Default
name str

The name of a DWBA model shape.

 required

Returns:

Type Description
DWBAorganism

An instance of DWBAorganism or None if there is no model with name.
Source code in src/echosms/dwbautils.py 
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
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 
181
182
183
def names(self):
    """Available DWBA model names."""
    return [*self.dwba_models]

DWBAorganism(rv_pos, a, g, h, name, source='', note='', rv_tan=None) dataclass

DWBA shape and property class to represent an organism.

Attributes:

Name Type Description
rv_pos iterable[ndarray]

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, as per the echoSMs coordinate system.
a ndarray

The radii [m] of each disc.
g ndarray

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

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

A link to the source of the data.
note str

Information about the data.
rv_tan ndarray

An interable 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.
plot()

Do a simple plot of the DWBA model data.

Source code in src/echosms/dwbautils.py 
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
def plot(self):
    """Do a simple plot of the DWBA model data."""
    import matplotlib.pyplot as plt
    fig, axs = plt.subplots(2, 1)
    x = self.rv_pos[:, 0]*1e3
    axs[0].plot(x, -self.rv_pos[:, 1]*1e3, '.-', c='C0')
    axs[0].plot(x, (-self.rv_pos[:, 1]+self.a/2)*1e3, c='C1')
    axs[0].plot(x, (-self.rv_pos[:, 1]-self.a/2)*1e3, c='C1')
    axs[0].set_title('Dorsal', loc='left', fontsize=8)
    axs[0].set_aspect('equal')

    axs[1].plot(x, -self.rv_pos[:, 2]*1e3, '.-', c='C0')
    axs[1].plot(x, (-self.rv_pos[:, 2]+self.a/2)*1e3, c='C1')
    axs[1].plot(x, (-self.rv_pos[:, 2]-self.a/2)*1e3, c='C1')
    axs[1].set_title('Lateral', loc='left', fontsize=8)
    axs[1].set_aspect('equal')

    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.

Parameters:

Name Type Description Default
radius float

The radius [m] of the cylinder.

 required
length float

The length [m] of the cylinder.

 required
spacing float

The spacing [m] between successive discs.

 0.0001

Returns:

Name Type Description
rv_pos iterable[ndarray]

An interable 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[ndarray]

An interable 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.
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 interable 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 interable 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_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).

Parameters:

Name Type Description Default
major_radius float

The major radius [m] of the spheroid.

 required
minor_radius float

The minor radius [m] of the spheroid.

 required
spacing float

The spacing [m] between successive discs.

 0.0001

Returns:

Name Type Description
rv_pos iterable[ndarray]

An interable 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[ndarray]

An interable 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.
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 interable 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 interable 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 = ['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.

Parameters:

Name Type Description Default
data Pandas DataFrame, Xarray DataArray or dict

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.
 required
multiprocess bool

Split the ts calculation across CPU cores. Multiprocessing is currently provided by mapply with little customisation. 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.

 False
expand bool

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.

 False
inplace bool

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.

 False
progress bool

If True, will produce a progress bar while running models

 False

Returns:

Type Description
None, list[float], Series, or 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
def calculate_ts(self, data, expand=False, inplace=False, multiprocess=False, progress=False):
    """Calculate the target strength (TS) for many parameters.

    Parameters
    ----------
    data : Pandas DataFrame, Xarray DataArray or dict
        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 : bool
        Split the ts calculation across CPU cores. Multiprocessing is currently provided by
        [mapply](https://github.com/ddelange/mapply) with little customisation. 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 : bool
        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 : bool
        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 : bool
        If `True`, will produce a progress bar while running models

    Returns
    -------
    : None, list[float], Series, or 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`.

    """
    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)
    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.

Parameters:

Name Type Description Default
medium_c float

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

 required
medium_rho float

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

 required
a float

Radius of the sphere [m].

 required
f float

Frequency to calculate the scattering at [Hz].

 required
target_longitudinal_c float

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

 required
target_transverse_c float

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

 required
target_rho float

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

 required
validate_parameters bool

Whether to validate the model parameters.

 True

Returns:

Type 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.

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, medium_rho, a, f,
                        target_longitudinal_c, target_transverse_c, target_rho,
                        validate_parameters=True,
                        **kwargs) -> float:
    """
    Calculate the backscatter from an elastic sphere for one set of parameters.

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

    Returns
    -------
    : 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.
    """
    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 = ['fluid filled', 'elastic', '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.

Parameters:

Name Type Description Default
data Pandas DataFrame, Xarray DataArray or dict

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.
 required
multiprocess bool

Split the ts calculation across CPU cores. Multiprocessing is currently provided by mapply with little customisation. 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.

 False
expand bool

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.

 False
inplace bool

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.

 False
progress bool

If True, will produce a progress bar while running models

 False

Returns:

Type Description
None, list[float], Series, or 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
def calculate_ts(self, data, expand=False, inplace=False, multiprocess=False, progress=False):
    """Calculate the target strength (TS) for many parameters.

    Parameters
    ----------
    data : Pandas DataFrame, Xarray DataArray or dict
        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 : bool
        Split the ts calculation across CPU cores. Multiprocessing is currently provided by
        [mapply](https://github.com/ddelange/mapply) with little customisation. 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 : bool
        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 : bool
        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 : bool
        If `True`, will produce a progress bar while running models

    Returns
    -------
    : None, list[float], Series, or 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`.

    """
    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)
    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.

Parameters:

Name Type Description Default
shape str

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

 required
medium_c float

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

 required
medium_rho float

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

 None
target_c float

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

 None
target_rho float

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

 None
a float

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

 required
f float

Frequency to calculate the scattering at [Hz].

 required
boundary_type str

The boundary type for the model.

 required
theta float

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

 None
L float

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.

 None
rho_c float

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

 None
irregular

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

 False
validate_parameters bool

Whether to validate the model parameters.

 True

Returns:

Type 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
194
def calculate_ts_single(self, 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) -> float:
    """
    Calculate the backscatter using the high pass model for one set of parameters.

    Parameters
    ----------
    shape : str
        The shape to model. Must be one of the shapes given in the `shapes` variable.
    medium_c : float
        Sound speed in the fluid medium surrounding the target [m/s].
    medium_rho : float
        Density of the fluid medium surrounding the target [kg/m³]. Not required when
        `boundary_type` is `fixed rigid`.
    target_c : float
        Longitudinal sound speed in the material inside the target [m/s]. Not required when
        `boundary_type` is `fixed rigid`.
    target_rho : float
        Density of the material inside the target [kg/m³]. Not required when
        `boundary_type` is `fixed rigid`.
    a : float
        Radius of the sphere, length of semi-minor axis of the prolate spheriod, or cylindrical
        radius of the straight or bent cylinder [m].
    f : float
        Frequency to calculate the scattering at [Hz].
    boundary_type : str
        The boundary type for the model.
    theta : float
        Pitch angle to calculate the scattering as per the echoSMs
        [coordinate system](https://ices-tools-dev.github.io/echoSMs/
        conventions/#coordinate-systems) [°]. Only required for the straight cylinder shape.
    L : float
        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 : float
        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 : bool
        Whether to validate the model parameters.

    Returns
    -------
    : 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](https://ices-tools-dev.github.io/echoSMs/
    conventions/#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 == '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 'fluid filled':
                        F = 40 * (k*a)**(-0.4)
                        G = 1-0.8*exp(-2.5*(k*a-2.25)**2)
                    case 'elastic':
                        F = 15 * (k*a)**(-1.9)
                    case 'rigid fixed':
                        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 'fluid filled':
                        F = 2.5 * (k*a)**(1.65)
                        G = 1-0.8*exp(-2.5*(k*a-2.3)**2)
                    case 'elastic':
                        F = 1.8 * (k*a)**(-0.4)
                    case 'rigid fixed':
                        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 'fluid filled':
                        F = 3 * (k*a)**(0.65)
                        G = 1-0.8*exp(-2.5*(k*a-2.0)**2)
                    case 'elastic':
                        F = 3.5 * (k*a)**(-1.0)
                    case 'rigid fixed':
                        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 'fluid filled':
                        F = 3.0 * (k*a)**(0.65)
                        G = 1-0.8*exp(-2.5*(k*a-2.0)**2)
                    case 'elastic':
                        F = 2.5 * (k*a)**(-1.0)
                    case 'rigid fixed':
                        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 
16
17
18
19
20
21
22
23
24
def __init__(self):
    super().__init__()
    self.long_name = 'Kirchhoff approximation'
    self.short_name = 'ka'
    self.analytical_type = 'approximate'
    self.boundary_types = ['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.

Parameters:

Name Type Description Default
data Pandas DataFrame, Xarray DataArray or dict

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.
 required
multiprocess bool

Split the ts calculation across CPU cores. Multiprocessing is currently provided by mapply with little customisation. 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.

 False
expand bool

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.

 False
inplace bool

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.

 False
progress bool

If True, will produce a progress bar while running models

 False

Returns:

Type Description
None, list[float], Series, or 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
def calculate_ts(self, data, expand=False, inplace=False, multiprocess=False, progress=False):
    """Calculate the target strength (TS) for many parameters.

    Parameters
    ----------
    data : Pandas DataFrame, Xarray DataArray or dict
        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 : bool
        Split the ts calculation across CPU cores. Multiprocessing is currently provided by
        [mapply](https://github.com/ddelange/mapply) with little customisation. 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 : bool
        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 : bool
        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 : bool
        If `True`, will produce a progress bar while running models

    Returns
    -------
    : None, list[float], Series, or 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`.

    """
    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)
    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.

Parameters:

Name Type Description Default
medium_c float

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

 required
theta float

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

 required
phi float

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

 required
f float

Frequency to calculate the scattering at [Hz].

 required
mesh Any

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.

 required
boundary_type str

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

 required
validate_parameters bool

Whether to validate the model parameters.

 True

Returns:

Type 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 
 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
def calculate_ts_single(self, medium_c, theta, phi, f, mesh,
                        boundary_type, validate_parameters=True, **kwargs) -> float:
    """
    Calculate the scatter using the ka model for one set of parameters.

    Parameters
    ----------
    medium_c : float
        Sound speed in the fluid medium surrounding the target [m/s].
    theta : float
        Pitch angle to calculate the scattering as per the echoSMs
        [coordinate system](https://ices-tools-dev.github.io/echoSMs/
        conventions/#coordinate-systems) [°].
    phi : float
        Roll angle to calculate the scattering as per the echoSMs
        [coordinate system](https://ices-tools-dev.github.io/echoSMs/
        conventions/#coordinate-systems) [°].
    f : float
        Frequency to calculate the scattering at [Hz].
    mesh : Any
        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).
    boundary_type : str
        The boundary type. Supported types are given in the `boundary_types` class variable.
    validate_parameters : bool
        Whether to validate the model parameters.

    Returns
    -------
    : 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>

    """
    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 
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_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 = ['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.

Parameters:

Name Type Description Default
data Pandas DataFrame, Xarray DataArray or dict

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.
 required
multiprocess bool

Split the ts calculation across CPU cores. Multiprocessing is currently provided by mapply with little customisation. 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.

 False
expand bool

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.

 False
inplace bool

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.

 False
progress bool

If True, will produce a progress bar while running models

 False

Returns:

Type Description
None, list[float], Series, or 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
def calculate_ts(self, data, expand=False, inplace=False, multiprocess=False, progress=False):
    """Calculate the target strength (TS) for many parameters.

    Parameters
    ----------
    data : Pandas DataFrame, Xarray DataArray or dict
        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 : bool
        Split the ts calculation across CPU cores. Multiprocessing is currently provided by
        [mapply](https://github.com/ddelange/mapply) with little customisation. 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 : bool
        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 : bool
        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 : bool
        If `True`, will produce a progress bar while running models

    Returns
    -------
    : None, list[float], Series, or 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`.

    """
    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)
    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, validate_parameters=True, **kwargs)

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

Warning

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

Parameters:

Name Type Description Default
medium_c float

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

 required
medium_rho float

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

 required
theta float

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

 required
f float

Frequency to calculate the scattering at [Hz].

 required
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.

 required
validate_parameters bool

Whether to validate the model parameters.

 True

Returns:

Type 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).

References

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

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

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
def calculate_ts_single(self, medium_c, medium_rho, theta, f, organism,
                        validate_parameters=True, **kwargs) -> float:
    """
    Calculate the scatter using the Kirchhoff ray mode model for one set of parameters.

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

    Parameters
    ----------
    medium_c : float
        Sound speed in the fluid medium surrounding the target [m/s].
    medium_rho : float
        Density in the fluid medium surrounding the target [kg/m³]
    theta : float
        Pitch angle to calculate the scattering at, as per the echoSMs
        [coordinate system](https://ices-tools-dev.github.io/echoSMs/
        conventions/#coordinate-systems) [°].
    f : float
        Frequency to calculate the scattering at [Hz].
    organism: KRMorganism
        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.
    validate_parameters : bool
        Whether to validate the model parameters.

    Returns
    -------
    : 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).

    References
    ----------
    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>

    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>
    """
    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
            sl.append(self._mode_solution(1/gp, 1/hp, k, a_e, incl.length(), theta))
        elif incl.boundary == 'soft':
            sl.append(self._soft_KA(incl, k, k_b, R_bc, TwbTbw, theta))
        elif incl.boundary == 'fluid':
            sl.append(self._fluid_KA(incl, k, k_b, 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 
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
def __init__(self):
    # Load in the NOAA KRM shapes data
    self.file = Path(__file__).parent/Path('resources')/Path('NOAA_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']:
        body = KRMshape('fluid', 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('soft', 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])

as_dict()

KRM model shapes as a dict.

Returns:

Type 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 
118
119
120
121
122
123
124
125
126
127
128
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.

Parameters:

Name Type Description Default
name str

The name of a KRM model shape.

 required

Returns:

Type Description
KRMorganism

An instance of KRMorganism or None if there is no model with name.
Source code in src/echosms/krmdata.py 
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
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 
114
115
116
def names(self):
    """Available KRM model names."""
    return [*self.krm_models]

ts(name) staticmethod

KRM model ts from model name.

Parameters:

Name Type Description Default
name str

The name of a KRM model shape.

 required

Returns:

Type 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 
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
@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) dataclass

KRM body and swimbladder model shape.

Attributes:

Name Type Description
name str

A name for the organism.
source str

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

The shape that represents the organism's body.
inclusions List[KRMshape]

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

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

KRM shape and property class.

Attributes:

Name Type Description
shape_type

The shape bounday condition - either soft or fluid.
x ndarray

The x-axis coordinates [m].
w ndarray

Width of the shape [m].
z_U ndarray

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

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

Sound speed in the shape [m/s].
rho float

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

length()

Length of the shape.

Returns:

Type Description
float

The length of the shape [m].
Source code in src/echosms/krmdata.py 
57
58
59
60
61
62
63
64
65
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:

Type Description
float

The volume of the shape [m³].
Source code in src/echosms/krmdata.py 
45
46
47
48
49
50
51
52
53
54
55
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 
19
20
21
22
23
24
25
26
27
28
def __init__(self):
    super().__init__()
    self.long_name = 'modal series solution'
    self.short_name = 'mss'
    self.analytical_type = 'exact'
    self.boundary_types = ['fixed rigid', 'pressure release', 'fluid filled',
                           'fluid shell fluid interior',
                           '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.

Parameters:

Name Type Description Default
data Pandas DataFrame, Xarray DataArray or dict

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.
 required
multiprocess bool

Split the ts calculation across CPU cores. Multiprocessing is currently provided by mapply with little customisation. 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.

 False
expand bool

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.

 False
inplace bool

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.

 False
progress bool

If True, will produce a progress bar while running models

 False

Returns:

Type Description
None, list[float], Series, or 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
def calculate_ts(self, data, expand=False, inplace=False, multiprocess=False, progress=False):
    """Calculate the target strength (TS) for many parameters.

    Parameters
    ----------
    data : Pandas DataFrame, Xarray DataArray or dict
        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 : bool
        Split the ts calculation across CPU cores. Multiprocessing is currently provided by
        [mapply](https://github.com/ddelange/mapply) with little customisation. 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 : bool
        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 : bool
        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 : bool
        If `True`, will produce a progress bar while running models

    Returns
    -------
    : None, list[float], Series, or 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`.

    """
    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)
    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.

Parameters:

Name Type Description Default
medium_c float

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

 required
medium_rho float

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

 required
a float

Radius of the spherical target [m].

 required
f float

Frequency to calculate the scattering at [Hz].

 required
boundary_type str

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

 required
target_c float

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

 None
target_rho float

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

 None
shell_c float

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

 None
shell_rho float

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

 None
shell_thickness float

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.

 None
validate_parameters bool

Whether to validate the model parameters.

 True

Returns:

Type 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 
 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
def calculate_ts_single(self, 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) -> float:
    """
    Calculate the scatter using the mss model for one set of parameters.

    Parameters
    ----------
    medium_c : float
        Sound speed in the fluid medium surrounding the target [m/s].
    medium_rho : float
        Density of the fluid medium surrounding the target [kg/m³].
    a : float
        Radius of the spherical target [m].
    f : float
        Frequency to calculate the scattering at [Hz].
    boundary_type : str
        The boundary type. Supported types are given in the `boundary_types` class variable.
    target_c : float, optional
        Sound speed in the fluid inside the sphere [m/s].
        Only required for `boundary_type` of ``fluid filled``.
    target_rho : float, optional
        Density of the fluid inside the sphere [kg/m³].
        Only required for `boundary_type` of ``fluid filled``.
    shell_c : float, optional
        Sound speed in the spherical shell [m/s].
        Only required for `boundary_type`s that include a fluid shell.
    shell_rho : float, optional
        Density in the spherical shell [kg/m³].
        Only required for `boundary_type`s that include a fluid shell.
    shell_thickness : float, optional
        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 : bool
        Whether to validate the model parameters.

    Returns
    -------
    : 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>
    """
    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 'fixed rigid':
            A = list(map(lambda x: -spherical_jn(x, ka, True) / h1(x, ka, True), n))
        case 'pressure release':
            A = list(map(lambda x: -spherical_jn(x, ka) / h1(x, ka), n))
        case '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 '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 '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 
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
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'])

    for bt in np.atleast_1d(p['boundary_type']):
        match bt:
            case 'fluid filled':
                super()._present_and_positive(p, ['target_c', 'target_rho'])
            case 'fluid shell fluid interior':
                super()._present_and_positive(p, ['target_c', 'target_rho', 'shell_c',
                                                  'shell_rho', 'shell_thickness'])
            case 'fluid shell pressure release interior':
                super()._present_and_positive(p, ['shell_c', 'shell_rho', 'shell_thickness'])

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 = ['fixed rigid', 'pressure release', '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.

Parameters:

Name Type Description Default
data Pandas DataFrame, Xarray DataArray or dict

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.
 required
multiprocess bool

Split the ts calculation across CPU cores. Multiprocessing is currently provided by mapply with little customisation. 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.

 False
expand bool

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.

 False
inplace bool

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.

 False
progress bool

If True, will produce a progress bar while running models

 False

Returns:

Type Description
None, list[float], Series, or 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
def calculate_ts(self, data, expand=False, inplace=False, multiprocess=False, progress=False):
    """Calculate the target strength (TS) for many parameters.

    Parameters
    ----------
    data : Pandas DataFrame, Xarray DataArray or dict
        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 : bool
        Split the ts calculation across CPU cores. Multiprocessing is currently provided by
        [mapply](https://github.com/ddelange/mapply) with little customisation. 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 : bool
        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 : bool
        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 : bool
        If `True`, will produce a progress bar while running models

    Returns
    -------
    : None, list[float], Series, or 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`.

    """
    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)
    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.

Parameters:

Name Type Description Default
medium_c float

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

 required
medium_rho float

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

 required
a float

Prolate spheroid major axis radius [m].

 required
b float

Prolate spheroid minor axis radius [m].

 required
theta float

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

 required
f float

Frequency to calculate the scattering at [Hz].

 required
boundary_type str

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

 required
target_c float

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

 None
target_rho float

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

 None
validate_parameters bool

Whether to validate the input parameters.

 True

Returns:

Type 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. 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.

Source code in src/echosms/psmsmodel.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
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
def calculate_ts_single(self, 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.

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

    Returns
    -------
    : 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.
    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.
    """
    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 == '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 == '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 '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 'pressure release':
                    Amn = -R1m/(R1m + 1j*R2m)
                case '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
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'])

    for bt in np.atleast_1d(p['boundary_type']):
        match bt:
            case 'fluid filled':
                super()._present_and_positive(p, ['target_c', 'target_rho'])

Utilities

echosms.BenchmarkData()

Convenient interface to the benchmark dataset.

This dataset contains the TS results from Jech et al. (2015).

Attributes:

Name Type Description
angle_dataset Pandas DataFrame

The angle dataset from the benchmark model runs.
freq_dataset Pandas DataFrame

The frequency dataset from the benchmark model runs.
Notes

The column names in the source benchmark files have been changed to be the same as those used in the [referencemodels] model definitions.

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/benchmarkdata.py 
59
60
61
62
63
64
65
66
67
68
69
70
71
def __init__(self):

    data_directory = Path(__file__).parent/Path('resources')/Path('BenchMark_Data')

    angle_data_file = data_directory/'Benchmark_Angle_TS.csv'
    freq_data_file = data_directory/'Benchmark_Frequency_TS.csv'

    self.angle_dataset = pd.read_csv(angle_data_file)
    self.freq_dataset = pd.read_csv(freq_data_file)

    # Change the column names to match the reference model names used in ReferenceModels
    self.angle_dataset.rename(columns=BenchmarkData.a_rename, inplace=True)
    self.freq_dataset.rename(columns=BenchmarkData.f_rename, inplace=True)

echosms.ReferenceModels()

Provide access to reference scattering model parameters.

Reference models are the models and parameters defined in Jech et al. (2015). The parameters are stored in a TOML-formatted file in the echoSMs repository and this class provides easy access to the data in that file. Additional reference models may be defined in the future and added to the TOML file (for example, entries have been added for all known calibration sphere sizes).

Attributes:

Name Type Description
definitions dict

A dict representation of the target definitions.toml file.

Raises:

Type Description
TOMLDecodeError

If the target definitions.toml file is not valid TOML.
KeyError

If the target definitions.toml file has multiple target entries with the same name.
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/referencemodels.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
def __init__(self):
    self.defs_filename = Path(__file__).parent/Path('resources')/Path('target definitions.toml')

    self.definitions = []

    with open(self.defs_filename, 'rb') as f:
        try:
            self.definitions = tomllib.load(f)
        except tomllib.TOMLDecodeError as e:
            raise SyntaxError(f'Error while parsing file "{self.defs_filename.name}"') from e

    # Flag duplicate target names
    pda = pd.Series(self.names())
    duplicates = list(pda[pda.duplicated()])
    if duplicates:
        raise KeyError(f'The "{self.defs_filename.name}" file has multiple targets '
                       f'with the same name: '+', '.join(duplicates))

    # Substitute parameters names in the target section by the values of those
    # parameters.
    for t in self.definitions['target']:
        for k, v in t.items():
            try:
                t[k] = self.definitions['parameters'][v]
            except (KeyError, TypeError):
                pass

names()

Names of all model definitions.

Returns:

Type Description
iterable of str

All model names in the target definitions.toml file.
Source code in src/echosms/referencemodels.py 
69
70
71
72
73
74
75
76
77
def names(self):
    """Names of all model definitions.

    Returns
    -------
    : iterable of str
        All model names in the ``target definitions.toml`` file.
    """
    return [n['name'] for n in self.definitions['target']]

parameters(name)

Model parameters for a particular model.

Model parameters are a subset of the model specification where the metadata items have been removed.

Parameters:

Name Type Description Default
name str

The name of a model in the target definitions.toml file.

 required

Returns:

Type Description
dict

The model parameters for the requested model or an empty set if no model with that name.
Source code in src/echosms/referencemodels.py 
 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
def parameters(self, name):
    """Model parameters for a particular model.

    Model parameters are a subset of the model specification where the metadata items have
    been removed.

    Parameters
    ----------
    name : str
        The name of a model in the ``target definitions.toml`` file.

    Returns
    -------
    : dict
        The model parameters for the requested model or an empty set if no model with that name.

    """
    s = self.specification(name)

    if not s:
        return []

    # Remove the entries that are not parameters
    p = s.copy()
    for k in ['name', 'shape', 'description', 'source', 'benchmark_model']:
        p.pop(k, None)
    return p

specification(name)

Model definitions for a particular model.

Parameters:

Name Type Description Default
name str

The name of a model in the target definitions.toml file.

 required

Returns:

Type Description
dict

The model definitions for the requested model or an empty set if no model with that name.
Source code in src/echosms/referencemodels.py 
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
def specification(self, name):
    """Model definitions for a particular model.

    Parameters
    ----------
    name : str
        The name of a model in the ``target definitions.toml`` file.

    Returns
    -------
    : dict
        The model definitions for the requested model or an empty set if no model
        with that name.
    """
    s = [t for t in self.definitions['target'] if t['name'] == name]
    if not s:
        return s

    return s[0]

echosms.JechEtAlData()

Simple access to all model results in Jech et al (2015).

Attributes:

Name Type Description
data dict[DataFrame]

One entry in the dict for each model results file.
data_directory Path

The directory containing the model results files.
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/jechetaldata.py 
25
26
27
def __init__(self):
    self.data_directory = Path(__file__).parent/Path('resources')/Path('Jechetal_allmodels')
    self.data = {f.stem: pd.read_csv(f) for f in self.data_directory.glob('*.csv')}

echosms.utils

Miscellaneous utility functions.

Neumann(m)

Neumann number.

Parameters:

Name Type Description Default
m int

The input integer.

 required

Returns:

Type Description
int

The Neumann number.
Source code in src/echosms/utils.py 
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
def Neumann(m: int) -> int:
    """Neumann number.

    Parameters
    ----------
    m :
        The input integer.

    Returns
    -------
    :
        The Neumann number.
    """
    if m == 0:
        return 1
    return 2

as_dataarray(params, no_expand=[])

Convert model parameters from dict form to a Xarray DataArray.

Parameters:

Name Type Description Default
params dict

The model parameters.

 required
no_expand list

Key values of the non-expandable model parameters in params.

 []

Returns:

Type Description
DataArray

Returns a multi-dimensional DataArray generated from the Cartesian product of all expandable items in the input dict. Non-expandable items are added to the DataArray attrs property. Expandable items are those that can be sensibly expanded into DataArray coordinates. Not all models have non-expandable items. The array is named ts, the values are initialised to nan, the dimension names are the dict keys, and the coordinate variables are the dict values.
Source code in src/echosms/utils.py 
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
def as_dataarray(params: dict, no_expand: list = []) -> xr.DataArray:
    """Convert model parameters from dict form to a Xarray DataArray.

    Parameters
    ----------
    params :
        The model parameters.

    no_expand :
        Key values of the non-expandable model parameters in `params`.

    Returns
    -------
    :
        Returns a multi-dimensional DataArray generated from the Cartesian product of all
        expandable items in the input dict. Non-expandable items are added to the DataArray
        attrs property. Expandable items are those that can be sensibly expanded into
        DataArray coordinates. Not all models have non-expandable items.
        The array is named `ts`, the values are initialised to `nan`, the
        dimension names are the dict keys, and the coordinate variables are the dict values.

    """
    expand, nexpand = split_dict(params, no_expand)

    # Convert scalars to iterables so xarray is happy
    for k, v in expand.items():
        if not isinstance(v, Iterable) or isinstance(v, str):
            expand[k] = [v]

    sz = [len(v) for k, v in expand.items()]
    return xr.DataArray(data=np.full(sz, np.nan), coords=expand, name='ts',
                        attrs={'units': 'dB', 'dB_reference': '1 m^2',
                               'parameters': nexpand})

as_dataframe(params, no_expand=[])

Convert model parameters from dict form to a Pandas DataFrame.

Parameters:

Name Type Description Default
params dict

The model parameters.

 required
no_expand list

Key values of the non-expandable model parameters in params.

 []

Returns:

Type Description
DataFrame

Returns a Pandas DataFrame generated from the Cartesian product of all expandable items in the input dict. DataFrame column names are obtained from the dict keys. Non-expandable items are added to the DataFrame attrs property. Expandable items are those that can be sensibly expanded into DataFrame columns. Not all models have non-expandable items.
Source code in src/echosms/utils.py 
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
def as_dataframe(params: dict, no_expand: list = []) -> pd.DataFrame:
    """Convert model parameters from dict form to a Pandas DataFrame.

    Parameters
    ----------
    params :
        The model parameters.

    no_expand :
        Key values of the non-expandable model parameters in `params`.

    Returns
    -------
    :
        Returns a Pandas DataFrame generated from the Cartesian product of all expandable
        items in the input dict. DataFrame column names are obtained from the dict keys.
        Non-expandable items are added to the DataFrame attrs property. Expandable items are
        those that can be sensibly expanded into DataFrame columns. Not all models have
        non-expandable items.

    """
    expand, nexpand = split_dict(params, no_expand)

    # Use meshgrid to do the Cartesian product then create a Pandas DataFrame from that, having
    # flattened the multidimensional arrays and using a dict to provide column names.
    # This preserves the differing dtypes in each column compared to other ways of
    # constructing the DataFrame).
    df = pd.DataFrame({k: t.flatten()
                       for k, t in zip(expand.keys(), np.meshgrid(*tuple(expand.values())))})
    df.attrs = {'parameters': nexpand}
    return df

as_dict(params)

Convert model parameters from DataFrame or DataArray to dict.

Parameters:

Name Type Description Default
params dict | DataFrame | DataArray

The model parameters

 required

Raises:

Type Description
TypeError:

If the input data type is not supported.

Returns:

Type Description
dict

A dict containing the model parameters.
Source code in src/echosms/utils.py 
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
def as_dict(params: dict | pd.DataFrame | xr.DataArray) -> dict:
    """Convert model parameters from DataFrame or DataArray to dict.

    Parameters
    ----------
    params:
        The model parameters

    Raises
    ------
    TypeError:
        If the input data type is not supported.

    Returns
    -------
    :
        A dict containing the model parameters.
    """
    if isinstance(params, dict):
        return params

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

    if isinstance(params, xr.DataArray):
        return dict(zip(params.coords, params.indexes.values())) | p
    elif isinstance(params, pd.DataFrame):
        # params.attrs = {}  # otherwise to_dict() exposes a bug in pandas?
        return params.to_dict(orient='series') | p

    raise TypeError('Only dict, DataFrame, or DataArray are accepted.')

h1(n, z, derivative=False)

Spherical Hankel function of the first kind or its' derivative.

Parameters:

Name Type Description Default
n int

Order (n ≥ 0).

 required
z float

Argument of the Hankel function.

 required
derivative

if True, the value of the derivative (rather than the function itself) is returned.

 False

Returns:

Type Description
complex

Value of the spherical Hankel function

Raises:

Type Description
ValueError

For negative n values.
Notes

The value of the Hankel function is calculated from spherical Bessel functions [1].

The derivative is computed from spherical Hankel functions [2].

References

[1] https://dlmf.nist.gov/10.47.E10

[2] https://dlmf.nist.gov/10.51.E2

Source code in src/echosms/utils.py 
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
def h1(n: int, z: float, derivative=False) -> complex:
    """Spherical Hankel function of the first kind or its' derivative.

    Parameters
    ----------
    n :
        Order (n ≥ 0).
    z :
        Argument of the Hankel function.
    derivative :
        if True, the value of the derivative (rather than the function itself) is returned.

    Returns
    -------
    :
        Value of the spherical Hankel function

    Raises
    ------
    ValueError
        For negative n values.

    Notes
    -----
    The value of the Hankel function is calculated from spherical Bessel functions [1].

    The derivative is computed from spherical Hankel functions [2].

    References
    ----------
    [1] <https://dlmf.nist.gov/10.47.E10>

    [2] <https://dlmf.nist.gov/10.51.E2>
    """
    if n < 0:
        raise ValueError('Negative n values are not supported for spherical Hankel functions.')

    if not derivative:
        return spherical_jn(n, z) + 1j*spherical_yn(n, z)
    return -h1(n+1, z) + (n/z) * h1(n, z)

pro_ang1(m, n, c, eta, norm=False)

Prolate spheroidal angular function of the first kind and derivative.

Calculates the prolate spheroidal angular function of the first kind and its' derivative with respect to eta.

Parameters:

Name Type Description Default
m int

The order parameter (≥ 0)

 required
n int

The degree parameter (≥ m).

 required
c float

The size parameter.

 required
eta float

The angular coordinate, η, where |η| ≤ 1.

 required
norm

If False, returned values are not normalised (i.e., the Meixner-Schäfke normlalisation scheme is used). For large m this norm becomes very large. If True the returned values are scaled by the square root of the normalisation of the corresponding Legendre function. This avoids the large values that occur when norm is False.

 False

Returns:

Type Description
tuple[float, float]

The value of the prolate spheroidal angular function and its' derivative.
Notes

This method uses the prolate spheroidal wave function code for non complex arguments (van Buren & Boisvert, 2002, and van Buren & Boisvert, 2024), available on github. This code is in Fortran90 and was interfaced to Python using numpy.f2py then wrapped with the current method to provide a similar calling convention as the Scipy function of the same name.

References

Van Buren, A. L., & Boisvert, J. E. (2002). Accurate calculation of prolate spheroidal radial functions of the first kind and their first derivatives. Quarterly of Applied Mathematics, 60(3), 589-599. https://doi.org/10.1090/qam/1914443

Van Buren, A. L., & Boisvert, J. E. (2004). Improved Calculation of Prolate Spheroidal Radial Functions of the Second Kind and Their First Derivatives. Quarterly of Applied Mathematics, 62(3), 493-507. https://doi.org/10.1090/qam/2086042

Source code in src/echosms/utils.py 
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
def pro_ang1(m: int, n: int, c: float, eta: float, norm=False) -> tuple[float, float]:
    """Prolate spheroidal angular function of the first kind and derivative.

    Calculates the prolate spheroidal angular function of the first kind and its'
    derivative with respect to `eta`.

    Parameters
    ----------
    m :
        The order parameter (≥ 0)
    n :
        The degree parameter (≥ `m`).
    c :
        The size parameter.
    eta :
        The angular coordinate, η, where |η| ≤ 1.
    norm :
        If `False`, returned values are not normalised (i.e., the Meixner-Schäfke normlalisation
        scheme is used). For large `m` this norm becomes very large. If `True` the returned values
        are scaled by the square root of the normalisation of the corresponding Legendre function.
        This avoids the large values that occur when `norm` is `False`.

    Returns
    -------
    :
        The value of the prolate spheroidal angular function and its' derivative.

    Notes
    -----
    This method uses the prolate spheroidal wave function code for non complex
    arguments (van Buren & Boisvert, 2002, and van Buren & Boisvert, 2024), available on
    [github](https://github.com/MathieuandSpheroidalWaveFunctions). This code is in Fortran90
    and was interfaced to Python using `numpy.f2py` then wrapped with the current method to
    provide a similar calling convention as the Scipy function of the same name.

    References
    ----------
    Van Buren, A. L., & Boisvert, J. E. (2002). Accurate calculation of prolate spheroidal
    radial functions of the first kind and their first derivatives. Quarterly of Applied
    Mathematics, 60(3), 589-599. <https://doi.org/10.1090/qam/1914443>

    Van Buren, A. L., & Boisvert, J. E. (2004). Improved Calculation of Prolate Spheroidal
    Radial Functions of the Second Kind and Their First Derivatives. Quarterly of Applied
    Mathematics, 62(3), 493-507. <https://doi.org/10.1090/qam/2086042>
    """
    if m < 0:
        raise ValueError('The m parameter must be positive.')
    if n < m:
        raise ValueError('The n parameter must be greater than or equal to the m parameter.')
    if abs(eta) > 1.0:
        raise ValueError('The eta parameter must be less than or equal to 1')

    a = prolate_swf.profcn(c=c, m=m, lnum=n-m+2, x1=0.0, ioprad=0, iopang=2,
                           iopnorm=int(norm), arg=[eta])
    p = swf_t._make(a)
    if np.isnan(p.s1c[n-m]) or np.isnan(p.s1dc[n-m]):
        # print('Root - trying again.')
        a = prolate_swf.profcn(c=c, m=m, lnum=n-m+2, x1=0.0, ioprad=0, iopang=2,
                               iopnorm=int(norm), arg=[eta+sys.float_info.epsilon])
        p = swf_t._make(a)

    s = p.s1c[n-m] * np.float_power(10.0, p.is1e[n-m])
    sp = p.s1dc[n-m] * np.float_power(10.0, p.is1de[n-m])

    return s[0], sp[0]

pro_rad1(m, n, c, xi)

Prolate spheroidal radial function of the first kind and derivative.

Calculates the prolate spheroidal radial function of the first kind and its' derivative with respect to xi.

Parameters:

Name Type Description Default
m int

The order parameter (≥ 0).

 required
n int

The degree parameter (≥ m).

 required
c float

The size parameter.

 required
xi float

The radial coordinate, ξ, where ξ ≥ 1.

 required

Returns:

Type Description
tuple[float, float]

The value of the prolate spheroidal radial function and its' derivative.
Notes

This method uses the prolate spheroidal wave function code for non complex arguments (van Buren & Boisvert, 2002, and van Buren & Boisvert, 2024), available on github. This code is in Fortran90 and was interfaced to Python using numpy.f2py then wrapped with the current method to provide a similar calling convention as the Scipy function of the same name.

References

Van Buren, A. L., & Boisvert, J. E. (2002). Accurate calculation of prolate spheroidal radial functions of the first kind and their first derivatives. Quarterly of Applied Mathematics, 60(3), 589-599. https://doi.org/10.1090/qam/1914443

Van Buren, A. L., & Boisvert, J. E. (2004). Improved Calculation of Prolate Spheroidal Radial Functions of the Second Kind and Their First Derivatives. Quarterly of Applied Mathematics, 62(3), 493-507. https://doi.org/10.1090/qam/2086042

Source code in src/echosms/utils.py 
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
def pro_rad1(m: int, n: int, c: float, xi: float) -> tuple[float, float]:
    """Prolate spheroidal radial function of the first kind and derivative.

    Calculates the prolate spheroidal radial function of the first kind and its'
    derivative with respect to `xi`.

    Parameters
    ----------
    m :
        The order parameter (≥ 0).
    n :
        The degree parameter (≥ `m`).
    c :
        The size parameter.
    xi :
        The radial coordinate, ξ, where ξ ≥ 1.

    Returns
    -------
    :
        The value of the prolate spheroidal radial function and its' derivative.

    Notes
    -----
    This method uses the prolate spheroidal wave function code for non complex
    arguments (van Buren & Boisvert, 2002, and van Buren & Boisvert, 2024), available on
    [github](https://github.com/MathieuandSpheroidalWaveFunctions). This code is in Fortran90
    and was interfaced to Python using `numpy.f2py` then wrapped with the current method to
    provide a similar calling convention as the Scipy function of the same name.

    References
    ----------
    Van Buren, A. L., & Boisvert, J. E. (2002). Accurate calculation of prolate spheroidal
    radial functions of the first kind and their first derivatives. Quarterly of Applied
    Mathematics, 60(3), 589-599. <https://doi.org/10.1090/qam/1914443>

    Van Buren, A. L., & Boisvert, J. E. (2004). Improved Calculation of Prolate Spheroidal
    Radial Functions of the Second Kind and Their First Derivatives. Quarterly of Applied
    Mathematics, 62(3), 493-507. <https://doi.org/10.1090/qam/2086042>
    """
    if m < 0:
        raise ValueError('The m parameter must be positive.')
    if n < m:
        raise ValueError('The n parameter must be greater than or equal to the m parameter.')
    if xi < 1.0:
        raise ValueError('The xi parameter must be greater than or equal to 1')

    a = prolate_swf.profcn(c=c, m=m, lnum=n-m+2, x1=xi-1.0, ioprad=1, iopang=0, iopnorm=0, arg=[0])
    p = swf_t._make(a)
    s = p.r1c * np.float_power(10.0, p.ir1e)
    sp = p.r1dc * np.float_power(10.0, p.ir1de)

    return s[n-m], sp[n-m]

pro_rad2(m, n, c, xi)

Prolate spheroidal radial function of the second kind and derivative.

Calculates the prolate spheroidal radial function of the second kind and its' derivative with respect to xi.

Parameters:

Name Type Description Default
m int

The order parameter (≥ 0).

 required
n int

The degree parameter (≥ m).

 required
c float

The size parameter.

 required
xi float

The radial coordinate, ξ, where ξ ≥ 1.

 required

Returns:

Type Description
tuple[float, float]

The value of the prolate spheroidal radial function and its' derivative.
Notes

This method uses the prolate spheroidal wave function code for non complex arguments (van Buren & Boisvert, 2002, and van Buren & Boisvert, 2024), available on github. This code is in Fortran90 and was interfaced to Python using numpy.f2py then wrapped with the current method to provide a similar calling convention as the Scipy function of the same name.

References

Van Buren, A. L., & Boisvert, J. E. (2002). Accurate calculation of prolate spheroidal radial functions of the first kind and their first derivatives. Quarterly of Applied Mathematics, 60(3), 589-599. https://doi.org/10.1090/qam/1914443

Van Buren, A. L., & Boisvert, J. E. (2004). Improved Calculation of Prolate Spheroidal Radial Functions of the Second Kind and Their First Derivatives. Quarterly of Applied Mathematics, 62(3), 493-507. https://doi.org/10.1090/qam/2086042

Source code in src/echosms/utils.py 
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
def pro_rad2(m: int, n: int, c: float, xi: float) -> tuple[float, float]:
    """Prolate spheroidal radial function of the second kind and derivative.

    Calculates the prolate spheroidal radial function of the second kind and its'
    derivative with respect to `xi`.

    Parameters
    ----------
    m :
        The order parameter (≥ 0).
    n :
        The degree parameter (≥ `m`).
    c :
        The size parameter.
    xi :
        The radial coordinate, ξ, where ξ ≥ 1.

    Returns
    -------
    :
        The value of the prolate spheroidal radial function and its' derivative.

    Notes
    -----
    This method uses the prolate spheroidal wave function code for non complex
    arguments (van Buren & Boisvert, 2002, and van Buren & Boisvert, 2024), available on
    [github](https://github.com/MathieuandSpheroidalWaveFunctions). This code is in Fortran90
    and was interfaced to Python using `numpy.f2py` then wrapped with the current method to
    provide a similar calling convention as the Scipy function of the same name.

    References
    ----------
    Van Buren, A. L., & Boisvert, J. E. (2002). Accurate calculation of prolate spheroidal
    radial functions of the first kind and their first derivatives. Quarterly of Applied
    Mathematics, 60(3), 589-599. <https://doi.org/10.1090/qam/1914443>

    Van Buren, A. L., & Boisvert, J. E. (2004). Improved Calculation of Prolate Spheroidal
    Radial Functions of the Second Kind and Their First Derivatives. Quarterly of Applied
    Mathematics, 62(3), 493-507. <https://doi.org/10.1090/qam/2086042>
    """
    if m < 0:
        raise ValueError('The m parameter must be positive.')
    if n < m:
        raise ValueError('The n parameter must be greater than or equal to the m parameter.')
    if xi < 1.0:
        raise ValueError('The xi parameter must be greater than or equal to 1')

    ioprad = 1 if xi-1.0 < 1e-10 else 2

    # Add +2 to lnum instead of +1 as it exposes a bug in the Fortran code - if n = 0, zeros
    # are returned instead of the correct value.
    a = prolate_swf.profcn(c=c, m=m, lnum=n-m+2, x1=xi-1.0,
                           ioprad=ioprad, iopang=0, iopnorm=0, arg=[0])
    p = swf_t._make(a)

    if ioprad == 1:
        s = np.inf
        sp = np.inf
    else:
        s = p.r2c * np.float_power(10.0, p.ir2e)
        sp = p.r2dc * np.float_power(10.0, p.ir2de)

    return s[n-m], sp[n-m]

spherical_jnpp(n, z)

Second derivative of the spherical Bessel function.

Parameters:

Name Type Description Default
n int

Order (n ≥ 0)

 required
z float

Argument of the Bessel function.

 required

Returns:

Type Description
float

The second derivative of the spherical Bessel function.
Source code in src/echosms/utils.py 
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
def spherical_jnpp(n: int, z: float) -> float:
    """Second derivative of the spherical Bessel function.

    Parameters
    ----------
    n :
        Order (n ≥ 0)
    z :
        Argument of the Bessel function.

    Returns
    -------
    :
        The second derivative of the spherical Bessel function.

    """
    return 1./z**2 * ((n**2-n-z**2)*spherical_jn(n, z) + 2.*z*spherical_jn(n+1, z))

split_dict(d, s)

Split a dict into two dicts based on a list of keys.

Parameters:

Name Type Description Default
d dict

Dict to be split.

 required
s list

List of dict keys to use for splitting d.

 required

Returns:

Type Description
tuple(dict, dict)

The input dict split into two dicts based on the keys in s. The first tuple item contains the items that do not have keys in s.
Source code in src/echosms/utils.py 
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
def split_dict(d: dict, s: list) -> tuple[dict, dict]:
    """Split a dict into two dicts based on a list of keys.

    Parameters
    ----------
    d : dict
        Dict to be split.

    s: list
        List of dict keys to use for splitting `d`.

    Returns
    -------
    : tuple(dict, dict)
        The `input` dict split into two dicts based on the keys in `s`. The first tuple item
        contains the items that do not have keys in `s`.
    """
    contains = {k: v for k, v in d.items() if k in s}
    ncontains = {k: v for k, v in d.items() if k not in s}
    return ncontains, contains

theoretical_Sa(ts, eba, r, nautical=False)

Theoretical area backscattering strength (Sₐ) of a target.

Parameters:

Name Type Description Default
ts float | ndarray

The target strength of the object [dB re 1 m²].

 required
eba float

The equivalent beam angle of the transducer [dB re 1 sr].

 required
r float

The range from the transducer to the target [m]. Used for acoustic beam spreading.

 required
nautical

If True, the nautical area scattering strength (SA) is returned instead of the area backscattering strength (Sₐ).

 False

Returns:

Type Description
float | ndarray

The theoretical Sₐ [dB re 1 m² m⁻²] or SA [dB re 1 m² nmi⁻²] of the input ts.

Raises:

Type Description
ValueError

If an input value is out of bounds.
Notes

While the calculation is valid for any target, the theoretical area strengths are most relevant when calibrating an echosounder using a sphere. The difference between the theoretical and measured can be used to calculate the calibration gain for an echosounder (when the sphere is on-axis).

Source code in src/echosms/utils.py 
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
def theoretical_Sa(ts: float | np.ndarray, eba: float, r: float, nautical=False)\
                   -> float | np.ndarray:
    """Theoretical area backscattering strength (Sₐ) of a target.

    Parameters
    ----------
    ts :
        The target strength of the object [dB re 1 m²].
    eba :
        The equivalent beam angle of the transducer [dB re 1 sr].
    r :
        The range from the transducer to the target [m]. Used for acoustic beam spreading.
    nautical :
        If `True`, the nautical area scattering strength (S<sub>A</sub>) is returned instead of the
        area backscattering strength (Sₐ).

    Returns
    -------
    :
        The theoretical Sₐ [dB re 1 m² m⁻²] or S<sub>A</sub> [dB re 1 m² nmi⁻²] of the input `ts`.

    Raises
    ------
    ValueError
        If an input value is out of bounds.

    Notes
    -----
    While the calculation is valid for any target, the theoretical area strengths are most relevant
    when calibrating an echosounder using a sphere. The difference between
    the theoretical and measured can be used to calculate the calibration gain for an
    echosounder (when the sphere is on-axis).
    """
    if eba > 0.0:
        raise ValueError('A positive eba value is not supported.')
    if r <= 0.0:
        raise ValueError('An r value less than or equal to 0 is not supported.')

    factor = 10.0*np.log10(4.0*np.pi*1852.0**2) if nautical else 0.0
    return ts - eba - 20.0*np.log10(r) + factor

wavelength(c, f)

Calculate the acoustic wavelength.

Parameters:

Name Type Description Default
c float

Sound speed [m/s]

 required
f float

Frequency [Hz]

 required

Returns:

Type Description
float

The acoustic wavelength [m].
Source code in src/echosms/utils.py 
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
def wavelength(c: float, f: float) -> float:
    """Calculate the acoustic wavelength.

    Parameters
    ----------
    c :
        Sound speed [m/s]

    f :
        Frequency [Hz]

    Returns
    -------
    :
        The acoustic wavelength [m].
    """
    return c/f

wavenumber(c, f)

Calculate the acoustic wavenumber.

Parameters:

Name Type Description Default
c float

Sound speed [m/s]

 required
f float

Frequency [Hz]

 required

Returns:

Type Description
float

The acoustic wavenumber [m⁻¹].
Source code in src/echosms/utils.py 
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
def wavenumber(c: float, f: float) -> float:
    """Calculate the acoustic wavenumber.

    Parameters
    ----------
    c :
        Sound speed [m/s]

    f :
        Frequency [Hz]

    Returns
    -------
    :
        The acoustic wavenumber [m⁻¹].
    """
    return 2*np.pi*f/c