diff --git a/PortalForecasts.R b/PortalForecasts.R
index d307b37872..9f1e3a16e5 100644
--- a/PortalForecasts.R
+++ b/PortalForecasts.R
@@ -12,14 +12,16 @@ forecast_date = Sys.Date()
 
 portalr::download_observations()
 moons = get_moon_data()
+#get dates of 12 new moons following newmoon of interest
+future_moons = get_future_moons(moons)
 
 #Beginning and end of the forecast timeperiod
 most_recent_newmoon = moons$newmoonnumber[which.max(moons$period)]
 first_forecast_newmoon=most_recent_newmoon+1
 last_forecast_newmoon=first_forecast_newmoon + 11
 forecast_newmoons = first_forecast_newmoon:last_forecast_newmoon
-forecast_months=month(forecast_date %m+% months(0:11))
-forecast_years=year(forecast_date %m+% months(0:11))
+forecast_months = future_moons$month[future_moons$newmoonnumber %in% forecast_newmoons]
+forecast_years = future_moons$year[future_moons$newmoonnumber %in% forecast_newmoons]
 
 rodent_data = get_rodent_data(moons, forecast_date, filename_suffix)
 weather_data = get_weather_data(moons, rodent_data$all, first_forecast_newmoon, last_forecast_newmoon)
diff --git a/PortalHindcasts.R b/PortalHindcasts.R
index 69da71ffec..8ac9bb9cea 100644
--- a/PortalHindcasts.R
+++ b/PortalHindcasts.R
@@ -40,6 +40,9 @@ for(this_newmoon in initial_time_newmoons){
 
   moons = get_moon_data() %>%
     filter(newmoonnumber<=this_newmoon)
+  #get dates of 12 new moons following newmoon of interest
+  future_moons = get_moon_data() %>%
+    filter(newmoonnumber>this_newmoon,newmoonnumber<=this_newmoon+12)
   
   #Beginning and end of the forecast timeperiod
   most_recent_newmoon = moons$newmoonnumber[which.max(moons$period)]
@@ -47,8 +50,8 @@ for(this_newmoon in initial_time_newmoons){
   first_forecast_newmoon=most_recent_newmoon+1
   last_forecast_newmoon=first_forecast_newmoon + 11
   forecast_newmoons = first_forecast_newmoon:last_forecast_newmoon
-  forecast_months=month(most_recent_newmoon_date %m+% months(1:12))
-  forecast_years=year(most_recent_newmoon_date %m+% months(1:12))
+  forecast_months = future_moons$month[future_moons$newmoonnumber %in% forecast_newmoons]
+  forecast_years = future_moons$year[future_moons$newmoonnumber %in% forecast_newmoons]
   
   rodent_data = get_rodent_data(moons, forecast_date, filename_suffix)
   rodent_data$all = rodent_data$all %>%
diff --git a/models/model_functions.R b/models/model_functions.R
index 270d539e6c..ac5f67a28d 100644
--- a/models/model_functions.R
+++ b/models/model_functions.R
@@ -14,6 +14,24 @@ get_moon_data <- function(){
   return(moons)
 }
 
+get_future_moons <- function(moons){
+  # Get dates of future new moons from navy website
+  # Returns data.frame of newmoons in the future in the same format as the output of get_moon_data() function
+  most_recent_year = tail(moons$year,1)
+  most_recent_month = tail(moons$month,1)+1
+  newmoondates = htmltab(doc=paste("http://aa.usno.navy.mil/cgi-bin/aa_phases.pl?year=",most_recent_year,"&month=",most_recent_month,"&day=1&nump=50&format=t", sep=""),which=1)
+  newmoondates = gsub('.{6}$', '', newmoondates$"Date and Time (Universal Time)"[newmoondates$"Moon Phase" == "New Moon"])
+  newmoondates = as.Date(ymd(newmoondates, format='%Y %m %d'))
+  #Set up dataframe for new moon dates to be added
+  newmoons = data.frame(newmoonnumber = max(moons$newmoonnumber)+1:length(newmoondates),
+                        newmoondate = as.Date(newmoondates),
+                        period = NA,
+                        censusdate = as.Date(NA),
+                        year = year(newmoondates),
+                        month = month(newmoondates))
+  return(newmoons)
+}
+
 get_rodent_data <- function(moons, forecast_date, filename_suffix){
   #Period 203/newmoonnumber 217 will be when the training data timeseries
   #begins. Corresponding to Jan 1995