R/FUNCTION-alodfromk_v2.R

In astrochron: A Computational Tool for Astrochronology

### This function is a component of astrochron: An R Package for Astrochronology
### Copyright (C) 2026 Stephen R. Meyers
###
###########################################################################
### function alodfromk - (SRM: September 12, 2025; July 16, 2026)
# Compute a (Earth radii), LOD (hrs) and their uncertainties for a given
# value of k (arcsec/yr) and its uncertainty. Also generate plots for
# figures. Based on the original equations of 
# Walker, J.C.G., Zahnle, K.J., 1986. Lunar nodal tide and distance to the
# Moon during the Precambrian. Nature 320, 600-602. 
# with adjusted constants Kconst and Aconst.
# alodfromk is an R translation of Alberto Malinverno's MATLAB function.
###########################################################################
alodfromk <- function(k,sdevk,genplot=TRUE,verbose=TRUE)
 {
  if(verbose) 
   {
    cat("\n----- alodfromk: Compute Earth-Moon distance and length of day -----\n")
    cat("    from precession frequency, as in Malinverno and Meyers (2024) \n\n")
   }

# error checking 
  if(is.null(k)) stop("\n**** ERROR: k must be defined. TERMINATING NOW!\n")
  if(is.null(sdevk)) stop("\n**** ERROR: sdevk must be defined. TERMINATING NOW!\n")

# uncertainties of K-curve and AM-curve
  sdevkcurve=0.005                              # assumed uncertainty of K-curve (0.5%)
  sdevAMcurve=0.005                             # assumed uncertainty of AM-curve (0.5%)

### define constants (these come from Alberto's astroconst MATLAB function)
# Values of astronomical constants (SI units). Values are from Yoder 1995 except where noted.
# Earth
  ac.reqearth=6378137                           # equatorial radius (m)
  ac.spinearth=7.292115e-5                      # mean spin rate (rad/s)
  ac.lod=24*2*pi/ac.spinearth/86400             # sidereal day (hrs)
  ac.k=50.475838                                # present precession freq. (arcsec/yr), eq. 4.19 of Laskar 2020
# Note: Yoder 1995 gives k0=50.290966.
# Moon
  ac.amoon=384400000                            # present semi-major axis of orbit (m)
# Adjusted constants in Walker and Zahnle (1986) equations
  ac.Kconst=0.358                               # originally 0.465
  ac.Aconst=4.81                                # originally 4.87
# Coefficients of polynomial that gives ar=a(t)/a(0) as a function of 
# log(kr) = log[k(t)/k(0)]
  ac.pcoeffar=c(0.00621404, -0.00060922, -0.217194, 1)
### END: astroconst MATLAB function

  a0=ac.amoon/ac.reqearth                       # Earth-Moon semimajor axis (Earth radii)
  k0=ac.k                                       # present precession frequency k ("/yr)
  lod0=ac.lod                                   # present sidereal LOD (hrs)
  Kconst=ac.Kconst                              # updated constants in Walker and Zahnle 1986 equations
  Aconst=ac.Aconst
# coefficients of polynomial that gives ar=a(t)/a(0) as a function of 
# log(kr) = log[k(t)/k(0)]
  pcoeffar=ac.pcoeffar

# compute a and LOD for given value of k
  kr=k/k0                                       # kr = k(t)/k(0)  
# this inspired by polyval in pracma            
  ar= outer(log(kr), 3:0, "^") %*% pcoeffar     # ar = a(t)/a(0)
  omegar=1+Aconst*(1-sqrt(ar))                  # omegar = omega(t)/omega(0)
  a=a0*ar                                       # a in Earth radii
  lod=lod0/omegar                               # LOD in hrs

# parameters to compute uncertainties of a and LOD
  sdevkr=sdevk/k0                               # uncertainty of kr = k(t)/k(0)
  sdevarlarge=1                                 # large uncertainty of ar for calculation of Ck and CAM
  slopek=kr*(1+Kconst)*3*ar^2/((Kconst*ar^3+1)^2)    # slope of k-curve at ar
  slopeAM=-Aconst/(2*sqrt(ar))                  # slope of AM-curve at ar

# compute covariance matrices Ck and CAM
  vararlarge=sdevarlarge^2
  for(i in 1:2)
   {
    if (i==1)
     {
      slope=slopek
      tausq=(sdevkr^2+sdevkcurve^2)            # conditional variance of k-curve
      }else{
      slope=slopeAM
      tausq=sdevAMcurve^2                      # conditional variance of AM-curve
     }
# covariance matrix for line
    varomegar=(slope^2)*vararlarge+tausq
    sdevomegar=sqrt(varomegar)
    rhosq=(varomegar-tausq)/varomegar
    rho=sign(slope)*sqrt(rhosq)
    C=matrix(data=c(vararlarge, rho*sdevarlarge*sdevomegar, rho*sdevarlarge*sdevomegar, varomegar),nrow=2,ncol=2,byrow=TRUE)
# save covariance matrix
    if(i==1)
     {
      Ck=C
     }else{ 
      CAM=C
     }
   }

# compute covariance matrix Cr for ar and omegar
# % Cr=inv(inv(Ck)+inv(CAM))

# compute covariance matrix C for ar and omegar analytically
  vark=(sdevkr^2+sdevkcurve^2)
  varAM=sdevAMcurve^2
  Cmult=1/(slopek-slopeAM)^2
  offdiag=slopek*varAM+slopeAM*vark
  C=as.double(Cmult)*matrix(data=c(vark+varAM, offdiag, offdiag, slopek^2*varAM+slopeAM^2*vark),nrow=2,ncol=2,byrow=TRUE)

# compute uncertainties of a and LOD
  sdevar=sqrt(C[1,1])
  sdeva=a*sdevar
  sdevomegar=sqrt(C[2,2])
  sdevlod=lod*sdevomegar

# ========= plot =========
  if(genplot)
   {
# plotting parameters
    armin=0.8
    armax=1.02
    omegarmin=0.8
    omegarmax=1.9
    nplot=300

# compute k-curve and AM-curve for plotting
    arplot=seq(from=armin,to=armax,length.out=nplot)
    omegark=kr*((1+Kconst)/(Kconst+arplot^(-3)))       # omegar as a function of kr
    omegarAM=1+Aconst-Aconst*arplot^(1/2)

# set plot bounds for zoomed-in figure
    arminzoom=ar-0.06
    armaxzoom=ar+0.06
    omegarminzoom=omegar-0.2
    omegarmaxzoom=omegar+0.2

    dev.new(height = 5.5, width = 5.5, units = "in")
    par(mfrow=c(1,1))
# plot
    plot(arplot,omegark,col="blue",type="l",lwd=2,xlim=c(armin,armax),ylim=c(omegarmin,omegarmax),xlab="",ylab="")
    lines(arplot,omegarAM,col="green",lwd=2)
    points(ar,omegar,col="red",cex=1,pch=16)
    lines(c(ar,ar),c(omegarmin,omegar),lty=2)
    lines(c(armin,ar),c(omegar,omegar),lty=2)
    text(ar+0.005,1.1*armin,sprintf('%.3f',ar),font=2)
    text(armin+0.005,1.05*omegar,sprintf('%.3f',omegar),font=2)
    mtext("Earth spin rate ratio w(t)/w(0)",side=2,line=2)
    mtext("Lunar distance ratio a(t)/a(0)",side=1,line=2)
    titlestr=sprintf('Earth-Moon Distance = %.2f +/- %.2f Earth radii (2 sigma)',a,2*sdeva)
    mtext(titlestr,side=3,line=2,font=2)
    titlestr=sprintf('Length of Day = %.2f +/- %.2f hrs (2 sigma)',lod,2*sdevlod)
    mtext(titlestr,side=3,line=0.75,font=2)
    legend(x="topright",legend=c('K-curve','AM-curve'),col=c("blue","green"),lty=c(1,1),lwd=c(2,2),cex=0.8,bty="n")

   } 

  out=cbind(a,sdeva,lod,sdevlod)
  colnames(out)=c("a","sdeva","lod","sdevlod")
  return(out)

# end function alodfromk
 }