util.cc

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

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

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 interpolate(const std::vector<double> &xvals, const std::vector<double> &yvals, double x, bool linear)
{
    // iterator to upper bound
    auto ub = std::upper_bound(xvals.begin(), xvals.end(), x);

    // below range, return first item
    if(ub == xvals.begin())
    {
        return yvals.front();
    }
    // above range, return last item
    else if(ub == xvals.end())
    {
        return yvals.back();
    }

    // indices to right/left side
    auto r = std::distance(xvals.begin(), ub);
    auto l = r - 1;

    // min/max vals
    auto xmin = xvals[l];
    auto xmax = xvals[r];
    auto ymin = yvals[l];
    auto ymax = yvals[r];

    // adjust for log interpolation
    if(!linear)
    {
        xmin = std::log(xmin);
        xmax = std::log(xmax);
        ymin = std::log(ymin);
        ymax = std::log(ymax);
        x = std::log(x);
    }

    // scaling ratio
    auto ratio = (x - xmin) / (xmax - xmin);

    // interpolation result
    auto result = yvals[l] + (ymax - ymin) * ratio;
    return result;
}

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