Implementing faster (?) Fortran version

.Rbuildignore

@@ -24,3 +24,4 @@
 ^\.httr-oauth$
 ^logo\..*$
 ^codemeta\.json$
+^CRAN-SUBMISSION$

.gitignore

@@ -5,3 +5,6 @@
 data-raw/
 data/
 docs/
+CRAN-SUBMISSION
+src/*.o
+src/*.mod

DESCRIPTION

@@ -15,14 +15,14 @@ Maintainer: Koen Hufkens <[email protected]>
 Description: A tool to calculate sky illuminance values (in lux) for both sun and moon. The
  model is a verbatim translation of the code by Janiczek and DeYoung (1987)
  <https://archive.org/details/DTIC_ADA182110>.
-URL: https://github.com/bluegreen-labs/skylight
+URL: https://github.com/bluegreen-labs/skylight, https://bluegreen-labs.github.io/skylight/
 BugReports: https://github.com/bluegreen-labs/skylight/issues
 Depends:
     R (>= 3.6)
 License: AGPL-3
 LazyData: false
 ByteCompile: true
-RoxygenNote: 7.2.1
+RoxygenNote: 7.3.1
 Suggests:
     tidyr,
     rnaturalearth,

NAMESPACE

@@ -1,3 +1,4 @@
 # Generated by roxygen2: do not edit by hand
 
 export(skylight)
+useDynLib(skylight, .registration=TRUE)

R/skylight.R

@@ -30,6 +30,7 @@
 #'  can be considered a scaling factor, substituting it with the (inverse)
 #'  slope parameter of an empirical fit should render more accurate results.
 #'  (this can be a single value or vector of values)
+#' @param fast fast processing (TRUE or FALSE)
 #'
 #' @return Sun and moon illuminance values (in lux), as well as their respective
 #' location in the sky (altitude, azimuth).
@@ -71,7 +72,8 @@ skylight <- function(
     longitude,
     latitude,
     date,
-    sky_condition = 1
+    sky_condition = 1,
+    fast = FALSE
     ){
 
   # pipe friendly function checks
@@ -115,131 +117,160 @@ skylight <- function(
   hour <- as.numeric(format(date,"%H"))
   minutes <- as.numeric(format(date, "%M"))
 
-  # calculate hours as a decimal number
-  hour_dec <- hour + minutes/60
-
-  # constant values
-  RD <- 57.29577951
-  DR <- 1 / RD
-  CE <- 0.91775
-  SE <- 0.39715
-
-  # convert latitude
-  latitude <-  latitude * DR
-
-  J <- 367 * year -
-    as.integer(7 * (year + as.integer((month + 9)/12))/4) +
-    as.integer(275 * month/9) +
-    day - 730531
-
-  E <- hour_dec/24
-  D <- J - 0.5 + E
-
-  #---- calculate solar parameters ----
-  solar_parameters <- sun(
-    D,
-    DR,
-    RD,
-    CE,
-    SE
-  )
+  if (fast) {
+
+    forcing = data.frame(
+        longitude = longitude,
+        latitude = latitude,
+        year = year,
+        month = month,
+        day = day,
+        hour = hour,
+        minutes = minutes,
+        sky_condition = sky_condition
+    )
 
-  # in place adjustments
-  solar_parameters$T <- solar_parameters$T + 360 * E + longitude
-  solar_parameters$H <- solar_parameters$T - solar_parameters$AS
-
-  # calculate celestial body
-  # parameters all these subroutines
-  # need proper clarifications as
-  # not provided in the original work
-  # and taken as is
-  altaz_parameters <- altaz(
-    solar_parameters$DS,
-    solar_parameters$H,
-    solar_parameters$SD,
-    cos(latitude),
-    sin(latitude),
-    DR,
-    RD
-  )
+    # C wrapper call
+    output <- .Call(
+      'c_skylight_f',
+      forcing = as.matrix(forcing),
+      n = as.integer(nrow(forcing))
+    )
 
-  H <- altaz_parameters$H
-  Z <- altaz_parameters$H * DR
-  solar_azimuth <- altaz_parameters$AZ
+    output <- matrix(output, nrow(forcing), 8, byrow = FALSE)
+    colnames(output) <- c("sun_azimuth","sun_altitude", "sun_illuminance",
+                       "moon_azimuth", "moon_altitude", "moon_illuminance",
+                       "moon_fraction","total_illuminance")
 
-  # solar altitude calculation
-  solar_altitude <- refr(
-    altaz_parameters$H,
-    DR
+    return(data.frame(output))
+  } else {
+
+    # calculate hours as a decimal number
+    hour_dec <- hour + minutes/60
+
+    # constant values
+    RD <- 57.29577951
+    DR <- 1 / RD
+    CE <- 0.91775
+    SE <- 0.39715
+
+    # convert latitude
+    latitude <-  latitude * DR
+
+    J <- 367 * year -
+      as.integer(7 * (year + as.integer((month + 9)/12))/4) +
+      as.integer(275 * month/9) +
+      day - 730531
+
+    E <- hour_dec / 24
+    D <- J - 0.5 + E
+
+    #---- calculate solar parameters ----
+    solar_parameters <- sun(
+      D,
+      DR,
+      RD,
+      CE,
+      SE
     )
 
-  # atmospheric calculations
-  # look up references
-  M <- atmos(
-    solar_altitude,
-    DR
-  )
+    # in place adjustments
+    solar_parameters$T <- solar_parameters$T + 360 * E + longitude
+    solar_parameters$H <- solar_parameters$T - solar_parameters$AS
+
+    # calculate celestial body
+    # parameters all these subroutines
+    # need proper clarifications as
+    # not provided in the original work
+    # and taken as is
+    altaz_parameters <- altaz(
+      solar_parameters$DS,
+      solar_parameters$H,
+      solar_parameters$SD,
+      cos(latitude),
+      sin(latitude),
+      DR,
+      RD
+    )
 
-  # Solar illuminance in lux, scaled using the value
-  # provided by sky_condition. The default does not
-  # scale the value, all other values > 1 scale the
-  # illuminance values
-  solar_illuminance <- 133775 * M / sky_condition
-
-  #---- calculate lunar parameters ----
-  lunar_parameters <- moon(
-    D,
-    solar_parameters$G,
-    CE,
-    SE,
-    RD,
-    DR
-  )
+    H <- altaz_parameters$H
+    Z <- altaz_parameters$H * DR
+    solar_azimuth <- altaz_parameters$AZ
+
+    # solar altitude calculation
+    solar_altitude <- refr(
+      altaz_parameters$H,
+      DR
+      )
+
+    # atmospheric calculations
+    # look up references
+    M <- atmos(
+      solar_altitude,
+      DR
+    )
+
+    # Solar illuminance in lux, scaled using the value
+    # provided by sky_condition. The default does not
+    # scale the value, all other values > 1 scale the
+    # illuminance values
+    solar_illuminance <- 133775 * M / sky_condition
+
+    #---- calculate lunar parameters ----
+    lunar_parameters <- moon(
+      D,
+      solar_parameters$G,
+      CE,
+      SE,
+      RD,
+      DR
+    )
 
-  lunar_parameters$H <- solar_parameters$T - lunar_parameters$AS
+    lunar_parameters$H <- solar_parameters$T - lunar_parameters$AS
 
-  altaz_parameters <- altaz(
-    lunar_parameters$DS,
-    lunar_parameters$H,
-    lunar_parameters$SD,
-    cos(latitude),
-    sin(latitude),
-    DR,
-    RD
-  )
+    altaz_parameters <- altaz(
+      lunar_parameters$DS,
+      lunar_parameters$H,
+      lunar_parameters$SD,
+      cos(latitude),
+      sin(latitude),
+      DR,
+      RD
+    )
 
-  # corrections?
-  Z <- altaz_parameters$H * DR
-  H <- altaz_parameters$H - 0.95 * cos(altaz_parameters$H * DR)
+    # corrections?
+    Z <- altaz_parameters$H * DR
+    H <- altaz_parameters$H - 0.95 * cos(altaz_parameters$H * DR)
 
-  # calculate lunar altitude
-  lunar_altitude <- refr(H, DR)
+    # calculate lunar altitude
+    lunar_altitude <- refr(H, DR)
 
-  # atmospheric conditions?
-  M <- atmos(lunar_altitude, DR)
+    # atmospheric conditions?
+    M <- atmos(lunar_altitude, DR)
 
-  E <- acos(cos(lunar_parameters$V - solar_parameters$LS) * lunar_parameters$CB)
-  P <- 0.892 * exp(-3.343/((tan(E/2.0))^0.632)) + 0.0344 * (sin(E) - E * cos(E))
-  P <- 0.418 * P/(1 - 0.005 * cos(E) - 0.03 * sin(Z))
+    E <- acos(cos(lunar_parameters$V - solar_parameters$LS) * lunar_parameters$CB)
+    P <- 0.892 * exp(-3.343/((tan(E/2.0))^0.632)) + 0.0344 * (sin(E) - E * cos(E))
+    P <- 0.418 * P/(1 - 0.005 * cos(E) - 0.03 * sin(Z))
 
-  # Lunar illuminance in lux, scaled using the value
-  # provided by sky_condition. The default does not
-  # scale the value, all other values > 1 scale the
-  # illuminance values
-  lunar_illuminance <- P * M / sky_condition
+    # Lunar illuminance in lux, scaled using the value
+    # provided by sky_condition. The default does not
+    # scale the value, all other values > 1 scale the
+    # illuminance values
+    lunar_illuminance <- P * M / sky_condition
 
-  # Lunar azimuth/altitude in degrees
-  # again forced to integers seems
-  # check if this requirement can be dropped
-  lunar_azimuth <- altaz_parameters$AZ
+    # Lunar azimuth/altitude in degrees
+    # again forced to integers seems
+    # check if this requirement can be dropped
+    lunar_azimuth <- altaz_parameters$AZ
 
-  # The percentage of the moon illuminated
-  lunar_fraction <- 50 * (1 - cos(E))
+    # The percentage of the moon illuminated
+    lunar_fraction <- 50 * (1 - cos(E))
 
-  # Total sky illuminance, this value is of importance when
-  # considering dusk/dawn conditions mostly, i.e. during hand-off
-  # between solar and lunar illumination conditions
-  total_illuminance <- solar_illuminance + lunar_illuminance + 0.0005 / sky_condition
+    # Total sky illuminance, this value is of importance when
+    # considering dusk/dawn conditions mostly, i.e. during hand-off
+    # between solar and lunar illumination conditions
+    total_illuminance <- solar_illuminance + lunar_illuminance + 0.0005 / sky_condition
+  }
 
   # format output data frame
   output <- data.frame(
@@ -263,3 +294,7 @@ skylight <- function(
     return(output)
   }
 }
+
+.onUnload <- function(libpath) {
+  library.dynam.unload("skylight", libpath)
+}