util.cc 6.86 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
160
Real hpaToAltitude(Real hpa)
{
161
  return 0.3048 * (1 - std::pow(hpa / 1013.25, 0.190284)) * 145366.45;
162
163
}

164
/** Real greatCircleArea(Real lat1, Real lon1, Real lat2, Real lon2, Real r1)
165
166
{

167
}*/
168

169
170
double gammaRayAbsorptionInAir(double energy, bool linear)
{
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
  // 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);
188
189
}

190
191
double gammaRayAttenuationInAir(double energy, bool linear)
{
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
  // 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);
209
}
210

211
212
double exponentialIntegral(double d)
{
213
214
  // based on polynomial approximation in document provided by Vince
  // specialization for d <= 1
215
216
  if (d <= 1)
  {
217
218
219
220
221
222
223
224
225
226
227
228
    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));
229
}
230
}  // namespace radix