sato_neutrons.RdCompute neutron flux using the analytical formulation of Sato et al. (2008)
sato_neutrons(h, Rc, s, w)atmospheric pressure (hPa)
cutoff rigidity (GV)
solar modulation potential
gravimetric fractional water content (non-dimensional)
Adaption of matlab code from Lifton et al. (2014)
Lifton, N., Sato, T., & Dunai, T. J. (2014). Scaling in situ cosmogenic nuclide production rates using analytical approximations to atmospheric cosmic-ray fluxes. Earth and Planetary Science Letters, 386, 149–160. https://doi.org/10.1016/j.epsl.2013.10.052
sato_neutrons
#> function (h, Rc, s, w)
#> {
#> x = h * 1.019716
#> E = 10^seq(0, 5.301, length.out = 200)
#> scaling = data.frame(matrix(NA, nrow = length(Rc), ncol = 6))
#> colnames(scaling) <- c("Rc", "nflux", "P3n", "P10n", "P14n",
#> "P26n")
#> scaling$Rc = Rc
#> PhiG_d = matrix(rep(0, length(Rc) * length(E)), nrow = length(Rc),
#> ncol = length(E))
#> Rc[Rc < 1] = 1
#> Et = 2.5e-08
#> smin = 400
#> smax = 1200
#> a6 = 0.00018882
#> a7 = 0.44791
#> a8 = 0.0014361
#> a12 = 0.014109
#> b11min = 25.702
#> b11max = -6.9221
#> b12min = -0.50931
#> b12max = 1.1336
#> b13min = 7.465
#> b13max = 26.961
#> b14min = 12.313
#> b14max = 11.746
#> b15min = 1.0498
#> b15max = 2.6171
#> b21min = 0.0065143
#> b21max = 0.0053211
#> b22min = 3.3511e-05
#> b22max = 8.4899e-05
#> b23min = 0.00094415
#> b23max = 0.0020704
#> b24min = 12.088
#> b24max = 11.714
#> b25min = 2.7782
#> b25max = 3.8051
#> b31min = 0.98364
#> b31max = 0.97536
#> b32min = 0.00014964
#> b32max = 0.00064733
#> b33min = -0.73249
#> b33max = -0.2275
#> b34min = -1.4381
#> b34max = 2.0713
#> b35min = 2.7448
#> b35max = 2.1689
#> b41min = 0.0088681
#> b41max = 0.0091435
#> b42min = -4.3322e-05
#> b42max = -6.4855e-05
#> b43min = 0.017293
#> b43max = 0.0058179
#> b44min = -1.0836
#> b44max = 1.0168
#> b45min = 2.6602
#> b45max = 2.4504
#> b121 = 0.931
#> b122 = 0.037
#> b123 = -2.02
#> b124 = 2.12
#> b125 = 5.34
#> b131 = 0.000667
#> b132 = -1.19e-05
#> b133 = 1e-04
#> b134 = 1.45
#> b135 = 4.29
#> b51 = 0.00097337
#> b52 = -9.6635e-05
#> b53 = 0.012121
#> b54 = 7.1726
#> b55 = 1.4601
#> b91 = 571.99
#> b92 = 7.1293
#> b93 = -107.03
#> b94 = 1.8538
#> b95 = 1.2142
#> b101 = 0.00068552
#> b102 = 2.7136e-05
#> b103 = 0.00057823
#> b104 = 8.8534
#> b105 = 3.6417
#> b111 = -0.508
#> b112 = 0.14783
#> b113 = 1.0068
#> b114 = 9.1556
#> b115 = 1.6369
#> c1 = 0.23555
#> c2 = 2.3779
#> c3 = 0.72597
#> c5 = 123.91
#> c6 = 2.2318
#> c7 = 0.0010791
#> c8 = 3.6435e-12
#> c9 = 1.6595
#> c10 = 8.4782e-08
#> c11 = 1.5054
#> h31 = -25.184
#> h32 = 2.7298
#> h33 = 0.071526
#> h51 = 0.3479
#> h52 = 3.3493
#> h53 = -1.5744
#> g1 = -0.023499
#> g2 = -0.012938
#> g3 = 10^(h31 + h32/(w + h33))
#> g4 = 0.96889
#> g5 = h51 + h52 * w + h53 * (w^2)
#> fG = 10^(g1 + g2 * log10(E/g3) * (1 - tanh(g4 * log10(E/g5))))
#> h61 = 0.118
#> h62 = 0.14438
#> h63 = 3.8733
#> h64 = 0.65298
#> h65 = 42.752
#> g6 = (h61 + h62 * exp(-h63 * w))/(1 + h64 * exp(-h65 * w))
#> PhiT = g6 * ((E/Et)^2) * exp(-E/Et)
#> PhiB = rep(0, length(E))
#> PhiG = rep(0, length(E))
#> PhiGMev = rep(0, length(E))
#> p10n = rep(0, length(E))
#> p14n = rep(0, length(E))
#> p26n = rep(0, length(E))
#> p3n = rep(0, length(E))
#> p36Can = rep(0, length(E))
#> p36Kn = rep(0, length(E))
#> p36Tin = rep(0, length(E))
#> p36Fen = rep(0, length(E))
#> for (a in 1:length(Rc)) {
#> a1min = b11min + b12min * Rc[a] + b13min/(1 + exp((Rc[a] -
#> b14min)/b15min))
#> a1max = b11max + b12max * Rc[a] + b13max/(1 + exp((Rc[a] -
#> b14max)/b15max))
#> a2min = b21min + b22min * Rc[a] + b23min/(1 + exp((Rc[a] -
#> b24min)/b25min))
#> a2max = b21max + b22max * Rc[a] + b23max/(1 + exp((Rc[a] -
#> b24max)/b25max))
#> a3min = b31min + b32min * Rc[a] + b33min/(1 + exp((Rc[a] -
#> b34min)/b35min))
#> a3max = b31max + b32max * Rc[a] + b33max/(1 + exp((Rc[a] -
#> b34max)/b35max))
#> a4min = b41min + b42min * Rc[a] + b43min/(1 + exp((Rc[a] -
#> b44min)/b45min))
#> a4max = b41max + b42max * Rc[a] + b43max/(1 + exp((Rc[a] -
#> b44max)/b45max))
#> a5 = b51 + b52 * Rc[a] + b53/(1 + exp((Rc[a] - b54)/b55))
#> a9 = b91 + b92 * Rc[a] + b93/(1 + exp((Rc[a] - b94)/b95))
#> a10 = b101 + b102 * Rc[a] + b103/(1 + exp((Rc[a] - b104)/b105))
#> a11 = b111 + b112 * Rc[a] + b113/(1 + exp((Rc[a] - b114)/b115))
#> b5 = b121 + b122 * Rc[a] + b123/(1 + exp((Rc[a] - b124)/b125))
#> b6 = b131 + b132 * Rc[a] + b133/(1 + exp((Rc[a] - b134)/b135))
#> c4 = a5 + a6 * x/(1 + a7 * exp(a8 * x))
#> c12 = a9 * (exp(-a10 * x) + a11 * exp(-a12 * x))
#> PhiLmin = a1min * (exp(-a2min * x) - a3min * exp(-a4min *
#> x))
#> PhiLmax = a1max * (exp(-a2max * x) - a3max * exp(-a4max *
#> x))
#> f3 = b5 + b6 * x
#> f2 = (PhiLmin - PhiLmax)/(smin^f3 - smax^f3)
#> f1 = PhiLmin - f2 * smin^f3
#> PhiL = f1 + f2 * s[a]^f3
#> PhiB = (c1 * (E/c2)^c3) * exp(-E/c2) + c4 * exp((-(log10(E) -
#> log10(c5))^2)/(2 * (log10(c6))^2)) + c7 * log10(E/c8) *
#> (1 + tanh(c9 * log10(E/c10))) * (1 - tanh(c11 * log10(E/c12)))
#> PhiG = PhiL * (PhiB * fG + PhiT)
#> PhiGMev = PhiG/E
#> PhiG_d[a, ] = PhiG
#> clipindex = (E >= 1)
#> scaling$P3n[a] = (pracma::trapz(E[clipindex], PhiGMev[clipindex] *
#> XSectsReedyAll$OnxHe3T[clipindex]) + pracma::trapz(E[clipindex],
#> PhiGMev[clipindex] * XSectsReedyAll$SinxHe3T[clipindex]/2)) *
#> XSectsReedyAll$Natoms3 * 1e-27 * 31536000
#> scaling$P10n[a] = (pracma::trapz(E[clipindex], PhiGMev[clipindex] *
#> XSectsReedyAll$O16nxBe10[clipindex]) + pracma::trapz(E[clipindex],
#> PhiGMev[clipindex] * XSectsReedyAll$SinxBe10[clipindex]/2)) *
#> XSectsReedyAll$Natoms10 * 1e-27 * 31536000
#> scaling$P14n[a] = (pracma::trapz(E[clipindex], PhiGMev[clipindex] *
#> XSectsReedyAll$O16nn2pC14[clipindex]) + pracma::trapz(E[clipindex],
#> PhiGMev[clipindex] * XSectsReedyAll$SinxC14[clipindex]/2)) *
#> XSectsReedyAll$Natoms14 * 1e-27 * 31536000
#> scaling$P26n[a] = pracma::trapz(E[clipindex], PhiGMev[clipindex] *
#> XSectsReedyAll$SinxAl26[clipindex]) * XSectsReedyAll$Natoms26 *
#> 1e-27 * 31536000
#> scaling$nflux[a] = pracma::trapz(E[clipindex], PhiGMev[clipindex])
#> }
#> res = list(E, PhiG_d, scaling)
#> names(res) <- c("E", "PhiG_d", "scaling")
#> return(res)
#> }
#> <bytecode: 0x56376ad9c9b0>
#> <environment: namespace:TCNtools>