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

#ifndef RADIX_RADIXMATH_UTIL_H
#define RADIX_RADIXMATH_UTIL_H

#include "radixbug/bug.hh"
#include "radixmath/constants.hh"

#include "radixcore/visibility.hh"

#include <vector>

namespace radix
{
/**
 * @brief hypsometric Calculates the thickness in meters between p1 and p2.
 * @param p1 z1 pressure kPa
 * @param p2 z2 pressure kPa
 * @param meanTemp mean temperature between z1 and z2.
 * @return thickness in meters
 */
template <typename type>
type RADIX_PUBLIC hypsometric(type p1, type p2, type meanTemp)
{
  type result = (SPECIFIC_GAS_CONSTANT * meanTemp) /
                GRAVITATIONAL_ACCELERATION * std::log(p1 / p2);
  // divide out the kg/mol using the molar mass of dry air and water
  result = result / MOLAR_MASS_DRY_AIR;
  return result;
}

/**
 * @brief hypsometricProfile Calculate thickness' for an entire profile
 * @param pressure std::vector<Real> pressures in kPa or millibars
 * @param temperature std::vector<Real>temperatures in Kelvin
 * @param sum whether to sum each thickness so it is the thickness between 0 &
 * the level.
 * @return std::vector<Real> thickness in meters
 */
template <typename type>
std::vector<type> RADIX_PUBLIC
hypsometricProfile(const std::vector<type> &pressure,
                   const std::vector<type> &temperature, bool sum = true)
{
  if (pressure.size() == 0) return std::vector<type>();
  if (pressure.size() != temperature.size())
    throw std::runtime_error(
        "hypsometricProfile requires pressure and temperature be the same "
        "size.");
  std::vector<type> altitude(pressure.size(), 0);
  type p1 = pressure[0];
  type t1 = temperature[0];
  for (size_t pi = 0; pi < pressure.size(); ++pi)
  {
    if (pi != 0)
    {
      p1 = pressure[pi - 1];
      t1 = (temperature[pi] + temperature[pi - 1]) / 2.0;
    }

    altitude[pi] = hypsometric(p1, pressure[pi], t1);
  }
  if (sum)
  {
    for (size_t ai = 1; ai < altitude.size(); ++ai)
    {
      altitude[ai] = altitude[ai] + altitude[ai - 1];
    }
  }
  return altitude;
}

/**
 * @brief absoluteTemperature Get the temperature of an air parcel given its
 * potential temperature and pressure
 * @param potentialTemp Potential temperature in Kelvin
 * @param pressure Pressure in mb
 * @return Absolute temperature (K)
 */
Real RADIX_PUBLIC absoluteTemperature(Real potentialTemp, Real pressure);

/**
 * @brief saturationVaporPressure Calculate the saturation pressure of an air
 * parcel given its absolute temperature and pressure
 * @param absTemp Absolute temperature in Kelvin
 * @param pressure Pressure in mb
 * @return Saturation vapour pressure (mb)
 */
Real RADIX_PUBLIC saturationVaporPressure(Real absTemp, Real pressure);

/**
 * @brief mixingRatioToRelativeHumidity Get the relative humidity of an air
 * parcel given a mixing ratio, temperature and pressure
 * @param mixingRatio Mixing ratio (g/g)
 * @param absTemp Absolute temperature (K)
 * @param pressure Pressure (mb)
 * @return Relative humidity (%)
 */
Real RADIX_PUBLIC mixingRatioToRelativeHumidity(Real mixingRatio, Real absTemp,
                                                Real pressure);

/**
 * @brief specificHumidityToRelativeHumidity Get the relative humidity of an air
 * parcel given a specific humidity, temperature and pressure
 * @param specificHumidity Specific humidity (g/kg)
 * @param absTemp Absolute temperature (K)
 * @param pressure Pressure (mb)
 * @return Relative humidity (%)
 */
Real RADIX_PUBLIC specificHumidityToRelativeHumidity(Real specificHumidity,
                                                     Real absTemp,
                                                     Real pressure);

/**
 * @brief cspanf returns a value in the interval (begin,end]
 * This function is used to reduce periodic variables to a standard range.
 * It adjusts for the behaviour of the mod function which provides positive
 * results for position input and negative results for negative input.
 * @param value number to be reduced to the span
 * @param begin first value of the span
 * @param end last value of the span
 * @return Real the reduced value
 */
Real RADIX_PUBLIC cspanf(Real value, Real begin, Real end);

/**
 * @brief haversine The Haversine formula for calculating great-circle
 * @param lat1
 * @param lon1
 * @param lat2
 * @param lon2
 * @return angle(radians) between two points
 */
Real RADIX_PUBLIC haversine(Real lat1, Real lon1, Real lat2, Real lon2);

/**
 * @brief greatCircleDistance Calculates the shortest distance between two
 * points on the surface of a sphere
 * @param lat1 Real point1 latitude
 * @param lon1 Real point1 longitude
 * @param lat2 Real point2 latitude
 * @param lon2 Real point2 longitude
 * @param radius
 * @return
 */
Real RADIX_PUBLIC greatCircleDistance(Real lat1, Real lon1, Real lat2,
                                      Real lon2, Real radius);

/**
 * @brief greatCircleVolume Calculates the volume giving latitude, longitude and
 * radius bounds.
 * @param lat1 Real min latitude
 * @param lon1 Real min longitude
 * @param lat2 Real max latitude
 * @param lon2 Real max longitude
 * @param r1 Real min radius
 * @param r2 Real max radius
 * @return Volume in cubic meters
 */
Real RADIX_PUBLIC greatCircleVolume(Real lat1, Real lon1, Real lat2, Real lon2,
                                    Real r1, Real r2);

/**
 * @brief greatCircleArea Calculates the area (m^2) given rectanglar lat/lon
 * corners at radius height
 * @param lat1
 * @param lon1
 * @param lat2
 * @param lon2
 * @param r1
 * @return Area in square meters
 */
// Real greatCircleArea(Real lat1, Real lon1, Real lat2, Real lon2, Real r1);

/**
 * @brief cylinderVolume Calculates the volume of a cylinder
 * @param r radius of cylinder
 * @param h height of cylinder
 * @return Real
 */
Real RADIX_PUBLIC cylinderVolume(Real r, Real h);

std::vector<Real> RADIX_PUBLIC
getEqualAreaRadii(const std::vector<Real> &radii,
                  const std::vector<unsigned short> &subrings);

inline Real RADIX_PUBLIC toDegrees(Real radians) { return radians * 180 / PI; }

inline Real RADIX_PUBLIC toRadians(Real degrees) { return degrees * PI / 180; }

// 360 degrees = 1 turn, and limit to being less than a full turn
inline Real RADIX_PUBLIC toTurns(Real degrees)
{
  // Convert to revolutions
  double turns = degrees / 360;
  // Remove extra turns, rotating into (-1, 1)
  turns -= std::floor(turns);
  // Convert negative turns to positive turns
  if (turns < 0) turns += 1.0;
  radix_ensure(turns >= 0 && turns < 1);
  return turns;
}

/**
 * @brief lineIntersect determines the intersection distance from p1 along v1 to
 * line described p2->p3
 * @param p - the origin of the line
 * @param v - the direction vector to determine intersection with.
 * @param sp - a start point on the line to intersect
 * @param ep - a end point on the line to intersect
 * @return  Real - the track length t from p along v to the line described by
 * sp,ep t == maxReal if no intersection occurs
 */
Real RADIX_PUBLIC lineIntersect(const class Point3D &p, const class Vector3D &v,
                                const class Point3D &sp,
                                const class Point3D &ep);

/**
 * @brief interpolate interpolates the given data vectors
 * @param xvals - a list of bin boundaries in the x-dimension
 * @param yvals - a list of y-dimension values at each x-value
 * @param x - an x-value at which to interpolate its associated y-value
 * @param linear - whether to interpolate in a linear fashion (false for log)
 */
template <typename T>
T RADIX_PUBLIC interpolate(const std::vector<T> &xvals,
                           const std::vector<T> &yvals, T x,
                           bool linear = false);

/**
 * @brief gammaRayAbsorptionInAir calculates the absorption of gamma rays as
 * they pass through air
 * @param energy - energy of the gamma ray in MeV
 * @param linear - whether to interpolate across energies in a linear fashion
 * (false for log) return double - the absorption of gamma rays of a given
 * energy as they pass through air
 */
double RADIX_PUBLIC gammaRayAbsorptionInAir(double energy, bool linear = false);

/**
 * @brief gammaRayAttenuationInAir calculates the attenuation of gamma rays as
 * they pass through air
 * @param energy - energy of the gamma ray in MeV
 * @param linear - whether to interpolate across energies in a linear fashion
 * (false for log) return double - the attentuation of gamma rays of a given
 * energy as they pass through air
 */
double RADIX_PUBLIC gammaRayAttenuationInAir(double energy,
                                             bool linear = false);

/**
 * @brief exponentialIntegral calculates the exponential integral of the given
 * value
 * @param d - value of which the exponential integral is calculated
 * @return double - the exponential integral of the given value
 */
double RADIX_PUBLIC exponentialIntegral(double d);
}  // namespace radix

#include "radixmath/util.i.hh"

#endif /* RADIX_RADIXMATH_UTIL_H */