moonlit: Predicting moonlight intensity for given time and location

Documented in calculateMoonlightIntensity

#' Calculate moonlight intensity
#'
#' This function predicts moonlight intensity on the ground for any given place
#' and time, based on the location, position of the moon and number of correction
#' factors
#'
#' @import suncalc
#' @import stats
#' @import graphics
#'
#' @param lat Latitude, numerical decimal
#' @param lon Longitude, numerical decimal
#' @param date Date time as POSIXct with the local time zone. If needed use
#'   as.POSIXct(date, tz=timezone)
#' @param e Extinction coefficient - a single numerical value depending on the
#'   altitude. Average extinction coefficients (magnitude per air mass) are as
#'   follows: At sea level: 0.28; at 500m asl: 0.24; at 1000m asl: 0.21; at
#'   2000m asl: 0.16
#'
#' @return A data frame with the following columns:
#'   \item{night}{Logical, TRUE when sun is below the horizon}
#'   \item{sunAltDegrees}{Solar altitude in degrees}
#'   \item{moonlightModel}{Predicted moonlight illumination, relative to an "average" full moon}
#'   \item{twilightModel}{Predicted twilight illumination in lux}
#'   \item{illumination}{Combined moon and twilight intensity, in lux}
#'   \item{moonPhase}{Lunar phase as a numerical value (\% of moon face illuminated)}
#'
#'
#' @export
#'
#' @examples
#' \dontrun{
#' lat <- 52.2297
#' lon <- 21.0122
#' date <- as.POSIXct("2023-06-15 22:00:00", tz = "Europe/Warsaw")
#' result <- calculateMoonlightIntensity(lat, lon, date, 0.21)
#' }
#'
#'


#library(suncalc)
#library(stats)
#library(graphics)

calculateMoonlightIntensity <- function(lat, lon, date, e)


{





  #extiction coefficient based on the user-assigned value
  extCoef <- e




  #suncalc library is needed to calculate sun and moon positions
  #library(suncalc)



  #creating a dataframe to calculate sun and moon positions
  d1 <- data.frame(date, lat, lon)


  #calculating values using suncalc library
  moonPosition <- suncalc::getMoonPosition(data = d1)
  moonIllum <- suncalc::getMoonIllumination(date = date)
  phaseAngle <- acos((2*moonIllum$fraction)-1)
  sunPosition <- suncalc::getSunlightPosition(data = d1)



  #using one of the outputs to start building a table for the results
  night<- moonIllum


  #integrating into the dataset


  #calculated values
  #night$time <- time
  night$moonAlt <- moonPosition$altitude
  night$sunAlt <- sunPosition$altitude
  night$moonAltDegrees <- (180/pi)*night$moonAlt
  night$sunAltDegrees <- (180/pi)*night$sunAlt
  night$moonIllum <- moonIllum$fraction
  night$moonIllumCorrected <- (night$moonAlt > 0) * night$moonIllum
  night$phase <- moonIllum$phase
  night$distance <- moonPosition$distance

#adding a "night" field - below 6 degrees equals to nautical twilight and darker
night$night <- night$sunAltDegrees < (0)

#subsetting the data for night only
#night <- subset(night, night == TRUE)

################################################################
### Modelling moonlight intensity on the ground ################
################################################################

# Model provides an iluminance value for a given location and time. Outcome of
# the model is a scalar quantity and represents a proportion of reference
# illumination value. The reference value, in this case, is luminance of full
# moon at 384400 km (average distance) in zenith. Model includes 4 variables,
# each providing a correction factor of 0-1: 1) visual extinction 2) distance to
# the moon relative to mean distance 3) phase effect 4) Angle of incidence of
# moonlight Final value is obtained by multiplying all 4 values.

############################## 1) Visual extinction of the stars depending on
###the angle of the atmosphere equations provided in: Green, Daniel W. E. 1992.
###"Magnitude Corrections for Atmospheric Extinction." International Comet
###Quarterly 14: 55-59 Average extinction coefficients (magnitude per air mass)
###are as follows:
#At sea level: 0.28
#At 500m asl: 0.24
#at 1000m asl: 0.21
#at 2000m asl: 0.16
#Default value in the model is for 1000m asl and, if needed,
###should be changed below: ############################

#extCoef <- 0.21

#Calculating thikness of the atmosphere (airmass value) at given elevation. Note
#that elevation is above the horizon, so we need to use 90-elevation from
#suncalc model. Might show false results for moon altitude <0, but these always
#have value 0 in end model so it doesn't matter. moon altitude must be in
#radians!

night$airMass <- 1/(cos(pi/2-night$moonAlt)+0.025*2.71827^(-11*cos(pi/2-night$moonAlt)))


#difference in magnitudes - in zenith the absortion reduces brightess by 0.24
#magnitude, or ~20%, a minimal value possible. We will calculate the difference
#between extinction at given angle and extinction at zenith and transform it to
#proportion

night$airMassMagnitude <- (night$airMass-1)*extCoef

#proportion
night$atmExtinction <- 10^(-0.4*night$airMassMagnitude)



#test plot angle vs magnitude and angle vs extinction
#plot(night$moonAltDegrees, night$airMassMagnitude, xlim = c(0,90),  ylim = c(0, 10))
#plot(night$moonAltDegrees, night$atmExtinction, xlim = c(0,90),  ylim = c(0, 1))


##########################################
####### 2) Distance to the moon ##########
##########################################

### the illumination is proportional to 1/(distance^2) (inverse-square law,
### Meeus 1991), therefore relative illumination:
night$distanceIllum = 1/(night$distance/384400)^2


###########################################
###### 3) phase effect ####################
###########################################

###Relationship between phase angle (Sun-Moon-observer) and normalised, disc-integrated brightness of the moon
###Data extracted from a plot and interpolation function created for further use in the model
###Source: Buratti, Bonnie J., John K. Hillier, and Michael Wang. 1996. "The
###Lunar Opposition Surge: Observations by Clementine." Icarus 124 (December):
###490-99. https://doi.org/10.1006/icar.1996.0225.

#Data points
#phase angle in degrees
phaseAngle <- c(0.367186581,0.561328356,0.755753776,1.017911051,1.343778816,1.638020082,1.9321973,2.201542706,2.446200869,2.690859033,2.960111111,3.304228275,3.798169858,4.566540124,5.659518624,6.926950784,8.294080819,9.760921815,11.2775818,12.86918679,14.51056309,16.15193939,17.79332234,19.43470363,21.07608991,22.71746788,24.55846678,26.59940044,28.6403448,30.68133202,32.72256159,34.76382063,36.805081,38.84640163,40.88817885,42.93016095,44.97243897,47.01472501,49.05701374,51.09936005,53.14188578,55.18442088,57.22695599,59.2695018,61.3123181,63.35526696,65.39821581,67.44116333,69.48411219,71.52711193,73.57035938,75.61364299,77.65692793,79.7002102,81.74354871,83.78698764,85.8304172,87.87385613,89.91732854,91.961129,94.00502721,96.04892543,98.09282096,100.1367821,102.1807767,104.2247874,106.2687927,108.3127874,110.3568141,112.4010042,114.4452425,116.4894862,118.5337607,120.5781851,122.6226792,124.6673768,126.7121119,128.7568778,130.8016611,132.8463627,134.8910482,136.9357699,138.9805559,141.0253124,143.0700395,145.1147611,147.1594494,149.2040787,151.2487227,153.2933882,155.3381112,157.382853,159.4275827,161.4723968,163.517239,165.5620973,167.6069529,169.6517964,171.6966668,173.7414943,175.786378,177.8311961,179.5020253
)

#converting phase angle to % illuminated
discIlluminated <- (1+ cos(phaseAngle*pi/180))/2

#disc-integrated brighness
phaseBrightness <- c(0.981495601,0.957317682,0.93441936,0.911719223,0.887147319,0.864543232,0.841650203,0.819234668,0.797948808,0.776662948,0.753826384,0.731224944,0.709518906,0.688329272,0.669028531,0.649231859,0.629197232,0.608983672,0.588516023,0.568644051,0.54830288,0.52796171,0.50765056,0.487331904,0.467035764,0.446702099,0.426908396,0.409070349,0.391280626,0.373684202,0.357181123,0.340810937,0.324446792,0.308354473,0.294321992,0.281213721,0.269440419,0.25770336,0.245978382,0.234513149,0.223857355,0.213243844,0.202630334,0.192065148,0.18272016,0.17397319,0.165226219,0.156473208,0.147726238,0.139208809,0.131808889,0.124572064,0.11734128,0.110098415,0.103109254,0.096573137,0.089994736,0.083458618,0.077073516,0.072168356,0.067704159,0.063239962,0.058763683,0.054583393,0.050554118,0.04659733,0.042616379,0.038587104,0.034702802,0.031555452,0.028625563,0.025719836,0.022953043,0.020862794,0.019086656,0.018228687,0.017539854,0.016989954,0.016518582,0.015678735,0.0147664,0.014017161,0.013557871,0.012965687,0.01224061,0.011491371,0.010591118,0.009425079,0.008325486,0.007322543,0.006579345,0.005920715,0.005207719,0.004875281,0.004669694,0.004536595,0.004391414,0.004191868,0.004113134,0.003841101,0.003822773,0.003508456,0.003536722
)

#plot(discIlluminated, phaseBrightness, type = "l")
#lines(discIlluminated, discIlluminated)


#spline interpolation function
brightness <- stats::splinefun(discIlluminated, phaseBrightness, method = "fmm", ties = mean)

#correction factor for model
night$phaseEffect <- brightness(night$moonIllum)

#plot(night$moonIllum, night$phaseEffect, type = "p")



########################################
###### 4) angle of incidence ###########
########################################

### correction factor is sine of the angle of moon's altitude (Austin et al. 1976)
night$incidence <- sin(night$moonAlt)

#test plot
#plot(night$date, night$incidence, ylim = c(0,1))


######################################################
##### total correction factor #######################

# assign values of 0 if moon not visible then multiply
# variables 1:4

# turns out there are some occasional NAs because atmExtinction is
# infinite, corrected for that by making a conditional version of
# the model output

night$illuminationModelComponents <- (night$moonAlt > 0 ) * night$atmExtinction * night$distanceIllum * night$phaseEffect * night$incidence
night$moonlightModel <- ifelse((night$moonAlt<=0), 0,  night$illuminationModelComponents)




#Alternative version of the model, assigning 0 when sun above the horizon. This
#is not really needed as there is a separate column night=TRUE/FALSE


# night$illuminationModelComponents <- (night$night > 0) * (night$moonAlt > 0 ) * night$atmExtinction * night$distanceIllum * night$phaseEffect * night$incidence
# night$moonlightModel <- ifelse((night$moonAlt<=0 | night$night <=0), 0,  night$illuminationModelComponents)






#############################################################
##############################################################
#### Twilight illumination from empirical data https://www.jstor.org/stable/44612241

####  Sun Position and Twilight Times for Driver Visibility Assessment; Duane D.
####  MacInnis, Peter B. Williamson and Geoffrey P. Nielsen; SAE Transactions;
####  Vol. 104, Section 6: JOURNAL OF PASSENGER CARS: Part 1 (1995), pp. 759-783

####  Values calulated from the paper for angles -18 to 5 degrees, by 0.25 degree, imported to R and spline function applied
##############################################################


solarAngle <- c(5,4.75,4.5,4.25,4,3.75,3.5,3.25,3,2.75,2.5,2.25,2,1.75,1.5,1.25,1,0.75,0.5,0.25,0,-0.25,-0.5,-0.75,-1,-1.25,-1.5,-1.75,-2,-2.25,-2.5,-2.75,-3,-3.25,-3.5,-3.75,-4,-4.25,-4.5,-4.75,-5,-5.25,-5.5,-5.75,-6,-6.25,-6.5,-6.75,-7,-7.25,-7.5,-7.75,-8,-8.25,-8.5,-8.75,-9,-9.25,-9.5,-9.75,-10,-10.25,-10.5,-10.75,-11,-11.25,-11.5,-11.75,-12,-12.25,-12.5,-12.75,-13,-13.25,-13.5,-13.75,-14,-14.25,-14.5,-14.75,-15,-15.25,-15.5,-15.75,-16,-16.25,-16.5,-16.75,-17,-17.25,-17.5,-17.75,-18)
twilightIllumination <- c(4499.368645,4251.166227,4010.398231,3776.215801,3547.994929,3325.316934,3107.949504,2895.827638,2689.034018,2487.778469,2292.376367,2103.225979,1920.784922,1745.546056,1578.01326,1418.677674,1267.99505,1126.364901,994.1121318,871.4718063,758.577575,655.4541802,562.0142693,478.0595544,404.556034,338.5269901,281.0455199,231.5470057,189.3612856,153.7599131,123.995336,99.33178815,79.06826738,62.55436811,49.1999767,38.47993901,29.93480689,23.16869331,17.84514167,13.68176809,10.44428227,7.911259933,5.993606008,4.527348005,3.411591605,2.566090598,1.927671576,1.447053311,1.086102387,0.815522737,0.612949604,0.461403805,0.348056383,0.263253638,0.199755997,0.152149406,0.116393841,0.089479462,0.069166322,0.053788354,0.04210643,0.033198579,0.026378238,0.02113351,0.017082103,0.01393796,0.01148655,0.009566583,0.008056477,0.006931022,0.00596833,0.00519104,0.004556611,0.004033257,0.003596965,0.003229411,0.002916468,0.002647147,0.002412827,0.002206692,0.002023321,0.001858385,0.001708419,0.001570657,0.001442901,0.001323428,0.001210913,0.001104369,0.001003101,0.000906659,0.000814804,0.000727464,0.000644703)

#plot(solarAngle, twilightIllumination)

twilightIlluminationFunction <- splinefun(x = solarAngle, y = twilightIllumination)



night$twilightModel <- twilightIlluminationFunction(night$sunAltDegrees)

#the model can only be interpolated to -18 degrees, so need to assign 0 to everything below -18 degrees
night$twilight <- night$sunAltDegrees > -18
night$twilightModel <- night$twilightModel*night$twilight


#plot(night$sunAltDegrees, night$twilightModel)

#summary(night$twilightModel)



##############################################################
### Combining twilight illumination with lunar illumination ###

night$illumination <- night$moonlightModel*0.32 + night$twilightModel

#plot(night$date, night$illuminationModel, type="l")
#plot(night$date, night$illumination, type="l")




#generating output data frame
d1$night <- night$night
d1$sunAltDegrees <- night$sunAltDegrees
d1$moonAltDegrees <- night$moonAltDegrees
d1$moonlightModel <- night$moonlightModel
d1$twilightModel <- night$twilightModel
d1$illumination <- night$illumination
d1$moonPhase <- night$moonIllum



#test plot
# night<- night[order(date),]
# plot(d1$date, d1$moonPhase)
# lines(d1$date, d1$moonlightModel)


return(d1)

}


#test <- calculateMoonlightIntensity(lat, lon, date, 0.26)



######################################################
###### Extinction coefficient - in development #######
######################################################

##### Calculating an extinction coefficient based on the observer's elevation above sea level
# extinction coefficient based on the
# Peak lunar irradiance is around 560 nm (Veilleux & Cummings, 2012  https://doi.org/10.1242/jeb.071415)
# Aerosol Optical Depth (AOD) is assumed to be 0.15, an average value

#' Calculate extinction coefficient based on elevation of the observer
#'
#' @param elev elevation in meters asl
#'
#' @return elevExtCoeff
#' @export
#'
#' @examples
#'\dontrun{
#'
#' elev <- 2228 # Elevation in m asl for Mount Kosciuszko
#' result <- elevExtCoeff(elev)
#' }
#'

elevExtCoeff <- function (elev) {

  e <- 0.1451*exp(-elev/7995)+0.136

  return(e)

}