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

reorganize the folders

parent d176bfa7
library(rvest)
url = "http://calibration.ice-d.org/cds/2"
cal = read_html(url)
cal2 <- cal %>%
html_node("pre") %>%
html_text()
cal3 <- read.csv(textConnection(cal2),header=FALSE)
# cal2 = html_node(cal,"title")
# cal3 = html_text(cal,cal2)
# cal3 = html_attr(cal,"title")
#l1 = "<!-- begin v3 --><pre>"
#l2 = "</pre><!-- end v3 -->"
# https://inbo.github.io/tutorials/tutorials/spatial_wfs_services/
library(httr)
wfs_octopus <- "http://earth.uow.edu.au/geoserver/wfs?"
properties_of_interest <- c("smpid1",
"smpid2",
"basin",
"be10np")
#
query <- list(service = "WFS",
request = "GetFeature",
version = "1.1.0",
typeName = "be10-denude:crn_int_basins",
propertyname = as.character(paste(properties_of_interest,collapse = ",")),
outputFormat = "csv",
CRS = "EPSG:4326")
result <- GET(wfs_octopus, query = query)
df <- read.csv(textConnection(content(result, 'text')))
library(TCNtools)
library("truncnorm")
#### 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)
rho=2.6
# solv_expo_eul <- function(Cobs,ero,C0,p,S,L){
# # mettre un garde fou pour les concentrations observées au dela de Cmax
# Cmax = solv_conc_euler(0,ero,4.55e9,C0,p,S,L)
# if (Cobs>Cmax){stop("Observed concentration higher than theoretical maximum")}
# a = uniroot(fun_solv_expo_eul,c(0,1e7),Cobs,ero,C0,p,S,L)
# return(a)
# }
#
# fun_solv_expo_eul <- function(t,Cobs,ero,C0,p,S,L){
# C = solv_conc_euler(0,ero,t,C0,p,S,L)
# return((C-Cobs)/Cobs)
# }
# solv_eros_eul <- function(Cobs,p,S,L,method="MC",n=10000){
# if (length(Cobs)==1){
# Cmax = solv_conc_euler(0,0,4.55e9,0,p,S,L)
# if (Cobs>Cmax){stop("Observed concentration higher than theoretical maximum")}
# res = uniroot(fun_solv_eros_eul,c(0,0.1),Cobs,p,S,L)
# } else {
# C = rnorm(n,Cobs[1],Cobs[2])
# ero = rep(NA,n)
# val = rep(NA,n)
# iter = rep(NA,n)
# prec = rep(NA,n)
# for (i in 1:n){
# a = uniroot(fun_solv_eros_eul,c(0,0.1),C[i],p,S,L,tol=1e-6)
# ero[i] = a$root
# val[i] = a$f.root
# iter[i] = a$iter
# prec[i] = a$estim.prec
# }
# res = data.frame(C,ero,val,iter,prec)
# }
# return(res)
# }
#
# fun_solv_eros_eul <- function(ero,Cobs,p,S,L){
# C = solv_conc_euler(0,ero,4.55e9,0,p,S,L)
# return((C-Cobs)/Cobs)
# }
# site and sample characteristics
lat = 20
lon = 120
alt = 200
P = atm_pressure(alt=alt,model="stone2000")
S = as.numeric(scaling_st(P,lat))
C = 0.2e5 # observed concentration at/g
C_e = C*0.1
#######exposure################################################
# t = seq(0,1e7,by=100)
# res = rep(0,length(t))
# for (i in 1:length(res)){
# res[i] = fun_solv_expo_eul(t[i],C,ero,C0,prm[,1],S,Lambda)
# }
#
# plot(t,res,type="l")
# abline(h=0)
# a = solv_expo_eul(C,ero,C0,prm[,1],S,Lambda)
# abline(v=a$root)
######erosion############################################################
iter=100e3
ero = seq(1,200,length.out = 100)*100/1e6*rho # m/Ma -> g/cm2/a
res = rep(0,length(ero))
for (i in 1:length(res)){
res[i] = fun_solv_ss_eros_eul(ero[i],C,prm[,1],S,Lambda)
}
plot(ero/100*1e6/rho,res,type="l",xlim=range(ero/100*1e6/rho))
abline(h=0)
#
a = solv_ss_eros_eul(C,prm[,1],S,Lambda)
err = a$root/100*1e6/rho
err_e = err*C_e/C
abline(v=err)
#
a = solv_ss_eros_eul(C,prm[,1],S,Lambda,C_e,method="MC",n=iter)
par(new=T)
plot(density(a$ero/100*1e6/rho,from=0,to=err*10),xlim=range(ero/100*1e6/rho))
#
a = solv_ss_eros_eul(C,prm[,1],S,Lambda,C_e,method="MCMC",n=iter)
lines(density(a$ero[(iter/5):iter]/100*1e6/rho),col="red")
#
x = seq(0,max(ero)/100*1e6/rho,length.out = 200)
lines(x,dnorm(x,mean=err,sd=err_e),col="green")
"letter" "ero10" "ero10_e" "ero26" "ero26_e" "lat" "lon" "alt" "C10" "C10_e" "C26" "C26_e"
"A" 4.8 0.43 5.41 0.731 -19.4024780935161 -43.4235417094978 1248.50483974316 882487 18933 5073030 154733
"B" 3.77 0.3447 4.03 0.549 -19.402396378158 -43.4175947411075 1276.71478887085 1124297 24458 6695314 188823
"C" 2.7 0.2648 3.01 0.454 -19.3965371753963 -43.4167807130148 1356.96503733099 1608682 38967 9041935 304007
"D" 1.26 0.1231 1.51 0.2371 -19.3891378019168 -43.4119850837232 1306.31758047264 3016943 52052 14965474 349306
"E" 3.08 0.2805 3.27 0.513 -19.3889474871341 -43.4083137042255 1228.35742686175 1304811 26794 7728252 298940
"F" 3.27 0.2938 3.95 0.559 -19.3874026842423 -43.4051680060566 1107.25281832585 1133214 22439 6043069 189478
"G" 1.34 0.1264 1.63 0.241 -19.3847890236156 -43.410518855339 1325.24786329589 2904057 46736 14322484 283208
"H" 2.3 0.2142 2.74 0.381 -19.3796942152289 -43.4097973968121 1319.86044814464 1807069 36266 9508044 233451
"I" 5.71 0.581 3.16 0.458 -19.3761026617733 -43.4083132811029 1184.49898199191 716624 21683 7706193 237219
"J" 1.8 0.1653 2.12 0.3347 -19.3936362173108 -43.4181629188829 1364.94565729397 2317167 40020 12071419 381933
"K" 2.64 0.2433 2.88 0.43 -19.3921852179763 -43.4180829518644 1339.30652595499 1621091 32629 9269024 297218
"letter" "ero10" "ero10_e" "ero26" "ero26_e" "lat" "lon" "alt" "C10" "C10_e" "C26" "C26_e" "cht0" "cht0_e" "cht_pred"
"M" 3.69 0.3303 3.76 0.542 -19.2636683333479 -43.5558033333332 1313 1176054 23955 7269896 235842 0.0265651963522562 0.00307606808833441 0.0321690865354726
"N" 1.84 0.1663 1.67 0.2624 -19.2613983333479 -43.5507833333332 1354 2254182 37199 14327173 368983 0.0206617868444067 0.00191000478387931 0.0151561736563099
"O" 2.37 0.2372 2.2 0.3399 -19.2611966666812 -43.550835 1353 1798477 44533 11588554 351978 0.0492564142134032 0.00143498964029475 0.0197332517405516
"P" 0.834 0.0876 0.789 0.1458 -19.2360750000144 -43.50575 1355 4436623 70831 24975976 429957 0.00346786787670033 0.000778915964449212 0.00700811106672107
"S" 0.614 0.0698 0.559 0.1207 -19.1594066666807 -43.5155716666665 1380 5837516 88727 32576587 529931 0.00742630714768265 0.000913374530932543 0.00506501188001468
"T" 0.35 0.04924 0.366 0.1013 -19.1582316666807 -43.5174 1379 9045913 132287 43197465 715748 0.00786339841143486 0.000440865366805912 0.00309168670595951
"U" 6 0.536 6.05 0.81 -19.220591666681 -43.4988733333332 1337 759349 16776 4867323 149041 0.050617843395582 0.00243856045563953 0.0520318782218047
"V" 2.32 0.2061 2.27 0.3351 -19.221751666681 -43.4967666666665 1301 1770416 30501 10910453 294062 0.0252178722202729 0.00230643051826881 0.0198196117044053
"W" 0.823 0.0901 0.742 0.1416 -19.2336700000144 -43.5050466666665 1357 4495359 82170 26155361 462037 0.015859809810652 0.000547990194318047 0.00675766717154558
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
P_spal = prm[1,1]*st$Nneutrons*exp(-1*z/Lambda[1])
col1a = "pink"
col1b = "blue"
col2 = "red"
#neg
plot(NA,xlim=range(muons1A$P_neg,muons1B$P_neg,muons2$P_neg),ylim=rev(range(z)),log="x",axes=F, xlab=NA, ylab=NA,xaxs="i",yaxs="i")
axis(side = 1,tcl=-0.25,mgp=c(1,0.25,0),cex.axis=0.8)
mtext(side = 1, line = 1.5, "Production rate (at/g/a)",cex=1.2)
axis(side = 2,tcl=-0.25,mgp=c(1,0.25,0),cex.axis=0.8)
mtext(side = 2, line = 1.5, "Depth (g/cm2)",cex=1.2)
mtext(side = 3, line = 1.5, "Negative muons",cex=1.2)
#
lines(muons1A$P_neg,muons1A$z,col=col1a)
lines(muons1B$P_neg,muons1B$z,col=col1b)
lines(muons2$P_neg,muons2$z,col=col2)
#
box()
legend("topleft",c("1A","1B","2"),lty=1,col=c(col1a,col1b,col2),cex=0.8)
#fast
plot(NA,xlim=range(muons1A$P_fast,muons1B$P_fast,muons2$P_fast),ylim=rev(range(z)),log="x",axes=F, xlab=NA, ylab=NA,xaxs="i",yaxs="i")
axis(side = 1,tcl=-0.25,mgp=c(1,0.25,0),cex.axis=0.8)
mtext(side = 1, line = 1.5, "Production rate (at/g/a)",cex=1.2)
axis(side = 2,tcl=-0.25,mgp=c(1,0.25,0),cex.axis=0.8)
mtext(side = 2, line = 1.5, "Depth (g/cm2)",cex=1.2)
mtext(side = 3, line = 1.5, "Fast muons",cex=1.2)
#
lines(muons1A$P_fast,muons1A$z,col=col1a)
lines(muons1B$P_fast,muons1B$z,col=col1b)
lines(muons2$P_fast,muons2$z,col=col2)
#
box()
legend("topleft",c("1A","1B","2"),lty=1,col=c(col1a,col1b,col2),cex=0.8)
#total
plot(NA,xlim=range(muons1A$P_tot,muons1B$P_tot,muons2$P_tot),ylim=rev(range(z)),log="x",axes=F, xlab=NA, ylab=NA,xaxs="i",yaxs="i")
axis(side = 1,tcl=-0.25,mgp=c(1,0.25,0),cex.axis=0.8)
mtext(side = 1, line = 1.5, "Production rate (at/g/a)",cex=1.2)
axis(side = 2,tcl=-0.25,mgp=c(1,0.25,0),cex.axis=0.8)
mtext(side = 2, line = 1.5, "Depth (g/cm2)",cex=1.2)
mtext(side = 3, line = 1.5, "Total muon production",cex=1.2)
#
lines(muons1A$P_tot,muons1A$z,col=col1a)
lines(muons1B$P_tot,muons1B$z,col=col1b)
lines(muons2$P_tot,muons2$z,col=col2)
lines(P_spal,z,col="grey",lty=2)
#
box()
legend("topleft",c("1A","1B","2"),lty=1,col=c(col1a,col1b,col2),cex=0.8)
# plot
plot(NA,xlim=range(muons1A$P_neg,muons1B$P_neg,muons2$P_neg,muons1A$P_fast,muons1B$P_fast,muons2$P_fast,muons1A$P_tot,muons1B$P_tot,muons2$P_tot),
ylim=rev(range(z)),log="x",axes=F, xlab=NA, ylab=NA,xaxs="i",yaxs="i")
axis(side = 1,tcl=-0.25,mgp=c(1,0.25,0),cex.axis=0.8)
mtext(side = 1, line = 1.5, "Production rate (at/g/a)",cex=1.2)
axis(side = 2,tcl=-0.25,mgp=c(1,0.25,0),cex.axis=0.8)
mtext(side = 2, line = 1.5, "Depth (g/cm2)",cex=1.2)
# neg
lines(muons1A$P_neg,muons1A$z,col=col1a,lty=2)
lines(muons1B$P_neg,muons1B$z,col=col1b,lty=2)
lines(muons2$P_neg,muons2$z,col=col2,lty=2)
# fast
lines(muons1A$P_fast,muons1A$z,col=col1a)
lines(muons1B$P_fast,muons1B$z,col=col1b)
lines(muons2$P_fast,muons2$z,col=col2)
# tot
lines(muons1A$P_tot,muons1A$z,col=col1a,lwd=3)
lines(muons1B$P_tot,muons1B$z,col=col1b,lwd=3)
lines(muons2$P_tot,muons2$z,col=col2,lwd=3)
#
box()
legend("topleft",c("1A","1B","2","neg","fast","total"),lty=c(1,1,1,2,1,1),lwd=c(1,1,1,1,1,3),col=c(col1a,col1b,col2,"black","black","black"),cex=0.8)
#
# plot(NA,xlim=range(muons2$P_neg,muons2$P_fast,muons2$P_tot),ylim=rev(range(z)),log="x",axes=F, xlab=NA, ylab=NA,xaxs="i",yaxs="i")
# axis(side = 1,tcl=-0.25,mgp=c(1,0.25,0),cex.axis=0.8)
# mtext(side = 1, line = 1.5, "Production rate (at/g/a)",cex=1.2)
# axis(side = 2,tcl=-0.25,mgp=c(1,0.25,0),cex.axis=0.8)
# mtext(side = 2, line = 1.5, "Depth (g/cm2)",cex=1.2)
# mtext(side = 3, line = 1.5, "Total muon production",cex=1.2)
# #
# lines(muons2$P_neg,muons1A$z,col=col1a)
# lines(muons2$P_fast,muons1B$z,col=col1b)
# lines(muons2$P_tot,muons2$z,col=col2)
# grid()
#
library(TCNtools)
#### test get_vdm ####
time=seq(0,4e6,length.out = 1000)
plot(NA,xlim=range(time),ylim=c(0,16),xlab="Time (a BP)",ylab="VDM (10^22 A.m^2)")
grid()
# 1 - Glopis
col1="coral1"
vdm=get_vdm(time,model="glopis")
lines(time,vdm/1e22,col=col1)
# 2 - Musch
col2 = "chartreuse"
vdm=get_vdm(time,model="musch")
lines(time,vdm/1e22,col=col2,lty=2)
# 3 - lsd
col3 = "cornflowerblue"
vdm=get_vdm(time,model="lsd")
lines(time,vdm/1e22,col=col3)
legend("topright",c("glopis","musch","lsd"),col=c(col1,col2,col3),lty=1)
#### test get_vdm ####
lat=0
lon=0
vdm=get_vdm(time,model="musch")
rc1 = vdm2rc(vdm,lat)
rc2 = vdm2rc(vdm,lat,model="lifton14")
rc3 = rc_ndp(time,lat,lon)
col1 = "blueviolet"
col2 = "darkorange"
plot(NA,xlim=range(time),ylim=range(rc1,rc2,na.omit(rc3)),xlab="Time (a BP)",ylab="Rc (GV)")
grid()
#
lines(time,rc3)
#
t = seq(7e3,max(time),length.out = 1000)
tmp = lsdrc(lat,lon=0,t)
lines(t,tmp,col="grey",lty=2,lwd=3)
#
lines(time,rc1,col=col1)
lines(time,rc2,col=col2)
legend("topright",c("elsasser54","lifton14","LSD>7ka","non dipolar for <7ka"),col=c(col1,col2,"grey","black"),lty=c(1,1,2,1),lwd=c(1,1,3,1))
# #### test lsdrc ####
# time=seq(0,4e4,length.out = 1000)
# lat=-30
# lon=-150
# Rc1 = lsdrc(lat,lon,time)
# vdm=get_vdm(time,model="lsd")
# Rc2 = vdm2rc(vdm,lat)
# #
# plot(NA,xlim=range(time),ylim=c(0,16),xlab="Time (a BP)",ylab="Rc (GV)")
# grid()
# lines(time,Rc1,col="deeppink2")
# lines(time,Rc2,col="dodgerblue3")
# legend("topright",c("LSD (long. variations, 6th order pol. fit)","LSD from VDM"),col=c("deeppink2","dodgerblue3"),lty=1)
# # les M/M0 sont identiques dans les 2 cas différence >7 ka due à l'utilisation de la formule 2 de Lifton et al (2014), vs formule 2 de Martin et al. (2017)
#
# plot(GMDB$lsd[1,]*1000,GMDB$lsd[2,]*1e22,type="l",xlim=c(0,40e3))
# lines(Pal_LSD$t_M,Pal_LSD$MM0*7.746e22,col="red")
borchers = as.data.frame(cbind(
c(3.92,4.01,4.09,4.00,3.69,3.70,4.06),
c(28.54,27.93,28.61,27.93,26.26,26.29,28.72),
c(114.55,118.20,118.64,117.23,122.47,122.75,131.32),
c(12.76,12.24,12.72,12.22,12.49,12.44,13.42)),
row.names = c("Sa","St","Sf","Lm","De","Du","Li"))
colnames(borchers) <- c("Be10","Al26","He3","C14")
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment