################################# BODY MASS TRAJECORIES - LAMERIS ET AL. BEHAV ECOL #################################
#####################################################################################################################

# load libraries

library(RNCEP)
library(plyr)
library(ggplot2)
library(geosphere)
library(fossil)
library(circular)
library(lme4)
library(MuMIn)
library(beepr)
library(data.table)
library(tidyr)

# Get data
setwd("E:/PAPERS/Fattening up/DATA")
gps_t<-read.csv("BMT_input_data_22122020.csv",sep=",")

gps_t$DateTimeString <- as.POSIXct(strptime(gps_t$date_time, "%d/%m/%Y %H:%M"))
gps_t$moon_date <- as.POSIXct(strptime(gps_t$moon_date, "%d/%m/%Y %H:%M"))
gps_t$DateTime=as.POSIXct(gps_t$DateTimeString,tz='UTC')
gps_t$lat <- as.numeric(gps_t$lat)
gps_t$lon <- as.numeric(gps_t$lon)
gps_t$ID <- as.factor(gps_t$device)
str(gps_t)
gps_t <- gps_t[!is.na(gps_t$device),]

## Select only GPS-fixes in NL
gps_t <- gps_t[gps_t$Area == "NL",]

## Set night-time behavior for BMT scenario "maximum foraging"

gps_t$behavior2 <- NA
gps_t$combined2 <- NA

j = 16
for (j in 1:length(unique(gps_t$device))){
  q <- gps_t[gps_t$device == unique(gps_t$device)[j],]
  p <- q[q$mooncycle == 0,]
  i = 1822
  for (i in 1:nrow(q)){
    q$behavior2[i] <- ifelse(q$daytime[i] == 0 & q$mooncycle[i] != 0 & q$mooncycle[i] != 5,
                               as.character(p$behavior[which(abs(p$DateTimeString-q$moon_date[i]) == min(abs(p$DateTimeString - q$moon_date[i])))]),
                               ifelse(q$mooncycle[i] == 5,
                                      as.character(q$behavior[i]),as.character(q$behavior[i])))
    q$combined2[i] <- ifelse(q$daytime[i] == 0 & q$mooncycle[i] != 0 & q$mooncycle[i] != 5,
                             as.character(p$combined[which(abs(p$DateTimeString-q$moon_date[i]) == min(abs(p$DateTimeString - q$moon_date[i])))]),
                             ifelse(q$mooncycle[i] == 5,
                                    as.character(q$combined[i]),as.character(q$combined[i])))
  }
  gps_t$behavior2[gps_t$device == unique(gps_t$device)[j]] <- q$behavior2
  gps_t$combined2[gps_t$device == unique(gps_t$device)[j]] <- q$combined2
  print(j)
}

## set night-time behavior for BMT scenario "no nigh-time foraging"

gps_t$behavior3 <- NA

gps_t$behavior3 <- ifelse(gps_t$daytime == 0, "inactive",as.character(gps_t$behavior))

# set BMR multipliers

behaviour <- c("inactive","grazing","active","flying")
multiplier <- c(1.5,1.6,1.9,0)
multipliers <- data.frame(behaviour, multiplier)
#multipliers$costs <- multipliers$multiplier*BMRNolet # in W -> J/s
#multipliers$costs_kJ_m <- multipliers$costs*60/1000

# Place BMR + multipliers in data
gps_t$BMR.multiplier <- multipliers$multiplier[match(gps_t$behavior, multipliers$behaviour)]
gps_t$BMR.multiplier_2 <- multipliers$multiplier[match(gps_t$behavior2, multipliers$behaviour)]
gps_t$BMR.multiplier_3 <- multipliers$multiplier[match(gps_t$behavior3, multipliers$behaviour)]

# To calculate thermoregulation costs for birds / GPS data
#
#   Using NCEP data build Thermoregulation model
### variables to take from data
#- bird weight (although here we take a fixed value, see below)
#- Temperature (degrees Celcius)
#- wind speed at 10m (m/s) -> Speed (m/s) = SQRT( u^2 + v^2) 
#- Rg = global radiation (Watt/m^2)-> nlwr + nswr

# calculate factors for to calculate windspeed at birds' level
# values taken from WMO 2008, lower saltmarsh from Wieringa 1993, higher salt marsh from De Vries 1999
habitats <- c("open sea","mud flats","lower salt marsh","open terrain / higher salt marsh","low crops")
zm <- c(0.0002,0.005,0.01,0.03,0.1)
WindSpeedFactors <- data.frame(habitats,zm)
factor=(log((.15)/zm)/(log((10)/zm)))
WindSpeedFactors$factor <- factor

# For now we use "open landscape / high salt marsh"
gps_t$windspeed.15m=gps_t$windspeed10m*WindSpeedFactors$factor[[4]]

#### Data on food intake
intake <- read.csv("food_intake_data_22122020.csv",sep=",")
intake <- intake[complete.cases(intake[12]),]
intake$ID <- paste(intake$area,intake$type)
#intake <- subset(intake,intake$year==2016)
intake$R_metab_3 <- intake$R_metab*(1/0.9)
intake$GDD2 <- intake$GDD^2
str(intake)

#NLA <- lm(R_metab_2~GDD+GDD2, data=subset(intake,intake$country == "NLA"))
summary(NLA <- lm(R_metab~GDD, data=subset(intake,intake$type == "agricultural")))
summary(NLS <- lm(R_metab~GDD, data=subset(intake,intake$type == "natural")))


#NLS1 <- mean(intake$R_metab_1[intake$type=="natural"], na.rm=TRUE)
NLS1 <- predict(NLS, newdata=data.frame(GDD=300)) #169
#NLA1 <- mean(intake$R_metab_1[intake$type=="agricultural"], na.rm=TRUE)
NLA1 <- predict(NLA, newdata=data.frame(GDD=300)) #169
#NLA1

col_gradient <- c( "schiermonnikoog" = "orange", "Zeeland" = "darkblue")

ggplot(data = intake)+
  geom_point(data = intake[intake$type == "natural",],
             aes(x=GDD,y=R_metab,colour=area),shape=21,fill=NA) +
  geom_point(data = intake[intake$type == "agricultural",],
             aes(x=GDD,y=R_metab,colour=area,fill=area),shape=21) +
  geom_segment(aes(x=300,xend=680,y=1.65,yend=2.37),linetype="dotted",colour="black") +
  geom_segment(aes(x=0,xend=300,y=1.65,yend=1.65),linetype="dotted",colour="black") +
geom_segment(aes(x=300,xend=680,y=1.75,yend=3.02),colour="black") +
  geom_segment(aes(x=0,xend=300,y=1.75,yend=1.75),colour="black") +
#  geom_smooth(aes(x=GDD,y=R_metab_1,linetype=type),
#              method=lm,formula=y~poly(x,1),se=FALSE,colour="black") +
  scale_colour_manual(name="Type",values=col_gradient) +
  scale_fill_manual(name="area",values=col_gradient) +
    ylab("Metabolisable energy intake") +
  theme(axis.title.y = element_text(size = rel(1.0), angle = 90)) +
  xlab ("GDD") +
  theme(axis.title.y = element_text(size = rel(1.5), angle = 90))  +
  theme(axis.text.y = element_text(size=15)) +
  #  theme(legend.position="none") +
  theme(axis.title.x = element_text(size = rel(1.5), angle = 0)) +
  theme(axis.text.x = element_text(size=15)) +
  theme(legend.position="none") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), 
        panel.background = element_blank(), axis.line = element_line(colour = "black")) +
  theme(axis.line.x = element_line(color="black", size = 1),
        axis.line.y = element_line(color="black", size = 1)) +
  coord_cartesian(xlim=c(150,700),)

gps_t$potential.intake <- NULL
gps_t$GDD2 <- gps_t$GDD^2


gps_t$potential.intake <- ifelse(gps_t$GDD < 301 & gps_t$combined == "NL Agricultural areas",
                                 (NLA1),
                                 ifelse(gps_t$GDD > 300 & gps_t$combined == "NL Agricultural areas",
                                        predict(NLA,newdata=gps_t),
                                        ifelse(gps_t$GDD < 301 & gps_t$combined == "NL Inland marshes" | gps_t$GDD < 301 & gps_t$combined == "NL Salt marshes",
                                               (NLS1),
                                               ifelse(gps_t$GDD > 300 & gps_t$combined == "NL Inland marshes" | gps_t$GDD > 300 & gps_t$combined == "NL Salt marshes",
                                                      predict(NLS,newdata=gps_t),NA))))

gps_t$potential.intake2 <- ifelse(gps_t$GDD < 301 & gps_t$combined2 == "NL Agricultural areas",
                                 (NLA1),
                                 ifelse(gps_t$GDD > 300 & gps_t$combined2 == "NL Agricultural areas",
                                        predict(NLA,newdata=gps_t),
                                        ifelse(gps_t$GDD < 301 & gps_t$combined2 == "NL Inland marshes" | gps_t$GDD < 301 & gps_t$combined == "NL Salt marshes",
                                               (NLS1),
                                               ifelse(gps_t$GDD > 300 & gps_t$combined2 == "NL Inland marshes" | gps_t$GDD > 300 & gps_t$combined == "NL Salt marshes",
                                                      predict(NLS,newdata=gps_t),NA))))


# airspeed

bearing <- sapply(1:(nrow(gps_t) - 1), function(i){
                  earth.bear(gps_t$lon[i], gps_t$lat[i],
                             gps_t$lon[i+1], gps_t$lat[i+1] )
                  })
bearing <- append(bearing,0)
gps_t$bearing <- bearing

# great circle distances between points (Spherical law of cosines)

gps_t$GCD <- 0
gps_t$groundspeed <- 0
#p <- subset(gps_t,gps_t$device == 6087)

for (j in 1:length(unique(gps_t$device))) {
  p <- gps_t[gps_t$device == unique(gps_t$device)[j],]
  for (i in 1:nrow(p))
  { p$GCD[i] <- distCosine(matrix(c(p$longitude[i], p$latitude[i]), ncol = 2), 
                           matrix(c(p$longitude[i+1], p$latitude[i+1]), ncol = 2)) # in meters
  p$groundspeed[i] <- ifelse(p$duration[i] > 0, (p$GCD[i]/1000)/(p$duration[i]/60), 0) # in km/h
  }
  
  gps_t$GCD[gps_t$device == unique(gps_t$device)[j]] <- p$GCD
  gps_t$groundspeed[gps_t$device == unique(gps_t$device)[j]] <- p$groundspeed
  print(j)
}

gps_t <- gps_t[!(gps_t$duration == 0),]

gps_t$GCD <- format(round(gps_t$GCD, 3), nsmall = 3)
gps_t$groundspeed <- format(round(gps_t$groundspeed, 3), nsmall = 3)
gps_t$groundspeed <- as.numeric(gps_t$groundspeed)
gps_t$GCD <- as.numeric(gps_t$GCD)

gps_t <- gps_t[!(gps_t$groundspeed > 150),]

gps_t$flight <- ifelse(gps_t$groundspeed > 1.000, 1, 0)

gps_t$time_flight <- 0

for (j in 1:length(unique(gps_t$device))) {
  p <- gps_t[gps_t$device == unique(gps_t$device)[j],]
  for (i in 1:nrow(p))
  { p$time_flight[i] <- ifelse(p$flight[i] == 1,
                               ifelse(p$flight[i-1] == 1 & p$flight[i+1] == 1,
                                      p$duration[i],
                                      ifelse(p$Area[i] == "NL",
                                             p$GCD[i]/(14.07*60),0)),0)
  }
  
  gps_t$time_flight[gps_t$device == unique(gps_t$device)[j]] <- p$time_flight
  
  print(j)
}

gps_t$time_flight <- ifelse(gps_t$groundspeed > 1.000,
                            ifelse(gps_t$Area == "NL",
                                   gps_t$GCD/(14.07*60),0),0) # 14.07 is average ground speed in m/s in NL

gps_t$duration <- gps_t$duration - gps_t$time_flight

gps_t$airspeed <- ifelse(gps_t$groundspeed > 1,
                         ifelse(gps_t$behavior == "flying",
                                sqrt((gps_t$gps_speed^2)+(gps_t$windspeed10m^2)-2*gps_t$gps_speed*gps_t$windspeed10m*cos(gps_t$bearing-gps_t$winddirection)),
                                sqrt((14.68^2)+(gps_t$windspeed10m^2)-2*14.68*gps_t$windspeed10m*cos(gps_t$bearing-gps_t$winddirection))),0)

gps_t$direction_wind_bird <- gps_t$bearing-gps_t$winddirection
gps_t$windsupport <- ifelse(gps_t$direction_wind_bird < 45 & gps_t$direction_wind_bird > -45, "head",
                            ifelse(gps_t$direction_wind_bird < -135 | gps_t$direction_wind_bird > 135, "tail", "side"))

#subset(gps_t,gps_t$index_clas == "flying" & gps_t$site == "0")
# Power is Power in watt
# Vmp is minimum power speed in meters per second
# Vmr is maxixum range speed in meters per second
# V: airspeed in m/s

# m: total body mass in kg
#m <- 1.69
# b: wingspan in meters
b <- 1.36
#b <- 1.45
# S_w: wing area in meters squared
S_w <- 0.229
#S_w <- 0.3
# rho: air density
rho <- 1.2
# g: acceleration due to gravity
g <- 9.81

#### ####


# Pennycuick (PEN) ####
# m, b, S_w, g, rho /// efficiency, C_db, S_b

# efficiency: Conversion efficiency, default value: 0.23
efficiency <- 0.23


PEN_Powerfun <- function(m,BMR,V) {
  # C_Db: body drag coefficient, default: 0.1. In other articles a value of 0.28 is used.
  C_Db <- 0.1
  # k: induced power factor, default: 1.2
  k <- 1.2
  # respiration: Respiration factor, default at 1.1
  respiration <- 1.1
  # S_b: frontal body area
  S_b <- 0.00813 * m^(2/3)
  # X_1: profile power ratio = profile power constant (8.4) divided by aspect ratio (wing span squared / wing area)
  X_1 <- 8.4 / (b^2 / S_w)
  # P_BMR: metabolic power
#  P_BMR <- (3.79 * ((m - (0))^0.723)) # replace 0 by additional weight in kg (e.g. logger)
  
  
  P_profile <- (1.05 * k^(3/4) * m^(3/2) * g^(3/2) * S_b^(1/4) * C_Db^(1/4) * rho^(-1/2) * b^(-3/2)) * X_1
  P_induced <- ((2 * k * (m * g)^2) / (V * pi * b^2 * rho))
  P_parasite <- ((rho * V^3 * S_b * C_Db)/2) # A = S_b * C_Db
  P_mech <- P_profile + P_induced + P_parasite
  P_chem <- ((P_mech / efficiency + BMR) * respiration)
  return(P_chem)
}

# thermoregulation

HMR <- function(temp,windspeed,Rg,weight) {

# input variables
# Example data
#T =5
#windspeed = 5
#Rg =100 
#weight=1300
#weight=1825 #for BMRNolet
AirSpecificHeat=1.01 #* in J/g/K *)
BodyT=40 #( specific for barnacle geese, F. Bokma unpublished)
rp=936 #(* ResistancePlumage [s/m] thesis Julia Stahl, barnacle geese, page 136 *);
#ResistancePlumage=867 #(* [s/m] thesis Julia Stahl, average barnacle / Brent Geese, page 136 *);
#rp=798 #(* [s/m] thesis Julia Stahl, value for Brent Geese, page 136 *);
#rp= 1200 # something random
EmissivityAnimal=0.98
EmissivitySurroundings=0.94
StefanBoltzmannConst=5.67*(10^-8) #* Stefan Boltzmann constant *);
AbsorptivityToRadiation=0.75
ARatio=0.29

#Thermoregulation formulas
BirdRadius=(((((485.6*weight)/1000)+592.83)/(4*pi))^0.5)/100
AirDensity=1292-(5*temp)+(0.01567*(temp^2)) #( Density of dry air [g/m^3] versus temperature in Celcius *)
#(* Heat flux [W/m^2] as a function of ambient temperature T, standard operative temperature Tes, equivalent outer resistance re *)
rr=AirDensity*AirSpecificHeat/(4*EmissivityAnimal*StefanBoltzmannConst*(temp+273)^3) # why the third power of the temperature and not the fourth? *);
rfo=307*sqrt((2*BirdRadius)/windspeed) # forced convection resistance [s/m] *)
BodySurfaceT=12.5+0.7*temp
rfr=820*(2*BirdRadius/(BodySurfaceT-temp))^.25 #* free convection resistance [s/m] *)
ra=((rfr*rfo)/(rfr+rfo)) # convection resistance [s/m] *)
re=((rr*ra)/(rr+ra)) # equivalent outer resistance
Rabs=AbsorptivityToRadiation*ARatio*Rg+EmissivitySurroundings*StefanBoltzmannConst*((temp+273)^4) #* radiation absorbed by the bird *)
Remi=EmissivityAnimal*StefanBoltzmannConst*((BodySurfaceT+273)^4) # radiation emitted by the bird *)
Te=temp+(Rabs-Remi)*(re/(AirDensity*AirSpecificHeat))
Tes=BodyT-(1+0.26*(sqrt(windspeed)))*(BodyT-Te) #* standard operative temperature *)
HeatFlux=(AirDensity*AirSpecificHeat)*((BodyT-Tes)/(rp+re))
HeatingMetabolicRate=(4*pi*(BirdRadius^2))*HeatFlux #energy flux (in J/s) needed to keep the birds body on temperature
HMR <- HeatingMetabolicRate*(60/1000)
return(HMR)
}

weights <- read.csv("biometrics_data_22122020.csv",sep=",")

weights <- subset(weights, weights$device == "6087" | weights$device == "6088" |  weights$device == "6090" | weights$device == "6091" |
                weights$device == "6097" | weights$device == "6099" |  weights$device == "6102" | weights$device == "6103" |
                weights$device == "6110" | weights$device == "6112" |  weights$device == "6115" | weights$device == "6122" |
                weights$device == "6129" | weights$device == "6131" |  weights$device == "6132" | weights$device == "6140" |
                weights$device == "6145" | weights$device == "6146")

weights$tarsus_head <- weights$tarsus*weights$head
weights$weight <- weights$tarsus_head*0.2186+254.43

#gps_t <- subset(gps_t, gps_t$DateTimeString > "2015-03-17")
#gps_t <- subset(gps_t, gps_t$DateTimeString > "2015-02-15")

gps_t <- subset(gps_t, gps_t$device == "6087" | gps_t$device == "6088" |  gps_t$device == "6090" | gps_t$device == "6091" |
                 gps_t$device == "6097" | gps_t$device == "6099" |  gps_t$device == "6102" | gps_t$device == "6103" |
                 gps_t$device == "6110" | gps_t$device == "6112" |  gps_t$device == "6115" | gps_t$device == "6122" |
                 gps_t$device == "6129" | gps_t$device == "6131" |  gps_t$device == "6132" | gps_t$device == "6140" |
               gps_t$device == "6145" | gps_t$device == "6146")

gps_t$device <- as.factor(gps_t$device)
gps_t$weight <- 0
gps_t$BMR <- 0
gps_t$activitycosts_n <- 0
gps_t$activitycosts <- 0
gps_t$flightcosts <- 0
gps_t$thermocosts <- 0
gps_t$costs <- 0
gps_t$intake <- 0
gps_t$energy <- 0
gps_t$energy2 <- 0
gps_t$eggs <- 0
gps_t$weight2 <- 0
gps_t$costs2 <- 0
gps_t$intake2 <- 0
gps_t$weight3 <- 0
gps_t$costs3 <- 0
gps_t$intake3 <- 0
gps_t$energy3 <- 0

i = 2
for (j in 1:length(unique(gps_t$device))) {
  p <- gps_t[gps_t$device == unique(gps_t$device)[j],]
  p$weight[[1]] <- weights$weight[j]
  for (i in 2:nrow(p))
  { 
    p$BMR[i] = (p$weight[i-1])*(5.54/1825)
    p$activitycosts_n[i] = (p$BMR.multiplier[i]*p$BMR[i]*p$duration[i])*60/1000
    p$flightcosts[i] = ifelse(p$time_flight[i]>0,
                              (PEN_Powerfun((p$weight[i-1]/1000),p$BMR[i],
                                            p$airspeed[i]))*(60/1000)*p$time_flight[i],
                              0)
    p$activitycosts[i] = p$activitycosts_n[i] + p$flightcosts[i]
    p$thermocosts[i] = HMR(p$temperature[i],p$windspeed.15m[i],p$radiation[i],p$weight[i-1])*p$duration[i]
    p$costs[i] = ifelse(p$activitycosts[i]>p$thermocosts[i],
                        p$activitycosts[i],
                        p$thermocosts[i])
    p$weight[i] <- p$weight[i-1]
    p$intake[i] = ifelse(p$behavior[i] == "grazing" |
                           p$behavior[i] == "active" & p$combined[i] == "NL Agricultural areas" |
                           p$behavior[i] == "active" & p$combined[i] == "NL Salt marshes" |
                           p$behavior[i] == "active" & p$combined[i] == "NL Inland marshes",
                         p$duration[i]*p$potential.intake[i],
                         0)
    p$energy[i] <- p$intake[i]-p$costs[i]
    p$weight[i] <- p$weight[i-1]+(ifelse(p$energy[i]<0,
                                         p$energy[i]*(1/29),      # 29kJ/g or 34.3 kJ/g
                                         p$energy[i]*0.8*(1/29)))   # efficiency 0.8 or 45.66kJ/g
  }
  
  gps_t$weight[gps_t$device == unique(gps_t$device)[j]] <- p$weight
  gps_t$costs[gps_t$device == unique(gps_t$device)[j]] <- p$costs 
  gps_t$activitycosts[gps_t$device == unique(gps_t$device)[j]] <- p$activitycosts 
  gps_t$flightcosts[gps_t$device == unique(gps_t$device)[j]] <- p$flightcosts
  gps_t$thermocosts[gps_t$device == unique(gps_t$device)[j]] <- p$thermocosts
  gps_t$eggs[gps_t$device == unique(gps_t$device)[j]] <- p$eggs
  gps_t$intake[gps_t$device == unique(gps_t$device)[j]] <- p$intake
  gps_t$energy[gps_t$device == unique(gps_t$device)[j]] <- p$energy
  gps_t$BMR[gps_t$device == unique(gps_t$device)[j]] <- p$BMR
  
  print(j)
}   
beep(4)

#### "maximum foraging" BMPs

i = 1659
j = 1
for (j in 1:length(unique(gps_t$device))) {
  p <- gps_t[gps_t$device == unique(gps_t$device)[j],]
  p$weight2[[1]] <- weights$weight[j]
  for (i in 2:nrow(p))
  { 
    p$BMR[i] = (p$weight2[i-1])*(5.54/1825)
    p$activitycosts_n[i] = (p$BMR.multiplier_2[i]*p$BMR[i]*p$duration[i])*60/1000
    p$flightcosts[i] = ifelse(p$time_flight[i]>0,
                              (PEN_Powerfun((p$weight2[i-1]/1000),p$BMR[i],
                                            p$airspeed[i]))*(60/1000)*p$time_flight[i],
                              0)
    p$activitycosts[i] = p$activitycosts_n[i] + p$flightcosts[i]
    p$thermocosts[i] = HMR(p$temperature[i],p$windspeed.15m[i],p$radiation[i],p$weight2[i-1])*p$duration[i]
    p$costs2[i] = ifelse(p$activitycosts[i]>p$thermocosts[i],
                        p$activitycosts[i],
                        p$thermocosts[i])
    p$weight2[i] <- p$weight2[i-1]
    p$intake2[i] = ifelse(p$behavior2[i] == "grazing" |
                           p$behavior2[i] == "active" & p$combined2[i] == "NL Agricultural areas" |
                           p$behavior2[i] == "active" & p$combined2[i] == "NL Salt marshes" |
                           p$behavior2[i] == "active" & p$combined2[i] == "NL Inland marshes",
                         p$duration[i]*p$potential.intake2[i],
                         0)
    p$energy2[i] <- p$intake2[i]-p$costs2[i]
    p$weight2[i] <- p$weight2[i-1]+(ifelse(p$energy2[i]<0,
                                         p$energy2[i]*(1/29),      # 29kJ/g or 34.3 kJ/g
                                         p$energy2[i]*0.8*(1/29)))   # efficiency 0.8 or 45.66kJ/g
  }
  
  gps_t$weight2[gps_t$device == unique(gps_t$device)[j]] <- p$weight2
  gps_t$costs2[gps_t$device == unique(gps_t$device)[j]] <- p$costs2 
  gps_t$activitycosts[gps_t$device == unique(gps_t$device)[j]] <- p$activitycosts 
  gps_t$flightcosts[gps_t$device == unique(gps_t$device)[j]] <- p$flightcosts
  gps_t$thermocosts[gps_t$device == unique(gps_t$device)[j]] <- p$thermocosts
  gps_t$eggs[gps_t$device == unique(gps_t$device)[j]] <- p$eggs
  gps_t$intake2[gps_t$device == unique(gps_t$device)[j]] <- p$intake2
  gps_t$energy2[gps_t$device == unique(gps_t$device)[j]] <- p$energy2
  gps_t$BMR[gps_t$device == unique(gps_t$device)[j]] <- p$BMR
  
  print(j)
}   
beep(4)

#### No night grazing

i = 2
for (j in 1:length(unique(gps_t$device))) {
  p <- gps_t[gps_t$device == unique(gps_t$device)[j],]
  p$weight3[[1]] <- weights$weight[j]
  for (i in 2:nrow(p))
  { 
    p$BMR[i] = (p$weight3[i-1])*(5.54/1825) # method Adriaan
    p$activitycosts_n[i] = (p$BMR.multiplier_3[i]*p$BMR[i]*p$duration[i])*60/1000
    p$flightcosts[i] = ifelse(p$time_flight[i]>0,
                              (PEN_Powerfun((p$weight3[i-1]/1000),p$BMR[i],
                                            p$airspeed[i]))*(60/1000)*p$time_flight[i],
                              0)
    p$activitycosts[i] = p$activitycosts_n[i] + p$flightcosts[i]
    p$thermocosts[i] = HMR(p$temperature[i],p$windspeed.15m[i],p$radiation[i],p$weight3[i-1])*p$duration[i]
    p$costs3[i] = ifelse(p$activitycosts[i]>p$thermocosts[i],
                         p$activitycosts[i],
                         p$thermocosts[i])
    p$weight3[i] <- p$weight3[i-1]
    p$intake3[i] = ifelse(p$behavior3[i] == "grazing" |
                            p$behavior3[i] == "active" & p$combined[i] == "NL Agricultural areas" |
                            p$behavior3[i] == "active" & p$combined[i] == "NL Salt marshes" |
                            p$behavior3[i] == "active" & p$combined[i] == "NL Inland marshes",
                          p$duration[i]*p$potential.intake[i],
                          0)
    p$energy3[i] <- p$intake3[i]-p$costs3[i]
    p$weight3[i] <- p$weight3[i-1]+(ifelse(p$energy3[i]<0,
                                           p$energy3[i]*(1/29),      # 29kJ/g or 34.3 kJ/g
                                           p$energy3[i]*0.8*(1/29)))   # efficiency 0.8 or 45.66kJ/g
  }
  
  gps_t$weight3[gps_t$device == unique(gps_t$device)[j]] <- p$weight3
  gps_t$costs3[gps_t$device == unique(gps_t$device)[j]] <- p$costs3
  gps_t$activitycosts[gps_t$device == unique(gps_t$device)[j]] <- p$activitycosts 
  gps_t$flightcosts[gps_t$device == unique(gps_t$device)[j]] <- p$flightcosts
  gps_t$thermocosts[gps_t$device == unique(gps_t$device)[j]] <- p$thermocosts
  gps_t$eggs[gps_t$device == unique(gps_t$device)[j]] <- p$eggs
  gps_t$intake3[gps_t$device == unique(gps_t$device)[j]] <- p$intake3
  gps_t$energy3[gps_t$device == unique(gps_t$device)[j]] <- p$energy3
  gps_t$BMR[gps_t$device == unique(gps_t$device)[j]] <- p$BMR
  
  print(j)
}   
beep(4)

library(dplyr)

gps3 <- gps_t
summary_gps2 <- ddply(gps3, c("device","date","day_number"), summarise,
                     thermo = sum(thermocosts, na.rm=TRUE),
                     activity = sum(activitycosts, na.rm=TRUE),
                     flight = sum(flightcosts, na.rm = TRUE),
                     total   = sum(costs, na.rm=TRUE),
                     intake = sum(intake, na.rm=TRUE),
                     energy = sum(energy, na.rm = TRUE),
                     weight = mean(weight, na.rm = TRUE),
                     weight2 = mean(weight2, na.rm = TRUE),
                     weight3 = mean(weight3, na.rm = TRUE))
summary_gps2 <- summary_gps2[order(as.Date(summary_gps2$date, format="%d-%m-%Y")),]
summary_gps2$date <- as.Date(summary_gps2$date, "%d-%m-%Y")
summary_gps2 <- summary_gps2[order(summary_gps2$device),]
summary_gps2$DateStringDevice <- paste(summary_gps2$date, summary_gps2$device)

summary_gps3 <- ddply(summary_gps2[summary_gps2$device != "6090",], c("date","day_number"), summarise,
                      thermo1 = mean(thermo, na.rm=TRUE),
                      activity1 = mean(activity, na.rm=TRUE),
                      flight1 = mean(flight, na.rm = TRUE),
                      total1   = mean(total, na.rm=TRUE),
                      intake1 = mean(intake, na.rm=TRUE),
                      energy1 = mean(energy, na.rm = TRUE),
                      weight1 = mean(weight, na.rm = TRUE),
                      weight2 = mean(weight2, na.rm = TRUE),
                      weight3 = mean(weight3, na.rm = TRUE))


ggplot()+
  geom_line(data=summary_gps3, aes(x=day_number, y = weight2),colour="darkgreen",size=2) +
  geom_line(data=summary_gps3, aes(x=day_number, y = weight1),colour="orange",size=2) +
  geom_line(data=summary_gps3, aes(x=day_number, y = weight3),colour="blue",size=2) +
  ylab("Weight (g)") +
  theme(axis.title.y = element_text(size = rel(1.0), angle = 90)) +
  xlab ("Julian date") +
  theme(axis.title.y = element_text(size = rel(1.5), angle = 90))  +
  theme(axis.text.y = element_text(size=15)) +
  #  theme(legend.position="none") +
  theme(axis.title.x = element_text(size = rel(1.5), angle = 0)) +
  theme(axis.text.x = element_text(size=15)) +
  #theme(legend.position="none") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), 
        panel.background = element_blank(), axis.line = element_line(colour = "black")) +
  theme(axis.line.x = element_line(color="black", size = 1),
        axis.line.y = element_line(color="black", size = 1))+
  #coord_cartesian(ylim=c(0,700),xlim=c(1,170))
  coord_cartesian(ylim=c(800,3250),xlim=c(1,170))