Commit ea97feae authored by Giovanni La Mura's avatar Giovanni La Mura
Browse files

Implement inertial definition of filling factor

parent 7f1f0781
Loading
Loading
Loading
Loading
+14 −8
Original line number Diff line number Diff line
@@ -25,6 +25,7 @@

import math
import numpy as np
#import pdb

from sys import argv

@@ -82,9 +83,9 @@ def main():
## \brief Compute the particle's filling factor.
#
#  In this application the filling factor is defined as the ratio of the
#  volume of a homogeneous ellipsoid filled with as much mass as the particle
#  having the same principal moments of inertia of the particle and the sum
#  of the volumes of the spheres that compose the particle.
#  sum of the volumes of the spheres that compose the particle with respect
#  to the volume of a homogeneous ellipsoid that has the same rotational
#  properties as the particle.
#
#  \param stypes: `dict` Dictionary of sphere types.
#  \param i_tot_cm: `numpy.ndarray` Inertia tensor of the particle with respect
@@ -98,10 +99,14 @@ def filling_factor(stypes, i_tot_cm, config):
    eigenval, eigenvec = np.linalg.eigh(i_tot_cm)
    print("INFO: ellipsoid axes:", eigenval)
    print("INFO: ellipsoid directions:", eigenvec)
    eigenval /= mtot
    gf = 5.0 / (2.0 * mtot)
    I1, I2, I3 = eigenval
    a = np.sqrt(max(0.0, gf * (I2 + I3 - I1)))
    b = np.sqrt(max(0.0, gf * (I1 + I3 - I2)))
    c = np.sqrt(max(0.0, gf * (I1 + I2 - I3)))
    # Volume of the ellipsoid that has the same moments of inertia as the particle
    vell = 4.0 * math.pi * math.sqrt(eigenval[0] * eigenval[1] * eigenval[2]) / 3.0
    ff = vell / vtot
    vell = 4.0 * math.pi * a * b * c / 3.0
    ff = vtot / vell
    return ff

## \brief Create a dictionary of sphere types.
@@ -355,6 +360,7 @@ def get_types(config):
            cur_type = int(split_line[ti])
            if (cur_type > read_spheres + ti):
                type_id += 1
            if (len(stypes['vec_types']) < config['nsph']):
                stypes['vec_types'].append(type_id)
        read_spheres += len(split_line)
    stypes['num_types'] = type_id