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

update shinyapp

parent b4b90094
......@@ -18,4 +18,5 @@ export(scaling_st)
export(solar_modulation)
export(solv_conc_eul)
export(solv_conc_lag)
export(solv_ss_eros_eul)
export(vdm2rc)
......@@ -36,7 +36,7 @@ We should the define the basic parameters we are going to use for the computatio
- neutrons for spallation reactions $\Lambda_{spal}$
- stopping muons $\Lambda_{stop}$
- fast muons $\Lambda_{fast}$
- a matrix with the SLHL production rates (in at/g/a), in this case for the *st* scaling scheme (@stone2000air), and decay constant $\lambda$ (in 1/a) for the nuclides of interest.
- a matrix with the SLHL production rates (in at/g/a), in this case for the *st* scaling scheme (@stone2000air), and decay constant $\lambda$ (in a$^{-1}$) for the nuclides of interest.
```{r}
# Attenuation lengths
......
This diff is collapsed.
---
title: 'R Notebook : introducing the Eulerian solver for concentration calculations'
output: html_notebook
bibliography: '`r system.file("tcntools.bib",package = "TCNtools") `'
output:
html_document:
df_print: paged
---
# R Notebook usage
......
This diff is collapsed.
library("shiny")
library("shinyWidgets")
library("TCNtools")
# Define UI for app that draws a histogram ----
ui <- fluidPage(
# img(src='amu.png', align = "right", width=100),
# App title ----
titlePanel("Transient response"),
# Sidebar layout with input and output definitions ----
sidebarLayout(
# Sidebar panel for inputs ----
sidebarPanel(
# input radio button
radioButtons(inputId = "N1",
label = "Nuclide 1:",
inline = TRUE,
choices = c("10Be" = "Be10",
"26Al" = "Al26",
"14C" = "C14"),
selected = "Be10"),
# input radio button
radioButtons(inputId = "N2",
label = "Nuclide 2:",
inline = TRUE,
choices = c("10Be" = "Be10",
"26Al" = "Al26",
"14C" = "C14"),
selected = "C14"),
sliderInput(inputId = "lat",
label = "Latitude (°):",
min = 0,
max = 90,
step = 1,
value = 40),
sliderInput(inputId = "alt",
label = "Altitude (m):",
min = 0,
max = 5000,
step = 100,
value = 2000),
sliderTextInput("tmax",
label = "Duration (ka):",
choices=c(10,20,50,100,200,500,1000,2000),
selected=50, grid = T),
sliderTextInput("ero",
label = "Initial denudation rate (m/Ma):",
choices=c(1,2,5,10,20,50,100,200,500),
selected=50, grid = T),
sliderInput(inputId = "fact",
label = "Factor of variation :",
min = 1,
max = 10,
step = 1,
value = 2),
radioButtons(inputId = "variation",
label = "Type of variation :",
inline = TRUE,
choices = c("Increase" = 1,
"Decrease" = -1))
),
# Main panel for displaying outputs ----
mainPanel(
# Output: ----
plotOutput(outputId = "distPlot"),
# plotOutput(outputId = "distPlot2")
)
)
)
# Define server logic required to draw a histogram ----
server <- function(input, output) {
output$distPlot <- renderPlot({
#unité longeurs: cm
#unité poids: g
#unité temps: an
# Attenuation lengths
Lambda = c(160,1500,4320) # g/cm2
names(Lambda) <- c("Lspal","Lstop","Lfast") # we just give names to the element of the vector
# Production and decay parameters
prm = matrix(c( 4.01 , 0.012 , 0.039 , log(2)/1.36e6,
27.93 , 0.84 , 0.081 , log(2)/0.717e6,
12.3 , 3.31 , 0 ,log(2)/5730),
nrow = 4,ncol=3 )
colnames(prm) <- c("Be10","Al26","C14") # we just give names to the columns of the matrix
rownames(prm) <- c("Pspal","Pstop","Pfast","lambda") # we just give names to the rows of the matrix
# material density
rho = 2.7
col1 = "blue"
col2 = "red"
#récupération des variables
N1 = input$N1
N2 = input$N2
altitude = input$alt
latitude = input$lat
tmax = input$tmax*1e3
ero = input$ero*100/1e6*rho # m/Ma -> g/cm2/a
fact = input$fact^as.numeric(input$variation)
time = seq(0,tmax,length.out = 100)
P = atm_pressure(alt=altitude,model="stone2000") # compute atmospheric pressure at site
st = scaling_st(P,latitude) # compute the scaling parameters according to Stone (2000)
layout(matrix(c(1,2), nrow = 1, ncol = 2, byrow = TRUE), widths =c(1,1))
C1 = solv_conc_eul(0,ero*fact,time,0,prm[,N1],st,Lambda,in_ero=ero) # compute concentration N1
C2 = solv_conc_eul(0,ero*fact,time,0,prm[,N2],st,Lambda,in_ero=ero) # compute concentration N2
#
plot(time/1e3,C1/1e3,ylim=range(C1,C2)/1000,type="l",col=col1,lwd=2,xlab="Time (ka)",ylab="Concentration (x1000 at/g)")
grid()
lines(time/1e3,C2/1e3,col=col2,lwd=2)
legend("topright",c(N1,N2),col=c(col1,col2),lty=1)
ero1 = rep(0,length(C1))
ero2 = rep(0,length(C2))
#
for (i in 1:length(C1)){
ero1[i] = solv_ss_eros_eul(C1[i],prm[,N1],st,Lambda)$root/100*1e6/rho
ero2[i] = solv_ss_eros_eul(C2[i],prm[,N2],st,Lambda)$root/100*1e6/rho
}
plot(time/1e3,ero1,ylim=range(ero1,ero2,c(ero,ero*fact)/100*1e6/rho),type="l",col=col1,lwd=2,xlab="Time (ka)",ylab="Denudation rate (m/Ma)")
grid()
abline(h=ero/100*1e6/rho)
abline(h=ero*fact/100*1e6/rho,lty=2)
lines(time/1e3,ero2,col=col2,lwd=2)
})
}
# Create Shiny app ----
shinyApp(ui = ui, server = server)
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/solvers.R
\name{fun_solv_ss_eros_eul}
\alias{fun_solv_ss_eros_eul}
\title{Function to be solved by solv_ss_eros_eul when method="single"}
\usage{
fun_solv_ss_eros_eul(ero, Cobs, p, S, L)
}
\arguments{
\item{L}{}
}
\description{
Function to be solved by solv_ss_eros_eul when method="single"
}
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/solvers.R
\name{likelihood1}
\alias{likelihood1}
\title{Likelihood function : P(obs|params) for solv_ss_eros_eul when method="MCMC"}
\usage{
likelihood1(ero, p, S, L, Cobs, Cobs_e)
}
\description{
Likelihood function : P(obs|params) for solv_ss_eros_eul when method="MCMC"
}
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/solvers.R
\name{posterior1}
\alias{posterior1}
\title{Un-normalized posterior : numerator of Bayes formula : P(obs|params)*P(params) for solv_ss_eros_eul when method="MCMC"}
\usage{
posterior1(ero, p, S, L, Cobs, Cobs_e, binf, bsup)
}
\description{
Un-normalized posterior : numerator of Bayes formula : P(obs|params)*P(params) for solv_ss_eros_eul when method="MCMC"
}
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/solvers.R
\name{prior1}
\alias{prior1}
\title{Prior distribution P(params) for solv_ss_eros_eul when method="MCMC"}
\usage{
prior1(ero, binf, bsup)
}
\description{
Prior distribution P(params) for solv_ss_eros_eul when method="MCMC"
}
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/solvers.R
\name{proposal1}
\alias{proposal1}
\title{Proposal function for solv_ss_eros_eul when method="MCMC"}
\usage{
proposal1(ero, binf, bsup)
}
\description{
Proposal function for solv_ss_eros_eul when method="MCMC"
}
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/solvers.R
\name{run_mcmc1}
\alias{run_mcmc1}
\title{RUN MCMC function for solv_ss_eros_eul when method="MCMC"}
\usage{
run_mcmc1(p, S, L, Cobs, Cobs_e, binf, bsup, n)
}
\description{
RUN MCMC function for solv_ss_eros_eul when method="MCMC"
}
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/solvers.R
\name{solv_conc_lag}
\alias{solv_conc_lag}
\title{Compute concentration evolution along a time-depth history
Lagrangian formulation}
\usage{
solv_conc_lag(t, z, C0, Psp0, Pmu0, lambda, S, final = FALSE)
}
\arguments{
\item{t}{time vector (a)}
\item{z}{depth vector (g/cm2), same length as t}
\item{C0}{initial concentration (at/g)}
\item{Psp0}{SLHL spallation production profile (at/g/a) at depth corresponding to z}
\item{Pmu0}{muonic production profil at depth corresponding to z}
\item{lambda}{radioactive decay constant (1/a)}
\item{S}{scaling parameters for PsP0 and Pmu0 (2 columns), at times corresponding to t (same number of rows)}
\item{final}{if TRUE only compute the final concentration (default=FALSE)}
}
\value{
}
\description{
Compute concentration evolution along a time-depth history
Lagrangian formulation
}
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/solvers.R
\name{solv_ss_eros_eul}
\alias{solv_ss_eros_eul}
\title{Compute steady state denudation
Eulerian formulation}
\usage{
solv_ss_eros_eul(Cobs, p, S, L, Cobs_e = 0, method = "single", n = 10000)
}
\arguments{
\item{Cobs}{Uncertainty on measured concentration (at/g), optional but required ffor method MC and MCMC}
\item{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)}
\item{S}{scaling factors (2 elements vector)
S[1] -> scaling factor for spallation
S[2] -> scaling factor for muons}
\item{L}{attenuation length (3 elements vector)
L[1] -> neutrons
L[2] -> stopped muons
L[3] -> fast muons}
\item{method}{One of single (default), MC, MCMC}
}
\value{
}
\description{
Compute steady state denudation
Eulerian formulation
}
......@@ -70,86 +70,6 @@ rho=2.6
solv_ss_eros_eul <- function(Cobs,p,S,L,Cobs_e=0,method="single",n=10000){
Cmax = solv_conc_euler(0,0,4.55e9,0,p,S,L)
if (Cobs>Cmax){stop("Observed concentration higher than theoretical maximum")}
if (!(method %in% c("single","MC","MCMC"))){stop("Parameter method should be one of single, MC or MCMC")}
if (method!="single" & Cobs_e==0){stop("Need to define the error on the observed concentration")}
if (method == "single") {
res = uniroot(fun_solv_ss_eros_eul,c(0,0.1),Cobs,p,S,L)
} else if (method == "MC") {
C0 = rtruncnorm(n,a=0,mean=Cobs,sd=Cobs_e)
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_ss_eros_eul,c(0,1),C0[i],p,S,L,tol=1e-6) # attention à la borne sup de l'interval que le taux d'érosion max soit assez rapide pour les cocentrations faibles
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)
} else if (method == "MCMC") {
ero0 = uniroot(fun_solv_ss_eros_eul,c(0,0.1),Cobs,p,S,L)$root
binf = max(0,ero0*(1-10*Cobs_e/Cobs))
bsup = ero0*(1+10*Cobs_e/Cobs)
res = run_mcmc1(p,S,L,Cobs,Cobs_e,binf,bsup,n)
}
return(res)
}
fun_solv_ss_eros_eul <- function(ero,Cobs,p,S,L){
C = solv_conc_euler(0,ero,4.55e9,0,p,S,L)
return((C-Cobs)/Cobs)
}
##### Likelihood function : P(obs|params) ##########
likelihood1 <- function(ero,p,S,L,Cobs,Cobs_e){
C = solv_conc_euler(0,ero,4.55e9,0,p,S,L)
chi2 = (C-Cobs)^2/Cobs_e^2
res = (-1*chi2/2) + log(1/(sqrt(2*pi)*Cobs_e))
return(res)
}
##### Prior distribution P(params) ##############
prior1 <- function(ero,binf,bsup){
return(dunif(ero, min=binf, max=bsup, log = T)) # log prior
}
###### Un-normalized posterior : numerator of Bayes formula : P(obs|params)*P(params) ####
##-> likelihood*prior (but we work with log)
posterior1 <- function(ero,p,S,L,Cobs,Cobs_e,binf,bsup){
return (likelihood1(ero,p,S,L,Cobs,Cobs_e) + prior1(ero,binf,bsup)) # log posterior
}
##### Proposal
proposal1 <- function(ero,binf,bsup){
return(rtruncnorm(1,a=binf,b=bsup,mean=ero,sd=(bsup-binf)/20))
}
###### RUN MCMC function #######################
run_mcmc1 <- function(p,S,L,Cobs,Cobs_e,binf,bsup,n){
# start chain
chain = as.data.frame(matrix(NA,nrow=n,ncol=3))
colnames(chain) <- c("ero","posterior","probab")
chain$ero[1] = proposal1(runif(1,binf,bsup),binf,bsup)
chain$posterior[1] = posterior1(chain$ero[1],p,S,L,Cobs,Cobs_e,binf,bsup)
for (i in 1:(n-1)){
# if (i%%(round(iter*1/100)+1)==0) {cat(round(i/iter*100),"% -")}
proposal = proposal1(chain$ero[i],binf,bsup)
pst = posterior1(proposal,p,S,L,Cobs,Cobs_e,binf,bsup)
probab = exp(pst - chain$posterior[i]) # ratio of posterior probability proposal/current
chain$probab[i]=probab
if (runif(1) < probab){
chain$ero[i+1] = proposal
chain$posterior[i+1] = pst
}else{
chain$ero[i+1] = chain$ero[i]
chain$posterior[i+1] = chain$posterior[i]
}
} # end loop
return(chain)
} # end function run_mcmc
......
......@@ -30,13 +30,13 @@ S=c(2.3,1.5)
# plot banana line ss erosion
ero = 10^seq(log10(0.0001),log10(2000),length.out = 100)/1.e6*100*rho # conversion en m/Ma -> g/cm2/a
C10Be = solv_conc_euler(0,ero,4.55e9,0,prm[,1],S,Lambda)
C26Al = solv_conc_euler(0,ero,4.55e9,0,prm[,2],S,Lambda)
C10Be = solv_conc_eul(0,ero,4.55e9,0,prm[,1],S,Lambda)
C26Al = solv_conc_eul(0,ero,4.55e9,0,prm[,2],S,Lambda)
ss_erosion = as.data.frame(cbind(ero,C10Be,C26Al))
# plot banana line ss exposure
time = 10^seq(log10(100),log10(200e6),length.out = 100) # a
C10Be = solv_conc_euler(0,0,time,0,prm[,1],S,Lambda)
C26Al = solv_conc_euler(0,0,time,0,prm[,2],S,Lambda)
C10Be = solv_conc_eul(0,0,time,0,prm[,1],S,Lambda)
C26Al = solv_conc_eul(0,0,time,0,prm[,2],S,Lambda)
ss_exposure = as.data.frame(cbind(time,C10Be,C26Al))
#
ss = cbind(c(ss_erosion$C10Be,(ss_exposure$C10Be))/1e6 , c(ss_erosion$C26Al/ss_erosion$C10Be,ss_exposure$C26Al/ss_exposure$C10Be))
......@@ -49,8 +49,8 @@ fact = 7
# eulerien
df_e = data.frame(t=seq(0,tmax,length.out = 10000))
df_e$C10=solv_conc_euler(0,ero*fact,df_e$t,0,prm[,1],S,Lambda,in_ero=ero)
df_e$C26=solv_conc_euler(0,ero*fact,df_e$t,0,prm[,2],S,Lambda,in_ero=ero)
df_e$C10=solv_conc_eul(0,ero*fact,df_e$t,0,prm[,1],S,Lambda,in_ero=ero)
df_e$C26=solv_conc_eul(0,ero*fact,df_e$t,0,prm[,2],S,Lambda,in_ero=ero)
# lagrangien
......@@ -58,8 +58,8 @@ n=10000 # nombre d'étapes
df_l = data.frame(t=seq(0,tmax,length.out = n))
df_l$z = ero*fact*df_l$t # g/cm2
df_l$z = max(df_l$z) - df_l$z
C10_0 = solv_conc_euler(max(df_l$z),ero,4.55e9,0,prm[,1],S,Lambda)
C26_0 = solv_conc_euler(max(df_l$z),ero,4.55e9,0,prm[,2],S,Lambda)
C10_0 = solv_conc_eul(max(df_l$z),ero,4.55e9,0,prm[,1],S,Lambda)
C26_0 = solv_conc_eul(max(df_l$z),ero,4.55e9,0,prm[,2],S,Lambda)
#C10_0 = 0
#C26_0 = 0
df_l$Ssp = rep(S[1],n)
......
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