PREM.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2019, The GENIE Collaboration
6
7  Author: Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
8  University of Liverpool & STFC Rutherford Appleton Lab
9
10  For documentation see the corresponding header file.
11
12  Important revisions after version 2.0.0 :
13  @ Aug 25, 2009 - CA
14  Was first added in the development version 2.5.1
15
16 */
17 //____________________________________________________________________________
18
19 #include <TMath.h>
20
23 #include "Framework/Utils/PREM.h"
24
25 //___________________________________________________________________________
27 {
28 // Return the Earth density according to the PREM model
29 // Inputs: r, Distance from the centre of the Earth (in std GENIE units)
30 // Outputs: rho, Earth density (in std GENIE units)
31 //
32
33  r = TMath::Max(0., r/units::km); // convert to km
34
35  double rE = constants::kREarth/units::km;
36  double rho = 0.;
37  double x = r / rE;
38
39  if (r <= 1221.5 )
40  {
41  rho = 13.0885 - 8.8381*x*x;
42  }
43  else if (r > 1221.5 && r <= 3480.0 )
44  {
45  rho = 12.5815 - 1.2638*x - 3.6426*x*x - 5.5281*x*x*x;
46  }
47  else if (r > 3480.0 && r <= 5701.0 )
48  {
49  rho = 7.9565 - 6.4761*x + 5.5283*x*x - 3.0807*x*x*x;
50  }
51  else if (r > 5701.0 && r <= 5771.0 )
52  {
53  rho = 5.3197 - 1.4836*x;
54  }
55  else if (r > 5771.0 && r <= 5971.0 )
56  {
57  rho = 11.2494 - 8.0298*x;
58  }
59  else if (r > 5971.0 && r <= 6151.0 )
60  {
61  rho = 7.1089 - 3.8045*x;
62  }
63  else if (r > 6151.0 && r <= 6346.6 )
64  {
65  rho = 2.691 + 0.6924*x;
66  }
67  else if (r > 6346.6 && r <= 6356.0 )
68  {
69  rho = 2.90;
70  }
71  else if (r > 6356.0 && r <= 6368.0 )
72  {
73  rho = 2.60;
74  }
75  else if (r > 6368.0 && r <= rE)
76  {
77  rho = 1.02;
78  }
79
80  rho = rho * units::g_cm3;
81
82  return rho;
83 }
84 //___________________________________________________________________________
static const double g_cm3
Definition: Units.h:154
static constexpr Double_t km
Definition: Munits.h:148
static const double kREarth
Definition: Constants.h:111
TRandom3 r(0)
double Density(double r)
Definition: PREM.cxx:26