util.hh 8.83 KB
 LEFEBVREJP email committed Jun 23, 2016 1 2 3 4 5 6 7 8 ``````/* * File: Util.h * Author: a16 * * Created on January 17, 2013, 9:27 AM */ #ifndef RADIX_RADIXMATH_UTIL_H `````` LEFEBVREJP email committed Jan 14, 2017 9 ``````#define RADIX_RADIXMATH_UTIL_H `````` LEFEBVREJP email committed Jun 23, 2016 10 11 12 13 `````` #include "radixbug/bug.hh" #include "radixmath/constants.hh" `````` LEFEBVREJP email committed Mar 27, 2018 14 ``````#include "radixcore/visibility.hh" `````` LEFEBVREJP email committed Nov 02, 2017 15 `````` `````` LEFEBVREJP email committed Jun 23, 2016 16 17 ``````#include `````` 18 19 ``````namespace radix { `````` LEFEBVREJP email committed Mar 19, 2018 20 ``````/** `````` LEFEBVREJP email committed Mar 20, 2018 21 22 23 24 25 `````` * @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 `````` LEFEBVREJP email committed Mar 19, 2018 26 `````` */ `````` LEFEBVREJP email committed Mar 20, 2018 27 28 29 30 31 32 33 34 35 ``````template 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; } `````` LEFEBVREJP email committed Mar 19, 2018 36 37 `````` /** `````` LEFEBVREJP email committed Mar 20, 2018 38 39 40 41 42 43 `````` * @brief hypsometricProfile Calculate thickness' for an entire profile * @param pressure std::vector pressures in kPa or millibars * @param temperature std::vectortemperatures in Kelvin * @param sum whether to sum each thickness so it is the thickness between 0 & * the level. * @return std::vector thickness in meters `````` LEFEBVREJP email committed Mar 19, 2018 44 `````` */ `````` LEFEBVREJP email committed Mar 20, 2018 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 ``````template std::vector RADIX_PUBLIC hypsometricProfile(const std::vector &pressure, const std::vector &temperature, bool sum = true) { if (pressure.size() == 0) return std::vector(); if (pressure.size() != temperature.size()) throw std::runtime_error( "hypsometricProfile requires pressure and temperature be the same " "size."); std::vector 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; } `````` LEFEBVREJP email committed Mar 19, 2018 65 `````` `````` LEFEBVREJP email committed Mar 20, 2018 66 67 68 69 70 71 72 73 74 75 76 `````` 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; } `````` Purves, Murray committed Mar 11, 2018 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 `````` /** * @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); `````` LEFEBVREJP email committed Jul 13, 2016 119 120 121 ``````/** * @brief cspanf returns a value in the interval (begin,end] * This function is used to reduce periodic variables to a standard range. `````` LEFEBVREJP email committed Dec 30, 2017 122 123 `````` * It adjusts for the behaviour of the mod function which provides positive * results for position input and negative results for negative input. `````` LEFEBVREJP email committed Jul 13, 2016 124 125 126 127 128 `````` * @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 */ `````` LEFEBVREJP email committed Nov 02, 2017 129 ``````Real RADIX_PUBLIC cspanf(Real value, Real begin, Real end); `````` LEFEBVREJP email committed Jul 13, 2016 130 `````` `````` LEFEBVREJP email committed Jun 25, 2016 131 132 133 134 135 136 137 138 ``````/** * @brief haversine The Haversine formula for calculating great-circle * @param lat1 * @param lon1 * @param lat2 * @param lon2 * @return angle(radians) between two points */ `````` LEFEBVREJP email committed Nov 02, 2017 139 ``````Real RADIX_PUBLIC haversine(Real lat1, Real lon1, Real lat2, Real lon2); `````` LEFEBVREJP email committed Jun 25, 2016 140 141 `````` /** `````` LEFEBVREJP email committed Dec 30, 2017 142 143 `````` * @brief greatCircleDistance Calculates the shortest distance between two * points on the surface of a sphere `````` LEFEBVREJP email committed Jul 01, 2016 144 145 146 147 `````` * @param lat1 Real point1 latitude * @param lon1 Real point1 longitude * @param lat2 Real point2 latitude * @param lon2 Real point2 longitude `````` LEFEBVREJP email committed Jun 25, 2016 148 149 150 `````` * @param radius * @return */ `````` LEFEBVREJP email committed Dec 30, 2017 151 152 ``````Real RADIX_PUBLIC greatCircleDistance(Real lat1, Real lon1, Real lat2, Real lon2, Real radius); `````` LEFEBVREJP email committed Jun 25, 2016 153 `````` `````` LEFEBVREJP email committed Jul 01, 2016 154 ``````/** `````` LEFEBVREJP email committed Dec 30, 2017 155 156 `````` * @brief greatCircleVolume Calculates the volume giving latitude, longitude and * radius bounds. `````` LEFEBVREJP email committed Jul 01, 2016 157 158 159 160 161 162 163 164 `````` * @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 */ `````` LEFEBVREJP email committed Dec 30, 2017 165 166 ``````Real RADIX_PUBLIC greatCircleVolume(Real lat1, Real lon1, Real lat2, Real lon2, Real r1, Real r2); `````` LEFEBVREJP email committed Jul 01, 2016 167 `````` `````` LEFEBVREJP email committed Nov 07, 2016 168 ``````/** `````` LEFEBVREJP email committed Dec 30, 2017 169 170 `````` * @brief greatCircleArea Calculates the area (m^2) given rectanglar lat/lon * corners at radius height `````` LEFEBVREJP email committed Nov 07, 2016 171 172 173 174 175 176 177 `````` * @param lat1 * @param lon1 * @param lat2 * @param lon2 * @param r1 * @return Area in square meters */ `````` LEFEBVREJP email committed Dec 30, 2017 178 ``````// Real greatCircleArea(Real lat1, Real lon1, Real lat2, Real lon2, Real r1); `````` LEFEBVREJP email committed Nov 07, 2016 179 `````` `````` LEFEBVREJP email committed Jul 01, 2016 180 181 182 183 184 185 ``````/** * @brief cylinderVolume Calculates the volume of a cylinder * @param r radius of cylinder * @param h height of cylinder * @return Real */ `````` LEFEBVREJP email committed Nov 02, 2017 186 ``````Real RADIX_PUBLIC cylinderVolume(Real r, Real h); `````` LEFEBVREJP email committed Jul 01, 2016 187 `````` `````` LEFEBVREJP email committed Dec 30, 2017 188 189 190 ``````std::vector RADIX_PUBLIC getEqualAreaRadii(const std::vector &radii, const std::vector &subrings); `````` LEFEBVREJP email committed Jun 23, 2016 191 `````` `````` LEFEBVREJP email committed Dec 30, 2017 192 ``````inline Real RADIX_PUBLIC toDegrees(Real radians) { return radians * 180 / PI; } `````` LEFEBVREJP email committed Jun 23, 2016 193 `````` `````` LEFEBVREJP email committed Dec 30, 2017 194 ``````inline Real RADIX_PUBLIC toRadians(Real degrees) { return degrees * PI / 180; } `````` LEFEBVREJP email committed Jun 23, 2016 195 196 `````` // 360 degrees = 1 turn, and limit to being less than a full turn `````` 197 198 ``````inline Real RADIX_PUBLIC toTurns(Real degrees) { `````` LEFEBVREJP email committed Dec 30, 2017 199 200 201 202 203 204 205 206 `````` // 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; `````` LEFEBVREJP email committed Jun 23, 2016 207 208 209 ``````} /** `````` LEFEBVREJP email committed Dec 30, 2017 210 211 `````` * @brief lineIntersect determines the intersection distance from p1 along v1 to * line described p2->p3 `````` LEFEBVREJP email committed Jun 23, 2016 212 213 214 215 `````` * @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 `````` LEFEBVREJP email committed Dec 30, 2017 216 217 `````` * @return Real - the track length t from p along v to the line described by * sp,ep t == maxReal if no intersection occurs `````` LEFEBVREJP email committed Jun 23, 2016 218 `````` */ `````` LEFEBVREJP email committed Dec 30, 2017 219 220 221 ``````Real RADIX_PUBLIC lineIntersect(const class Point3D &p, const class Vector3D &v, const class Point3D &sp, const class Point3D &ep); `````` 222 223 224 225 226 227 228 229 `````` /** * @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) */ `````` LEFEBVREJP email committed Dec 30, 2017 230 231 232 233 ``````template T RADIX_PUBLIC interpolate(const std::vector &xvals, const std::vector &yvals, T x, bool linear = false); `````` 234 235 `````` /** `````` LEFEBVREJP email committed Dec 30, 2017 236 237 `````` * @brief gammaRayAbsorptionInAir calculates the absorption of gamma rays as * they pass through air `````` 238 `````` * @param energy - energy of the gamma ray in MeV `````` LEFEBVREJP email committed Dec 30, 2017 239 240 241 `````` * @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 `````` 242 `````` */ `````` LEFEBVREJP email committed Nov 02, 2017 243 ``````double RADIX_PUBLIC gammaRayAbsorptionInAir(double energy, bool linear = false); `````` 244 245 `````` /** `````` LEFEBVREJP email committed Dec 30, 2017 246 247 `````` * @brief gammaRayAttenuationInAir calculates the attenuation of gamma rays as * they pass through air `````` 248 `````` * @param energy - energy of the gamma ray in MeV `````` LEFEBVREJP email committed Dec 30, 2017 249 250 251 `````` * @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 `````` 252 `````` */ `````` LEFEBVREJP email committed Dec 30, 2017 253 254 ``````double RADIX_PUBLIC gammaRayAttenuationInAir(double energy, bool linear = false); `````` Thompson, Adam committed Aug 24, 2016 255 256 `````` /** `````` LEFEBVREJP email committed Dec 30, 2017 257 258 `````` * @brief exponentialIntegral calculates the exponential integral of the given * value `````` Thompson, Adam committed Aug 24, 2016 259 260 261 `````` * @param d - value of which the exponential integral is calculated * @return double - the exponential integral of the given value */ `````` LEFEBVREJP email committed Nov 02, 2017 262 ``````double RADIX_PUBLIC exponentialIntegral(double d); `````` LEFEBVREJP email committed Dec 30, 2017 263 ``````} // namespace radix `````` 264 265 `````` #include "radixmath/util.i.hh" `````` LEFEBVREJP email committed Jun 23, 2016 266 `````` `````` LEFEBVREJP email committed Dec 30, 2017 267 ``#endif /* RADIX_RADIXMATH_UTIL_H */``