Commit 773e52c3 authored by GODARD Vincent's avatar GODARD Vincent
Browse files

LSD scaling and corresponding tests

parent 30e28449
......@@ -11,6 +11,7 @@ export(sato_neutrons)
export(sato_neutrons_low)
export(sato_protons)
export(scaling_lm)
export(scaling_lsd)
export(scaling_st)
export(solar_modulation)
export(vdm2rc)
......@@ -65,3 +65,122 @@ lsdrc<-function(lat,lon,time){
##############################################################################
#' scaling_lsd
#'
#' Implements the Lifton-Sato-Dunai scaling scheme from Lifton et al. (2014)
#'
#' Compute the scaling paramters for both the Sf (flux based) and Sa (flux ad cross-section based, nuclide dependent) schemes.
#'
#' Lifton, N., Sato, T., & Dunai, T. J. (2014).
#' Scaling in situ cosmogenic nuclide production rates using analytical approximations to atmospheric cosmic-ray fluxes.
#' Earth and Planetary Science Letters, 386, 149–160.
#' https://doi.org/10.1016/j.epsl.2013.10.052
#'
#' Initial Matlab code from Lifton et al. (2014)
#'
#' @param h atmospheric pressure (hPa)
#' @param Rc cutoff rigidity (GV)
#' @param SPhi solar modulation potential (Phi, see source paper)
#' @param w fractional water content of ground (nondimensional)
#'
#' @return
#' @export
#'
#' @examples
scaling_lsd<- function(h,Rc,SPhi,w){
# Written by Nat Lifton 2013, Purdue University
# nlifton@purdue.edu
# Based on code by Greg Balco -- Berkeley Geochronology Lab
# balcs@bgc.org
# April, 2007
# Part of the CRONUS-Earth online calculators:
# http://hess.ess.washington.edu/math
#
# Copyright 2001-2013, University of Washington, Purdue University
# All rights reserved
# Developed in part with funding from the National Science Foundation.
#
# This program is free software you can redistribute it and/or modify
# it under the terms of the GNU General Public License, version 3,
# as published by the Free Software Foundation (www.fsf.org).
# define structures to store results
scaling=data.frame(matrix(NA, nrow=length(Rc), ncol=13))
colnames(scaling)<-c("Rc","sp","He","Be","C","Al","eth","th","muTotal","mn","mp","mnabs","mpabs")
scaling$Rc=Rc
#
Site=list()
if (w < 0) {w = 0.066} # default gravimetric water content for Sato et al. (2008)
# reference muon flux
mfluxRef = Ref_LSD$mfluxRef
muRef = (unlist(mfluxRef$neg) + unlist(mfluxRef$pos))
# reference values for nuclide of interest or flux
HeRef = Ref_LSD$P3nRef + Ref_LSD$P3pRef
BeRef = Ref_LSD$P10nRef + Ref_LSD$P10pRef
CRef = Ref_LSD$P14nRef + Ref_LSD$P14pRef
AlRef = Ref_LSD$P26nRef + Ref_LSD$P26pRef
SpRef = Ref_LSD$nfluxRef + Ref_LSD$pfluxRef # Sato et al. (2008) Reference hadron flux integral >1 MeV
EthRef = Ref_LSD$ethfluxRef
ThRef = Ref_LSD$thfluxRef
# Site nucleon fluxes
NSite = sato_neutrons(h,Rc,SPhi,w)
NlowSite = sato_neutrons_low(h,Rc,SPhi,w)
PSite = sato_protons(h,Rc,SPhi)
# Site omnidirectional muon flux
mflux = sato_muons(h,Rc,SPhi) #Generates muon flux at site from Sato et al. (2008) model
muSite = (mflux$flux_diff_neg + mflux$flux_diff_pos)
#Nuclide-specific scaling factors as f(Rc)
scaling$He = (NSite$scaling$P3n + PSite$scaling$P3p)/HeRef
scaling$Be = (NSite$scaling$P10n + PSite$scaling$P10p)/BeRef
scaling$C = (NSite$scaling$P14n + PSite$scaling$P14p)/CRef
scaling$Al = (NSite$scaling$P26n + PSite$scaling$P26p)/AlRef
scaling$sp = ((NSite$scaling$nflux + PSite$scaling$pflux))/SpRef # Sato et al. (2008) Reference hadron flux integral >1 MeV
#
scaling$eth = NlowSite$flux$ethflux/EthRef #Epithermal neutron flux scaling factor as f(Rc)
scaling$th = NlowSite$flux$thflux/ThRef#Thermal neutron flux scaling factor as f(Rc)
#
#
Site$E = NSite$E #Nucleon flux energy bins
#Differential muon flux scaling factors as f(Energy, Rc)
Site$muE = mflux$E #Muon flux energy bins (in MeV)
Site$mup = mflux$p #Muon flux momentum bins (in MeV/c)
muSF = matrix(rep(0,length(Rc)*length(Site$muE)),nrow=length(Rc),ncol=length(Site$muE))
for (i in 1:length(Rc)){
muSF[i,] = muSite[i,]/muRef
}
Site$muSF=muSF
#Integral muon flux scaling factors as f(Rc)
scaling$muTotal = mflux$flux_int$total/unlist(mfluxRef$total) #Integral total muon flux scaling factor
scaling$mn = mflux$flux_int$neg/unlist(mfluxRef$nint) #Integral neg muon flux scaling factor
scaling$mp = mflux$flux_int$pos/unlist(mfluxRef$pint) #Integral pos muon flux scaling factor
scaling$mnabs = mflux$flux_int$neg #Integral neg muon flux
scaling$mpabs = mflux$flux_int$pos #Integral pos muon flux
Site$scaling=scaling
return(Site)
}
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/scaling_lsd.R
\name{scaling_lsd}
\alias{scaling_lsd}
\title{scaling_lsd}
\usage{
scaling_lsd(h, Rc, SPhi, w)
}
\arguments{
\item{h}{atmospheric pressure (hPa)}
\item{Rc}{cutoff rigidity (GV)}
\item{SPhi}{solar modulation potential (Phi, see source paper)}
\item{w}{fractional water content of ground (nondimensional)}
}
\description{
Implements the Lifton-Sato-Dunai scaling scheme from Lifton et al. (2014)
}
\details{
Compute the scaling paramters for both the Sf (flux based) and Sa (flux ad cross-section based, nuclide dependent) schemes.
Lifton, N., Sato, T., & Dunai, T. J. (2014).
Scaling in situ cosmogenic nuclide production rates using analytical approximations to atmospheric cosmic-ray fluxes.
Earth and Planetary Science Letters, 386, 149–160.
https://doi.org/10.1016/j.epsl.2013.10.052
Initial Matlab code from Lifton et al. (2014)
}
library(TCNtools)
# test st
z = c(1000,2000,3000)
lat = c(-30,20,60)
#
time=seq(0,5e4,length.out = 1000)
lat=30
lon=-100
z=1000
P = atm_pressure(alt=z,model="stone2000")
# test st
#
scal = scaling_st(P,lat)
st = scaling_st(P,lat)
# test lm
time=seq(0,1e5,length.out = 1000)
lat=40
z=2000
P = atm_pressure(alt=z,model="stone2000")
vdm=get_vdm(time,model="glopis")
vdm=get_vdm(time,model="musch")
Rc = vdm2rc(vdm,lat)
lm = scaling_lm(P,Rc)
plot(NA,xlim=range(time),ylim=range(lm),xlab="Time (a BP)",ylab="lm scaling factor")
grid()
lines(time,lm)
# test lsd
S = solar_modulation(time)
Rc = lsdrc(lat,lon,time)
w=0.066
lsd = scaling_lsd(P,Rc,S,w)
avg=10e3
col_st = "black"
col_lm = "gold"
col_sf = "brown2"
col_sa = "deepskyblue2"
# plot
plot(NA,xlim=range(time),ylim=range(lm*4.01,lsd$scaling$sp*4.09,lsd$scaling$Be*3.92),xlab="Time (a BP)",ylab="10Be production rate (at/g/a)")
grid()
abline(h=st[1]*4.01,col=col_st)
lines(time,lsd$scaling$sp*4.09,col=col_sf)
lines(time,lsd$scaling$Be*3.92,col=col_sa)
lines(time,lm*4,col=col_lm)
legend("bottomright",c("st","lm","sf","sa"),lty=1,col=c(col_st,col_lm,col_sf,col_sa))
points(0,mean(lm[time<avg])*4.00,pch=21,bg=col_lm )
points(0,mean(lsd$scaling$Be[time<avg])*3.92,pch=21,bg=col_sa )
points(0,mean(lsd$scaling$sp[time<avg])*4.09,pch=21,bg=col_sf )
lines(c(0,avg),rep(min(lm*4.01,lsd$scaling$sp*4.09,lsd$scaling$Be*3.92),2),lty=2 )
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