diff --git a/.Rbuildignore b/.Rbuildignore index c773db6..a70489f 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -24,3 +24,4 @@ ^\.httr-oauth$ ^logo\..*$ ^codemeta\.json$ +^CRAN-SUBMISSION$ diff --git a/.gitignore b/.gitignore index 8cde772..2323c33 100644 --- a/.gitignore +++ b/.gitignore @@ -5,3 +5,6 @@ data-raw/ data/ docs/ +CRAN-SUBMISSION +src/*.o +src/*.mod diff --git a/DESCRIPTION b/DESCRIPTION index 79b969f..4847365 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -15,14 +15,14 @@ Maintainer: Koen Hufkens 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) . -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, diff --git a/NAMESPACE b/NAMESPACE index 5496f4b..3f967b4 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,3 +1,4 @@ # Generated by roxygen2: do not edit by hand export(skylight) +useDynLib(skylight, .registration=TRUE) diff --git a/R/skylight.R b/R/skylight.R index 75bfc33..84952ca 100644 --- a/R/skylight.R +++ b/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) +} diff --git a/analysis/scraps b/analysis/scraps new file mode 100644 index 0000000..c32455a --- /dev/null +++ b/analysis/scraps @@ -0,0 +1,73 @@ + ! In-place adjustments (UTTER LLM BULLSHIT, double ckeck everything) + output%sun_azimuth = output%sun_azimuth + 360 * E + output%longitude + output%sun_altitude = output%sun_azimuth - output%sun_illuminance + + ! Calculate celestial body parameters + call altaz( & + output%sun_altitude, & + output%sun_azimuth, & + output%sun_illuminance, & + cos(latitude), & + sin(latitude), & + DR, & + RD, & + Z, & + output%moon_azimuth & + ) + + ! Solar altitude calculation + call refr(output%sun_altitude, DR, output%sun_altitude) + + ! Atmospheric calculations + call atmos(output%sun_altitude, DR, M) + + ! Solar illuminance + output%sun_illuminance = 133775 * M / input%sky_condition + + ! Calculate lunar parameters + call moon( & + D, & + output%sun_azimuth, & + CE, & + SE, & + RD, & + DR, & + output%moon_azimuth, & + output%moon_altitude, & + output%moon_illuminance, & + output%moon_fraction & + ) + + ! Lunar altazimuth + call altaz( & + output%moon_altitude, & + output%moon_azimuth, & + output%moon_illuminance, & + cos(latitude), & + sin(latitude), & + DR, & + RD, & + Z, & + output%moon_azimuth & + ) + + ! Lunar altitude + call refr( & + output%moon_altitude, & + DR, & + output%moon_altitude & + ) + + ! Atmospheric conditions + call atmos( & + output%moon_altitude, & + DR, & + M & + ) + + ! Lunar illuminance + output%moon_illuminance = P * M / input%sky_condition + + +---- model + diff --git a/analysis/test.R b/analysis/test.R new file mode 100755 index 0000000..f48f04c --- /dev/null +++ b/analysis/test.R @@ -0,0 +1,70 @@ +#!/usr/bin/env Rscript + +#R CMD INSTALL --preclean --no-multiarch --with-keep.source skylight + +# load the library +library(skylight) +library(tidyr) +library(dplyr) +library(microbenchmark) + +input <- data.frame( + longitude = 0, + latitude = 50, + date = as.POSIXct("2020-06-18 00:00:00", tz = "GMT") + seq(0, 60*24*3600, 900), + sky_condition = 1 +) + +input_quick <- input[1:20,] + +# calculate sky illuminance values for +# a single date/time and location +df1 <- skylight( + input_quick +) |> + select( + (starts_with("sun") | starts_with("moon")) + ) + +print(head(df1)) + +# calculate sky illuminance values for +# a single date/time and location +df2 <- skylight( + input_quick, + fast = TRUE +) |> + select( + (starts_with("sun") | starts_with("moon")) + ) + +print(head(df2)) + +b <- microbenchmark( + "R" = {skylight( + input + )}, + "Fortran" = {skylight( + input, + fast = TRUE + )}, + times = 100 +) + +par(mfrow = c(3,1)) +boxplot(b) + +R <- b$time[b$expr == "R"] +hist(log(R)) + +F <- b$time[b$expr != "R"] +hist(log(F)) + +# library(profvis) +# +# profvis::profvis({ +# skylight( +# input, +# fast = TRUE +# ) +# }) diff --git a/man/skylight.Rd b/man/skylight.Rd index b631408..19e0d67 100644 --- a/man/skylight.Rd +++ b/man/skylight.Rd @@ -4,7 +4,7 @@ \alias{skylight} \title{Sky illuminance values for the sun and moon} \usage{ -skylight(.data, longitude, latitude, date, sky_condition = 1) +skylight(.data, longitude, latitude, date, sky_condition = 1, fast = FALSE) } \arguments{ \item{.data}{A data frame or data frame extension (e.g. a tibble) with @@ -23,6 +23,8 @@ illuminance values (1 = cloud cover < 30%, 2 = thin veiled clouds 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)} + +\item{fast}{fast processing (TRUE or FALSE)} } \value{ Sun and moon illuminance values (in lux), as well as their respective @@ -49,13 +51,33 @@ The original code has been vectorized, as such vectors of location, time and/or sky conditions can be provided. } \examples{ + + # run the function on standard + # input variables (single values or vectors of equal size) df <- skylight( longitude = -135.8, latitude = -23.4, date = as.POSIXct("1986-12-18 21:00:00", tz = "GMT"), sky_condition = 1 -) + ) + + print(df) + + # create data frame of input variables + input <- data.frame( + longitude = 0, + latitude = 50, + date = as.POSIXct("2020-06-18 00:00:00", tz = "GMT") + seq(0, 1*24*3600, 1800), + sky_condition = 1 + ) + + # calculate on data frame + df <- skylight(input) + + print(df) -print(df) + # the above statement can also be used + # in a piped fashion in R >= 4.2 + # input |> skylight() } diff --git a/src/Makevars b/src/Makevars new file mode 100644 index 0000000..35833a1 --- /dev/null +++ b/src/Makevars @@ -0,0 +1,25 @@ +# C objects +C_OBJS = wrappersc.o + +# Fortran objects: refer to file names , +# order reflects dependency structure +FT_OBJS = subroutines.o skylight.o + +all: $(SHLIB) clean + +$(SHLIB): $(FT_OBJS) $(C_OBJS) + +# Dependency of objects (?) +# : +skylight.o: subroutines.o subroutines.mod + +# Source (object) of Fortran modules +# : +skylight_r_mod.mod: skylight.o +subroutines.mod: subroutines.o + +# Dependency of the C wrapper +wrappersc.o: skylight_r_mod.mod + +clean: + @rm -rf *.mod *.o diff --git a/src/skylight.f90 b/src/skylight.f90 new file mode 100644 index 0000000..5fe553c --- /dev/null +++ b/src/skylight.f90 @@ -0,0 +1,130 @@ +module skylight_r_mod + use, intrinsic :: iso_c_binding + implicit none + + private + public :: skylight_f + +contains + + subroutine skylight_f( & + forcing, & + n, & + output & + ) bind(C, name = "skylight_f_") + + use subroutines + implicit none + + integer(kind = c_int), intent(in) :: n + real(kind = c_double), intent(in), dimension(n, 8) :: forcing + real(kind = c_double), intent(out), dimension(n, 8) :: output + + ! internal state variables (i.e. subroutine output) + real(kind = c_double), dimension(n) :: T, G, LS, AS, SD, DS, V, & + CB, H, CI, SI, HA, D, E, J, AZ, Z, M, P + + real(kind = c_double), dimension(n) :: latitude, solar_azimuth, solar_altitude, & + hour_dec, solar_illuminance, lunar_illuminance, lunar_altitude, & + lunar_fraction, lunar_azimuth, total_illuminance + + !longitude, latitude, & + !hour, minutes, sky_conditions, + + ! assign constants + real(kind = c_double) :: RD, DR, CE, SE + + ! Constant values + RD = 57.29577951 + DR = 1.0 / RD + CE = 0.91775 + SE = 0.39715 + + ! directly referencing these details + ! keep for clarity + !longitude = forcing(:, 1) + !latitude = forcing(:, 2) + !year = forcing(:, 3) + !month = forcing(:, 4) + !day = forcing(:, 5) + !hour = forcing(:, 6) + !minutes = forcing(:, 7) + !sky_conditions = forcing(:, 8) + + ! calculate decimal hours + hour_dec = forcing(:, 6) + (forcing(:, 7) / 60.0) + + ! Convert latitude + latitude = forcing(:, 2) * DR + + ! Julian day calculation + J = 367 * int(forcing(:, 3)) - int(7 * (forcing(:, 3) + int((forcing(:, 4) + 9) / 12)) / 4) + & + int(275 * forcing(:, 4) / 9) + forcing(:, 5) - 730531 + + E = hour_dec / 24.0 + D = J - 0.5 + E + + !------------- SUN routines ------------------------ + + ! Calculate solar parameters returning second line values + call sun(D, DR, RD, CE, SE, T, G, LS, AS, SD, DS) + + T = T + (360 * E) + forcing(:,1) + H = T - AS + + call altaz(DS, H, SD, cos(latitude), sin(latitude), DR, RD, AZ) + Z = H * DR + + ! corrections + call refr(H, DR, HA) + solar_altitude = HA + + call atmos(HA, DR, M) + + ! readable output + solar_azimuth = AZ + solar_illuminance = 133775 * M / forcing(:, 8) + + !------------- MOON routines ------------------------ + call moon(D, G, CE, SE, RD, DR, & + V, SD, AS, DS, CB & + ) + + ! reclycling parameters NOT SAFE + H = T - AS + + call altaz(DS, H, SD, cos(latitude), sin(latitude), DR, RD, AZ) + Z = H * DR + H = H - (0.95 * cos(Z)) + + call refr(H, DR, HA) + lunar_altitude = HA + + call atmos(HA, DR, M) + + E = acos(cos(V - LS) * 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 = P * M / forcing(:, 8) + lunar_azimuth = AZ + lunar_fraction = 50.0 * (1 - cos(E)) + + !------------- OUTPUT routines ------------------------ + + ! Total illuminance + total_illuminance = solar_illuminance + lunar_illuminance + 0.0005 / forcing(:, 8) + + ! assign T + output(:,1) = solar_azimuth + output(:,2) = solar_altitude + output(:,3) = solar_illuminance + output(:,4) = lunar_azimuth + output(:,5) = lunar_altitude + output(:,6) = lunar_illuminance + output(:,7) = lunar_fraction + output(:,8) = total_illuminance + + end subroutine skylight_f + +end module skylight_r_mod diff --git a/src/subroutines.f90 b/src/subroutines.f90 new file mode 100644 index 0000000..35433a3 --- /dev/null +++ b/src/subroutines.f90 @@ -0,0 +1,192 @@ +module subroutines + + use, intrinsic :: iso_c_binding + + ! basically no implicit assingments, need to make them + ! using array size (n) per subroutine or can I use (:)? + + contains + + ! see SUN subroutine on p.21 of + ! Computer Programs for Sun and Moon Illuminance + ! With Contingent Tables and Diagrams of + ! Janiczek and DeYoung, US Naval observatory circular + ! nr. 171, 1987 + + subroutine sun(& + D, DR, RD, CE, SE, & + T, G, LS, AS, SD, DS & + ) + + !implicit none + real(kind = c_double), intent(in) :: D(:) + real(kind = c_double), intent(in) :: DR, RD, CE, SE + real(kind = c_double), intent(out), dimension(size(D,1)) :: T, G, LS, AS, SD, DS + real(kind = c_double), dimension(size(D,1)) :: Y + + T = 280.46 + 0.98565 * D + T = T - (int(T/360.0) * 360.0) + + where (T < 0.0) + T = T + 360.0 + end where + + G = (357.5 + 0.98560 * D) * DR + LS = (T + 1.91 * sin(G)) * DR + AS = atan(CE * tan(LS)) * RD + Y = cos(LS) + + where (Y < 0.0) + AS = AS + 180.0 + end where + + SD = SE * sin(LS) + DS = asin(SD) + T = T - 180.0 + + end subroutine sun + + ! see MOON subroutine on p.21 of + ! Computer Programs for Sun and Moon Illuminance + ! With Contingent Tables and Diagrams of + ! Janiczek and DeYoung, US Naval observatory circular + ! nr. 171, 1987 + + subroutine moon( & + D, G, CE, SE, RD, DR, & + V, SD, AS, DS, CB & + ) + + implicit none + real(kind = c_double), intent(in) :: D(:), G(:) + real(kind = c_double), intent(in) :: CE, SE, RD, DR + real(kind = c_double), intent(out), dimension(size(D,1)) :: V, SD, AS, DS, CB + real(kind = c_double), dimension(size(D,1)) :: Y, X, W, CD, O, P, Q, S, SB, SV + + V = 218.32 + 13.1764 * D + V = V - int(V / 360.0) * 360.0 + + where (V < 0.0) + V = V + 360.0 + end where + + Y = (134.96 + 13.06499 * D) * DR + O = (93.27 + 13.22935 * D) * DR + W = (235.7 + 24.38150 * D) * DR + SB = sin(Y) + CB = cos(Y) + X = sin(O) + S = cos(O) + SD = sin(W) + CD = cos(W) + V = (V + (6.29-1.27 * CD + 0.43 * CB) * SB + (0.66 + 1.27 * CB) * SD - & + 0.19 * sin(G) - 0.23 * X * S) * DR + + Y = ((5.13 - 0.17 * CD) * X + (0.56 * SB + 0.17 * SD) * S) * DR + + SV = sin(V) + SB = sin(Y) + CB = cos(Y) + + Q = CB * cos(V) + P = CE * SV * CB - SE * SB + SD = SE * SV * CB + CE * SB + + DS = asin(SD) + AS = atan(P/Q) * RD + + where (Q < 0.0) + AS = AS + 180.0 + end where + + end subroutine moon + + ! see ALTAZ subroutine on p.22 of + ! Computer Programs for Sun and Moon Illuminance + ! With Contingent Tables and Diagrams of + ! Janiczek and DeYoung, US Naval observatory circular + ! nr. 171, 1987 + + subroutine altaz( & + DS, H, SD, CI, SI, DR, RD, AZ & + ) + + implicit none + + real(kind = c_double), intent(inout) :: H(:) + real(kind = c_double), intent(in) :: DS(:), SD(:), CI(:), SI(:) + real(kind = c_double), intent(in) :: DR, RD + real(kind = c_double), intent(out), dimension(size(H, 1)) :: AZ + real(kind = c_double), dimension(size(H, 1)) :: CS, Q, P + + real(kind = c_double), dimension(size(DS,1)) :: CD + + CD = cos(DS) + CS = cos(H * DR) + Q = (SD * CI) - (CD * SI * CS) + P = -CD * sin(H * DR) + AZ = atan(P/Q) * RD + + where (Q < 0.0) + AZ = AZ + 180.0 + end where + + where (AZ < 0.0) + AZ = AZ + 360.0 + end where + + AZ = int(AZ + 0.5) + H = asin((SD * SI) + (CD * CI * CS)) * RD + + end subroutine altaz + + ! see REFR subroutine on p.22 of + ! Computer Programs for Sun and Moon Illuminance + ! With Contingent Tables and Diagrams of + ! Janiczek and DeYoung, US Naval observatory circular + ! nr. 171, 1987 + + subroutine refr(H, DR, HA) + implicit none + + real(kind = c_double), intent(in) :: H(:) + real(kind = c_double), intent(in) :: DR + real(kind = c_double), intent(out), dimension(size(H,1)) :: HA + + where (H < -5.0 / 6.0) + HA = H + elsewhere + HA = H + 1.0 / tan((H + 8.6 / (H + 4.42)) * DR) / 60.0 + end where + + end subroutine refr + + ! see ATMOS subroutine on p.22 of + ! Computer Programs for Sun and Moon Illuminance + ! With Contingent Tables and Diagrams of + ! Janiczek and DeYoung, US Naval observatory circular + ! nr. 171, 1987 + + + subroutine atmos(HA, DR, M) + implicit none + + real(kind = c_double), intent(in) :: HA(:) + real(kind = c_double), intent(in) :: DR + real(kind = c_double), intent(out), dimension(size(HA, 1)) :: M + + ! internal variables + real(kind = c_double), dimension(size(HA, 1)) :: U, S + + ! constant + real(kind = c_double) :: X + + U = sin(HA * DR) + X = 753.66156 + S = asin(X * cos(HA * DR) / (X + 1.0)) + M = X * (cos(S) - U) + cos(S) + M = exp(-0.21 * M) * U + 0.0289 * exp(-0.042 * M) * & + (1.0 + (HA + 90.0) * U / 57.29577951) + + end subroutine atmos +end module subroutines diff --git a/src/wrappersc.c b/src/wrappersc.c new file mode 100644 index 0000000..a812f5d --- /dev/null +++ b/src/wrappersc.c @@ -0,0 +1,49 @@ +#include +#include +#include +#include +#include + +void F77_NAME(skylight_f)( + double *input, + int *n, + double *output + ); + +// C wrapper function, order of arguments is fixed +extern SEXP c_skylight_f( + SEXP input, + SEXP n + ){ + + // nr rows + const int nt = INTEGER(n)[0]; + + // Specify output + // 2nd argument to allocMatrix is number of rows, + // 3rd is number of columns + SEXP output = PROTECT( allocMatrix(REALSXP, nt, 8)); + + // Fortran subroutine call + F77_CALL(skylight_f)( + REAL(input), + INTEGER(n), + REAL(output) + ); + + UNPROTECT(1); + return(output); +}; + +// Specify number of arguments to C wrapper as the last number here +static const R_CallMethodDef CallEntries[] = { + {"c_skylight_f", (DL_FUNC) &c_skylight_f, 2}, + {NULL,NULL,0} +}; + +void R_init_skylight(DllInfo *dll) +{ + R_registerRoutines(dll, NULL, CallEntries, NULL, NULL); + R_useDynamicSymbols(dll, FALSE); + R_RegisterCCallable("skylight", "c_skylight_f", (DL_FUNC) &c_skylight_f); +}