util.hh 8.83 KB
Newer Older
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
9
#define RADIX_RADIXMATH_UTIL_H
10
11
12
13

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

14
#include "radixcore/visibility.hh"
15

16
17
#include <vector>

18
19
namespace radix
{
20
/**
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
26
 */
27
28
29
30
31
32
33
34
35
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;
}
36
37

/**
38
39
40
41
42
43
 * @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
44
 */
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
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;
    }
65

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;
}
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);

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.
122
123
 * It adjusts for the behaviour of the mod function which provides positive
 * results for position input and negative results for negative input.
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
 */
129
Real RADIX_PUBLIC cspanf(Real value, Real begin, Real end);
130

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
 */
139
Real RADIX_PUBLIC haversine(Real lat1, Real lon1, Real lat2, Real lon2);
140
141

/**
142
143
 * @brief greatCircleDistance Calculates the shortest distance between two
 * points on the surface of a sphere
144
145
146
147
 * @param lat1 Real point1 latitude
 * @param lon1 Real point1 longitude
 * @param lat2 Real point2 latitude
 * @param lon2 Real point2 longitude
148
149
150
 * @param radius
 * @return
 */
151
152
Real RADIX_PUBLIC greatCircleDistance(Real lat1, Real lon1, Real lat2,
                                      Real lon2, Real radius);
153

154
/**
155
156
 * @brief greatCircleVolume Calculates the volume giving latitude, longitude and
 * radius bounds.
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
 */
165
166
Real RADIX_PUBLIC greatCircleVolume(Real lat1, Real lon1, Real lat2, Real lon2,
                                    Real r1, Real r2);
167

168
/**
169
170
 * @brief greatCircleArea Calculates the area (m^2) given rectanglar lat/lon
 * corners at radius height
171
172
173
174
175
176
177
 * @param lat1
 * @param lon1
 * @param lat2
 * @param lon2
 * @param r1
 * @return Area in square meters
 */
178
// Real greatCircleArea(Real lat1, Real lon1, Real lat2, Real lon2, Real r1);
179

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
 */
186
Real RADIX_PUBLIC cylinderVolume(Real r, Real h);
187

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

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

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

// 360 degrees = 1 turn, and limit to being less than a full turn
197
198
inline Real RADIX_PUBLIC toTurns(Real degrees)
{
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;
207
208
209
}

/**
210
211
 * @brief lineIntersect determines the intersection distance from p1 along v1 to
 * line described p2->p3
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
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
218
 */
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)
 */
230
231
232
233
template <typename T>
T RADIX_PUBLIC interpolate(const std::vector<T> &xvals,
                           const std::vector<T> &yvals, T x,
                           bool linear = false);
234
235

/**
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
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
 */
243
double RADIX_PUBLIC gammaRayAbsorptionInAir(double energy, bool linear = false);
244
245

/**
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
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
 */
253
254
double RADIX_PUBLIC gammaRayAttenuationInAir(double energy,
                                             bool linear = false);
255
256

/**
257
258
 * @brief exponentialIntegral calculates the exponential integral of the given
 * value
259
260
261
 * @param d - value of which the exponential integral is calculated
 * @return double - the exponential integral of the given value
 */
262
double RADIX_PUBLIC exponentialIntegral(double d);
263
}  // namespace radix
264
265

#include "radixmath/util.i.hh"
266

267
#endif /* RADIX_RADIXMATH_UTIL_H */