Commit 82c05a96 authored by Giovanni La Mura's avatar Giovanni La Mura
Browse files

Define particle center as the average of sphere centers weighted on sphere radii

parent 038fe855

src/libnptm/Configuration.cpp

+11 −8
Original line number Diff line number Diff line
@@ -1150,16 +1150,19 @@ double ScattererConfiguration::get_particle_radius(GeometryConfiguration *gc) { 
  } else {

    double dist2, max_dist;

    double max_dist2 = 0.0;

    double avgX = 0.0, avgY = 0.0, avgZ = 0.0;

#pragma omp parallel for reduction(+: avgX, avgY, avgZ)

    double avgX = 0.0, avgY = 0.0, avgZ = 0.0, avg_radius = 0.0;

#pragma omp parallel for reduction(+: avgX, avgY, avgZ, avg_radius)

    for (int si = 0; si < _number_of_spheres; si++) {

      avgX += gc->get_sph_x(si);

      avgY += gc->get_sph_y(si);

      avgZ += gc->get_sph_z(si);

    }

    avgX /= _number_of_spheres;

    avgY /= _number_of_spheres;

    avgZ /= _number_of_spheres;

      double sph_radius = get_radius(si);

      avgX += (sph_radius * gc->get_sph_x(si));

      avgY += (sph_radius * gc->get_sph_y(si));

      avgZ += (sph_radius * gc->get_sph_z(si));

      avg_radius += sph_radius;

    }

    avg_radius /= _number_of_spheres;

    avgX /= (_number_of_spheres * avg_radius);

    avgY /= (_number_of_spheres * avg_radius);

    avgZ /= (_number_of_spheres * avg_radius);

    // int configuration = 0;

    for (int si = 0; si < _number_of_spheres; si++) {

      // int iogi = _iog_vec[si];

src/scripts/model_maker.py

+10 −7
Original line number Diff line number Diff line
@@ -628,6 +628,7 @@ def print_model_summary(scatterer, geometry): 
    avgX = 0.0

    avgY = 0.0

    avgZ = 0.0

    avgROS = 0.0

    Rmin = 0.0

    Rmax = 0.0

    Reqm = 0.0
@@ -635,20 +636,22 @@ def print_model_summary(scatterer, geometry): 
    Rcirc = 0.0

    square_farthest = 0.0

    for i in range(scatterer['nsph']):

        avgX += geometry['vec_sph_x'][i]

        avgY += geometry['vec_sph_y'][i]

        avgZ += geometry['vec_sph_z'][i]

        sph_type_index = scatterer['vec_types'][i] - 1

        ros = scatterer['ros'][sph_type_index]

        avgX += (ros * geometry['vec_sph_x'][i])

        avgY += (ros * geometry['vec_sph_y'][i])

        avgZ += (ros * geometry['vec_sph_z'][i])

        avgROS += ros

        sph_type_index = scatterer['vec_types'][i] - 1

        R3tot += math.pow(ros, 3.0)

        if (ros > Rmax):

            Rmax = ros

        if (Rmin == 0.0 or ros < Rmin):

            Rmin = ros

    Reqm = math.pow(R3tot, 1.0 / 3.0)

    avgX /= scatterer['nsph']

    avgY /= scatterer['nsph']

    avgZ /= scatterer['nsph']

    avgROS /= scatterer['nsph']

    avgX /= (avgROS * scatterer['nsph'])

    avgY /= (avgROS * scatterer['nsph'])

    avgZ /= (avgROS * scatterer['nsph'])

    for i in range(scatterer['nsph']):

        sph_type_index = scatterer['vec_types'][i] - 1

        ros = scatterer['ros'][sph_type_index]