Skip to content

GitLab

Commit 221c5d38 authored by GODARD Vincent's avatar GODARD Vincent
Browse files

Travail sur les tutos

parent e7b00bea
Expand all Hide whitespace changes
Inline Side-by-side
DESCRIPTION
View file @ 221c5d38
Package: TCNtools
Title: What the Package Does (One Line, Title Case)
Version: 0.0.0.9000
Version: 0.1
Authors@R:
person(given = "Vincent",
family = "Godard",
......@@ -11,4 +11,8 @@ Description: What the package does (one paragraph).
License: What license it uses
Encoding: UTF-8
LazyData: true
RoxygenNote: 6.1.1
RoxygenNote: 7.1.1
Imports:
ellipse,
pracma,
shiny
NAMESPACE
View file @ 221c5d38
# Generated by roxygen2: do not edit by hand
export(angdist)
export(atm_pressure)
export(d2r)
export(ellipse_2nuclides)
export(get_vdm)
export(lsdrc)
export(r2d)
export(mu_model1a)
export(mu_model1b)
export(rc_ndp)
export(sapp_launch)
export(sato_muons)
export(sato_neutrons)
export(sato_neutrons_low)
......@@ -14,4 +16,6 @@ export(scaling_lm)
export(scaling_lsd)
export(scaling_st)
export(solar_modulation)
export(solv_conc_eul)
export(solv_conc_lag)
export(vdm2rc)
R/balco2017production.R 0 → 100644
View file @ 221c5d38
#'
#' Muons production profile (model 1A)
#'
#' Compute muons production profile as a function of depth below the surface z (g/cm2)
#' and site atmospheric pressure h (hPa). According to model 1A from Balco (2017)
#'
#' Adapted from G. Balco matlab code
#'
#' Balco, G. (\strong{2017}). Production rate calculations for cosmic-ray-muon-produced 10Be and 26Al benchmarked against geological calibration data.
#' \emph{Quaternary Geochronology}, 39, 150–173.
#' https://doi.org/10.1016/j.quageo.2017.02.001
#'
#' cst is a list containing nuclide-specific constants, as follows:
#' cst$Natoms - atom number density of the target atom (atoms/g)
#' cst$k_neg - summary cross-section for negative muon capture (atoms/muon)
#' cst$sigma190 - 190 GeV x-section for fast muon production (cm2)
#'
#' @param z depth below the surface z (g/cm2)
#' @param h site atmospheric pressure h (hPa)
#' @param cst list of nuclide specific parameters (see below)
#' @keywords
#' @export
#' @examples
#'
mu_model1a<-function(z,h,cst){
# TODO case alpha = 1
# TODO use for function phi_v_slhl(z) ?
H = (1013.25 - h)*1.019716 # figure the atmospheric depth in g/cm2
phi_vert_slhl = phi_v_slhl(z) # find the vertical flux at SLHL
R_vert_slhl = Rv0(z) # find the stopping rate of vertical muons at SLHL
R_vert_site = R_vert_slhl*exp(H/LZ(z)) # find the stopping rate of vertical muons at site
# find the flux of vertical muons at site
fi<-function(x){return(Rv0(x)*exp(H/LZ(x)))}
phi_vert_site=rep(0,length(z))
for(a in 1:length(z)) {
# integration ends at 200001 g/cm2 to avoid being asked for an zero
# range of integration --
# get integration tolerance -- want relative tolerance around 1 part in 10^4
tol = phi_vert_slhl[a] * 1e-4
temp = integrate(fi,z[a],2e5+1,abs.tol=tol)
phi_vert_site[a] = temp$value
}
# invariant flux at 2e5 g/cm2 depth - constant of integration
# calculated using commented-out formula above
a = 258.5*(100^2.66) # à supprimer en generalisant la formule plus haut
b = 75*(100^1.66)
phi_200k = (a/((2e5+21000)*(((2e5+1000)^1.66) + b)))*exp(-5.5e-6 * 2e5)
phi_vert_site = phi_vert_site + phi_200k
# find the total flux of muons at site
# angular distribution exponent
nofz = 3.21 - 0.297*log((z+H)/100 + 42) + 1.21e-5*(z+H)
# derivative of same
dndz = (-0.297/100)/((z+H)/100 + 42) + 1.21e-5
phi_temp = phi_vert_site * 2 * pi / (nofz+1)
# that was in muons/cm2/s
# convert to muons/cm2/yr
phi = phi_temp*60*60*24*365
#find the total stopping rate of muons at site
R_temp = (2*pi/(nofz+1))*R_vert_site - phi_vert_site*(-2*pi*((nofz+1)^-2))*dndz
# that was in total muons/g/s
# convert to negative muons/g/yr
R = R_temp*0.44*60*60*24*365
#Now calculate the production rates.
#Depth-dependent parts of the fast muon reaction cross-section
Beta = 0.846 - 0.015 * log((z/100)+1) + 0.003139 * (log((z/100)+1)^2)
Ebar = 7.6 + 321.7*(1 - exp(-8.059e-6*z)) + 50.7*(1-exp(-5.05e-7*z))
#internally defined constants
aalpha = 0.75
# attention prendre en compte le cas sigma0
sigma0 = cst$sigma190/(190^aalpha)
if ( "sigma190" %in% names(cst)) {
sigma0 = cst$sigma190/(190^aalpha)
} else if ( "sigma0" %in% names(cst)) {
sigma0=cst$sigma0
} else {
stop("ERROR : Undefined cross section (either sigma0 or sigma190)")
}
# fast muon production
P_fast = phi*Beta*(Ebar^aalpha)*sigma0*cst$Natoms
# negative muon capture
P_neg = R*cst$k_neg
#return
out=data.frame(z,P_fast,P_neg,phi_vert_slhl,R_vert_slhl,phi_vert_site,R_vert_site,phi,R,Beta,Ebar)
return(out)
} # end of mu_model1a
#'
#' Muons production profile (model 1B)
#'
#' Compute muons production profile as a function of depth below the surface z (g/cm2)
#' and site atmospheric pressure h (hPa). According to model 1B from Balco (2017)
#'
#' Adapted from G. Balco matlab code
#'
#' Balco, G. (\strong{2017}). Production rate calculations for cosmic-ray-muon-produced 10Be and 26Al benchmarked against geological calibration data.
#' \emph{Quaternary Geochronology}, 39, 150–173.
#' https://doi.org/10.1016/j.quageo.2017.02.001
#'
#' cst is a list containing nuclide-specific constants, as follows:
#' cst$Natoms - atom number density of the target atom (atoms/g)
#' cst$k_neg - summary cross-section for negative muon capture (atoms/muon)
#' cst$sigma190 - 190 GeV x-section for fast muon production (cm2)
#'
#' @param z depth below the surface z (g/cm2)
#' @param h site atmospheric pressure h (hPa)
#' @param Rc cutoff rigidity (GV)
#' @param Sphi solar modulation parameter
#' @param cst list of nuclide specific parameters (see below)
#' @keywords
#' @export
#' @examples
#'
mu_model1b<-function(z,h,Rc,Sphi,cst){
# figure the atmospheric depth in g/cm2
H = (1013.25 - h)*1.019716
Href = 1013.25
# find the omnidirectional flux at the site
mflux = sato_muons(h,Rc,Sphi) #Generates muon flux at site from Sato et al. (2008) model
mfluxRef = sato_muons(Href,0,462.04)
# plot(tmp[["consts"]][["mfluxRef"]][["E"]][[1]],tmp[["consts"]][["mfluxRef"]][["pos"]][[1]],type="l",log="xy")
# lines(mfluxRef[["E"]],mfluxRef[["flux_diff_pos"]],col="red")
# plot(tmp[["consts"]][["mfluxRef"]][["E"]][[1]],tmp[["consts"]][["mfluxRef"]][["neg"]][[1]],type="l",log="xy")
# lines(mfluxRef[["E"]],mfluxRef[["flux_diff_neg"]],col="red")
phi_site = (mflux$flux_diff_neg + mflux$flux_diff_pos)
phiRef = (mfluxRef$flux_diff_neg + mfluxRef$flux_diff_pos)
# find the vertical flux at SLHL
a = 258.5*(100^2.66)
b = 75*(100^1.66)
phi_vert_slhl = (a/((z+21000)*(((z+1000)^1.66) + b)))*exp(-5.5e-6 * z)
# Convert E to Range
Temp = E2R(mflux$E)
RTemp = Temp$R
#Set upper limit to stopping range to test comparability with measurements
StopLimit = 10
# find the stopping rate of vertical muons at site
# find all ranges <10 g/cm2
stopindex = RTemp<StopLimit
SFmu = phi_site/phiRef
SFmuslow = sum(phi_site[stopindex]/sum(phiRef[stopindex]))
# Prevent depths less than the minimum range in E2R to be used below
z[z < min(RTemp)] = min(RTemp)
# Find scaling factors appropriate for energies associated with stopping
# muons at depths z
Rz = approx(RTemp,SFmu,z)$y
Rz[Rz>SFmuslow] = SFmuslow
RzSpline = splinefun(RTemp, SFmu)
# find the stopping rate of vertical muons at the site, scaled from SLHL
# this is done in a subfunction Rv0, because it gets integrated below.
R_vert_slhl = Rv0(z)
R_vert_site = R_vert_slhl*Rz
# find the flux of vertical muons at site
fi<-function(x){return(Rv0(x)*RzSpline(x))}
phi_vert_site=rep(0,length(z))
for(a in 1:length(z)) {
# integration ends at 200001 g/cm2 to avoid being asked for an zero
# range of integration --
# get integration tolerance -- want relative tolerance around 1 part in 10^4
tol = phi_vert_slhl[a] * 1e-4
temp = integrate(fi,z[a],2e5+1,abs.tol=tol)
phi_vert_site[a] = temp$value
}
# invariant flux at 2e5 g/cm2 depth - constant of integration
# calculated using commented-out formula above
a = 258.5*(100^2.66)
b = 75*(100^1.66)
phi_200k = (a/((2e5+21000)*(((2e5+1000)^1.66) + b)))*exp(-5.5e-6 * 2e5)
phi_vert_site = phi_vert_site + phi_200k
# find the total flux of muons at site
# angular distribution exponent
nofz = 3.21 - 0.297*log((z+H)/100 + 42) + 1.21e-5*(z+H)
# derivative of same
dndz = (-0.297/100)/((z+H)/100 + 42) + 1.21e-5
phi_temp = phi_vert_site*2*pi/(nofz+1)
# that was in muons/cm2/s
# convert to muons/cm2/yr
phi = phi_temp*60*60*24*365
R_temp = (2*pi/(nofz+1))*R_vert_site - phi_vert_site*(-2*pi*((nofz+1)^-2))*dndz
# that was in total muons/g/s
# convert to negative muons/g/yr
R = R_temp*0.44*60*60*24*365
# Attenuation lengths
LambdaMu =(Href-h)/(log(phi_site)-log(phiRef))
# Now calculate the production rates.
# Depth-dependent parts of the fast muon reaction cross-section
# Per John Stone, personal communication 2011 - see text
aalpha = 1.0
Beta = 1.0
Ebar = 7.6 + 321.7*(1 - exp(-8.059e-6*z)) + 50.7*(1-exp(-5.05e-7*z))
# internally defined constants
# GB modified here to pass consts as does v 1.2 of P_mu_total.m
Natoms = cst$Natoms
k_neg = cst$k_neg
# attention prendre en compte le cas sigma0
sigma0 = cst$sigma190/(190^aalpha)
if ( "sigma190" %in% names(cst)) {
sigma0 = cst$sigma190/(190^aalpha)
} else if ( "sigma0" %in% names(cst)) {
sigma0=cst$sigma0
} else {
stop("ERROR : Undefined cross section (either sigma0 or sigma190)")
}
# fast muon production
P_fast = phi*Beta*(Ebar^aalpha)*sigma0*Natoms
# negative muon capture
P_neg = R*k_neg
#return
out=data.frame(z,P_fast,P_neg,phi_vert_slhl,R_vert_slhl,phi_vert_site,R_vert_site,phi,R,Beta,Ebar)
return(out)
} # end of mu_model1b
#########################################################################
LZ<-function(z){
# Balco 2017 matlab code
# this subfunction returns the effective atmospheric attenuation length for
# muons of range Z
# define range/momentum relation
# table for muons in standard rock in Groom and others 2001
# col 1 : momentum (MeV/c)
# col 2 : range (g/cm2)
data = as.data.frame(matrix(c(4.704e1, 8.516e-1,
5.616e1, 1.542e0,
6.802e1, 2.866e0,
8.509e1, 5.698e0,
1.003e2, 9.145e0,
1.527e2, 2.676e1,
1.764e2, 3.696e1,
2.218e2, 5.879e1,
2.868e2, 9.332e1,
3.917e2, 1.524e2,
4.945e2, 2.115e2,
8.995e2, 4.418e2,
1.101e3, 5.534e2,
1.502e3, 7.712e2,
2.103e3, 1.088e3,
3.104e3, 1.599e3,
4.104e3, 2.095e3,
8.105e3, 3.998e3,
1.011e4, 4.920e3,
1.411e4, 6.724e3,
2.011e4, 9.360e3,
3.011e4, 1.362e4,
4.011e4, 1.776e4,
8.011e4, 3.343e4,
1.001e5, 4.084e4,
1.401e5, 5.495e4,
2.001e5, 7.459e4,
3.001e5, 1.040e5,
4.001e5, 1.302e5,
8.001e5, 2.129e5),ncol=2,byrow=TRUE))
colnames(data) <- c("momentum","range")
# obtain momenta using log-linear interpolation
P_MeVc = exp(approx(log(data$range),log(data$momentum),log(z),rule=2)$y)
# obtain attenuation lengths
out = 263 + 150 * (P_MeVc/1000)
return(out)
}
#########################################################################
Rv0<-function(z){
# Balco 2017 matlab code
# this subfunction returns the stopping rate of vertically traveling muons
# as a function of depth z at sea level and high latitude.
a = exp(-5.5e-6*z)
b = z + 21000
c = (z + 1000)^1.66 + 1.567e5
dadz = -5.5e-6 * exp(-5.5e-6*z)
dbdz = 1
dcdz = 1.66*(z + 1000)^0.66
out = -5.401e7 * (b*c*dadz - a*(c*dbdz + b*dcdz))/(b^2 * c^2)
return(out)
}
###################################################
E2R<-function(x){
# this subfunction returns the range and energy loss values for
# muons of energy E in MeV
# define range/energy/energy loss relation
# table for muons in standard rock
# http://pdg.lbl.gov/2010/AtomicNuclearProperties/ Table 281
data =matrix(c(1.0e1, 8.400e-1, 6.619,
1.4e1, 1.530e0, 5.180,
2.0e1, 2.854e0, 4.057,
3.0e1, 5.687e0, 3.157,
4.0e1, 9.133e0, 2.702,
8.0e1, 2.675e1, 2.029,
1.0e2, 3.695e1, 1.904,
1.4e2, 5.878e1, 1.779,
2.0e2, 9.331e1, 1.710,
3.0e2, 1.523e2, 1.688,
4.0e2, 2.114e2, 1.698,
8.0e2, 4.418e2, 1.775,
1.0e3, 5.534e2, 1.808,
1.4e3, 7.712e2, 1.862,
2.0e3, 1.088e3, 1.922,
3.0e3, 1.599e3, 1.990,
4.0e3, 2.095e3, 2.038,
8.0e3, 3.998e3, 2.152,
1.0e4, 4.920e3, 2.188,
1.4e4, 6.724e3, 2.244,
2.0e4, 9.360e3, 2.306,
3.0e4, 1.362e4, 2.383,
4.0e4, 1.776e4, 2.447,
8.0e4, 3.343e4, 2.654,
1.0e5, 4.084e4, 2.747,
1.4e5, 5.495e4, 2.925,
2.0e5, 7.459e4, 3.187,
3.0e5, 1.040e5, 3.611,
4.0e5, 1.302e5, 4.037,
8.0e5, 2.129e5, 5.748),ncol=3,byrow=TRUE)
# units are range in g cm-2 (column 2)
# energy in MeV (column 1)
# Total energy loss/g/cm2 in MeV cm2/g(column 3)
# deal with zero situation
x[x<10] = 1
# obtain ranges
# use log-linear interpolation
R = exp(approx(log(data[,1]),log(data[,2]),log(x),rule=2)$y)
Eloss = exp(approx(log(data[,1]),log(data[,3]),log(x),rule=2)$y)
res=data.frame(cbind(R,Eloss))
colnames(res)<-c("R","Eloss")
return(res)
}
##########################################################################
# Heisinger et al. (2002a) eqs 1 and 2 input in g/cm2
phi_v_slhl<-function(z){
a = 258.5*(100^2.66)
b = 75*(100^1.66)
out = (a/((z+21000)*(((z+1000)^1.66) + b)))*exp(-5.5e-6 * z)
return(out)
}
R/solvers.R 0 → 100644
View file @ 221c5d38
###############################################################################################
#' Prediction of concentration evolution over a time interval on an Eulerian grid
#' Exponential attenuation models for neutrons and muons
#'
#' @param z depth coordinate of the profile (g/cm2)
#' @param ero erosion rate (g/cm2/a)
#' @param t time interval (a)
#' @param C0 inherited concentration constant with depth (at/g)
#' @param p production and decay parameters (4 elements vector)
#' p[1] -> unscaled spallation production rate (at/g/a)
#' p[2] -> unscaled stopped muons production rate (at/g/a)
#' p[3] -> unscaled fast muons production rate (at/g/a)
#' p[4] -> decay constant (1/a)
#' @param S scaling factors (2 elements vector)
#' S[1] -> scaling factor for spallation
#' S[2] -> scaling factor for muons
#' @param L attenuation length (3 elements vector)
#' L[1] -> neutrons
#' L[2] -> stopped muons
#' L[3] -> fast muons
#' @param in_ero initial erosion rate yielding a steady state inherited concentration (g/cm2/a), C0 not used in this case
#' @keywords
#' @export
#' @examples
solv_conc_eul <- function(z,ero,t,C0,p,S,L,in_ero=NULL){
p = as.numeric(p)
L = as.numeric(L)
S = as.numeric(S)
# concentration acquired over the time increment
Cspal = (S[1]*p[1])/((ero/L[1])+p[4])*exp(-1*z/L[1])*(1-exp(-1*(p[4]+(ero/L[1]))*t))
Cstop = (S[2]*p[2])/((ero/L[2])+p[4])*exp(-1*z/L[2])*(1-exp(-1*(p[4]+(ero/L[2]))*t))
Cfast = (S[2]*p[3])/((ero/L[3])+p[4])*exp(-1*z/L[3])*(1-exp(-1*(p[4]+(ero/L[3]))*t))
C_produced = Cspal + Cstop + Cfast
if (is.null(in_ero)) { # inheritance constant with depth, defined by the C0 value
C_inherited = C0*exp(-1*p[4]*t)
}
else{
Cspal_ss = (S[1]*p[1])/((in_ero/L[1])+p[4])*exp(-1*ero*t/L[1])
Cstop_ss = (S[2]*p[2])/((in_ero/L[2])+p[4])*exp(-1*ero*t/L[2])
Cfast_ss = (S[2]*p[3])/((in_ero/L[3])+p[4])*exp(-1*ero*t/L[3])
C_inherited = (Cspal_ss + Cstop_ss + Cfast_ss)*exp(-1*p[4]*t)
}
return(C_produced + C_inherited)
# # evolution of preexisting concentration
# # total = (S[1]*p[1] + S[2]*p[2] + S[2]*p[3])
# # f1 = (S[1]*p[1]) / total # fraction of production attributed to spallation
# # f2 = (S[2]*p[2]) / total
# # f3 = (S[2]*p[3]) / total
# if (C0 == 0) {
# Cspal_in = 0
# Cstop_in = 0
# Cfast_in = 0
# } else {
# total = (S[1]*p[1]*exp(-1*(p[4]+(ero/L[1]))*t) + S[2]*p[2]*exp(-1*(p[4]+(ero/L[2]))*t) + S[2]*p[3]*exp(-1*(p[4]+(ero/L[3]))*t))
# f1 = (S[1]*p[1]*exp(-1*(p[4]+(ero/L[1]))*t)) / total # fraction of production attributed to spallation
# f2 = (S[2]*p[2]*exp(-1*(p[4]+(ero/L[2]))*t)) / total
# f3 = (S[2]*p[3]*exp(-1*(p[4]+(ero/L[3]))*t)) / total
# Cspal_in = C0*f1*exp(-1*(p[4]+(ero/L[1]))*t)
# Cstop_in = C0*f2*exp(-1*(p[4]+(ero/L[2]))*t)
# Cfast_in = C0*f3*exp(-1*(p[4]+(ero/L[3]))*t)
# }
# #
# C = + Cstop_in + Cfast_in
}
#######################################################################
#' Compute concentration evolution along a time-depth history
#' Lagrangian formulation
#'
#' @param t time vector (a)
#' @param z depth vector (g/cm2), same length as t
#' @param C0 initial concentration (at/g)
#' @param Psp0 SLHL spallation production profile (at/g/a) at depth corresponding to z
#' @param Pmu0 muonic production profil at depth corresponding to z
#' @param lambda radioactive decay constant (1/a)
#' @param S scaling parameters for PsP0 and Pmu0 (2 columns), at times corresponding to t (same number of rows)
#' @param final if TRUE only compute the final concentration (default=FALSE)
#'
#' @return
#' @export
#'
#' @examples
solv_conc_lag <- function(t,z,C0,Psp0,Pmu0,lambda,S,final=FALSE){
nt=length(t)
if (!final){
C=rep(0,nt)
C[1]=C0
for(i in 2:nt) {
dt=abs(t[i]-t[i-1])
P = ( Psp0[i-1]*S[i-1,1] + Pmu0[i-1]*S[i-1,2] + Psp0[i]*S[i,1] + Pmu0[i]*S[i,2] )/2
C[i] = C0 + P*dt - (C0+P*dt/2)*lambda*dt
C0 = C[i]
}
}
else{
# C =pracma::trapz(t, (Psp0*S[,1] + Pmu0*S[,2])*exp(-1*lambda*(max(t)-t)) ) + C0*exp(-1*lambda*max(t))
C = abs(pracma::trapz(t, (Psp0*S[,1] + Pmu0*S[,2])*exp(-1*lambda*abs(t[nt]-t)) )) + C0*exp(-1*lambda*(max(t)-min(t)))
}
return(C)
}
R/utils.R
View file @ 221c5d38
......@@ -60,24 +60,30 @@ r2d<-function(a){
#' Compute atmospheric pressure at site
#'
#' Returns site pressure in hPa
#'
#'
#' Choice of models :
#'
#' stone2000 (default) : equation 1 of Stone (2000) JGR paper