util.cc

/*
 * File:   util.cc
 * Author: a16
 *
 * Created on January 17, 2013, 9:27 AM
 */

#include "radixmath/util.hh"
#include <algorithm>
#include <cstddef>
#include "radixmath/point3d.hh"
#include "radixmath/vector3d.hh"

namespace radix
{
std::vector<Real> getEqualAreaRadii(const std::vector<Real> &radii,
                                    const std::vector<unsigned short> &subrings)
{
  std::vector<Real> newRadii;
  Real prevRadius;

  prevRadius = 0;

  for (size_t index = 0; index < radii.size(); index++)
  {
    Real _prevRadius;
    int _subrings;
    Real area;
    Real radius;
    Real subarea;

    radius      = radii[index];
    _prevRadius = prevRadius;
    _subrings   = subrings[index];

    area    = PI * radius * radius - PI * prevRadius * prevRadius;
    subarea = area / _subrings;

    for (int count = 1; count < _subrings; count++)
    {
      Real _radius;

      _radius     = std::sqrt((subarea + PI * _prevRadius * _prevRadius) / PI);
      _prevRadius = _radius;

      newRadii.push_back(_radius);
    }

    newRadii.push_back(radius);

    prevRadius = radius;
  }

  return newRadii;
}

Real lineIntersect(const Point3D &p, const Vector3D &v, const Point3D &sp,
                   const Point3D &ep)
{
  Real t = maxReal;
  Point3D p2(p + v);  // second point for line intersect
  Real x1 = p.x;
  Real y1 = p.y;
  Real x2 = p2.x;
  Real y2 = p2.y;

  Real x3 = sp.x;
  Real y3 = sp.y;
  Real x4 = ep.x;
  Real y4 = ep.y;

  Real denominator = (y1 - y2) * (x3 - x4) - (x1 - x2) * (y3 - y4);

  if (isWithin(0.0, denominator)) return t;

  Real n1 = (x1 - x2) * (x3 * y4 - y3 * x4) -
            (x1 * y2 - y1 * x2) * (x3 - x4);  // px numerator
  Real n2 = (y1 - y2) * (x3 * y4 - y3 * x4) -
            (x1 * y2 - y1 * x2) * (y3 - y4);  // py numerator

  Real px = n1 / denominator;
  Real py = n2 / denominator;

  Real tx     = px - p.x;
  Real ty     = py - p.y;
  Real xymag  = sqrt(v.x * v.x + v.y * v.y);
  Real zratio = xymag == 0 ? 0 : v.z / xymag;

  t       = std::sqrt(tx * tx + ty * ty);
  Real tz = t * zratio;
  t       = std::sqrt(tx * tx + ty * ty + tz * tz);
  t *= denominator < 0 ? -1 : 1;
  return t;
}

Real haversine(Real lat1, Real lon1, Real lat2, Real lon2)
{
  // Haversine formula
  Real dlat           = toRadians(lat2 - lat1);
  Real dlon           = toRadians(lon2 - lon1);
  Real haversine_dlat = std::sin(dlat / 2.0);
  haversine_dlat *= haversine_dlat;
  Real haversine_dlon = std::sin(dlon / 2.0);
  haversine_dlon *= haversine_dlon;
  Real delta_sigma = haversine_dlat + std::cos(toRadians(lat1)) *
                                          std::cos(toRadians(lat2)) *
                                          haversine_dlon;
  delta_sigma = 2.0 * std::asin(std::sqrt(delta_sigma));
  return delta_sigma;
}

Real greatCircleDistance(Real lat1, Real lon1, Real lat2, Real lon2,
                         Real radius)
{
  return radius * haversine(lat1, lon1, lat2, lon2);
}

Real greatCircleVolume(Real lat1, Real lon1, Real lat2, Real lon2, Real r1,
                       Real r2)
{
  Real r11 = r1, r22 = r2;
  r1 = std::min(r11, r22);
  r2 = std::max(r11, r22);

  double northa = radix::greatCircleDistance(lat1, lon1, lat1, lon2, r1);
  double southa = radix::greatCircleDistance(lat2, lon1, lat2, lon2, r1);
  double westa  = radix::greatCircleDistance(lat1, lon1, lat2, lon1, r1);

  double northb = radix::greatCircleDistance(lat1, lon1, lat1, lon2, r2);
  double southb = radix::greatCircleDistance(lat2, lon1, lat2, lon2, r2);
  double westb  = radix::greatCircleDistance(lat1, lon1, lat2, lon1, r2);
  //
  // approximate with a square
  double north  = (northa + northb + southa + southb) / 4.0;
  double west   = (westa + westb) / 2.0;
  double height = r2 - r1;

  return (north * west * height);
}

Real cylinderVolume(Real r, Real h) { return (PI * std::pow(r, 2) * h); }

Real cspanf(Real value, Real begin, Real end)
{
  Real first = 0.0, last = 0.0;
  first = std::min(begin, end);
  last  = std::max(begin, end);
  value = std::fmod(value - first, last - first);
  if (value <= 0)
  {
    return value + last;
  }
  else
  {
    return value + first;
  }
}

Real hpaToAltitude(Real hpa)
{
  return 0.3048 * (1 - std::pow(hpa / 1013.25, 0.190284)) * 145366.45;
}

/** Real greatCircleArea(Real lat1, Real lon1, Real lat2, Real lon2, Real r1)
{

}*/

double gammaRayAbsorptionInAir(double energy, bool linear)
{
  // supported energies in MeV
  static const std::vector<double> energies = {
      0.001, 0.0015, 0.002, 0.003, 0.004, 0.005, 0.006, 0.008, 0.01,
      0.015, 0.02,   0.03,  0.04,  0.05,  0.06,  0.08,  0.1,   0.15,
      0.2,   0.3,    0.4,   0.5,   0.6,   0.8,   1,     1.25,  1.5,
      2,     3,      4,     5,     6,     8,     10,    15,    20};

  // absorption coefficients
  const std::vector<double> absorption = {
      3599,    1188,    526.2,   161.4,   76.36,   39.31,   22.7,    9.446,
      4.742,   1.334,   0.5389,  0.1537,  0.6833,  0.04098, 0.03041, 0.02407,
      0.02325, 0.02496, 0.02672, 0.02872, 0.02949, 0.02966, 0.02953, 0.02882,
      0.02789, 0.02666, 0.02547, 0.02345, 0.02057, 0.0187,  0.0174,  0.01647,
      0.01525, 0.0145,  0.01353, 0.01311};

  // interpolate to determine the absorption for the given energy
  return interpolate(energies, absorption, energy, linear);
}

double gammaRayAttenuationInAir(double energy, bool linear)
{
  // supported energies in MeV
  static const std::vector<double> energies = {
      0.001, 0.0015, 0.002, 0.003, 0.004, 0.005, 0.006, 0.008, 0.01,
      0.015, 0.02,   0.03,  0.04,  0.05,  0.06,  0.08,  0.1,   0.15,
      0.2,   0.3,    0.4,   0.5,   0.6,   0.8,   1,     1.25,  1.5,
      2,     3,      4,     5,     6,     8,     10,    15,    20};

  // attenuation coefficients
  static const std::vector<double> attenuation = {
      3606,    1191,    527.9,   162.5,   77.88,   40.27,   23.41,   9.921,
      5.12,    1.614,   0.7779,  0.3538,  0.2485,  0.208,   0.1875,  0.1662,
      0.1541,  0.1356,  0.1233,  0.1067,  0.09549, 0.08712, 0.08055, 0.07074,
      0.06358, 0.05687, 0.05175, 0.04447, 0.03581, 0.03079, 0.02751, 0.02522,
      0.02225, 0.02045, 0.0181,  0.01705};

  // interpolate to determine the attenuation for the given energy
  return interpolate(energies, attenuation, energy, linear);
}

double exponentialIntegral(double d)
{
  // based on polynomial approximation in document provided by Vince
  // specialization for d <= 1
  if (d <= 1)
  {
    return 0.00107857 * (d * d * d * d * d) - 0.00976004 * (d * d * d * d) +
           0.05519968 * (d * d * d) - 0.24991055 * (d * d) + 0.99999193 * d -
           0.57721566 - std::log(d);
  }

  // specialization for d > 1
  auto numer = 1 * (d * d * d * d) + 8.5733287401 * (d * d * d) +
               18.0590169730 * (d * d) + 8.6347608925 * d + 0.2677737343;
  auto denom = 1 * (d * d * d * d) + 9.5733223454 * (d * d * d) +
               25.6329561486 * (d * d) + 21.0996530827 * d + 3.9584969228;

  return numer / denom / (d * std::exp(d));
}
}  // namespace radix