sato_protons.RdCompute protons flux using the analytical formulation of Sato et al. (2008)
sato_protons(h, Rc, s)atmospheric pressure (hPa)
cutoff rigidity (GV)
solar modulation potential
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_protons
#> function (h, Rc, s)
#> {
#> 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", "pflux", "P3p", "P10p", "P14p",
#> "P26p")
#> scaling$Rc = Rc
#> phiPtot_d = matrix(rep(0, length(Rc) * length(E)), nrow = length(Rc),
#> ncol = length(E))
#> A = 1
#> Z = 1
#> Ep = 938.27
#> U = (4 - 1.675) * pi * A/Z * 1e-07
#> Rc[Rc < 1] = 1
#> smin = 400
#> smax = 1200
#> a1 = 2.1153
#> a2 = 0.44511
#> a3 = 0.010064
#> a4 = 0.039564
#> a5 = 2.9236
#> a6 = 2.7076
#> a7 = 12663
#> a8 = 4828.8
#> a9 = 32822
#> a10 = 7437.8
#> a11 = 3.4643
#> a12 = 1.6752
#> a13 = 1.3691
#> a14 = 2.0665
#> a15 = 108.33
#> a16 = 2301.3
#> Etoa = E + a1 * x
#> Rtoa = 0.001 * sqrt((A * Etoa)^2 + 2 * A * Ep * Etoa)/Z
#> Elis = rep(0, length(E))
#> Beta = rep(0, length(E))
#> Rlis = rep(0, length(E))
#> phiTOA = rep(0, length(E))
#> phiLIS = rep(0, length(E))
#> phiSec = rep(0, length(E))
#> phiPtot = rep(0, length(E))
#> p10p = rep(0, length(E))
#> c11 = 1.256
#> c12 = 0.003226
#> c13 = -4.8077e-06
#> c14 = 2.2825e-09
#> c21 = 0.43783
#> c22 = -0.00055759
#> c23 = 7.8388e-07
#> c24 = -3.8671e-10
#> c31 = 0.00018102
#> c32 = -5.1754e-07
#> c33 = 7.5876e-10
#> c34 = -3.822e-13
#> c41 = 1.7065
#> c42 = 0.00071608
#> c43 = -9.322e-07
#> c44 = 5.2665e-10
#> b1 = c11 + c12 * x + c13 * x^2 + c14 * x^3
#> b2 = c21 + c22 * x + c23 * x^2 + c24 * x^3
#> b3 = c31 + c32 * x + c33 * x^2 + c34 * x^3
#> b4 = c41 + c42 * x + c43 * x^2 + c44 * x^3
#> h11min = 0.0024354
#> h11max = 0.002545
#> h12min = -6.0339e-05
#> h12max = -7.1807e-05
#> h13min = 0.0021951
#> h13max = 0.001458
#> h14min = 6.6767
#> h14max = 6.915
#> h15min = 0.93228
#> h15max = 0.99366
#> h21min = 0.0077872
#> h21max = 0.0076828
#> h22min = -9.5771e-06
#> h22max = -2.4119e-06
#> h23min = 0.00062229
#> h23max = 0.00066411
#> h24min = 7.7842
#> h24max = 7.7461
#> h25min = 1.8502
#> h25max = 1.9431
#> h31min = 0.9634
#> h31max = 0.97353
#> h32min = 0.0015974
#> h32max = 0.0010577
#> h33min = -0.071179
#> h33max = -0.021383
#> h34min = 2.232
#> h34max = 3.0058
#> h35min = 0.788
#> h35max = 0.91845
#> h41min = 0.0078132
#> h41max = 0.0073482
#> h42min = 9.7085e-11
#> h42max = 2.5598e-05
#> h43min = 0.00082392
#> h43max = 0.0012457
#> h44min = 8.5138
#> h44max = 8.1896
#> h45min = 2.3125
#> h45max = 2.9368
#> h51 = 0.191
#> h52 = 0.0703
#> h53 = -0.645
#> h54 = 2.03
#> h55 = 1.3
#> h61 = 0.000571
#> h62 = 6.13e-06
#> h63 = 0.000547
#> h64 = 1.11
#> h65 = 0.837
#> for (a in 1:length(Rc)) {
#> Elis = Etoa + s[a] * Z/A
#> Beta = sqrt(1 - (Ep/(Ep + Elis * A))^2)
#> Rlis = 0.001 * sqrt((A * Elis)^2 + 2 * A * Ep * Elis)/Z
#> C = a7 + a8/(1 + exp((Elis - a9)/a10))
#> phiTOA = (C * (Beta^a5)/(Rlis^a6)) * (Rtoa/Rlis)^2
#> phiPri = (U/Beta) * phiTOA * (a2 * exp(-a3 * x) + (1 -
#> a2) * exp(-a4 * x))
#> g1min = h11min + h12min * Rc[a] + h13min/(1 + exp((Rc[a] -
#> h14min)/h15min))
#> g1max = h11max + h12max * Rc[a] + h13max/(1 + exp((Rc[a] -
#> h14max)/h15max))
#> g2min = h21min + h22min * Rc[a] + h23min/(1 + exp((Rc[a] -
#> h24min)/h25min))
#> g2max = h21max + h22max * Rc[a] + h23max/(1 + exp((Rc[a] -
#> h24max)/h25max))
#> g3min = h31min + h32min * Rc[a] + h33min/(1 + exp((Rc[a] -
#> h34min)/h35min))
#> g3max = h31max + h32max * Rc[a] + h33max/(1 + exp((Rc[a] -
#> h34max)/h35max))
#> g4min = h41min + h42min * Rc[a] + h43min/(1 + exp((Rc[a] -
#> h44min)/h45min))
#> g4max = h41max + h42max * Rc[a] + h43max/(1 + exp((Rc[a] -
#> h44max)/h45max))
#> phiPmin = g1min * (exp(-g2min * x) - g3min * exp(-g4min *
#> x))
#> phiPmax = g1max * (exp(-g2max * x) - g3max * exp(-g4max *
#> x))
#> g5 = h51 + h52 * Rc[a] + h53/(1 + exp((Rc[a] - h54)/h55))
#> g6 = h61 + h62 * Rc[a] + h63/(1 + exp((Rc[a] - h64)/h65))
#> f3 = g5 + g6 * x
#> f2 = (phiPmin - phiPmax)/(smin^f3 - smax^f3)
#> f1 = phiPmin - f2 * smin^f3
#> phiP = f1 + f2 * s[a]^f3
#> phiSec = (phiP * b1 * E^b2)/(1 + b3 * E^b4)
#> Ec = (sqrt((1000 * Rc[a] * Z)^2 + Ep^2) - Ep)/A
#> Es = a13 * (Ec - a14 * x)
#> Es1 = max(a15, Es)
#> Es2 = max(a16, Es)
#> phiPtot = phiPri * (tanh(a11 * (E/Es1 - 1)) + 1)/2 +
#> phiSec * (tanh(a12 * (1 - E/Es2)) + 1)/2
#> phiPtot_d[a, ] = phiPtot
#> clipindex = (E >= 1)
#> scaling$P3p[a] = (pracma::trapz(E, phiPtot * XSectsReedyAll$OpxHe3T) +
#> pracma::trapz(E, phiPtot * XSectsReedyAll$SipxHe3T/2)) *
#> XSectsReedyAll$Natoms3 * 1e-27 * 31536000
#> scaling$P10p[a] = (pracma::trapz(E, phiPtot * XSectsReedyAll$O16pxBe10) +
#> pracma::trapz(E, phiPtot * XSectsReedyAll$SipxBe10/2)) *
#> XSectsReedyAll$Natoms10 * 1e-27 * 31536000
#> scaling$P14p[a] = (pracma::trapz(E, phiPtot * XSectsReedyAll$O16pxC14) +
#> pracma::trapz(E, phiPtot * XSectsReedyAll$SipxC14/2)) *
#> XSectsReedyAll$Natoms14 * 1e-27 * 31536000
#> scaling$P26p[a] = pracma::trapz(E, phiPtot * XSectsReedyAll$SipxAl26) *
#> XSectsReedyAll$Natoms26 * 1e-27 * 31536000
#> scaling$pflux[a] = pracma::trapz(E[clipindex], phiPtot[clipindex])
#> }
#> res = list(E, phiPtot_d, scaling)
#> names(res) <- c("E", "phiPtot_d", "scaling")
#> return(res)
#> }
#> <bytecode: 0x56376eb12c70>
#> <environment: namespace:TCNtools>