Commit e7b00bea authored by GODARD Vincent's avatar GODARD Vincent
Browse files

mise à jour des fonctions des base

parent 773e52c3
...@@ -20,9 +20,6 @@ solar_modulation<-function(t){ ...@@ -20,9 +20,6 @@ solar_modulation<-function(t){
#' Compute cutoff rigidities time series for the calculation of scaling factors according to Lifton et al. (2014) #' Compute cutoff rigidities time series for the calculation of scaling factors according to Lifton et al. (2014)
#' 'LSD' scaling scheme #' 'LSD' scaling scheme
#' #'
......
...@@ -70,13 +70,15 @@ scaling_lm<-function(h,Rc){ ...@@ -70,13 +70,15 @@ scaling_lm<-function(h,Rc){
# ought to be linear above 10 GV. Note that this is speculative, but # ought to be linear above 10 GV. Note that this is speculative, but
# relatively unimportant, as the approximation is pretty much only used # relatively unimportant, as the approximation is pretty much only used
# for low latitudes for a short time in the Holocene. # for low latitudes for a short time in the Holocene.
fits=lm(log(sf[1:3])~log(iRcs[1:3])) fits=lm(log10(sf[1:3])~log10(iRcs[1:3]))
add_iRcs=c(21,20,19,18,17,16) add_iRcs=c(30,28,26,24,22,21,20,19,18,17,16)
add_sf = exp( log(sf[1]) + fits$coefficients[2]*( log(add_iRcs) - log(iRcs[1]) ) ) add_sf = exp( log(sf[1]) + fits$coefficients[2]*( log10(add_iRcs) - log10(iRcs[1]) ) )
sf = c(add_sf,sf) sf = c(add_sf,sf)
iRcs = c(add_iRcs,iRcs) iRcs = c(add_iRcs,iRcs)
tmp = data.frame(iRcs,sf)
tmp =tmp[order(iRcs),]
# interpolate # interpolate
return(approx(iRcs,sf,Rc)$y) return(approx(tmp[,1],tmp[,2],Rc)$y)
} }
......
#' Compute atmospheric pressure at site
#'
#' Returns site pressure in hPa
#'
#' Choice of models :
#'
#' stone2000 (default) : equation 1 of Stone (2000) JGR paper
#'
#' era40 : Looks up mean sea level pressure and mean 1000 mb temp from ERA-40 reanalysis
#' and calculates site atmospheric pressures using these as inputs to the standard atmosphere equation
#'
#'
#' @param alt Altitude (m)
#' @param lon Longitude (m) optional (only for "era40" model)
#' @param lat Latitude (m) optional (only for "era40" model)
#' @param model Model used for computation, one of "stone2000" or "era40"
#' @keywords
#' @export
#' @examples
#' atm_pressure
atm_pressure<-function(alt,lon=NULL,lat=NULL,model="stone2000"){
if (model == "stone2000"){
Ps=1013.25 # sea level pressure (hPa)
Ts=288.15 # sea level temperature (K)
ksi=0.0065 # adiabatic lapse rate (K/m)
gMR=0.03417 # g*M/R (K/m)
#hydrostatic pressure (equation 1 of Stone 2000)
P=Ps*exp(-1*gMR/ksi*(log(Ts) - log(Ts - ksi* alt)))
}
else if (model == "era40"){
if(is.null(lat) | is.null(lon)) {stop("must specify lat and lon for model era40")}
z = alt
# get longitude from 0 to 360
lon_c=lon+(lon<0)*360
# Some constants and definitions
gmr = -0.03417
lr = c(-6.1517E-03,-3.1831E-06,-1.5014E-07,1.8097E-09,1.1791E-10,-6.5359E-14,-9.5209E-15)
# interpolate
# ERA40 packaged in internal data (sysdata.rda)
site_slp=pracma::interp2(ERA40$ERA40lon,ERA40$ERA40lat,ERA40$meanP, lon_c, lat)
site_T=pracma::interp2(ERA40$ERA40lon,ERA40$ERA40lat,ERA40$meanT, lon_c, lat)
# Lifton Lapse Rate Fit to COSPAR CIRA-86 <10 km altitude
dtdz = lr[1] + lr[2]*lat + lr[3]*lat^2 + lr[4]*lat^3 + lr[5]*lat^4 + lr[6]*lat^5 + lr[7]*lat^6
dtdz = -dtdz
P = site_slp * exp( (gmr/dtdz) * ( log(site_T) - log(site_T - (z*dtdz)) ) )
}
else {
stop("model must be stone2000 or era40")
}
return(P)
}
############################################################################################### ###############################################################################################
...@@ -68,7 +10,6 @@ atm_pressure<-function(alt,lon=NULL,lat=NULL,model="stone2000"){ ...@@ -68,7 +10,6 @@ atm_pressure<-function(alt,lon=NULL,lat=NULL,model="stone2000"){
#' @param lat2 latitude of point 2 (degrees) #' @param lat2 latitude of point 2 (degrees)
#' @param lon2 longitude of point 2 (degrees) #' @param lon2 longitude of point 2 (degrees)
#' @keywords #' @keywords
#' @export
#' @examples #' @examples
#' angdist #' angdist
angdist<-function(lat1,lon1,lat2,lon2){ angdist<-function(lat1,lon1,lat2,lon2){
...@@ -89,7 +30,6 @@ angdist<-function(lat1,lon1,lat2,lon2){ ...@@ -89,7 +30,6 @@ angdist<-function(lat1,lon1,lat2,lon2){
#' #'
#' @param a angle (degrees) #' @param a angle (degrees)
#' @keywords #' @keywords
#' @export
#' @examples #' @examples
#' d2r #' d2r
d2r<-function(a){ d2r<-function(a){
...@@ -104,7 +44,6 @@ d2r<-function(a){ ...@@ -104,7 +44,6 @@ d2r<-function(a){
#' #'
#' @param a angle (radians) #' @param a angle (radians)
#' @keywords #' @keywords
#' @export
#' @examples #' @examples
#' r2d #' r2d
r2d<-function(a){ r2d<-function(a){
...@@ -113,7 +52,79 @@ r2d<-function(a){ ...@@ -113,7 +52,79 @@ r2d<-function(a){
# TODO case t>2 Ma
#' Compute atmospheric pressure at site
#'
#' Returns site pressure in hPa
#'
#' Choice of models :
#'
#' stone2000 (default) : equation 1 of Stone (2000) JGR paper
#'
#' era40 : Looks up mean sea level pressure and mean 1000 mb temp from ERA-40 reanalysis
#' and calculates site atmospheric pressures using these as inputs to the standard atmosphere equation
#'
#'
#' @param alt Altitude (m)
#' @param lon Longitude (m) optional (only for "era40" model)
#' @param lat Latitude (m) optional (only for "era40" model)
#' @param model Model used for computation, one of "stone2000" or "era40"
#' @keywords
#' @export
#' @examples
#' atm_pressure
atm_pressure<-function(alt,lon=NULL,lat=NULL,model="stone2000"){
if (model == "stone2000"){
Ps=1013.25 # sea level pressure (hPa)
Ts=288.15 # sea level temperature (K)
ksi=0.0065 # adiabatic lapse rate (K/m)
gMR=0.03417 # g*M/R (K/m)
#hydrostatic pressure (equation 1 of Stone 2000)
P=Ps*exp(-1*gMR/ksi*(log(Ts) - log(Ts - ksi* alt)))
}
else if (model == "era40"){
if(is.null(lat) | is.null(lon)) {stop("must specify lat and lon for model era40")}
z = alt
# get longitude from 0 to 360
lon_c=lon+(lon<0)*360
# Some constants and definitions
gmr = -0.03417
lr = c(-6.1517E-03,-3.1831E-06,-1.5014E-07,1.8097E-09,1.1791E-10,-6.5359E-14,-9.5209E-15)
# interpolate
# ERA40 packaged in internal data (sysdata.rda)
site_slp=pracma::interp2(ERA40$ERA40lon,ERA40$ERA40lat,ERA40$meanP, lon_c, lat)
site_T=pracma::interp2(ERA40$ERA40lon,ERA40$ERA40lat,ERA40$meanT, lon_c, lat)
# Lifton Lapse Rate Fit to COSPAR CIRA-86 <10 km altitude
dtdz = lr[1] + lr[2]*lat + lr[3]*lat^2 + lr[4]*lat^3 + lr[5]*lat^4 + lr[6]*lat^5 + lr[7]*lat^6
dtdz = -dtdz
P = site_slp * exp( (gmr/dtdz) * ( log(site_T) - log(site_T - (z*dtdz)) ) )
}
else {
stop("model must be stone2000 or era40")
}
return(P)
}
####################################################################### #######################################################################
#' Get Virtual Dipole Moment time series #' Get Virtual Dipole Moment time series
#' #'
...@@ -129,15 +140,19 @@ r2d<-function(a){ ...@@ -129,15 +140,19 @@ r2d<-function(a){
#' get_vdm #' get_vdm
get_vdm <- function(time,model){ get_vdm <- function(time,model){
if (!(model %in% c("musch","glopis","lsd") )) {stop("Argument model must be one of musch, glopis or lsd")} if (!(model %in% c("musch","glopis","lsd") )) {stop("Argument model must be one of musch, glopis or lsd")}
val=0
if (model == "musch") if (model == "musch")
{ {
vdm = approx(GMDB$musch[1,]*1000,GMDB$musch[2,]*1e22,time)$y if(max(time)>max(GMDB$musch[1,]*1000)){val = mean(GMDB$musch[2,GMDB$musch[1,]>1800])*1e22}
vdm = approx(GMDB$musch[1,]*1000,GMDB$musch[2,]*1e22,time,yright = val)$y
} else if (model == "glopis") } else if (model == "glopis")
{ {
vdm = approx(GMDB$glopis[1,]*1000,GMDB$glopis[2,]*1e22,time)$y if(max(time)>max(GMDB$glopis[1,]*1000)){val = mean(GMDB$glopis[2,GMDB$glopis[1,]>1800])*1e22}
vdm = approx(GMDB$glopis[1,]*1000,GMDB$glopis[2,]*1e22,time,yright = val)$y
} else } else
{ {
vdm = approx(GMDB$lsd[1,]*1000,GMDB$lsd[2,]*1e22,time)$y if(max(time)>max(GMDB$lsd[1,]*1000)){val = mean(GMDB$lsd[2,GMDB$lsd[1,]>1800])*1e22}
vdm = approx(GMDB$lsd[1,]*1000,GMDB$lsd[2,]*1e22,time,yright = val)$y
} }
return(vdm) return(vdm)
} }
...@@ -145,28 +160,200 @@ get_vdm <- function(time,model){ ...@@ -145,28 +160,200 @@ get_vdm <- function(time,model){
####################################################################### #######################################################################
#' calculate cutoff rigidity from Virtual Dipole MOment #' calculate cutoff rigidity from Virtual Dipole Moment time series
#' #'
#' #'
#' #'
#' @param vdm Virtual Dipole Moment (A.m^2) #' @param vdm Virtual Dipole Moment (A.m^2)
#' @param lat Latitude (deg) #' @param lat Latitude (deg)
#' @param model Model used to compute Rc from the virtual dipole moment value (one of "elsasser54" or "lifton14")
#' #'
#' @return #' @return
#' @export #' @export
#' #'
#' @examples #' @examples
#' vdm2rc #' vdm2rc
vdm2rc <- function(vdm,lat){ vdm2rc <- function(vdm,lat,model="elsasser54"){
M0 = 7.746e22 M0 = 7.746e22 # A/m2 reference dipole moment (2010)
# mu0 = 4*pi*10^-7 # mu0 = 4*pi*10^-7
# c = 3.0e8 # c = 3.0e8
# r = 6.3712e6 # r = 6.3712e6
# Rc = (mu0*c)/(16*pi*10^9*r^2) * vdm * (cos(d2r(lat)))^4 if (model == "elsasser54"){
Rc = 14.31187 * vdm/M0 * (cos(d2r(lat)))^4 # Rc = (mu0*c)/(16*pi*10^9*r^2) * vdm * (cos(d2r(lat)))^4
Rc = 14.31187 * vdm/M0 * (cos(d2r(lat)))^4
} else if (model == "lifton14") {
dd = c(6.89901,-103.241,522.061,-1152.15,1189.18,-448.004)
Rc = vdm/M0 * (dd[1]*cos(d2r(lat)) +
dd[2]*(cos(d2r(lat)))^2 +
dd[3]*(cos(d2r(lat)))^3 +
dd[4]*(cos(d2r(lat)))^4 +
dd[5]*(cos(d2r(lat)))^5 +
dd[6]*(cos(d2r(lat)))^6)
} else {
stop("model must be one of elsasser54 or lifton14")
}
#
return(Rc)
}
# TODO check flipud
#######################################################################
#' calculate cutoff rigidity from non-dipolar gridded model
#'
#'
#'
#' @param time time (a)
#' @param lat Latitude (deg)
#' @param lon Longitude (deg)
#'
#' @return
#' @export
#'
#' @examples
#' vdm2rc
rc_ndp <- function(time,lat,lon){
# get longitude from 0 to 360
lon_c=lon+(lon<0)*360
tempRc=rep(0,dim(Pal_LSD$TTRc)[3])
for(i in 1:dim(Pal_LSD$TTRc)[3]){
tempRc[i]=pracma::interp2(Pal_LSD$lon_Rc,rev(Pal_LSD$lat_Rc),pracma::flipud(Pal_LSD$TTRc[,,i]), lon_c, lat)
}
Rc=approx(Pal_LSD$t_Rc,tempRc,time,rule=1)$y
return(Rc) return(Rc)
} }
#######################################################################
#' 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)
}
###############################################################################################
#' 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 (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
#' @keywords
#' @export
#' @examples
solv_conc_euler <- function(z,ero,t,C0,p,S,L){
# 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))
# evolution of preexisting concentration
# total = (S[1]*p[1]*exp(-1*z/L[1]) + S[2]*p[2]*exp(-1*z/L[2]) + S[2]*p[3]*exp(-1*z/L[3]))
# f1 = (S[1]*p[1]*exp(-1*z/L[1])) / total # fraction of production attributed to spallation
# f2 = (S[2]*p[2]*exp(-1*z/L[2])) / total
# f3 = (S[2]*p[3]*exp(-1*z/L[3])) / total
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
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 = Cspal + Cstop + Cfast + Cspal_in + Cstop_in + Cfast_in
return(C)
}
#######################################################################
#' Compute uncertainty ellipses for two-nuclides plot
#'
#' Two nuclides A and B plotted as B/A on the Y-axis and A on the X-axis
#' Return a list with given confidence level ellipses coordinates for each pair of concentrations
#'
#' @param Ca Concentration of first nuclide (usually 10Be) (at)
#' @param Ca_e 1-sigma uncertainty on the concentration of first nuclide (at)
#' @param Cb Concentration of second nuclide (usually 26Al) (at)
#' @param Cb_e 1-sigma uncertainty on the concentration of second nuclide (at)
#' @param confidence confidence level (default 0.98)
#'
#' @return
#' @export
#'
#' @examples
ellipse_2nuclides <- function(Ca, Ca_e, Cb, Cb_e,confidence=0.98) {
if (!( (length(Ca)==(length(Ca_e))) & (length(Cb)==(length(Cb_e))) & (length(Ca)==(length(Cb))) )) {stop("All inputs should be of same size")}
n = 10000
ll = list()
for (i in 1:length(Ca)) {
a = rnorm(n, Ca[i], Ca_e[i])
b = rnorm(n, Cb[i], Cb_e[i])
ratio = b / a
cc = cor(ratio, a)
corr = matrix(c(1, cc, cc, 1), nrow = 2)
#
r_e = sqrt((Cb_e[i] / Cb[i])^2 + (Ca_e[i] / Ca[i])^2) * (Cb[i] / Ca[i])
ll[[i]] <-
ellipse::ellipse(
corr,
level = confidence,
scale = c(Ca_e[i], r_e),
centre = c(Ca[i], Cb[i] / Ca[i])
)
}
return(ll)
}
## Version 0.1
### main features
Basic functionalities to compute concentration under varying erosion conditions
### todo
- robust documentation for all functions (with examples)
- a clear test for each function
- hide unecessary functions -> done
- simple function to calculation erosion rate from concentration
- simple function for two nuclides plots
- update eulerian solver -> problem de l'héritage
## Version 0.2
Focus on teaching
### main features
- Shiny apps
- scalings
- muons
- Sato fluxes
- profiles
- make vignette of notebook with examples
## Version 0.5
### main features
36Cl
library(TCNtools)
data1 = read.table("tests/data/data_hilltop.dat",header=T)
data2 = read.table("tests/data/data_basin.dat",header=T)
plot(NA, xlim=c(0,10) , ylim=c(3,7),xaxs="i",yaxs="i",axes=F, xlab=NA, ylab=NA)
axis(side = 1,tcl=-0.25,mgp=c(1,0.25,0),cex.axis=0.8)
mtext(side = 1, line = 1.5,expression('['^10*'Be]'~'('*10^6~'at/g)'),cex=1)
axis(side = 2,tcl=-0.25,mgp=c(1,0.25,0),cex.axis=0.8)
mtext(side = 2, line = 1.5,expression(''^26*'Al/'^10*'Be'),cex=1)
grid()
box()
#
col1 = "red"
el = ellipse_2nuclides(data1$C10/1e6, data1$C10_e/1e6, data1$C26/1e6, data1$C26_e/1e6,confidence=0.68)
for (i in 1:length(el)) {polygon(el[[i]],col=col1,border="grey")}
# error bars
f=1
data1$ratio = data1$C26/data1$C10
data1$ratio_e = sqrt((data1$C10_e/data1$C10)^2 + (data1$C26_e/data1$C26)^2)*data1$ratio
arrows( (data1$C10-f*data1$C10_e)/1e6 ,data1$ratio,(data1$C10+f*data1$C10_e)/1e6,data1$ratio,length=0)
arrows(data1$C10/1e6,data1$ratio-f*data1$ratio_e,data1$C10/1e6,data1$ratio+f*data1$ratio_e,length=0)
#
col1 = "lightblue"
el = ellipse_2nuclides(data2$C10/1e6, data2$C10_e/1e6, data2$C26/1e6, data2$C26_e/1e6,confidence=0.68)
for (i in 1:length(el)) {polygon(el[[i]],col=col1,border="grey")}
library(TCNtools)
Natoms = 2.006e22
# 1A
sigma0 = 0.739e-30
k_neg = 0.704*0.1828*0.00157
consts_1A = as.list(c(Natoms,sigma0,k_neg))
names(consts_1A) <- c("Natoms","sigma0","k_neg")
# 1B
sigma0 = 0.237e-30
k_neg = 0.704*0.1828*0.00192
consts_1B = as.list(c(Natoms,sigma0,k_neg))
names(consts_1B) <- c("Natoms","sigma0","k_neg")
Sphi=445 # TODO verif
#
z = seq(1,1e4,by=100) # g/cm2
lat=0
alt=1000
P = atm_pressure(alt,model="stone2000")
Rc = 14.31187 * 1 * (cos(d2r(lat)))^4
st = scaling_st(P,lat)
# modele 1A
muons1A = mu_model1a(z,P,consts_1A)
muons1A$P_tot =muons1A$P_neg + muons1A$P_fast
# modele 1B
muons1B = mu_model1b(z,P,Rc,Sphi,consts_1B)
muons1B$P_tot =muons1B$P_neg + muons1B$P_fast
#### cosmo parameters definition ####
# standard attenuation parameters for initial eulerian approximation
Lspal=160. #spalation (g/cm2)
Lslow=1500. #stopping muons (g/cm2)
Lfast=4320. #fast muons (g/cm2)
Lambda=c(Lspal,Lslow,Lfast)
#params production and decay 10Be, 26Al (St scaling from Borchers et al., muons from braucher et al. 2011)
prm=data.frame(matrix(NA, nrow=4, ncol=2))
# 10Be
prm[1,1]=4.01 # production spalation slhl (at/g)
prm[2,1]=0.012 # production stopping muons slhl (at/g)
prm[3,1]=0.039 # production fast muons slhl (at/g)
prm[4,1]=log(2)/1.36e6 #decay constant (1/an)
# 26Al
prm[1,2]=27.93 # production spalation slhl (at/g)
prm[2,2]=0.84 # production stopping muons slhl (at/g)
prm[3,2]=0.081 # production fast muons slhl (at/g)
prm[4,2]=log(2)/0.717e6 #decay constant (1/an)
# model 2
P_neg = prm[2,1]*st$Nmuons*exp(-1*z/Lambda[2])
P_fast = prm[3,1]*st$Nmuons*exp(-1*z/Lambda[3])
muons2 = data.frame(z,P_neg,P_fast)
muons2$P_tot =muons2$P_neg + muons2$P_fast