Nothing
### This function is a component of astrochron: An R Package for Astrochronology
### Copyright (C) 2016 Stephen R. Meyers
###
###########################################################################
### etp: generate combined eccentricity-tilt-precession curve
### (SRM: March 23, 2012; May 3, 2012; October 11, 2012; April 23, 2013; April 29, 2013;
### May 20, 2013; July 10, 2013; August 15, 2013; April 7, 2014; January 21, 2015;
### July 28, 2016)
###
###########################################################################
etp <- function (tmin=NULL,tmax=NULL,dt=1,eWt=1,oWt=1,pWt=1,esinw=T,solution=NULL,standardize=T,genplot=T,verbose=T)
{
if(verbose) cat("\n----- GENERATING ECCENTRICITY-TILT-PRECESSION MODEL -----\n")
if(is.null(solution))
{
if(is.null(tmin)) tmin=0
if(is.null(tmax)) tmax=1000
if(tmax<=249000)
{
if(verbose) cat("\n * Downloading Laskar et al. (2004) astronomical solution\n")
tempLaskar <- tempfile()
download.file("http://www.geology.wisc.edu/~smeyers/astrochron/la04.txt.bz2",tempLaskar)
if(verbose) cat(" * Decompressing Laskar et al. (2004) astronomical solution\n")
la <- read.table(bzfile(tempLaskar),header=T)
unlink(tempLaskar)
la=data.frame(cbind((1:249001)-1,la))
}
if(tmax>249000)
{
cat("\n**** ERROR: Laskar et al. (2004) solution is not available > 249000 ka.\n")
stop(" TERMINATING NOW!")
}
}
if(!is.null(solution))
{
la=solution
if(is.null(tmin)) tmin=min(la[,1])
if(is.null(tmax)) tmax=max(la[,1])
}
if(verbose) cat("\n * Creating model for interval from",tmin,"to",tmax,"ka \n")
# replace preccession angle with sin(precession angle)
la[2] <- sin(la[2])
# replace precession angle with esin(precession angle)
if(esinw) la[2] <- la[4]*la[2]
# isolate portion of record that is needed
la <- subset(la, (la[1] >= (tmin-1)) & (la[1] <= (tmax+1)) )
# calculate mean values for eccentricity, obliquity, precession
laMean <- colMeans(la)
# calculate std deviations for eccentricity, obliquity, precession
laStdev <- sapply(la,sd)
# standardize eccentricity, tilt, precession
if(standardize)
{
prec <- ( la[2]-laMean[2] ) / laStdev[2]
obl <- ( la[3]-laMean[3] ) / laStdev[3]
ecc <- ( la[4]-laMean[4] ) / laStdev[4]
}
if(!standardize)
{
prec <- la[2]
obl <- la[3]
ecc <- la[4]
}
# now weight each term as desired
prec <- prec * pWt
obl <- obl * oWt
ecc <- ecc * eWt
# combine eccentricity, tilt and precession
combined <- ecc + obl + prec
out <- data.frame( cbind(la[1],combined) )
### then interpolate as needed
xout <- seq(tmin,tmax,by=dt)
### redefine npts
npts <- length(xout)
interp <- approx(out[,1],out[,2],xout,method="linear",n=npts)
interp <- data.frame(interp)
if (genplot)
{
### plots
par(mfrow=c(2,2))
plot(la[,1],la[,4],cex=0.5,xlab="Time (ka BP)",ylab="Eccentricity",main="Eccentricity",type="l")
plot(la[,1],la[,3],cex=0.5,xlab="Time (ka BP)",ylab="Obliquity (radians)",main="Obliquity",type="l")
if(esinw) plot(la[,1],la[,2],cex=0.5,xlab="Time (ka BP)",ylab="Eccentricity*sin(angle)",main="Precession",type="l")
if(!esinw) plot(la[,1],la[,2],cex=0.5,xlab="Time (ka BP)",ylab="sin(angle)",main="Precession",type="l")
plot(interp,cex=0.25,xlab="Time (ka BP)",ylab="Value",main="ETP");lines(interp,col="red")
}
return(interp)
### END function etp
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.