<filename>tephrange/density.py<gh_stars>0
# -*- coding: utf-8 -*-
"""Functions for calculating the density of an ash particle based upon the
size and composition."""

import numpy as np


def bp2003(diameter, rho_pumice=440, rho_glass=2300):
    """
    Return the density of particles of given diameter (m) based on the model
    of Bonadonna and Phillips (2003), which assumes a linear decrease (on
    the phi scale) from lithic at 8 microns (6 phi) to pumice at 2
    millimeters (-1 phi).

    Example values are:

    Rhyolite (default):
        Eruption: Askja 1875 Pumice: 440 kg/m3, Lithic: 2300 kg/m3
    Dacite:
        Eruption: Quizapú 1932  Pumice: 655 kg/m3, Lithic: 2400 kg/m3
    Andesite:
        Eruption: Hudson 1991  Pumice: 1000 kg/m3, Lithic: 2600 kg/m3

    :param diameter: Numpy array of particle diameters
    :return density: Numpy array of particle density in kg/m3.
    """
    diam_lithic = 7
    diam_pumice = -1

    diam_phi = -np.log2(diameter * 1000)
    if diam_phi > diam_lithic:
        return rho_glass
    elif diam_phi < diam_pumice:
        return rho_pumice
    else:
        diff_diam = diam_pumice - diam_lithic
        diff_rho = abs(rho_pumice - rho_glass)
        density = rho_glass - (((diam_phi - diam_lithic) /
                                diff_diam) * diff_rho)
    return density


def convert_to_solidity(density, rho_glass=2300, rho_atm=1.225):
    """
    Calculate the solidity of a particle given the bulk density (kg/m3)
    and the lithic density.  Default rho_glass corresponds to Askja
    rhyolite glass.

    :density: np.array of density values
    :return solidity: Proportion of particle volume occupied by solid.
    """
    # density = X*rho_glass + (1-X)*rho_air
    # X = (density - rho_air) / (rho_glass - rho_air)
    result = np.array(
        [single_solidity(x, rho_glass=rho_glass, rho_atm=rho_atm)
         for x in density])

    return result


def single_solidity(density, rho_glass=2300, rho_atm=1.225):
    if density > rho_glass:
        raise ValueError('Particle density ({}) is more than glass '
                         'density ({})'.format(density, rho_glass))
    if density < rho_atm:
        raise ValueError('Particle density ({}) is less than '
                         'atmosphere density ({})'.format(density,
                                                          rho_atm))
    return (density - rho_atm) / (rho_glass - rho_atm)

<filename>tephrange/density.py<gh_stars>0
# -*- coding: utf-8 -*-
"""Functions for calculating the density of an ash particle based upon the
size and composition."""

import numpy as np


def bp2003(diameter, rho_pumice=440, rho_glass=2300):
    """
    Return the density of particles of given diameter (m) based on the model
    of Bonadonna and Phillips (2003), which assumes a linear decrease (on
    the phi scale) from lithic at 8 microns (6 phi) to pumice at 2
    millimeters (-1 phi).

    Example values are:

    Rhyolite (default):
        Eruption: Askja 1875 Pumice: 440 kg/m3, Lithic: 2300 kg/m3
    Dacite:
        Eruption: Quizapú 1932  Pumice: 655 kg/m3, Lithic: 2400 kg/m3
    Andesite:
        Eruption: Hudson 1991  Pumice: 1000 kg/m3, Lithic: 2600 kg/m3

    :param diameter: Numpy array of particle diameters
    :return density: Numpy array of particle density in kg/m3.
    """
    diam_lithic = 7
    diam_pumice = -1

    diam_phi = -np.log2(diameter * 1000)
    if diam_phi > diam_lithic:
        return rho_glass
    elif diam_phi < diam_pumice:
        return rho_pumice
    else:
        diff_diam = diam_pumice - diam_lithic
        diff_rho = abs(rho_pumice - rho_glass)
        density = rho_glass - (((diam_phi - diam_lithic) /
                                diff_diam) * diff_rho)
    return density


def convert_to_solidity(density, rho_glass=2300, rho_atm=1.225):
    """
    Calculate the solidity of a particle given the bulk density (kg/m3)
    and the lithic density.  Default rho_glass corresponds to Askja
    rhyolite glass.

    :density: np.array of density values
    :return solidity: Proportion of particle volume occupied by solid.
    """
    # density = X*rho_glass + (1-X)*rho_air
    # X = (density - rho_air) / (rho_glass - rho_air)
    result = np.array(
        [single_solidity(x, rho_glass=rho_glass, rho_atm=rho_atm)
         for x in density])

    return result


def single_solidity(density, rho_glass=2300, rho_atm=1.225):
    if density > rho_glass:
        raise ValueError('Particle density ({}) is more than glass '
                         'density ({})'.format(density, rho_glass))
    if density < rho_atm:
        raise ValueError('Particle density ({}) is less than '
                         'atmosphere density ({})'.format(density,
                                                          rho_atm))
    return (density - rho_atm) / (rho_glass - rho_atm)

# -*- coding: utf-8 -*- Functions for calculating the density of an ash particle based upon the size and composition. Return the density of particles of given diameter (m) based on the model     of Bonadonna and Phillips (2003), which assumes a linear decrease (on     the phi scale) from lithic at 8 microns (6 phi) to pumice at 2     millimeters (-1 phi).      Example values are:      Rhyolite (default):         Eruption: Askja 1875 Pumice: 440 kg/m3, Lithic: 2300 kg/m3     Dacite:         Eruption: Quizapú 1932  Pumice: 655 kg/m3, Lithic: 2400 kg/m3     Andesite:         Eruption: Hudson 1991  Pumice: 1000 kg/m3, Lithic: 2600 kg/m3      :param diameter: Numpy array of particle diameters     :return density: Numpy array of particle density in kg/m3. Calculate the solidity of a particle given the bulk density (kg/m3)     and the lithic density.  Default rho_glass corresponds to Askja     rhyolite glass.      :density: np.array of density values     :return solidity: Proportion of particle volume occupied by solid. # density = X*rho_glass + (1-X)*rho_air # X = (density - rho_air) / (rho_glass - rho_air)