scaling_lm.RdAdapted from G. Balco matlab code (http://hess.ess.washington.edu/math)
scaling_lm(h, Rc)atmospheric pressure (hPa)
cutoff rigidity (GV)
scaling_lm
#> function (h, Rc)
#> {
#> if (length(h) > 1) {
#> stop("Function not vectorized for atmospheric pressure")
#> }
#> ilats_d = c(0, 10, 20, 30, 40, 50, 60)
#> ilats_r = d2r(ilats_d)
#> iRcs = 14.9 * ((cos(ilats_r))^4)
#> a = c(31.8518, 34.3699, 40.3153, 42.0983, 56.7733, 69.072,
#> 71.8733)
#> b = c(250.3193, 258.4759, 308.9894, 512.6857, 649.1343, 832.4566,
#> 863.1927)
#> c = c(-0.083393, -0.089807, -0.106248, -0.120551, -0.160859,
#> -0.199252, -0.207069)
#> d = c(7.426e-05, 7.9457e-05, 9.4508e-05, 0.00011752, 0.00015463,
#> 0.00019391, 0.00020127)
#> e = c(-2.2397e-08, -2.3697e-08, -2.8234e-08, -3.8809e-08,
#> -5.033e-08, -6.3653e-08, -6.6043e-08)
#> sf = a + (b * exp(h/(-150))) + (c * h) + (d * (h^2)) + (e *
#> (h^3))
#> iRcs[8] = 0
#> sf[8] = sf[7]
#> fits = lm(log10(sf[1:3]) ~ log10(iRcs[1:3]))
#> add_iRcs = c(30, 28, 26, 24, 22, 21, 20, 19, 18, 17, 16)
#> add_sf = exp(log(sf[1]) + fits$coefficients[2] * (log10(add_iRcs) -
#> log10(iRcs[1])))
#> sf = c(add_sf, sf)
#> iRcs = c(add_iRcs, iRcs)
#> tmp = data.frame(iRcs, sf)
#> tmp = tmp[order(iRcs), ]
#> return(approx(tmp[, 1], tmp[, 2], Rc)$y)
#> }
#> <bytecode: 0x56376d006fd0>
#> <environment: namespace:TCNtools>