###Optimization of weather parameters and flight speed for Monarch Migration####
###based on mark-recapture data                                             ####  

#Packages
library(tidyverse)
library(reshape2)
library(readxl)
library(geosphere)
library(sp)

# Can't use total met data, closest station can be temporally way off! Need to use seperate met file for each butterfly (some same met data)
# NAs removed

metfilenameformonarch <- read_excel("Monarch_25HourlyPrunedMetData.xlsx",col_types = c("numeric", "numeric", "numeric","date", "numeric", "numeric", "numeric","numeric"))
attach(metfilenameformonarch)
#View(metfilenameformonarch)

# load butterfly data
DaPLmigdata<-read_excel("Tagged Data Full.xlsx",col_types = c("text","text","numeric","numeric","date","date","numeric","numeric","numeric"))
attach(DaPLmigdata)
#View(DaPLmigdata)

out<-data.frame()
Best<-data.frame()
SimHours<-data.frame()
FlightPath<-data.frame()

# 25 recoveries 

# Butterfly loop (need to run each butterfly [i] individually, check right met file)

for (i in 25:25) {

mtime=StartDate[i]   
pp<-as.numeric(EndDate[i]-StartDate[i],unit="hours")
SimHours<-pp
bCords<-matrix(0,SimHours,2)   

# Flight speed loop  
  
#Speeds:
#18 km/hr (Urquhart 1960 - book: the monarch butterfly)
#7 km/hr (Article: Notes on a large dragonfly and butterfly migration; NE Nat.; Davis and Howard / someone else with 4 similar names)
#14 km/hr (Article: Davis and Garland 2002)
#3.9 km/hr (Davis et al. 2012; laboratory setting)  

 for(ii in 1:8) {
   if (ii==1) {FlightSpeed=2.5}
   if (ii==2) {FlightSpeed=5.0}
   if (ii==3) {FlightSpeed=7.5}
   if (ii==4) {FlightSpeed=10.0}
   if (ii==5) {FlightSpeed=12.5}
   if (ii==6) {FlightSpeed=15.0}
   if (ii==7) {FlightSpeed=17.5}
   if (ii==8) {FlightSpeed=20.0}

# Temperature min. and max. based on Masters et al. 1988; Kammer 1970
  
  for(iii in 1:5){
   if (iii==1) {Tmax=70.0}
   if (iii==2) {Tmax=75.0}
   if (iii==3) {Tmax=80.0}
   if (iii==4) {Tmax=85.0}
   if (iii==5) {Tmax=90.0}
  
   for(iiii in 1:6){
    if (iiii==1) {Tmin=40.0}
    if (iiii==2) {Tmin=45.0}
    if (iiii==3) {Tmin=50.0}
    if (iiii==4) {Tmin=55.0}
    if (iiii==5) {Tmin=60.0}
    if (iiii==6) {Tmin=65.0}
     
# loop based on number of hours, start and end at noon
  
       for (j in 2:SimHours) {
       
       bCords[1,1]<-StartLon[i]
       bCords[1,2]<-StartLat[i]
       
       days<-as.integer(j/24)
       daytime<-((j+12)-(24*(days))) 
       
# Flight only between between 9 and 5 
# Model 5 on 7 to 7 ??

if (daytime<7|daytime>19) { 
   bCords[j,]<-bCords[j-1,]
   mtime=mtime+(60*60)### need time to increment
      }          
else {

## Minimizing distance to find closest met station  

bmatrix<-(c(bCords[j-1,1],bCords[j-1,2]))  
Smatrix<-matrix(c(SLongitude,SLatitude),ncol=2)

xx<-which((distHaversine(bmatrix,Smatrix))==min(distHaversine(bmatrix,Smatrix)))

##For that met station find closest time
xxx<-which(abs(mtime-StationDate[min(xx):max(xx)])==min(abs(mtime-StationDate[min(xx):max(xx)])))
   
### int is line in met data corresponding to closest time and date
int<-min(xx)+xxx-1

### rain or temp extremes =no fly that hour      
ifelse(Precipitation[int]>0.0|Temp[int]>Tmax|Temp[int]<Tmin,bCords[j,]<-bCords[j-1,],bCords[j,]<-geosphere::destPoint(bCords[j-1,1:2], b=Bearing[i], d=(FlightSpeed*1000)))
  
 
# Movement based on speed and bearing; move to next hour 

      mtime=mtime+(60*60) 
               }
    } #SimHours[i] 

      mtime=StartDate[i]   
#      SimHours<-(EndDate[i]-StartDate[i])
      
      FlightPath<-data.frame(i,FlightSpeed,Tmax,Tmin,bCords)
#      Flightpath[i]<-as.numeric(flightpath)
#     print(FlightPath)
      out<-rbind(out,FlightPath) 

    }#Tmin
  } #Tmax
 } # Flight Speed 

names(out)[5]<-paste("FLon")
names(out)[6]<-paste("FLat")
attach(out)

#pull every last point from each simulation
out.new<-as.data.frame(out[seq(SimHours, nrow(out), SimHours), ])
attach(out.new)

## Find closest simulated end to actual end lat and lon

bbmatrix<-matrix(c(EndLon[i],EndLat[i]),ncol=2)  
Fmatrix<-matrix(c(FLon,FLat),ncol=2)

yy<-which((distHaversine(bbmatrix,Fmatrix))==min(distHaversine(bbmatrix,Fmatrix)))

Best<-rbind((out.new[yy,])) #

} # Butterfly
Best