# NOT RUN {
### Example 1: Determine age/stratigraphic height of single objects
##define deposition rate
my.xdep=seq(0,12,length.out=100)
my.ydep=splinefunH(x=c(0,2,4,6,8,10,12),y=c(1,5,6,1,0.5,1,6),m=c(0,1.5,-0.5,-0.5,0,0.5,0))(my.xdep)
#unit of deposition rate is sediment per time unit (default setting)
usedunit='sediment per time' #unit of deposition rate is sediment per time unit (default setting)
#Plot deposition rate (in stratigraphic height)
plot(my.xdep,my.ydep,type='l',main='Deposition Rate',xlab='Stratigraphic Height',ylab=usedunit,
ylim=c(0,max(my.ydep)))
##at what time was the object found at stratigraphic height 9 deposited?
#using default setting for depositionmodel (depositionmodel = 'piecewise linear deposition rate')
pointtransform(points=9,xdep=my.xdep,ydep=my.ydep,direction='height to time', unit=usedunit)
#change unit used
usedunit='time per sediment'
pointtransform(points=9,xdep=my.xdep,ydep=my.ydep,direction='height to time', unit=usedunit)
#note how different the results are!
##Now, take deposition rate as deposition rate in time
##at what stratigraphic height will an object appear that was deposited at time 5?
pointtransform(points=5,xdep=my.xdep,ydep=my.ydep,direction='time to height')
#The option "unit" is unused when transforming from time to height
### Example 2: Create Age models based on a deposition rate
##create an age model. Essentially transform many points, which then approximate the age model
stratheights=seq(min(my.xdep),max(my.xdep),length.out=1000) #many points to approx. age model
usedunit='sediment per time'
reslist=pointtransform(points=stratheights,xdep=my.xdep,ydep=my.ydep,
direction='height to time',unit=usedunit)
reslist$report
agemodelage=reslist$time
agemodelheight=reslist$height
#plot age model
plot(agemodelage,agemodelheight,xlab='Time',ylab='Stratigraphic Height',
main=paste('Age model based on deposition rate \n with unit',usedunit))
#create age model but with other units for sedimentn input
usedunit='time per sediment'
reslist=pointtransform(points=stratheights,xdep=my.xdep,ydep=my.ydep,
direction='height to time',unit=usedunit)
reslist$report
agemodelage=reslist$time
agemodelheight=reslist$height
#plot age model (note the difference the setting of unit makes in terms of time
#required to deposit the section!)
plot(agemodelage,agemodelheight,xlab='Time',ylab='Stratigraphic Height',
main=paste('Age model based on deposition rate \n with unit',usedunit))
##create age model with a hiatus 1: height to time
stratigraphicheight=5 #strat. height of the hiatus
duration=10 #duration of the hiatus
my.hiatuslist=list(c(stratigraphicheight,duration)) #required input format for hiatuses
reslist=pointtransform(points=stratheights,xdep=my.xdep,ydep=my.ydep,
direction='height to time',hiatuslist=my.hiatuslist)
reslist$report
agemodelage=reslist$time
agemodelheight=reslist$height
#!using default setting for unit (sediment per time) again!
plot(agemodelage,agemodelheight,xlab='Time',ylab='stratigraphic height')
#the gap corresponds to the hiatus
##create age model with a hiatus 2: time to height
my.xdep2=c(0,6,8,12)
my.ydep2=c(0,6,5,12)
plot(my.xdep2,my.ydep2,type='l',main='Age Model, not eroded',xlab='time',ylab='height')
reslist=pointtransform(points=stratheights,xdep=my.xdep2,ydep=my.ydep2,
direction='time to height',depositionmodel='age model')
reslist$report
agemodelage=reslist$time
agemodelheight=reslist$height
plot(agemodelage,agemodelheight,xlab='Time',ylab='stratigraphic height'
,main='Age model, eroded (with hiatus)')
###Example 3: Transform (isotope) ratios
##define deposition rate
my.xdep3=seq(0,12,length.out=100)
my.ydep3=splinefunH(x=c(0,2,4,6,8,10,12),y=c(1,5,6,1,0.5,1,6),
m=c(0,1.5,-0.5,-0.5,0,0.5,0))(my.xdep3)
#create fake (oxygen) isotope curves
samplelocation=seq(0,12,length.out=20) #where the samples are taken
isotoperatio=(-1)^(0:19) +0.2*0:19
plot(my.xdep3,my.ydep3,type='l',ylim=c(0,7),xlab='Stratigraphic Height',ylab='')
lines(samplelocation,isotoperatio,type='l',lwd=4)
legend('topleft',lwd=c(1,4),legend=c('Deposition rate','(Isotope) ratio'))
#transform only (!) sample locations, NOT values
#again using the default setting for unit
reslist=pointtransform(points=samplelocation,xdep=my.xdep3,ydep=my.ydep3,
direction='height to time')
#Isotope ratios in time
plot(reslist$time,isotoperatio,type='l',xlab='Time',ylab='Isotope Ratio',lwd=4)
# }
Run the code above in your browser using DataLab