main.cc 2.31 KB
Newer Older
Thompson, Adam's avatar
Thompson, Adam committed
1
2
#include "radixmath/util.hh"

3
#include "radixcore/commandline.hh"
4
5

#include <cstdlib>
Thompson, Adam's avatar
Thompson, Adam committed
6
#include <iomanip>
7
8
9
#include <iostream>

#include <cmath>
Thompson, Adam's avatar
Thompson, Adam committed
10

11
using namespace radix;
12
13
int main(int argc, char **argv)
{
14
15
16
17
18
19
20
21
  CommandLine command(argc, argv);
  command.declareArgument(
      "altitude",
      "Altitude at which to calcuate dose rate (0 for sea level) in m");
  command.declareArgument("activity", "Activity of the area in Ci/m^2");
  command.declareArgument("energy", "Gamma ray energy in MeV");

  // expecting three arguments
22
23
  if (!command.validate(std::cerr))
  {
24
    command.help(std::cerr);
Thompson, Adam's avatar
Thompson, Adam committed
25
    return 0;
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
  }

  // altitude above sea level - default to 0
  double z = 0.0;
  z        = command.arg<double>(0);

  // activity in Ci/m^2 of the area of interest - default to 1.767e3
  double a = 1.767e3;
  a        = command.arg<double>(1);

  // gamma-ray energy in MeV - default to 7e-1
  double e = 7e-1;
  e        = command.arg<double>(2);

  // radius of area of interest (assumed to be altitude), distance
  //    double r = z;
  //    double s = std::sqrt(z*z + r*r);

  // sigma scalar; constant at sea level
  double sigma = 1.225;

  // not at sea level, update sigma
48
49
  if (std::fabs(z) > 1e-16)
  {
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
    sigma = 1.911515e-16 * (z * z * z * z) - 9.857582e-13 * (z * z * z) +
            5.399174e-9 * (z * z) - 1.180686e-4 * z + 1.22507;
  }

  // absorption/attenuation of gamma ray as it passes through air
  //    auto mua = radix::gammaRayAbsorptionInAir(e);
  auto mut = radix::gammaRayAttenuationInAir(e);

  // mu scalars
  auto mutg = mut * 1.225;
  //    auto muag = -mua * 1.225;
  auto mutz = mut * sigma;
  //    auto muaz = -mua * sigma;

  // exponential integral values and their ratio
  auto EIg     = 1.0 / radix::exponentialIntegral(mutg);
  auto EIz     = 1.0 / radix::exponentialIntegral(mutz * z);
  auto EIRatio = EIz / EIg;

  // dose rate at sea level/altitude
  auto DRg = 9.62 * a;
  auto DRz = DRg * EIRatio;

  // report dose rate at sea level
74
75
  if (std::fabs(z) < 1e-16)
  {
76
77
78
79
    std::cout << std::scientific << DRg
              << " is the Dose Rate at Ground Sea Level in R/hr" << std::endl;
  }
  // report dose rate at altitude above sea level
80
81
  else
  {
82
83
84
85
86
87
88
    std::cout
        << std::scientific << DRz
        << " is the Dose Rate at your Altitude above Ground Sea Level in R/hr"
        << std::endl;
  }

  return 0;
Thompson, Adam's avatar
Thompson, Adam committed
89
}