UPDATE: Add OudinFortran in PE_Oudin and add tests

R/PE_Oudin.R

 PE_Oudin <- function(JD, Temp,
                      Lat, LatUnit = c("rad", "deg"),
-                     TimeStepIn = "daily", TimeStepOut = "daily") {
+                     TimeStepIn = "daily", TimeStepOut = "daily",
+                     RunFortran = F) {
   
   
   ## ---------- check arguments
@@ -18,22 +19,26 @@ PE_Oudin <- function(JD, Temp,
     stop("'Temp' must be of class 'numeric'")
   }
   if (length(JD) != length(Temp)) {
-    stop("'Temp' and 'LatUnit' must have the same length")
+    stop("'JD' and 'Temp' must have the same length")
   }
-  LatUnit <- match.arg(LatUnit, choices = c("rad", "deg"))
-  if (!inherits(Lat, "numeric") | length(Lat) != 1) {
+
+  if (!RunFortran & (!inherits(Lat, "numeric") | length(Lat) != 1)) {
     stop("'Lat' must be a 'numeric' of length one")
   }
-  if (LatUnit[1L] == "rad" & ((Lat >= pi/2) | (Lat <= -pi/2))) {
+  if (RunFortran & (!inherits(Lat, "numeric") | (!length(Lat) %in% c(1, length(Temp))))) {
+    stop("'Lat' must be a 'numeric' of length one or of the same length as 'Temp'")
+  }
+  LatUnit <- match.arg(LatUnit, choices = c("rad", "deg"))
+  if (LatUnit[1L] == "rad" & (all(Lat >= pi/2) | all(Lat <= -pi/2))) {
     stop("'Lat' must be comprised between -pi/2 and +pi/2 degrees")
   }
-  if (LatUnit[1L] == "deg" & ((Lat >= 90) | (Lat <= -90))) {
+  if (LatUnit[1L] == "deg" & (all(Lat >= 90) | all(Lat <= -90))) {
     stop("'Lat' must be  comprised between -90 and +90 degrees")
   }
-  if (LatUnit[1L] == "rad") {
+  if (!RunFortran & LatUnit[1L] == "rad") {
     FI <- Lat
   }
-  if (LatUnit[1L] == "deg") {
+  if (!RunFortran & LatUnit[1L] == "deg") {
     FI <- Lat / (180 / pi)
   }
   if (any(JD < 0) | any(JD > 366)) {
@@ -62,6 +67,31 @@ PE_Oudin <- function(JD, Temp,
   
   ## ---------- Oudin's formula
   
+  if (RunFortran)
+  {
+    if (LatUnit[1L] == "rad") {
+      Lat = Lat * 180 / pi
+    }
+    
+    LInputs = as.integer(length(Temp))
+    
+    if (length(Lat) == 1) {
+      Lat = rep(Lat, LInputs)
+    }
+    
+    RESULTS <- .Fortran("frun_etp_oudin", PACKAGE="airGR",
+                        ##inputs
+                        LInputs = LInputs,
+                        InputsLAT = as.double(Lat),
+                        InputsTT = as.double(Temp),
+                        InputsJJ = as.double(JD),
+                        ##outputs                                                               
+                        PE_Oudin_D = rep(as.double(-999.999), LInputs)
+    )
+    PE_Oudin_D = RESULTS$PE_Oudin_D
+    
+  } else {
+    
   PE_Oudin_D <- rep(NA, length(Temp))
   COSFI <- cos(FI)
   AFI <- abs(FI / 42) 
@@ -111,6 +141,7 @@ PE_Oudin <- function(JD, Temp,
     
   }
   
+  }
   
   ## ---------- disaggregate PE from daily to hourly
   

tests/testthat/test-evap.R

 context("Test evaporation")
 
-test_that("PEdaily_Oudin works", {
+comp_evap = function(BasinObs,
+                     Lat, LatUnit,
+                     TimeStepIn = "daily",
+                     TimeStepOut = "daily")
+{
+  PotEvap <- PE_Oudin(JD = as.POSIXlt(BasinObs$DatesR)$yday + 1,
+                      Temp = BasinObs$T,
+                      Lat, LatUnit, 
+                      TimeStepIn, TimeStepOut)
+  PotEvapFor <- PE_Oudin(JD = as.POSIXlt(BasinObs$DatesR)$yday + 1,
+                         Temp = BasinObs$T,
+                         Lat, LatUnit,
+                         TimeStepIn, TimeStepOut,
+                         RunFortran = T)
+  all(range(PotEvap - PotEvapFor) < 0.000001)
+}
+
+test_that("PE_Oudin works", {
   skip_on_cran()
   rm(list = ls())
   
-  data(L0123001)
-  PotEvap <- PEdaily_Oudin(JD = as.POSIXlt(BasinObs$DatesR)$yday + 1,
-                           Temp = BasinObs$T,
-                           Lat = 0.8, LatUnit = "rad")
-  PotEvapFor <- PEdaily_Oudin(JD = as.POSIXlt(BasinObs$DatesR)$yday + 1,
-                              Temp = BasinObs$T,
-                              Lat = 0.8, LatUnit = "rad", run_fortran = T)
-  expect_true(all(range(PotEvap - PotEvapFor) < 0.001))
+  data(L0123001); BasinObs_L0123001 <- BasinObs
+  data(L0123002); BasinObs_L0123002 <- BasinObs
+
+  expect_true(comp_evap(BasinObs_L0123001,
+                        0.8, "rad",
+                        "daily", "daily"))
+  expect_true(comp_evap(BasinObs_L0123001,
+                        0.8, "rad",
+                        "daily", "hourly"))
+  expect_true(comp_evap(BasinObs_L0123002,
+                        0.9, "rad", 
+                        "daily", "daily"))
+  expect_true(comp_evap(BasinObs_L0123002,
+                        0.9, "rad",
+                        "daily", "hourly"))
+  
+  ## check with several catchments using different values for Lat
+  
+  ## one by one
+  PotEvapFor1 <- PE_Oudin(JD = as.POSIXlt(BasinObs_L0123001$DatesR)$yday + 1,
+                          Temp = BasinObs_L0123001$T,
+                          Lat = 0.8, LatUnit = "rad",
+                          RunFortran = T)
+  PotEvapFor2 <- PE_Oudin(JD = as.POSIXlt(BasinObs_L0123002$DatesR)$yday + 1,
+                          Temp = BasinObs_L0123002$T,
+                          Lat = 0.9, LatUnit = "rad",
+                          RunFortran = T)
+  
+  ## all in one
+  BasinObs_L0123001$Lat <- 0.8
+  BasinObs_L0123002$Lat <- 0.9
+  BasinObs <- rbind(BasinObs_L0123001, BasinObs_L0123002)
+  PotEvapFor <- PE_Oudin(JD = as.POSIXlt(BasinObs$DatesR)$yday + 1,
+                         Temp = BasinObs$T,
+                         Lat = BasinObs$Lat, LatUnit = "rad",
+                         RunFortran = T)
+  
+  expect_equal(PotEvapFor, c(PotEvapFor1, PotEvapFor2))
+  
 })