util.cc 9.92 KB
Newer Older
1
2
3
4
5
6
7
8
/*
 * File:   util.cc
 * Author: a16
 *
 * Created on January 17, 2013, 9:27 AM
 */

#include "radixmath/util.hh"
9
#include <algorithm>
10
#include <cstddef>
11
12
#include "radixmath/point3d.hh"
#include "radixmath/vector3d.hh"
13

14
15
16
17
18
namespace radix
{
std::vector<Real> getEqualAreaRadii(const std::vector<Real> &radii,
                                    const std::vector<unsigned short> &subrings)
{
19
20
  std::vector<Real> newRadii;
  Real prevRadius;
21

22
  prevRadius = 0;
23

24
25
  for (size_t index = 0; index < radii.size(); index++)
  {
26
27
28
29
30
    Real _prevRadius;
    int _subrings;
    Real area;
    Real radius;
    Real subarea;
31

32
33
34
    radius      = radii[index];
    _prevRadius = prevRadius;
    _subrings   = subrings[index];
35

36
37
    area    = PI * radius * radius - PI * prevRadius * prevRadius;
    subarea = area / _subrings;
38

39
40
    for (int count = 1; count < _subrings; count++)
    {
41
      Real _radius;
42

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

46
      newRadii.push_back(_radius);
47
48
    }

49
    newRadii.push_back(radius);
50

51
52
53
54
    prevRadius = radius;
  }

  return newRadii;
55
}
56

57
Real lineIntersect(const Point3D &p, const Vector3D &v, const Point3D &sp,
58
59
                   const Point3D &ep)
{
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
  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;
94
95
}

96
97
Real haversine(Real lat1, Real lon1, Real lat2, Real lon2)
{
98
99
100
101
102
103
104
105
106
107
108
109
  // 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;
110
111
}

112
Real greatCircleDistance(Real lat1, Real lon1, Real lat2, Real lon2,
113
114
                         Real radius)
{
115
  return radius * haversine(lat1, lon1, lat2, lon2);
116
117
}

118
Real greatCircleVolume(Real lat1, Real lon1, Real lat2, Real lon2, Real r1,
119
120
                       Real r2)
{
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
  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);
139
140
}

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

143
144
Real cspanf(Real value, Real begin, Real end)
{
145
146
147
148
  Real first = 0.0, last = 0.0;
  first = std::min(begin, end);
  last  = std::max(begin, end);
  value = std::fmod(value - first, last - first);
149
150
  if (value <= 0)
  {
151
    return value + last;
152
153
154
  }
  else
  {
155
156
    return value + first;
  }
157
158
}

159
Real hpaToAltitude(Real hpa, Real relHum, Real temperature, Real msle)
160
{
161
162
163
  // convert temperature to virtual temperature
  temperature = virtualTemperature(temperature, hpa, relHum);

164
165
  Real result = (SPECIFIC_GAS_CONSTANT * temperature) /
                GRAVITATIONAL_ACCELERATION * std::log(msle / hpa);
166
167
  // divide out the kg/mol using the molar mass of dry air and water
  result = result / MOLAR_MASS_DRY_AIR;
168
169
170
  radix_tagged_line(result << " = hpaToAltitude(" << hpa << ", " << temperature
                           << ", " << msle << ")");
  return result;
171
172
}

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

176
}*/
177

178
179
double gammaRayAbsorptionInAir(double energy, bool linear)
{
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
  // 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);
197
198
}

199
200
double gammaRayAttenuationInAir(double energy, bool linear)
{
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
  // 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);
218
}
219

220
221
double exponentialIntegral(double d)
{
222
223
  // based on polynomial approximation in document provided by Vince
  // specialization for d <= 1
224
225
  if (d <= 1)
  {
226
227
228
229
230
231
232
233
234
235
236
237
    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));
238
}
LEFEBVREJP email's avatar
LEFEBVREJP email committed
239

240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
Real absoluteTemperature(Real potentialTemp, Real pressure)
{
  // Gas constant / specific heat capacity at constant pressure for air
  double rDivCp = 0.286;

  // Reference pressure (mb)
  double refPressure = 1000.0;

  double expComponent = pow((refPressure / pressure), rDivCp);

  double temp = potentialTemp / expComponent;

  return temp;
}

Real saturationVaporPressure(Real absTemp, Real pressure)
{
  // Convert temperature to C and limit range to -45 < temp < +60
  double tempC = absTemp + ABS_ZERO_CELSIUS;
  tempC        = std::max(-45.0, std::min(60.0, tempC));

  // Calculate saturation vapour pressure
  double satVapPress = 6.112 * exp(17.62 * tempC / (243.12 + tempC));

  // Apply correction for moist air
  double moistCorrection =
      1.0 + (4.5 + 0.0006 * pow(tempC, 2)) * pressure * 1.0E-6;
  satVapPress = satVapPress * moistCorrection;

  return satVapPress;
}

Real mixingRatioToRelativeHumidity(Real mixingRatio, Real absTemp,
                                   Real pressure)
{
  // Calculate the saturation vapour pressure
  double satVapPress = saturationVaporPressure(absTemp, pressure);

  // Calculate the vapour pressure
  double vapPress =
      mixingRatio / (MASS_RATIO_WATER_VAPOR_DRY_AIR + mixingRatio) * pressure;

  // Calculate the relative humidity - ensure it is 0 < humidity < 100
  double relHum = 100.0 * vapPress / satVapPress;
  relHum        = std::max(0.0, std::min(100.0, relHum));

  return relHum;
}

Real specificHumidityToRelativeHumidity(Real specificHumidity, Real absTemp,
                                        Real pressure)
{
  // Convert input specific humidity from g/kg to kg/kg
  specificHumidity = specificHumidity / 1000.0;

  // Calculate the saturation vapour pressure
  double satVapPress = saturationVaporPressure(absTemp, pressure);

  // Calculate the vapour pressure
  double vapPress = specificHumidity /
                    (MASS_RATIO_WATER_VAPOR_DRY_AIR +
                     (1 - MASS_RATIO_WATER_VAPOR_DRY_AIR) * specificHumidity) *
                    pressure;

  // Calculate the relative humidity - ensure it is 0 < humidity < 100
  double relHum = 100.0 * vapPress / satVapPress;
  relHum        = std::max(0.0, std::min(100.0, relHum));

  return relHum;
}

311
312
313
314
315
316
317
318
319
320
321
322
323
Real virtualTemperature(Real temperature, Real pressure, Real relHum)
{
  Real mR = mixingRatio(temperature, pressure, relHum);
  return temperature * (1. + (mR / MASS_RATIO_WATER_VAPOR_DRY_AIR)) / (1. + mR);
}

Real mixingRatio(Real temperature, Real pressure, Real relHum)
{
  Real e_s = saturationVaporPressure(temperature, pressure);
  Real e_a = e_s * (relHum / 100.);
  return (0.622 * e_a) / (pressure - e_a);
}

324
}  // namespace radix