From 9e9569fbaf1588d7714ff00861892a1676c14fe6 Mon Sep 17 00:00:00 2001 From: KH Date: Sun, 1 Sep 2024 10:40:51 +0200 Subject: [PATCH 01/15] cran sub --- .Rbuildignore | 1 + 1 file changed, 1 insertion(+) 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$ From eb9b8b93bed95e82a7f5261cdb42804b4241b4e0 Mon Sep 17 00:00:00 2001 From: KH Date: Sun, 1 Sep 2024 10:42:18 +0200 Subject: [PATCH 02/15] fortran version --- DESCRIPTION | 2 +- src/skylight.f90 | 92 ++++++++++++++++++++++++ src/subroutines.f90 | 168 ++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 261 insertions(+), 1 deletion(-) create mode 100644 src/skylight.f90 create mode 100644 src/subroutines.f90 diff --git a/DESCRIPTION b/DESCRIPTION index 79b969f..72e001c 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -15,7 +15,7 @@ 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) diff --git a/src/skylight.f90 b/src/skylight.f90 new file mode 100644 index 0000000..566bd11 --- /dev/null +++ b/src/skylight.f90 @@ -0,0 +1,92 @@ +! Defines the main program setup, which will call upon +! different subroutines for flexibility. + +module skylight_mod + + implicit none + + ! derive type for data management + type skylight_data + real :: longitude, latitude, year, month, day, hour, minutes, sky_condition + real :: sun_azimuth, sun_altitude, sun_illuminance + real :: moon_azimuth, moon_altitude, moon_illuminance, moon_fraction + real :: total_illuminance + end type skylight_data + + subroutine skylight(data, n_data, output_data) bind(C, name = "skylight_f_") + implicit none + + type(skylight_data), intent(in) :: data(:) + integer, intent(in) :: n_data + type(skylight_data), intent(out) :: output_data(:) + + integer :: year, month, day, hour, minutes + real :: hour_dec, RD, DR, CE, SE, latitude, J, E, solar_altitude, lunar_altitude, M, P, Z + + ! Constant values + RD = 57.29577951 + DR = 1.0 / RD + CE = 0.91775 + SE = 0.39715 + + ! Parameter conversions + do i = 1, n_data + + hour_dec = hour + minutes / 60.0 + + ! Convert latitude + latitude = data(i)%latitude * DR + + ! Julian day calculation + J = 367 * year - int(7 * (year + int((month + 9) / 12)) / 4) + & + int(275 * month / 9) + day - 730531 + + E = hour_dec / 24.0 + D = J - 0.5 + E + + ! Calculate solar parameters + call sun(D, DR, RD, CE, SE, output_data(i)%sun_azimuth, output_data(i)%sun_altitude, & + output_data(i)%sun_illuminance) + + ! In-place adjustments + output_data(i)%sun_azimuth = output_data(i)%sun_azimuth + 360 * E + data(i)%longitude + output_data(i)%sun_altitude = output_data(i)%sun_azimuth - output_data(i)%sun_illuminance + + ! Calculate celestial body parameters + call altaz(output_data(i)%sun_altitude, output_data(i)%sun_azimuth, output_data(i)%sun_illuminance, & + cos(latitude), sin(latitude), DR, RD, Z, output_data(i)%moon_azimuth) + + ! Solar altitude calculation + call refr(output_data(i)%sun_altitude, DR, output_data(i)%sun_altitude) + + ! Atmospheric calculations + call atmos(output_data(i)%sun_altitude, DR, M) + + ! Solar illuminance + output_data(i)%sun_illuminance = 133775 * M / data(i)%sky_condition + + ! Calculate lunar parameters + call moon(D, output_data(i)%sun_azimuth, CE, SE, RD, DR, & + output_data(i)%moon_azimuth, output_data(i)%moon_altitude, & + output_data(i)%moon_illuminance, output_data(i)%moon_fraction) + + ! Lunar altazimuth + call altaz(output_data(i)%moon_altitude, output_data(i)%moon_azimuth, output_data(i)%moon_illuminance, & + cos(latitude), sin(latitude), DR, RD, Z, output_data(i)%moon_azimuth) + + ! Lunar altitude + call refr(output_data(i)%moon_altitude, DR, output_data(i)%moon_altitude) + + ! Atmospheric conditions + call atmos(output_data(i)%moon_altitude, DR, M) + + ! Lunar illuminance + output_data(i)%moon_illuminance = P * M / data(i)%sky_condition + + ! Total illuminance + output_data(i)%total_illuminance = output_data(i)%sun_illuminance + output_data(i)%moon_illuminance + 0.0005 / data(i)%sky_condition + end do + + end subroutine skylight + +end module skylight_module diff --git a/src/subroutines.f90 b/src/subroutines.f90 new file mode 100644 index 0000000..d2dd8ed --- /dev/null +++ b/src/subroutines.f90 @@ -0,0 +1,168 @@ +module subroutines + +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, intent(in) :: D, DR, RD, CE, SE + real, intent(out) :: T, G, LS, AS, SD, DS + + real :: T_temp + + T = 280.46 + 0.98565 * D + T_temp = T / 360.0 + T = T - int(T_temp) * 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, intent(in) :: D, G, CE, SE, RD, DR + real, intent(out) :: V, SD, AS, DS, CB + + real :: V_temp + + V = 218.32 + 13.1764 * D + V_temp = V / 360.0 + V = V - int(V_temp) * 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, H_out, AZ_out) + implicit none + + real, intent(in) :: DS, H, SD, CI, SI, DR, RD + real, intent(out) :: H_out, AZ_out + + real :: CD, CS, Q, P + + CD = cos(DS) + CS = cos(H * DR) + Q = SD * CI - CD * SI * CS + P = -CD * sin(H * DR) + AZ_out = atan(P/Q) * RD + + where (Q < 0.0) + AZ_out = AZ_out + 180.0 + end where + + where (AZ_out < 0.0) + AZ_out = AZ_out + 360.0 + end where + + AZ_out = int(AZ_out + 0.5) + H_out = 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, intent(in) :: H, DR + real, intent(out) :: HA + + where (H < -5.0 / 6.0) + HA = H + else + 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, intent(in) :: HA, DR + real, intent(out) :: M + + real :: U, X, S + + 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 From 84f8923d44a7b33f6f674225512bed2d8deb5864 Mon Sep 17 00:00:00 2001 From: KH Date: Sun, 1 Sep 2024 10:47:41 +0200 Subject: [PATCH 03/15] comments --- src/subroutines.f90 | 50 ++++++++++++++++++++++----------------------- 1 file changed, 25 insertions(+), 25 deletions(-) diff --git a/src/subroutines.f90 b/src/subroutines.f90 index d2dd8ed..703cd3a 100644 --- a/src/subroutines.f90 +++ b/src/subroutines.f90 @@ -2,11 +2,11 @@ module subroutines 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 +! 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 @@ -39,11 +39,11 @@ subroutine sun(D, DR, RD, CE, SE, T, G, LS, AS, SD, DS) 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 +! 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 @@ -90,11 +90,11 @@ subroutine moon(D, G, CE, SE, RD, DR, V, SD, AS, DS, CB) 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 +! 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, H_out, AZ_out) implicit none @@ -123,11 +123,11 @@ subroutine altaz(DS, H, SD, CI, SI, DR, RD, H_out, AZ_out) 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 +! 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 @@ -143,11 +143,11 @@ subroutine refr(H, DR, HA) 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 +! 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) From ebe7810b2ee061a580c3ddaf70f70a400ffcb832 Mon Sep 17 00:00:00 2001 From: KH Date: Sun, 1 Sep 2024 10:48:03 +0200 Subject: [PATCH 04/15] git ignore --- .gitignore | 1 + 1 file changed, 1 insertion(+) diff --git a/.gitignore b/.gitignore index 8cde772..376bff2 100644 --- a/.gitignore +++ b/.gitignore @@ -5,3 +5,4 @@ data-raw/ data/ docs/ +CRAN-SUBMISSION From ed77667820989d7f6dfae30fc6577efb85cee82e Mon Sep 17 00:00:00 2001 From: KH Date: Sun, 1 Sep 2024 11:10:12 +0200 Subject: [PATCH 05/15] coding C wrapper --- src/C_wrapper.c | 58 +++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 58 insertions(+) create mode 100644 src/C_wrapper.c diff --git a/src/C_wrapper.c b/src/C_wrapper.c new file mode 100644 index 0000000..2e04fc6 --- /dev/null +++ b/src/C_wrapper.c @@ -0,0 +1,58 @@ +#include +#include +#include +#include +#include + +void F77_NAME(skylight_f)( + double *longitude, + double *latitude, + double *altitude, + double *whc, + double *par, + double *forcing, + double *output + ); + +// C wrapper function, order of arguments is fixed +extern SEXP skylight_C( + SEXP longitude, + SEXP latitude, + SEXP altitude, + SEXP n + ){ + + // Number of time steps (same in forcing and output) + 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, 19) ); + + // Fortran subroutine call + F77_CALL(pmodel_f)( + REAL(longitude), + REAL(latitude), + REAL(altitude), + 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[] = { + {"skylight_C", (DL_FUNC) &skylight_C, 23}, + {NULL,NULL,0} +}; + +void R_init_rsofun(DllInfo *dll) +{ + R_registerRoutines(dll, NULL, CallEntries, NULL, NULL); + R_useDynamicSymbols(dll, FALSE); + R_RegisterCCallable("skylight_C", "skylight_C", (DL_FUNC) &skylight_C); +} From 499fc0d0e4133853089d9d7b02281276d1f534ed Mon Sep 17 00:00:00 2001 From: KH Date: Sun, 1 Sep 2024 11:10:33 +0200 Subject: [PATCH 06/15] indentation --- src/subroutines.f90 | 261 ++++++++++++++++++++++---------------------- 1 file changed, 131 insertions(+), 130 deletions(-) diff --git a/src/subroutines.f90 b/src/subroutines.f90 index 703cd3a..cefcb50 100644 --- a/src/subroutines.f90 +++ b/src/subroutines.f90 @@ -1,168 +1,169 @@ module subroutines -contains + 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 + ! 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, intent(in) :: D, DR, RD, CE, SE - real, intent(out) :: T, G, LS, AS, SD, DS + real, intent(in) :: D, DR, RD, CE, SE + real, intent(out) :: T, G, LS, AS, SD, DS - real :: T_temp - - T = 280.46 + 0.98565 * D - T_temp = T / 360.0 - T = T - int(T_temp) * 360.0 - - where (T < 0.0) - T = T + 360.0 - end where + real :: T_temp + + T = 280.46 + 0.98565 * D + T_temp = T / 360.0 + T = T - int(T_temp) * 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) + 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 + where (Y < 0.0) + AS = AS + 180.0 + end where - SD = SE * sin(LS) - DS = asin(SD) - T = T - 180.0 + SD = SE * sin(LS) + DS = asin(SD) + T = T - 180.0 -end subroutine sun + 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 + ! 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, intent(in) :: D, G, CE, SE, RD, DR - real, intent(out) :: V, SD, AS, DS, CB + subroutine moon(D, G, CE, SE, RD, DR, V, SD, AS, DS, CB) + implicit none + + real, intent(in) :: D, G, CE, SE, RD, DR + real, intent(out) :: V, SD, AS, DS, CB - real :: V_temp + real :: V_temp - V = 218.32 + 13.1764 * D - V_temp = V / 360.0 - V = V - int(V_temp) * 360.0 + V = 218.32 + 13.1764 * D + V_temp = V / 360.0 + V = V - int(V_temp) * 360.0 - where (V < 0.0) - V = V + 360.0 - end where + 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 + 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 + DS = asin(SD) + AS = atan(P/Q) * RD - where (Q < 0.0) - AS = AS + 180.0 - end where + where (Q < 0.0) + AS = AS + 180.0 + end where -end subroutine moon + 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 + ! 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, H_out, AZ_out) - implicit none + subroutine altaz(DS, H, SD, CI, SI, DR, RD, H_out, AZ_out) + implicit none - real, intent(in) :: DS, H, SD, CI, SI, DR, RD - real, intent(out) :: H_out, AZ_out + real, intent(in) :: DS, H, SD, CI, SI, DR, RD + real, intent(out) :: H_out, AZ_out - real :: CD, CS, Q, P + real :: CD, CS, Q, P - CD = cos(DS) - CS = cos(H * DR) - Q = SD * CI - CD * SI * CS - P = -CD * sin(H * DR) - AZ_out = atan(P/Q) * RD + CD = cos(DS) + CS = cos(H * DR) + Q = SD * CI - CD * SI * CS + P = -CD * sin(H * DR) + AZ_out = atan(P/Q) * RD - where (Q < 0.0) - AZ_out = AZ_out + 180.0 - end where + where (Q < 0.0) + AZ_out = AZ_out + 180.0 + end where - where (AZ_out < 0.0) - AZ_out = AZ_out + 360.0 - end where + where (AZ_out < 0.0) + AZ_out = AZ_out + 360.0 + end where - AZ_out = int(AZ_out + 0.5) - H_out = asin(SD * SI + CD * CI * CS) * RD + AZ_out = int(AZ_out + 0.5) + H_out = asin(SD * SI + CD * CI * CS) * RD -end subroutine altaz + 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 + ! 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 + subroutine refr(H, DR, HA) + implicit none - real, intent(in) :: H, DR - real, intent(out) :: HA + real, intent(in) :: H, DR + real, intent(out) :: HA - where (H < -5.0 / 6.0) - HA = H - else - HA = H + 1.0 / tan((H + 8.6 / (H + 4.42)) * DR) / 60.0 - end where + where (H < -5.0 / 6.0) + HA = H + else + HA = H + 1.0 / tan((H + 8.6 / (H + 4.42)) * DR) / 60.0 + end where -end subroutine refr + 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 + ! 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 + subroutine atmos(HA, DR, M) + implicit none - real, intent(in) :: HA, DR - real, intent(out) :: M + real, intent(in) :: HA, DR + real, intent(out) :: M + + real :: U, X, S + + 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) - real :: U, X, S - - 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 subroutine atmos +end module subroutine From f119061f305affa8396842f416c7991ea4480604 Mon Sep 17 00:00:00 2001 From: KH Date: Sun, 1 Sep 2024 11:10:43 +0200 Subject: [PATCH 07/15] main routine conversion --- src/skylight.f90 | 77 ++++++++++++++++++++++-------------------------- 1 file changed, 36 insertions(+), 41 deletions(-) diff --git a/src/skylight.f90 b/src/skylight.f90 index 566bd11..82277a9 100644 --- a/src/skylight.f90 +++ b/src/skylight.f90 @@ -17,7 +17,6 @@ subroutine skylight(data, n_data, output_data) bind(C, name = "skylight_f_") implicit none type(skylight_data), intent(in) :: data(:) - integer, intent(in) :: n_data type(skylight_data), intent(out) :: output_data(:) integer :: year, month, day, hour, minutes @@ -29,63 +28,59 @@ subroutine skylight(data, n_data, output_data) bind(C, name = "skylight_f_") CE = 0.91775 SE = 0.39715 - ! Parameter conversions - do i = 1, n_data + hour_dec = hour + minutes / 60.0 - hour_dec = hour + minutes / 60.0 + ! Convert latitude + latitude = data(i)%latitude * DR - ! Convert latitude - latitude = data(i)%latitude * DR + ! Julian day calculation + J = 367 * year - int(7 * (year + int((month + 9) / 12)) / 4) + & + int(275 * month / 9) + day - 730531 - ! Julian day calculation - J = 367 * year - int(7 * (year + int((month + 9) / 12)) / 4) + & - int(275 * month / 9) + day - 730531 + E = hour_dec / 24.0 + D = J - 0.5 + E - E = hour_dec / 24.0 - D = J - 0.5 + E - - ! Calculate solar parameters - call sun(D, DR, RD, CE, SE, output_data(i)%sun_azimuth, output_data(i)%sun_altitude, & + ! Calculate solar parameters + call sun(D, DR, RD, CE, SE, output_data(i)%sun_azimuth, output_data(i)%sun_altitude, & output_data(i)%sun_illuminance) - ! In-place adjustments - output_data(i)%sun_azimuth = output_data(i)%sun_azimuth + 360 * E + data(i)%longitude - output_data(i)%sun_altitude = output_data(i)%sun_azimuth - output_data(i)%sun_illuminance + ! In-place adjustments + output_data(i)%sun_azimuth = output_data(i)%sun_azimuth + 360 * E + data(i)%longitude + output_data(i)%sun_altitude = output_data(i)%sun_azimuth - output_data(i)%sun_illuminance - ! Calculate celestial body parameters - call altaz(output_data(i)%sun_altitude, output_data(i)%sun_azimuth, output_data(i)%sun_illuminance, & - cos(latitude), sin(latitude), DR, RD, Z, output_data(i)%moon_azimuth) + ! Calculate celestial body parameters + call altaz(output_data(i)%sun_altitude, output_data(i)%sun_azimuth, output_data(i)%sun_illuminance, & + cos(latitude), sin(latitude), DR, RD, Z, output_data(i)%moon_azimuth) - ! Solar altitude calculation - call refr(output_data(i)%sun_altitude, DR, output_data(i)%sun_altitude) + ! Solar altitude calculation + call refr(output_data(i)%sun_altitude, DR, output_data(i)%sun_altitude) - ! Atmospheric calculations - call atmos(output_data(i)%sun_altitude, DR, M) + ! Atmospheric calculations + call atmos(output_data(i)%sun_altitude, DR, M) - ! Solar illuminance - output_data(i)%sun_illuminance = 133775 * M / data(i)%sky_condition + ! Solar illuminance + output_data(i)%sun_illuminance = 133775 * M / data(i)%sky_condition - ! Calculate lunar parameters - call moon(D, output_data(i)%sun_azimuth, CE, SE, RD, DR, & - output_data(i)%moon_azimuth, output_data(i)%moon_altitude, & - output_data(i)%moon_illuminance, output_data(i)%moon_fraction) + ! Calculate lunar parameters + call moon(D, output_data(i)%sun_azimuth, CE, SE, RD, DR, & + output_data(i)%moon_azimuth, output_data(i)%moon_altitude, & + output_data(i)%moon_illuminance, output_data(i)%moon_fraction) - ! Lunar altazimuth - call altaz(output_data(i)%moon_altitude, output_data(i)%moon_azimuth, output_data(i)%moon_illuminance, & + ! Lunar altazimuth + call altaz(output_data(i)%moon_altitude, output_data(i)%moon_azimuth, output_data(i)%moon_illuminance, & cos(latitude), sin(latitude), DR, RD, Z, output_data(i)%moon_azimuth) - ! Lunar altitude - call refr(output_data(i)%moon_altitude, DR, output_data(i)%moon_altitude) + ! Lunar altitude + call refr(output_data(i)%moon_altitude, DR, output_data(i)%moon_altitude) - ! Atmospheric conditions - call atmos(output_data(i)%moon_altitude, DR, M) + ! Atmospheric conditions + call atmos(output_data(i)%moon_altitude, DR, M) - ! Lunar illuminance - output_data(i)%moon_illuminance = P * M / data(i)%sky_condition + ! Lunar illuminance + output_data(i)%moon_illuminance = P * M / data(i)%sky_condition - ! Total illuminance - output_data(i)%total_illuminance = output_data(i)%sun_illuminance + output_data(i)%moon_illuminance + 0.0005 / data(i)%sky_condition - end do + ! Total illuminance + output_data(i)%total_illuminance = output_data(i)%sun_illuminance + output_data(i)%moon_illuminance + 0.0005 / data(i)%sky_condition end subroutine skylight From c840a5ade0bbde3aa0e1969a2da4d23848d1494c Mon Sep 17 00:00:00 2001 From: KH Date: Sun, 1 Sep 2024 20:07:08 +0200 Subject: [PATCH 08/15] restructuring logic --- R/skylight.R | 237 +++++++++++++++++-------------- src/Makevars.in | 25 ++++ src/skylight.f90 | 159 +++++++++++++++------ src/subroutines.f90 | 50 +++---- src/{C_wrapper.c => wrappersc.c} | 27 +--- 5 files changed, 302 insertions(+), 196 deletions(-) create mode 100644 src/Makevars.in rename src/{C_wrapper.c => wrappersc.c} (64%) diff --git a/R/skylight.R b/R/skylight.R index 75bfc33..047ed0d 100644 --- a/R/skylight.R +++ b/R/skylight.R @@ -115,131 +115,148 @@ 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 - ) - - # 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 - ) + if (fast) { + # C wrapper call + output <- .Call( + 'skylight_f_C', + n = n, + forcing = data.frame( + .data, + year = year, + month = month, + day = day, + hour = hour, + minutes = minutes + ) + ) + } else { - H <- altaz_parameters$H - Z <- altaz_parameters$H * DR - solar_azimuth <- altaz_parameters$AZ + # 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 + ) - # solar altitude calculation - solar_altitude <- refr( - altaz_parameters$H, - 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 ) - # atmospheric calculations - # look up references - M <- atmos( - solar_altitude, - 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 - ) + # 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( diff --git a/src/Makevars.in b/src/Makevars.in new file mode 100644 index 0000000..e1def93 --- /dev/null +++ b/src/Makevars.in @@ -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 +# : +subroutines.mod: subroutines.o +skylight_r_mod.mod: skylight.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 index 82277a9..acc0281 100644 --- a/src/skylight.f90 +++ b/src/skylight.f90 @@ -1,23 +1,57 @@ -! Defines the main program setup, which will call upon -! different subroutines for flexibility. - -module skylight_mod +module skylight_r_mod + use, intrinsic :: iso_c_binding, only: c_double, c_int, c_char, c_bool implicit none - ! derive type for data management - type skylight_data - real :: longitude, latitude, year, month, day, hour, minutes, sky_condition - real :: sun_azimuth, sun_altitude, sun_illuminance - real :: moon_azimuth, moon_altitude, moon_illuminance, moon_fraction - real :: total_illuminance - end type skylight_data - - subroutine skylight(data, n_data, output_data) bind(C, name = "skylight_f_") - implicit none - - type(skylight_data), intent(in) :: data(:) - type(skylight_data), intent(out) :: output_data(:) + private + public :: skylight_f + + ! set input output data types + type skylight_forcing + real(kind=sp) :: longitude + real(kind=sp) :: latitude + real(kind=sp) :: date + real(kind=sp) :: year + real(kind=sp) :: month + real(kind=sp) :: day + real(kind=sp) :: hour + real(kind=sp) :: minutes + real(kind=sp) :: sky_conditions + end type skylight_forcing + + type skylight_output + real(kind=sp) :: solar_azimuth + real(kind=sp) :: solar_altitude + real(kind=sp) :: solar_illuminance + real(kind=sp) :: lunar_azimuth + real(kind=sp) :: lunar_altitude + real(kind=sp) :: lunar_illuminance + real(kind=sp) :: lunar_fraction + real(kind=sp) :: total_illuminance + end type skylight_output + + subroutine skylight_f( & + forcing, & + output & + ) bind(C, name = "skylight_f_") + + use subroutines + + integer :: n + real(kind=dp), dimension(n, 9), intent(in) :: forcing + type(skylight_forcing), dimension(n, 9) :: input + type(skylight_output), dimension(n, 8), intent(out) :: output + + ! warning: column indices in forcing array are hard coded + input(:)%longitude = real(forcing(:, 1)) + input(:)%latitude = real(forcing(:, 2)) + input(:)%date = real(forcing(:, 3)) + input(:)%year = real(forcing(:, 4)) + input(:)%month = real(forcing(:, 5)) + input(:)%day = real(forcing(:, 6)) + input(:)%hour = real(forcing(:, 7)) + input(:)%minutes = real(forcing(:, 8)) + input(:)%sky_conditions = real(forcing(:, 9)) integer :: year, month, day, hour, minutes real :: hour_dec, RD, DR, CE, SE, latitude, J, E, solar_altitude, lunar_altitude, M, P, Z @@ -31,57 +65,98 @@ subroutine skylight(data, n_data, output_data) bind(C, name = "skylight_f_") hour_dec = hour + minutes / 60.0 ! Convert latitude - latitude = data(i)%latitude * DR + latitude = input%latitude * DR ! Julian day calculation - J = 367 * year - int(7 * (year + int((month + 9) / 12)) / 4) + & - int(275 * month / 9) + day - 730531 + J = 367 * int(input%year) - int(7 * (input%year + int((input%month + 9) / 12)) / 4) + & + int(275 * input%month / 9) + day - 730531 E = hour_dec / 24.0 D = J - 0.5 + E ! Calculate solar parameters - call sun(D, DR, RD, CE, SE, output_data(i)%sun_azimuth, output_data(i)%sun_altitude, & - output_data(i)%sun_illuminance) - - ! In-place adjustments - output_data(i)%sun_azimuth = output_data(i)%sun_azimuth + 360 * E + data(i)%longitude - output_data(i)%sun_altitude = output_data(i)%sun_azimuth - output_data(i)%sun_illuminance + call sun( & + D, & + DR, & + RD, & + CE, & + SE, & + output + ) + + ! 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_data(i)%sun_altitude, output_data(i)%sun_azimuth, output_data(i)%sun_illuminance, & - cos(latitude), sin(latitude), DR, RD, Z, output_data(i)%moon_azimuth) + 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_data(i)%sun_altitude, DR, output_data(i)%sun_altitude) + call refr(output%sun_altitude, DR, output%sun_altitude) ! Atmospheric calculations - call atmos(output_data(i)%sun_altitude, DR, M) + call atmos(output%sun_altitude, DR, M) ! Solar illuminance - output_data(i)%sun_illuminance = 133775 * M / data(i)%sky_condition + output%sun_illuminance = 133775 * M / input%sky_condition ! Calculate lunar parameters - call moon(D, output_data(i)%sun_azimuth, CE, SE, RD, DR, & - output_data(i)%moon_azimuth, output_data(i)%moon_altitude, & - output_data(i)%moon_illuminance, output_data(i)%moon_fraction) + 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_data(i)%moon_altitude, output_data(i)%moon_azimuth, output_data(i)%moon_illuminance, & - cos(latitude), sin(latitude), DR, RD, Z, output_data(i)%moon_azimuth) + 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_data(i)%moon_altitude, DR, output_data(i)%moon_altitude) + call refr( & + output%moon_altitude, & + DR, & + output%moon_altitude & + ) ! Atmospheric conditions - call atmos(output_data(i)%moon_altitude, DR, M) + call atmos( & + output%moon_altitude, & + DR, & + M & + ) ! Lunar illuminance - output_data(i)%moon_illuminance = P * M / data(i)%sky_condition + output%moon_illuminance = P * M / input%sky_condition ! Total illuminance - output_data(i)%total_illuminance = output_data(i)%sun_illuminance + output_data(i)%moon_illuminance + 0.0005 / data(i)%sky_condition + output%total_illuminance = output%sun_illuminance + output%moon_illuminance + 0.0005 / input%sky_condition - end subroutine skylight + end subroutine skylight_f -end module skylight_module +end module skylight_r_mod diff --git a/src/subroutines.f90 b/src/subroutines.f90 index cefcb50..c4b2af9 100644 --- a/src/subroutines.f90 +++ b/src/subroutines.f90 @@ -15,23 +15,24 @@ subroutine sun(D, DR, RD, CE, SE, T, G, LS, AS, SD, DS) real, intent(out) :: T, G, LS, AS, SD, DS real :: T_temp + real :: Y T = 280.46 + 0.98565 * D T_temp = T / 360.0 T = T - int(T_temp) * 360.0 - where (T < 0.0) - T = T + 360.0 - end where + !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 + !where (Y < 0.0) + ! AS = AS + 180.0 + !end where SD = SE * sin(LS) DS = asin(SD) @@ -52,14 +53,15 @@ subroutine moon(D, G, CE, SE, RD, DR, V, SD, AS, DS, CB) real, intent(out) :: V, SD, AS, DS, CB real :: V_temp + real :: Y, X, W, CD, O, P, Q, S, SB, SV V = 218.32 + 13.1764 * D V_temp = V / 360.0 V = V - int(V_temp) * 360.0 - where (V < 0.0) - V = V + 360.0 - end where + !where (V < 0.0) + ! V = V + 360.0 + !end where Y = (134.96 + 13.06499 * D) * DR O = (93.27 + 13.22935 * D) * DR @@ -83,9 +85,9 @@ subroutine moon(D, G, CE, SE, RD, DR, V, SD, AS, DS, CB) DS = asin(SD) AS = atan(P/Q) * RD - where (Q < 0.0) - AS = AS + 180.0 - end where + !where (Q < 0.0) + ! AS = AS + 180.0 + !end where end subroutine moon @@ -110,13 +112,13 @@ subroutine altaz(DS, H, SD, CI, SI, DR, RD, H_out, AZ_out) P = -CD * sin(H * DR) AZ_out = atan(P/Q) * RD - where (Q < 0.0) - AZ_out = AZ_out + 180.0 - end where + !where (Q < 0.0) + ! AZ_out = AZ_out + 180.0 + !end where - where (AZ_out < 0.0) - AZ_out = AZ_out + 360.0 - end where + !where (AZ_out < 0.0) + ! AZ_out = AZ_out + 360.0 + !end where AZ_out = int(AZ_out + 0.5) H_out = asin(SD * SI + CD * CI * CS) * RD @@ -135,11 +137,11 @@ subroutine refr(H, DR, HA) real, intent(in) :: H, DR real, intent(out) :: HA - where (H < -5.0 / 6.0) - HA = H - else - HA = H + 1.0 / tan((H + 8.6 / (H + 4.42)) * DR) / 60.0 - end where + !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 @@ -166,4 +168,4 @@ subroutine atmos(HA, DR, M) (1.0 + (HA + 90.0) * U / 57.29577951) end subroutine atmos -end module subroutine +end module subroutines diff --git a/src/C_wrapper.c b/src/wrappersc.c similarity index 64% rename from src/C_wrapper.c rename to src/wrappersc.c index 2e04fc6..921cdf1 100644 --- a/src/C_wrapper.c +++ b/src/wrappersc.c @@ -5,42 +5,29 @@ #include void F77_NAME(skylight_f)( - double *longitude, - double *latitude, - double *altitude, - double *whc, - double *par, - double *forcing, + double *input, double *output ); // C wrapper function, order of arguments is fixed extern SEXP skylight_C( - SEXP longitude, - SEXP latitude, - SEXP altitude, - SEXP n + SEXP input, + SEXP output + //SEXP n ){ - // Number of time steps (same in forcing and output) - 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, 19) ); + //SEXP output = PROTECT( allocMatrix(REALSXP, nt, 19) ); // Fortran subroutine call - F77_CALL(pmodel_f)( - REAL(longitude), - REAL(latitude), - REAL(altitude), - INTEGER(n), + F77_CALL(skylight_f)( + REAL(input), REAL(output) ); UNPROTECT(1); - return output; }; From 6801df73a9000f1e1aba33ae2f2d9fc5d8a29b6a Mon Sep 17 00:00:00 2001 From: KH Date: Wed, 4 Sep 2024 19:01:11 +0200 Subject: [PATCH 09/15] updating code --- DESCRIPTION | 2 +- NAMESPACE | 1 + R/skylight.R | 25 +++-- man/skylight.Rd | 29 +++++- src/{Makevars.in => Makevars} | 2 +- src/skylight.f90 | 173 ++++++++++------------------------ src/subroutines.f90 | 126 ++++++++++++++----------- src/wrappersc.c | 20 ++-- 8 files changed, 178 insertions(+), 200 deletions(-) rename src/{Makevars.in => Makevars} (100%) diff --git a/DESCRIPTION b/DESCRIPTION index 72e001c..4847365 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -22,7 +22,7 @@ Depends: 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 047ed0d..15b3003 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 #' #' @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 @@ -116,19 +118,26 @@ skylight <- function( minutes <- as.numeric(format(date, "%M")) if (fast) { - # C wrapper call - output <- .Call( - 'skylight_f_C', - n = n, + forcing = data.frame( - .data, + longitude = longitude, + latitude = latitude, year = year, month = month, day = day, hour = hour, - minutes = minutes - ) + minutes = minutes, + sky_condition = sky_condition + ) + + # C wrapper call + output <- .Call( + 'c_skylight_f', + forcing = as.matrix(forcing), + n = as.integer(nrow(forcing)) ) + + return(output) } else { # calculate hours as a decimal number diff --git a/man/skylight.Rd b/man/skylight.Rd index b631408..ce00b4b 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 @@ -22,7 +22,8 @@ illuminance values (1 = cloud cover < 30%, 2 = thin veiled clouds 3 = average clouds, 10 = dark stratus clouds). By and large this 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)} +(this can be a single value or vector of values) +@param fast fast processing} } \value{ Sun and moon illuminance values (in lux), as well as their respective @@ -49,13 +50,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.in b/src/Makevars similarity index 100% rename from src/Makevars.in rename to src/Makevars index e1def93..35833a1 100644 --- a/src/Makevars.in +++ b/src/Makevars @@ -15,8 +15,8 @@ skylight.o: subroutines.o subroutines.mod # Source (object) of Fortran modules # : -subroutines.mod: subroutines.o skylight_r_mod.mod: skylight.o +subroutines.mod: subroutines.o # Dependency of the C wrapper wrappersc.o: skylight_r_mod.mod diff --git a/src/skylight.f90 b/src/skylight.f90 index acc0281..2a81d3d 100644 --- a/src/skylight.f90 +++ b/src/skylight.f90 @@ -1,60 +1,35 @@ module skylight_r_mod + use, intrinsic :: iso_c_binding - use, intrinsic :: iso_c_binding, only: c_double, c_int, c_char, c_bool implicit none private public :: skylight_f - ! set input output data types - type skylight_forcing - real(kind=sp) :: longitude - real(kind=sp) :: latitude - real(kind=sp) :: date - real(kind=sp) :: year - real(kind=sp) :: month - real(kind=sp) :: day - real(kind=sp) :: hour - real(kind=sp) :: minutes - real(kind=sp) :: sky_conditions - end type skylight_forcing - - type skylight_output - real(kind=sp) :: solar_azimuth - real(kind=sp) :: solar_altitude - real(kind=sp) :: solar_illuminance - real(kind=sp) :: lunar_azimuth - real(kind=sp) :: lunar_altitude - real(kind=sp) :: lunar_illuminance - real(kind=sp) :: lunar_fraction - real(kind=sp) :: total_illuminance - end type skylight_output +contains subroutine skylight_f( & forcing, & + n, & output & ) bind(C, name = "skylight_f_") use subroutines + implicit none - integer :: n - real(kind=dp), dimension(n, 9), intent(in) :: forcing - type(skylight_forcing), dimension(n, 9) :: input - type(skylight_output), dimension(n, 8), intent(out) :: output - - ! warning: column indices in forcing array are hard coded - input(:)%longitude = real(forcing(:, 1)) - input(:)%latitude = real(forcing(:, 2)) - input(:)%date = real(forcing(:, 3)) - input(:)%year = real(forcing(:, 4)) - input(:)%month = real(forcing(:, 5)) - input(:)%day = real(forcing(:, 6)) - input(:)%hour = real(forcing(:, 7)) - input(:)%minutes = real(forcing(:, 8)) - input(:)%sky_conditions = real(forcing(:, 9)) - - integer :: year, month, day, hour, minutes - real :: hour_dec, RD, DR, CE, SE, latitude, J, E, solar_altitude, lunar_altitude, M, P, Z + 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, dimension(n) :: T, G, LS, AS, SD, DS, V, & + CB, H, CI, SI, HA, D, E, J, hour_dec + + real, dimension(n) :: longitude, latitude, year, month, day, & + hour, minutes, sky_conditions + + ! assign constants + real :: RD, DR, CE, SE ! Constant values RD = 57.29577951 @@ -62,100 +37,50 @@ subroutine skylight_f( & CE = 0.91775 SE = 0.39715 - hour_dec = hour + minutes / 60.0 + longitude = real(forcing(:, 1)) + latitude = real(forcing(:, 2)) + year = real(forcing(:, 3)) + month = real(forcing(:, 4)) + day = real(forcing(:, 5)) + hour = real(forcing(:, 6)) + minutes = real(forcing(:, 7)) + sky_conditions = real(forcing(:, 8)) + + ! calculate decimal hours + hour_dec = (hour + minutes) / 60.0 ! Convert latitude - latitude = input%latitude * DR + !latitude = latitude * DR ! Julian day calculation - J = 367 * int(input%year) - int(7 * (input%year + int((input%month + 9) / 12)) / 4) + & - int(275 * input%month / 9) + day - 730531 + J = 367 * int(year) - int(7 * (year + int((month + 9) / 12)) / 4) + & + int(275 * month / 9) + int(day) - 730531 E = hour_dec / 24.0 D = J - 0.5 + E - ! Calculate solar parameters - call sun( & - D, & - DR, & - RD, & - CE, & - SE, & - output - ) + !------------- SUN routines ------------------------ + + ! Calculate solar parameters returning second line values + call sun(D, DR, RD, CE, SE, T, G, LS, AS, SD, DS) ! 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 + T = T + 360 * E + longitude + H = T - AS + + + !------------- MOON routines ------------------------ + + + + + !------------- OUTPUT routines ------------------------ ! Total illuminance - output%total_illuminance = output%sun_illuminance + output%moon_illuminance + 0.0005 / input%sky_condition + !total_illuminance = sun_illuminance + moon_illuminance + 0.0005 / sky_condition + + ! assign T + output(:,1) = longitude end subroutine skylight_f diff --git a/src/subroutines.f90 b/src/subroutines.f90 index c4b2af9..3035cec 100644 --- a/src/subroutines.f90 +++ b/src/subroutines.f90 @@ -1,5 +1,8 @@ module subroutines + ! 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 @@ -8,31 +11,33 @@ module subroutines ! 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, intent(in) :: D, DR, RD, CE, SE - real, intent(out) :: T, G, LS, AS, SD, DS + subroutine sun(& + D, DR, RD, CE, SE, & + T, G, LS, AS, SD, DS & + ) - real :: T_temp - real :: Y + !implicit none + real, intent(in) :: D(:) + real, intent(in) :: DR, RD, CE, SE + real, intent(out), dimension(size(D,1)) :: T, G, LS, AS, SD, DS + real, dimension(size(D,1)) :: Y T = 280.46 + 0.98565 * D - T_temp = T / 360.0 - T = T - int(T_temp) * 360.0 + T = T / 360.0 + T = T - int(T) * 360.0 - !where (T < 0.0) - ! T = T + 360.0 - !end where + 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 + where (Y < 0.0) + AS = AS + 180.0 + end where SD = SE * sin(LS) DS = asin(SD) @@ -46,22 +51,24 @@ end subroutine sun ! 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, intent(in) :: D, G, CE, SE, RD, DR - real, intent(out) :: V, SD, AS, DS, CB + subroutine moon( & + D, G, CE, SE, RD, DR, & + V, SD, AS, DS, CB & + ) - real :: V_temp - real :: Y, X, W, CD, O, P, Q, S, SB, SV + implicit none + real, intent(in) :: D(:), G(:) + real, intent(in) :: CE, SE, RD, DR + real, intent(out), dimension(size(D,1)) :: V, SD, AS, DS, CB + real, dimension(size(D,1)) :: Y, X, W, CD, O, P, Q, S, SB, SV V = 218.32 + 13.1764 * D - V_temp = V / 360.0 - V = V - int(V_temp) * 360.0 + V = V / 360.0 + V = V - int(V) * 360.0 - !where (V < 0.0) - ! V = V + 360.0 - !end where + where (V < 0.0) + V = V + 360.0 + end where Y = (134.96 + 13.06499 * D) * DR O = (93.27 + 13.22935 * D) * DR @@ -85,43 +92,48 @@ subroutine moon(D, G, CE, SE, RD, DR, V, SD, AS, DS, CB) DS = asin(SD) AS = atan(P/Q) * RD - !where (Q < 0.0) - ! AS = AS + 180.0 - !end where + 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, H_out, AZ_out) + subroutine altaz( & + DS, H, SD, CI, SI, DR, RD, AZ & + ) + implicit none - real, intent(in) :: DS, H, SD, CI, SI, DR, RD - real, intent(out) :: H_out, AZ_out + real, intent(inout) :: H(:) + real, intent(in) :: DS, SD, CI, SI + real, intent(in) :: DR, RD + real, intent(out), dimension(size(H, 1)) :: AZ + real, dimension(size(H, 1)) :: CS, Q, P - real :: CD, CS, Q, P + real :: CD CD = cos(DS) CS = cos(H * DR) Q = SD * CI - CD * SI * CS P = -CD * sin(H * DR) - AZ_out = atan(P/Q) * RD + AZ = atan(P/Q) * RD - !where (Q < 0.0) - ! AZ_out = AZ_out + 180.0 - !end where + where (Q < 0.0) + AZ = AZ + 180.0 + end where - !where (AZ_out < 0.0) - ! AZ_out = AZ_out + 360.0 - !end where + where (AZ < 0.0) + AZ = AZ + 360.0 + end where - AZ_out = int(AZ_out + 0.5) - H_out = asin(SD * SI + CD * CI * CS) * RD + AZ = int(AZ + 0.5) + H = asin(SD * SI + CD * CI * CS) * RD end subroutine altaz @@ -134,14 +146,15 @@ end subroutine altaz subroutine refr(H, DR, HA) implicit none - real, intent(in) :: H, DR - real, intent(out) :: HA + real, intent(in) :: H(:) + real, intent(in) :: DR + real, 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 + 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 @@ -155,10 +168,15 @@ end subroutine refr subroutine atmos(HA, DR, M) implicit none - real, intent(in) :: HA, DR - real, intent(out) :: M + real, intent(in) :: HA(:) + real, intent(in) :: DR + real, intent(out), dimension(size(HA, 1)) :: M + + ! internal variables + real, dimension(size(HA, 1)) :: U, S - real :: U, X, S + ! constant + real :: X U = sin(HA * DR) X = 753.66156 diff --git a/src/wrappersc.c b/src/wrappersc.c index 921cdf1..a812f5d 100644 --- a/src/wrappersc.c +++ b/src/wrappersc.c @@ -6,40 +6,44 @@ void F77_NAME(skylight_f)( double *input, + int *n, double *output ); // C wrapper function, order of arguments is fixed -extern SEXP skylight_C( +extern SEXP c_skylight_f( SEXP input, - SEXP output - //SEXP n + 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, 19) ); + SEXP output = PROTECT( allocMatrix(REALSXP, nt, 8)); // Fortran subroutine call F77_CALL(skylight_f)( REAL(input), + INTEGER(n), REAL(output) ); UNPROTECT(1); - return output; + return(output); }; // Specify number of arguments to C wrapper as the last number here static const R_CallMethodDef CallEntries[] = { - {"skylight_C", (DL_FUNC) &skylight_C, 23}, + {"c_skylight_f", (DL_FUNC) &c_skylight_f, 2}, {NULL,NULL,0} }; -void R_init_rsofun(DllInfo *dll) +void R_init_skylight(DllInfo *dll) { R_registerRoutines(dll, NULL, CallEntries, NULL, NULL); R_useDynamicSymbols(dll, FALSE); - R_RegisterCCallable("skylight_C", "skylight_C", (DL_FUNC) &skylight_C); + R_RegisterCCallable("skylight", "c_skylight_f", (DL_FUNC) &c_skylight_f); } From 2e63edb5a5728b5e41d0f26dfd86224cda1ec399 Mon Sep 17 00:00:00 2001 From: KH Date: Thu, 5 Sep 2024 12:14:55 +0200 Subject: [PATCH 10/15] indentation etc --- src/skylight.f90 | 51 +++++++++++++++++++++++++++++++++++++++------ src/subroutines.f90 | 17 ++++++++------- 2 files changed, 54 insertions(+), 14 deletions(-) diff --git a/src/skylight.f90 b/src/skylight.f90 index 2a81d3d..2b1a1e1 100644 --- a/src/skylight.f90 +++ b/src/skylight.f90 @@ -23,10 +23,12 @@ subroutine skylight_f( & ! internal state variables (i.e. subroutine output) real, dimension(n) :: T, G, LS, AS, SD, DS, V, & - CB, H, CI, SI, HA, D, E, J, hour_dec + CB, H, CI, SI, HA, D, E, J, AZ, Z, M, P real, dimension(n) :: longitude, latitude, year, month, day, & - hour, minutes, sky_conditions + hour, minutes, sky_conditions, solar_azimuth, solar_altitude, & + hour_dec, solar_illuminance, lunar_illuminance, lunar_altitude, & + lunar_fraction, lunar_azimuth, total_illuminance ! assign constants real :: RD, DR, CE, SE @@ -50,7 +52,7 @@ subroutine skylight_f( & hour_dec = (hour + minutes) / 60.0 ! Convert latitude - !latitude = latitude * DR + latitude = latitude * DR ! Julian day calculation J = 367 * int(year) - int(7 * (year + int((month + 9) / 12)) / 4) + & @@ -65,22 +67,59 @@ subroutine skylight_f( & call sun(D, DR, RD, CE, SE, T, G, LS, AS, SD, DS) ! In-place adjustments (UTTER LLM BULLSHIT, double ckeck everything) - T = T + 360 * E + longitude + T = T + (360 * E) + longitude 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) / sky_conditions !------------- 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 / sky_conditions + lunar_azimuth = AZ + lunar_fraction = 50 * (1 - cos(E)) !------------- OUTPUT routines ------------------------ ! Total illuminance - !total_illuminance = sun_illuminance + moon_illuminance + 0.0005 / sky_condition + total_illuminance = solar_illuminance + lunar_illuminance + 0.0005 / sky_conditions ! assign T - output(:,1) = longitude + 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 diff --git a/src/subroutines.f90 b/src/subroutines.f90 index 3035cec..e68a7b9 100644 --- a/src/subroutines.f90 +++ b/src/subroutines.f90 @@ -23,8 +23,7 @@ subroutine sun(& real, dimension(size(D,1)) :: Y T = 280.46 + 0.98565 * D - T = T / 360.0 - T = T - int(T) * 360.0 + T = T - (int(T/360.0) * 360.0) where (T < 0.0) T = T + 360.0 @@ -63,8 +62,7 @@ subroutine moon( & real, dimension(size(D,1)) :: Y, X, W, CD, O, P, Q, S, SB, SV V = 218.32 + 13.1764 * D - V = V / 360.0 - V = V - int(V) * 360.0 + V = V - int(V / 360.0) * 360.0 where (V < 0.0) V = V + 360.0 @@ -81,10 +79,13 @@ subroutine moon( & 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 @@ -111,16 +112,16 @@ subroutine altaz( & implicit none real, intent(inout) :: H(:) - real, intent(in) :: DS, SD, CI, SI + real, intent(in) :: DS(:), SD(:), CI(:), SI(:) real, intent(in) :: DR, RD real, intent(out), dimension(size(H, 1)) :: AZ real, dimension(size(H, 1)) :: CS, Q, P - real :: CD + real, dimension(size(DS,1)) :: CD CD = cos(DS) CS = cos(H * DR) - Q = SD * CI - CD * SI * CS + Q = (SD * CI) - (CD * SI * CS) P = -CD * sin(H * DR) AZ = atan(P/Q) * RD @@ -133,7 +134,7 @@ subroutine altaz( & end where AZ = int(AZ + 0.5) - H = asin(SD * SI + CD * CI * CS) * RD + H = asin((SD * SI) + (CD * CI * CS)) * RD end subroutine altaz From 2f12892adc7b4c75d926cc74c4023072dde03c19 Mon Sep 17 00:00:00 2001 From: KH Date: Thu, 5 Sep 2024 12:15:12 +0200 Subject: [PATCH 11/15] process benchmarking --- R/skylight.R | 11 ++++++-- analysis/scraps | 73 +++++++++++++++++++++++++++++++++++++++++++++++++ analysis/test.R | 58 +++++++++++++++++++++++++++++++++++++++ 3 files changed, 139 insertions(+), 3 deletions(-) create mode 100644 analysis/scraps create mode 100755 analysis/test.R diff --git a/R/skylight.R b/R/skylight.R index 15b3003..5a309a0 100644 --- a/R/skylight.R +++ b/R/skylight.R @@ -119,7 +119,7 @@ skylight <- function( if (fast) { - forcing = data.frame( + forcing = data.frame( longitude = longitude, latitude = latitude, year = year, @@ -128,7 +128,7 @@ skylight <- function( hour = hour, minutes = minutes, sky_condition = sky_condition - ) + ) # C wrapper call output <- .Call( @@ -137,7 +137,12 @@ skylight <- function( n = as.integer(nrow(forcing)) ) - return(output) + 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") + + return(data.frame(output)) } else { # calculate hours as a decimal number 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..f5bdbc2 --- /dev/null +++ b/analysis/test.R @@ -0,0 +1,58 @@ +#!/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 +) + +microbenchmark( + "R" = {skylight( + input + )}, + "Fortran" = {skylight( + input, + fast = TRUE + )} +) + +# calculate sky illuminance values for +# a single date/time and location +df <- skylight( + input +) |> + select( + (starts_with("sun") | starts_with("moon")) + ) + +print(head(df)) + +# calculate sky illuminance values for +# a single date/time and location +df <- skylight( + input, + fast = TRUE +) |> + select( + (starts_with("sun") | starts_with("moon")) + ) + +print(head(df)) + +library(profvis) + +profvis::profvis({ + skylight( + input, + fast = TRUE + ) +}) From 2f22dd0fa3a1be3fe7b17b0f7fda10ffa3d57735 Mon Sep 17 00:00:00 2001 From: KH Date: Thu, 5 Sep 2024 14:06:29 +0200 Subject: [PATCH 12/15] add parameter switch --- R/skylight.R | 6 +++++- man/skylight.Rd | 5 +++-- 2 files changed, 8 insertions(+), 3 deletions(-) diff --git a/R/skylight.R b/R/skylight.R index 5a309a0..2357e77 100644 --- a/R/skylight.R +++ b/R/skylight.R @@ -30,7 +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 +#' @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). @@ -294,3 +294,7 @@ skylight <- function( return(output) } } + +.onUnload <- function(libpath) { + library.dynam.unload("skylight", libpath) +} diff --git a/man/skylight.Rd b/man/skylight.Rd index ce00b4b..19e0d67 100644 --- a/man/skylight.Rd +++ b/man/skylight.Rd @@ -22,8 +22,9 @@ illuminance values (1 = cloud cover < 30%, 2 = thin veiled clouds 3 = average clouds, 10 = dark stratus clouds). By and large this 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} +(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 From 91deb1b1e9c00f3ff1b256bdead35a2ead232089 Mon Sep 17 00:00:00 2001 From: KH Date: Thu, 5 Sep 2024 14:19:36 +0200 Subject: [PATCH 13/15] ignore compiled chunks --- .gitignore | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.gitignore b/.gitignore index 376bff2..2323c33 100644 --- a/.gitignore +++ b/.gitignore @@ -6,3 +6,5 @@ data-raw/ data/ docs/ CRAN-SUBMISSION +src/*.o +src/*.mod From b680568b68b3beeb774ca8fffb5946193b9652cd Mon Sep 17 00:00:00 2001 From: KH Date: Sat, 16 Nov 2024 18:41:34 +0100 Subject: [PATCH 14/15] updated decimal time + consistent use of real variables --- R/skylight.R | 2 +- src/skylight.f90 | 52 ++++++++++++++++++++++++--------------------- src/subroutines.f90 | 46 ++++++++++++++++++++------------------- 3 files changed, 53 insertions(+), 47 deletions(-) diff --git a/R/skylight.R b/R/skylight.R index 2357e77..84952ca 100644 --- a/R/skylight.R +++ b/R/skylight.R @@ -162,7 +162,7 @@ skylight <- function( as.integer(275 * month/9) + day - 730531 - E <- hour_dec/24 + E <- hour_dec / 24 D <- J - 0.5 + E #---- calculate solar parameters ---- diff --git a/src/skylight.f90 b/src/skylight.f90 index 2b1a1e1..5fe553c 100644 --- a/src/skylight.f90 +++ b/src/skylight.f90 @@ -1,6 +1,5 @@ module skylight_r_mod use, intrinsic :: iso_c_binding - implicit none private @@ -22,16 +21,18 @@ subroutine skylight_f( & real(kind = c_double), intent(out), dimension(n, 8) :: output ! internal state variables (i.e. subroutine output) - real, dimension(n) :: T, G, LS, AS, SD, DS, V, & + 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, dimension(n) :: longitude, latitude, year, month, day, & - hour, minutes, sky_conditions, solar_azimuth, solar_altitude, & + 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 :: RD, DR, CE, SE + real(kind = c_double) :: RD, DR, CE, SE ! Constant values RD = 57.29577951 @@ -39,24 +40,26 @@ subroutine skylight_f( & CE = 0.91775 SE = 0.39715 - longitude = real(forcing(:, 1)) - latitude = real(forcing(:, 2)) - year = real(forcing(:, 3)) - month = real(forcing(:, 4)) - day = real(forcing(:, 5)) - hour = real(forcing(:, 6)) - minutes = real(forcing(:, 7)) - sky_conditions = real(forcing(:, 8)) + ! 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 = (hour + minutes) / 60.0 + hour_dec = forcing(:, 6) + (forcing(:, 7) / 60.0) ! Convert latitude - latitude = latitude * DR + latitude = forcing(:, 2) * DR ! Julian day calculation - J = 367 * int(year) - int(7 * (year + int((month + 9) / 12)) / 4) + & - int(275 * month / 9) + int(day) - 730531 + 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 @@ -66,8 +69,7 @@ subroutine skylight_f( & ! Calculate solar parameters returning second line values call sun(D, DR, RD, CE, SE, T, G, LS, AS, SD, DS) - ! In-place adjustments (UTTER LLM BULLSHIT, double ckeck everything) - T = T + (360 * E) + longitude + T = T + (360 * E) + forcing(:,1) H = T - AS call altaz(DS, H, SD, cos(latitude), sin(latitude), DR, RD, AZ) @@ -76,11 +78,12 @@ subroutine skylight_f( & ! corrections call refr(H, DR, HA) solar_altitude = HA + call atmos(HA, DR, M) ! readable output solar_azimuth = AZ - solar_illuminance = (133775 * M) / sky_conditions + solar_illuminance = 133775 * M / forcing(:, 8) !------------- MOON routines ------------------------ call moon(D, G, CE, SE, RD, DR, & @@ -92,24 +95,25 @@ subroutine skylight_f( & call altaz(DS, H, SD, cos(latitude), sin(latitude), DR, RD, AZ) Z = H * DR - H = H - 0.95 * cos(Z) + 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 / sky_conditions + lunar_illuminance = P * M / forcing(:, 8) lunar_azimuth = AZ - lunar_fraction = 50 * (1 - cos(E)) + lunar_fraction = 50.0 * (1 - cos(E)) !------------- OUTPUT routines ------------------------ ! Total illuminance - total_illuminance = solar_illuminance + lunar_illuminance + 0.0005 / sky_conditions + total_illuminance = solar_illuminance + lunar_illuminance + 0.0005 / forcing(:, 8) ! assign T output(:,1) = solar_azimuth diff --git a/src/subroutines.f90 b/src/subroutines.f90 index e68a7b9..35433a3 100644 --- a/src/subroutines.f90 +++ b/src/subroutines.f90 @@ -1,5 +1,7 @@ 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 (:)? @@ -17,10 +19,10 @@ subroutine sun(& ) !implicit none - real, intent(in) :: D(:) - real, intent(in) :: DR, RD, CE, SE - real, intent(out), dimension(size(D,1)) :: T, G, LS, AS, SD, DS - real, dimension(size(D,1)) :: Y + 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) @@ -56,10 +58,10 @@ subroutine moon( & ) implicit none - real, intent(in) :: D(:), G(:) - real, intent(in) :: CE, SE, RD, DR - real, intent(out), dimension(size(D,1)) :: V, SD, AS, DS, CB - real, dimension(size(D,1)) :: Y, X, W, CD, O, P, Q, S, SB, SV + 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 @@ -111,13 +113,13 @@ subroutine altaz( & implicit none - real, intent(inout) :: H(:) - real, intent(in) :: DS(:), SD(:), CI(:), SI(:) - real, intent(in) :: DR, RD - real, intent(out), dimension(size(H, 1)) :: AZ - real, dimension(size(H, 1)) :: CS, Q, P + 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, dimension(size(DS,1)) :: CD + real(kind = c_double), dimension(size(DS,1)) :: CD CD = cos(DS) CS = cos(H * DR) @@ -147,9 +149,9 @@ end subroutine altaz subroutine refr(H, DR, HA) implicit none - real, intent(in) :: H(:) - real, intent(in) :: DR - real, intent(out), dimension(size(H,1)) :: HA + 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 @@ -169,15 +171,15 @@ end subroutine refr subroutine atmos(HA, DR, M) implicit none - real, intent(in) :: HA(:) - real, intent(in) :: DR - real, intent(out), dimension(size(HA, 1)) :: M + 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, dimension(size(HA, 1)) :: U, S + real(kind = c_double), dimension(size(HA, 1)) :: U, S ! constant - real :: X + real(kind = c_double) :: X U = sin(HA * DR) X = 753.66156 From a6356348dcf8df41c9a2f42cee527de9287959ed Mon Sep 17 00:00:00 2001 From: KH Date: Sat, 16 Nov 2024 18:41:57 +0100 Subject: [PATCH 15/15] updated test routines --- analysis/test.R | 56 ++++++++++++++++++++++++++++++------------------- 1 file changed, 34 insertions(+), 22 deletions(-) diff --git a/analysis/test.R b/analysis/test.R index f5bdbc2..f48f04c 100755 --- a/analysis/test.R +++ b/analysis/test.R @@ -1,6 +1,6 @@ #!/usr/bin/env Rscript -# R CMD INSTALL --preclean --no-multiarch --with-keep.source skylight +#R CMD INSTALL --preclean --no-multiarch --with-keep.source skylight # load the library library(skylight) @@ -15,44 +15,56 @@ input <- data.frame( sky_condition = 1 ) -microbenchmark( - "R" = {skylight( - input - )}, - "Fortran" = {skylight( - input, - fast = TRUE - )} -) +input_quick <- input[1:20,] # calculate sky illuminance values for # a single date/time and location -df <- skylight( - input +df1 <- skylight( + input_quick ) |> select( (starts_with("sun") | starts_with("moon")) ) -print(head(df)) +print(head(df1)) # calculate sky illuminance values for # a single date/time and location -df <- skylight( - input, +df2 <- skylight( + input_quick, fast = TRUE ) |> select( (starts_with("sun") | starts_with("moon")) ) -print(head(df)) - -library(profvis) +print(head(df2)) -profvis::profvis({ - skylight( +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 +# ) +# })