---
title: "An evolvable celestial homing mechanism - Code"
author: "Yuval Werber"
date: "2026-07-26"
output: html_document
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
knitr::opts_chunk$set(root.dir="")
setwd("")
library(pracma)
library(suncalc)
library(tidyverse)
library(dplyr)
library(hms)
library(geosphere)
library(circular)
library(RcppAlgos)
library(ggpubr)
library(cowplot)
library(ggnewscale)

```

Preparing simulation datasets

```{r}
#Randomly sample each vector along the relevant tange and then join them manually
##Lats and Lons
Latitude_displacement<-round(runif(n = 1000000, min = 0, max = 90),digits=2)
Longitude_displacement<-round(runif(n = 1000000, min = 0, max = 180),digits=2)
Distance<-round(runif(n = 1000000, min = 0, max = 100000),digits=0)#in meters
azimuth<-round(runif(n = 1000000, min = 0, max = 359.99),digits=2)
##Date_time
# 1. Define start and end times
start_time <- as.POSIXct("2023-01-01 00:00:00", tz = "UTC")
end_time <- as.POSIXct("2024-01-01 00:00:00", tz = "UTC")

# 2. Define the number of random timestamps you need
n_timestamps <- 1000000

# 3. Set a seed for reproducibility (optional)
set.seed(123)

# 4. Generate 'n' random numbers (seconds) between the integer representation of 
#    start and end times using runif()
random_seconds <- runif(n = n_timestamps, 
                        min = as.numeric(start_time), 
                        max = as.numeric(end_time))

# 5. Convert the random seconds back to POSIXct format
random_timestamps <- as.POSIXct(random_seconds, origin = "1970-01-01", tz = "UTC")

# 6. Create a data frame with the new column
df <- data.frame(
  id = 1:n_timestamps,
  timestamp = random_timestamps
)
Date_time<-df$timestamp

Time_shift<-round(runif(n = 1000000, min = -6, max = 6),digits=0)
##Raw table-using expand grid not scaleable because makes enormous tables
#This seems to break at wider ranges of longitue, test whiy this is 
##First trial - 0-90 Lat, 0-40 lon
Table<-data.frame(Latitude_displacement=Latitude_displacement,Longitude_displacement=Longitude_displacement,Distance=Distance,azimuth=azimuth,Date_time=Date_time,Time_shift=Time_shift)
#Table<-expand.grid(Latitude_displacement = seq(0,90,10), Longitude_displacement = seq(0, 180,10),
#            Date_time=seq(from = start_time, to = end_time, by = "1 hour"),Time_shift=seq(from =0, to=6,by=1))

#set.seed(10)
#Table$Latitude_home<-sample(Table$Latitude_displacement)
#Table$Longitude_home<-sample(Table$Longitude_displacement)


Table$time_epoch<-as.numeric(Table$Date_time)
##Make a date column
Table$Date<-as.Date(Table$Date_time)


# Convert to an hms object (only stores time internally, uses a dummy date)
Table$Hour<-  as_hms(Table$Date_time)

# Calculate destination points
destination_points <- destPoint(
  p = cbind(Table$Longitude_displacement, Table$Latitude_displacement),
  b = Table$azimuth,
  d = Table$Distance
)
destination_points<-as.data.frame(destination_points)
Table$Latitude_home<-destination_points$lat
Table$Longitude_home<-destination_points$lon
Table<-Table[c(10,11,1:9)]
Table$Distance_km<-Table$Distance/1000

##Get sunrise and sunset tHour##Get sunrise and sunset times
lat<-Table$Latitude_displacement
lon<-Table$Longitude_displacement
date<-Table$Date

Sun_phase_prep_table <- data.frame(cbind(date,lat,lon))
Sun_phase_prep_table$date <- as.Date(Sun_phase_prep_table$date, origin="1970-01-01")


Sun_phase_table <- getSunlightTimes(data =
Sun_phase_prep_table, tz = "UTC", keep = c("sunrise", "sunset"))
Table$Sunrise<-Sun_phase_table$sunrise
Table$sunrise_epoch<-as.numeric(Table$Sunrise)
Table$Sunset<-Sun_phase_table$sunset
Table$sunset_epoch<-as.numeric(Table$Sunset)
##Calculate day length
Table$Day_length<-difftime(Table$Sunset, Table$Sunrise, units = "hours")
Table$Day_length<-as.numeric(Table$Day_length)
##Calculate Day/night
Table$Daynight<-ifelse(Table$time_epoch<Table$sunrise_epoch|Table$time_epoch>Table$sunset_epoch,"Night","Day")
Table<-Table%>%filter(Daynight=="Day")


##Calculate time from sunrise
Table$Time_from_sunrise<-difftime(Table$Date_time, Table$Sunrise, units = "hours")
Table$Time_from_sunrise<-as.numeric(Table$Time_from_sunrise)

##Calculate the observed hour angle based on the Expected formula (which uses time) and the difference in longitude between displacement and home.
Table$Observed_hour_angle<-Table$Time_from_sunrise*180/Table$Day_length+(Table$Longitude_displacement-Table$Longitude_home)
##Filter only home-displacement combinations that produce obsered hour angles between 0-180, below which and above which the sun is not in the sky
##at dispalcement. Expected hour angle doesnt need filtering becuase the forula restricts it to 0-180
Table<-Table%>%filter(as.numeric(Observed_hour_angle)>0)
Table<-Table%>%filter(as.numeric(Observed_hour_angle)<180)


Table$time_epoch<-NULL
Table$sunrise_epoch<-NULL
Table$sunset_epoch<-NULL
Table$Daynight<-NULL
#Table<-Table[,c(1,2,5,6,3,7,8,4,9:13)]
rm(Sun_phase_prep_table,Sun_phase_table,date,end_time,start_time,lat,lon,Latitude_displacement,Latitude_home,Longitude_displacement, Longitude_home,df,n_timestamps,Date_time,random_seconds,random_timestamps,Time_shift,destination_points,azimuth,Distance)

##Add distance in km and 3 azimuth types using the geosphere package as reference to compare the output of our functions



# 1. Calculate the distance (in meters by default)
# Use distGeo for the most accurate ellipsoidal distance
library(geosphere)
library(dplyr)
#For northern Hemisphere
# Calculate the distance for each row
Table <- Table %>%
  dplyr::rowwise() %>% # perform the operation row by row
  dplyr::mutate(distance_km_geo = as.numeric(distGeo(p1 = c(Longitude_displacement, Latitude_displacement), 
                                         p2 = c(Longitude_home, Latitude_home))/1000)) %>%
  ungroup() # Return to a standard data frame

# Calculate the distance for each row
Table <- Table %>%
  dplyr::rowwise() %>% # perform the operation row by row
  dplyr::mutate(distance_km_rhumb_geo = as.numeric(distRhumb(p1 = c(Longitude_displacement, Latitude_displacement), 
                                         p2 = c(Longitude_home, Latitude_home))/1000)) %>%
  ungroup() # Return to a standard data frame


# 2. Calculate the initial azimuth (bearing in degrees from North)
# The bearing function returns the initial bearing
Table <- Table %>%
  dplyr::rowwise() %>% # perform the operation row by row
  dplyr::mutate(Azimuth_initial_geo = ?bearing(p1 = c(Longitude_displacement, Latitude_displacement), 
                                         p2 = c(Longitude_home, Latitude_home))) %>%
  ungroup() # Return to a standard data frame

Table$Azimuth_initial_geo <- (Table$Azimuth_initial_geo + 360) %% 360


# 3. Calculate the final azimuth (bearing when arriving at the destination)
Table <- Table %>%
  dplyr::rowwise() %>% # perform the operation row by row
  dplyr::mutate(Azimuth_final_geo = finalBearing(p1 = c(Longitude_displacement, Latitude_displacement), 
                                         p2 = c(Longitude_home, Latitude_home))) %>%
  ungroup() # Return to a standard data frame

Table$Azimuth_final_geo <- (Table$Azimuth_final_geo + 360) %% 360


# 4. Calculate the Rhump azimuth (the bearing of the rhumb line- a line that will bring the animal home in a fixed azimuth)
Table <- Table %>%
  dplyr::rowwise() %>% # perform the operation row by row
  dplyr::mutate(Azimuth_rhumb_geo = bearingRhumb(p1 = c(Longitude_displacement, Latitude_displacement), 
                                         p2 = c(Longitude_home, Latitude_home))) %>%
  ungroup() # Return to a standard data frame

Table$Azimuth_rhumb_geo <- (Table$Azimuth_rhumb_geo + 360) %% 360
Table_100<-Table
Table_100$range<-rep(100,nrow(Table_100))

#1000 km 






#Think of a way to creare a table without producing all the combinations:
#Randomly sample each vector along the relevant tange and then join them manually
##Lats and Lons
Latitude_displacement<-round(runif(n = 1000000, min = 0, max = 90),digits=2)
Longitude_displacement<-round(runif(n = 1000000, min = 0, max = 180),digits=2)
Distance<-round(runif(n = 1000000, min = 100000, max = 1000000),digits=0)#in meters
azimuth<-round(runif(n = 1000000, min = 0, max = 359.99),digits=2)
##Date_time
# 1. Define start and end times
start_time <- as.POSIXct("2023-01-01 00:00:00", tz = "UTC")
end_time <- as.POSIXct("2024-01-01 00:00:00", tz = "UTC")

# 2. Define the number of random timestamps you need
n_timestamps <- 1000000

# 3. Set a seed for reproducibility (optional)
set.seed(123)

# 4. Generate 'n' random numbers (seconds) between the integer representation of 
#    start and end times using runif()
random_seconds <- runif(n = n_timestamps, 
                        min = as.numeric(start_time), 
                        max = as.numeric(end_time))

# 5. Convert the random seconds back to POSIXct format
random_timestamps <- as.POSIXct(random_seconds, origin = "1970-01-01", tz = "UTC")

# 6. Create a data frame with the new column
df <- data.frame(
  id = 1:n_timestamps,
  timestamp = random_timestamps
)
Date_time<-df$timestamp

Time_shift<-round(runif(n = 1000000, min = -6, max = 6),digits=0)
##Raw table-using expand grid not scaleable because makes enormous tables
#This seems to break at wider ranges of longitue, test whiy this is 
##First trial - 0-90 Lat, 0-40 lon
Table<-data.frame(Latitude_displacement=Latitude_displacement,Longitude_displacement=Longitude_displacement,Distance=Distance,azimuth=azimuth,Date_time=Date_time,Time_shift=Time_shift)
#Table<-expand.grid(Latitude_displacement = seq(0,90,10), Longitude_displacement = seq(0, 180,10),
#            Date_time=seq(from = start_time, to = end_time, by = "1 hour"),Time_shift=seq(from =0, to=6,by=1))

#set.seed(10)
#Table$Latitude_home<-sample(Table$Latitude_displacement)
#Table$Longitude_home<-sample(Table$Longitude_displacement)


Table$time_epoch<-as.numeric(Table$Date_time)
##Make a date column
Table$Date<-as.Date(Table$Date_time)


# Convert to an hms object (only stores time internally, uses a dummy date)
Table$Hour<-  as_hms(Table$Date_time)

# Calculate destination points
destination_points <- destPoint(
  p = cbind(Table$Longitude_displacement, Table$Latitude_displacement),
  b = Table$azimuth,
  d = Table$Distance
)
destination_points<-as.data.frame(destination_points)
Table$Latitude_home<-destination_points$lat
Table$Longitude_home<-destination_points$lon
Table<-Table[c(10,11,1:9)]
Table$Distance_km<-Table$Distance/1000

##Get sunrise and sunset tHour##Get sunrise and sunset times
lat<-Table$Latitude_displacement
lon<-Table$Longitude_displacement
date<-Table$Date

Sun_phase_prep_table <- data.frame(cbind(date,lat,lon))
Sun_phase_prep_table$date <- as.Date(Sun_phase_prep_table$date, origin="1970-01-01")


Sun_phase_table <- getSunlightTimes(data =
Sun_phase_prep_table, tz = "UTC", keep = c("sunrise", "sunset"))
Table$Sunrise<-Sun_phase_table$sunrise
Table$sunrise_epoch<-as.numeric(Table$Sunrise)
Table$Sunset<-Sun_phase_table$sunset
Table$sunset_epoch<-as.numeric(Table$Sunset)
##Calculate day length
Table$Day_length<-difftime(Table$Sunset, Table$Sunrise, units = "hours")
Table$Day_length<-as.numeric(Table$Day_length)
##Calculate Day/night
Table$Daynight<-ifelse(Table$time_epoch<Table$sunrise_epoch|Table$time_epoch>Table$sunset_epoch,"Night","Day")
Table<-Table%>%filter(Daynight=="Day")


##Calculate time from sunrise
Table$Time_from_sunrise<-difftime(Table$Date_time, Table$Sunrise, units = "hours")
Table$Time_from_sunrise<-as.numeric(Table$Time_from_sunrise)

##Calculate the observed hour angle based on the Expected formula (which uses time) and the difference in longitude between displacement and home.
Table$Observed_hour_angle<-Table$Time_from_sunrise*180/Table$Day_length+(Table$Longitude_displacement-Table$Longitude_home)
##Filter only home-displacement combinations that produce obsered hour angles between 0-180, below which and above which the sun is not in the sky
##at dispalcement. Expected hour angle doesnt need filtering becuase the forula restricts it to 0-180
Table<-Table%>%filter(as.numeric(Observed_hour_angle)>0)
Table<-Table%>%filter(as.numeric(Observed_hour_angle)<180)


Table$time_epoch<-NULL
Table$sunrise_epoch<-NULL
Table$sunset_epoch<-NULL
Table$Daynight<-NULL
#Table<-Table[,c(1,2,5,6,3,7,8,4,9:13)]
rm(Sun_phase_prep_table,Sun_phase_table,date,end_time,start_time,lat,lon,Latitude_displacement,Latitude_home,Longitude_displacement, Longitude_home,df,n_timestamps,Date_time,random_seconds,random_timestamps,Time_shift,destination_points,azimuth,Distance)

##Add distance in km and 3 azimuth types using the geosphere package as reference to compare the output of our functions



# 1. Calculate the distance (in meters by default)
# Use distGeo for the most accurate ellipsoidal distance
library(geosphere)
library(dplyr)
# Calculate the distance for each row
Table <- Table %>%
  dplyr::rowwise() %>% # perform the operation row by row
  dplyr::mutate(distance_km_geo = as.numeric(distGeo(p1 = c(Longitude_displacement, Latitude_displacement), 
                                         p2 = c(Longitude_home, Latitude_home))/1000)) %>%
  ungroup() # Return to a standard data frame

# Calculate the distance for each row
Table <- Table %>%
  dplyr::rowwise() %>% # perform the operation row by row
  dplyr::mutate(distance_km_rhumb_geo = as.numeric(distRhumb(p1 = c(Longitude_displacement, Latitude_displacement), 
                                         p2 = c(Longitude_home, Latitude_home))/1000)) %>%
  ungroup() # Return to a standard data frame


# 2. Calculate the initial azimuth (bearing in degrees from North)
# The bearing function returns the initial bearing
Table <- Table %>%
  dplyr::rowwise() %>% # perform the operation row by row
  dplyr::mutate(Azimuth_initial_geo = bearing(p1 = c(Longitude_displacement, Latitude_displacement), 
                                         p2 = c(Longitude_home, Latitude_home))) %>%
  ungroup() # Return to a standard data frame

Table$Azimuth_initial_geo <- (Table$Azimuth_initial_geo + 360) %% 360


# 3. Calculate the final azimuth (bearing when arriving at the destination)
Table <- Table %>%
  dplyr::rowwise() %>% # perform the operation row by row
  dplyr::mutate(Azimuth_final_geo = finalBearing(p1 = c(Longitude_displacement, Latitude_displacement), 
                                         p2 = c(Longitude_home, Latitude_home))) %>%
  ungroup() # Return to a standard data frame

Table$Azimuth_final_geo <- (Table$Azimuth_final_geo + 360) %% 360


# 4. Calculate the Rhump azimuth (the bearing of the rhumb line- a line that will bring the animal home in a fixed azimuth)
Table <- Table %>%
  dplyr::rowwise() %>% # perform the operation row by row
  dplyr::mutate(Azimuth_rhumb_geo = bearingRhumb(p1 = c(Longitude_displacement, Latitude_displacement), 
                                         p2 = c(Longitude_home, Latitude_home))) %>%
  ungroup() # Return to a standard data frame

Table$Azimuth_rhumb_geo <- (Table$Azimuth_rhumb_geo + 360) %% 360

Table_1000<-Table

Table_1000$range<-rep(1000,nrow(Table_1000))

##20000 km

#Think of a way to creare a table without producing all the combinations:
#Randomly sample each vector along the relevant tange and then join them manually
##Lats and Lons
Latitude_displacement<-round(runif(n = 1000000, min = 0, max = 90),digits=2)
Longitude_displacement<-round(runif(n = 1000000, min = 0, max = 180),digits=2)
Distance<-round(runif(n = 1000000, min = 1000000, max = 20000000),digits=0)#in meters
azimuth<-round(runif(n = 1000000, min = 0, max = 359.99),digits=2)
##Date_time
# 1. Define start and end times
start_time <- as.POSIXct("2023-01-01 00:00:00", tz = "UTC")
end_time <- as.POSIXct("2024-01-01 00:00:00", tz = "UTC")

# 2. Define the number of random timestamps you need
n_timestamps <- 1000000

# 3. Set a seed for reproducibility (optional)
set.seed(123)

# 4. Generate 'n' random numbers (seconds) between the integer representation of 
#    start and end times using runif()
random_seconds <- runif(n = n_timestamps, 
                        min = as.numeric(start_time), 
                        max = as.numeric(end_time))

# 5. Convert the random seconds back to POSIXct format
random_timestamps <- as.POSIXct(random_seconds, origin = "1970-01-01", tz = "UTC")

# 6. Create a data frame with the new column
df <- data.frame(
  id = 1:n_timestamps,
  timestamp = random_timestamps
)
Date_time<-df$timestamp

Time_shift<-round(runif(n = 1000000, min = -6, max = 6),digits=0)
##Raw table-using expand grid not scaleable because makes enormous tables
#This seems to break at wider ranges of longitue, test whiy this is 
##First trial - 0-90 Lat, 0-40 lon
Table<-data.frame(Latitude_displacement=Latitude_displacement,Longitude_displacement=Longitude_displacement,Distance=Distance,azimuth=azimuth,Date_time=Date_time,Time_shift=Time_shift)
#Table<-expand.grid(Latitude_displacement = seq(0,90,10), Longitude_displacement = seq(0, 180,10),
#            Date_time=seq(from = start_time, to = end_time, by = "1 hour"),Time_shift=seq(from =0, to=6,by=1))

#set.seed(10)
#Table$Latitude_home<-sample(Table$Latitude_displacement)
#Table$Longitude_home<-sample(Table$Longitude_displacement)


Table$time_epoch<-as.numeric(Table$Date_time)
##Make a date column
Table$Date<-as.Date(Table$Date_time)


# Convert to an hms object (only stores time internally, uses a dummy date)
Table$Hour<-  as_hms(Table$Date_time)

# Calculate destination points
destination_points <- destPoint(
  p = cbind(Table$Longitude_displacement, Table$Latitude_displacement),
  b = Table$azimuth,
  d = Table$Distance
)
destination_points<-as.data.frame(destination_points)
Table$Latitude_home<-destination_points$lat
Table$Longitude_home<-destination_points$lon
Table<-Table[c(10,11,1:9)]
Table$Distance_km<-Table$Distance/1000

##Get sunrise and sunset tHour##Get sunrise and sunset times
lat<-Table$Latitude_displacement
lon<-Table$Longitude_displacement
date<-Table$Date

Sun_phase_prep_table <- data.frame(cbind(date,lat,lon))
Sun_phase_prep_table$date <- as.Date(Sun_phase_prep_table$date, origin="1970-01-01")


Sun_phase_table <- getSunlightTimes(data =
Sun_phase_prep_table, tz = "UTC", keep = c("sunrise", "sunset"))
Table$Sunrise<-Sun_phase_table$sunrise
Table$sunrise_epoch<-as.numeric(Table$Sunrise)
Table$Sunset<-Sun_phase_table$sunset
Table$sunset_epoch<-as.numeric(Table$Sunset)
##Calculate day length
Table$Day_length<-difftime(Table$Sunset, Table$Sunrise, units = "hours")
Table$Day_length<-as.numeric(Table$Day_length)
##Calculate Day/night
Table$Daynight<-ifelse(Table$time_epoch<Table$sunrise_epoch|Table$time_epoch>Table$sunset_epoch,"Night","Day")
Table<-Table%>%filter(Daynight=="Day")


##Calculate time from sunrise
Table$Time_from_sunrise<-difftime(Table$Date_time, Table$Sunrise, units = "hours")
Table$Time_from_sunrise<-as.numeric(Table$Time_from_sunrise)

##Calculate the observed hour angle based on the Expected formula (which uses time) and the difference in longitude between displacement and home.
Table$Observed_hour_angle<-Table$Time_from_sunrise*180/Table$Day_length+(Table$Longitude_displacement-Table$Longitude_home)
##Filter only home-displacement combinations that produce obsered hour angles between 0-180, below which and above which the sun is not in the sky
##at dispalcement. Expected hour angle doesnt need filtering becuase the forula restricts it to 0-180
Table<-Table%>%filter(as.numeric(Observed_hour_angle)>0)
Table<-Table%>%filter(as.numeric(Observed_hour_angle)<180)


Table$time_epoch<-NULL
Table$sunrise_epoch<-NULL
Table$sunset_epoch<-NULL
Table$Daynight<-NULL
#Table<-Table[,c(1,2,5,6,3,7,8,4,9:13)]
rm(Sun_phase_prep_table,Sun_phase_table,date,end_time,start_time,lat,lon,Latitude_displacement,Latitude_home,Longitude_displacement, Longitude_home,df,n_timestamps,Date_time,random_seconds,random_timestamps,Time_shift,destination_points,azimuth,Distance)

##Add distance in km and 3 azimuth types using the geosphere package as reference to compare the output of our functions

# 1. Calculate the distance (in meters by default)
# Use distGeo for the most accurate ellipsoidal distance
library(geosphere)
library(dplyr)
# Calculate the distance for each row
Table <- Table %>%
  dplyr::rowwise() %>% # perform the operation row by row
  dplyr::mutate(distance_km_geo = as.numeric(distGeo(p1 = c(Longitude_displacement, Latitude_displacement), 
                                         p2 = c(Longitude_home, Latitude_home))/1000)) %>%
  ungroup() # Return to a standard data frame

# Calculate the distance for each row
Table <- Table %>%
  dplyr::rowwise() %>% # perform the operation row by row
  dplyr::mutate(distance_km_rhumb_geo = as.numeric(distRhumb(p1 = c(Longitude_displacement, Latitude_displacement), 
                                         p2 = c(Longitude_home, Latitude_home))/1000)) %>%
  ungroup() # Return to a standard data frame


# 2. Calculate the initial azimuth (bearing in degrees from North)
# The bearing function returns the initial bearing
Table <- Table %>%
  dplyr::rowwise() %>% # perform the operation row by row
  dplyr::mutate(Azimuth_initial_geo = bearing(p1 = c(Longitude_displacement, Latitude_displacement), 
                                         p2 = c(Longitude_home, Latitude_home))) %>%
  ungroup() # Return to a standard data frame

Table$Azimuth_initial_geo <- (Table$Azimuth_initial_geo + 360) %% 360


# 3. Calculate the final azimuth (bearing when arriving at the destination)
Table <- Table %>%
  dplyr::rowwise() %>% # perform the operation row by row
  dplyr::mutate(Azimuth_final_geo = finalBearing(p1 = c(Longitude_displacement, Latitude_displacement), 
                                         p2 = c(Longitude_home, Latitude_home))) %>%
  ungroup() # Return to a standard data frame

Table$Azimuth_final_geo <- (Table$Azimuth_final_geo + 360) %% 360


# 4. Calculate the Rhump azimuth (the bearing of the rhumb line- a line that will bring the animal home in a fixed azimuth)
Table <- Table %>%
  dplyr::rowwise() %>% # perform the operation row by row
  dplyr::mutate(Azimuth_rhumb_geo = bearingRhumb(p1 = c(Longitude_displacement, Latitude_displacement), 
                                         p2 = c(Longitude_home, Latitude_home))) %>%
  ungroup() # Return to a standard data frame

Table$Azimuth_rhumb_geo <- (Table$Azimuth_rhumb_geo + 360) %% 360

Table_20000<-Table

Table_20000$range<-rep(20000,nrow(Table_20000))

#Table<-rbind(Table_100,Table_1000,Table_20000)
#rm(Table)
```

##Now apply the functions for 100,000 simulations of each range separately, other wise it will crash to computer.


A chunck that defines all the functoin we want to test 
For explanation about hte function go to script v4 and look by name


```{r}
#Full Binary function - No time shift:
########################################


#Parameter definition

#NCP_altitude_home<-as.circular(40, type = "angles", units ="degrees", template = "geographics")
#NCP_altitude_displacement<-as.circular(40, type = "angles", units ="degrees", template = "geographics")
#Day_length<-12#in hours
#Time_from_sunrise_at_displacement<-2#in hours
#Observed_hour_angle_displacement<-as.circular(30, type = "angles", units ="degrees", template = "geographics")
#Observed_hour_angle_displacement_no_shift<-as.circular(30, type = "angles", units ="degrees", template = "geographics")
#Time_shift<-1#in hours

#Function:

N_Binary_No_Shift<- function(NCP_altitude_displacement,NCP_altitude_home,Time_from_sunrise_at_displacement,Day_length,Observed_hour_angle_displacement) {
  
    NS_correction<-as.circular(ifelse(NCP_altitude_displacement>NCP_altitude_home,180,ifelse(NCP_altitude_displacement<NCP_altitude_home,0,NA)) , type = "angles", units = "degrees", template = "geographics")#NS correction, NA marks NCP equivalence, meaning no correction needee
    
    Expected_hour_angle_displacement<-as.circular(Time_from_sunrise_at_displacement*180/Day_length, type = "angles", units ="degrees", template = "geographics")#Expected hor angle calculation
    
    EW_correction<-as.circular(ifelse(Observed_hour_angle_displacement<Expected_hour_angle_displacement,90,ifelse(Observed_hour_angle_displacement>Expected_hour_angle_displacement,270,NA)), type = "angles", units = "degrees", template = "geographics")# EW correction, NA marks Obsereved-Expected equivalence, meaning no correction needed
    
    t<-mean.circular(as.circular(na.omit(c(NS_correction,EW_correction)),  type = "angles", units = "degrees", template = "geographics"))#na.omit remove NAs from the mean calculation - because they signify ne correction needed on the relevant axis

  return(ifelse(all(is.na(t)),"Home",ifelse( t>0,t,360+t))) # If both axis get NA, it means you are home, otherwise, it correctes for the negative output automatically produced by mean.circular by adding 360  
  
}


#Full Binary function - with time shift
########################################
#Time shift is integrated in the observed hour angle

N_Binary_Shifted <- function(NCP_altitude_displacement,NCP_altitude_home,Time_from_sunrise_at_displacement,Day_length,Observed_hour_angle_displacement_no_shift, Time_shift) {
  
    NS_correction<-as.circular(ifelse(NCP_altitude_displacement>NCP_altitude_home,180,ifelse(NCP_altitude_displacement<NCP_altitude_home,0,NA)) , type = "angles", units = "degrees", template = "geographics")#NS correction, NA marks NCP equivalence, meaning no correction needee
    
    Expected_hour_angle_displacement<-as.circular(Time_from_sunrise_at_displacement*180/Day_length, type = "angles", units ="degrees", template = "geographics")#Expected hor angle calculation
    
    Observed_hour_angle_displacement<-Observed_hour_angle_displacement_no_shift+Time_shift*180/Day_length#There might be a better way to do this becasue the observed that the function inputs is meant to contain the time shift
    
    EW_correction<-as.circular(ifelse(Observed_hour_angle_displacement<Expected_hour_angle_displacement,90,ifelse(Observed_hour_angle_displacement>Expected_hour_angle_displacement,270,NA)), type = "angles", units = "degrees", template = "geographics")# EW correction, NA marks Obsereved-Expected equivalence, meaning no correction needed
    
    t<-mean.circular(as.circular(na.omit(c(NS_correction,EW_correction)),  type = "angles", units = "degrees", template = "geographics"))#na.omit remove NAs from the mean calculation - because they signify ne correction needed on the relevant axis

  return(ifelse(all(is.na(t)),"Home",ifelse( t>0,t,360+t))) # If both axis get NA, it means you are home, otherwise, it correctes for the negative output automatically produced by mean.circular by adding 360  
  
}





#Separate functions for North-South and East-West components
############################################################
#North-South component
######################
#Time shift is irrelevant for this compoonent

N_Binary_NS <- function(NCP_altitude_displacement,NCP_altitude_home) {
  
  NS_correction<-as.circular(ifelse(NCP_altitude_displacement>NCP_altitude_home,180,ifelse(NCP_altitude_displacement<NCP_altitude_home,0,NA)) , type = "angles", units = "degrees", template = "geographics")#NS correction, NA marks NCP equivalence, meaning no correction needee
  
  return(NS_correction)
  
}


#East-West component
####################

#No time shift
###############

N_Binary_EW_No_Shift <- function(Time_from_sunrise_at_displacement,Day_length,Observed_hour_angle_displacement) {
  
    Expected_hour_angle_displacement<-as.circular(Time_from_sunrise_at_displacement*180/Day_length)

        EW_correction<-as.circular(ifelse(Observed_hour_angle_displacement<Expected_hour_angle_displacement,90,ifelse(Observed_hour_angle_displacement>Expected_hour_angle_displacement,270,NA)), type = "angles", units = "degrees", template = "geographics")# EW correction, NA marks Obsereved-Expected equivalence, meaning no correction needed
    
  return(EW_correction)
  
}



#With time shift
###############

N_Binary_EW_Shifted <- function(Time_from_sunrise_at_displacement,Day_length,Observed_hour_angle_displacement_no_shift, Time_shift) {
  
    Expected_hour_angle_displacement<-as.circular(Time_from_sunrise_at_displacement*180/Day_length)
    
        Observed_hour_angle_displacement<-Observed_hour_angle_displacement_no_shift+Time_shift*180/Day_length#There might be a better way to do this becasue the observed that the function inputs is meant to contain the time shift
    

        EW_correction<-as.circular(ifelse(Observed_hour_angle_displacement<Expected_hour_angle_displacement,90,ifelse(Observed_hour_angle_displacement>Expected_hour_angle_displacement,270,NA)), type = "angles", units = "degrees", template = "geographics")# EW correction, NA marks Obsereved-Expected equivalence, meaning no correction needed
    
  return(EW_correction)
  
}



###############################################################
###############################################################
#VECTOR HOMING FUNCTIONS
###############################################################
###############################################################

#No time shift 
#############

#String function
###############

#Returns string: "Distance ------, Azimuth-------" For use in testing and isolated cases

N_Vector_String_No_Shift <- function(NCP_altitude_displacement,NCP_altitude_home,Time_from_sunrise_at_displacement,Day_length,Observed_hour_angle_displacement) {
  
NS_distance<-as.numeric(abs(NCP_altitude_displacement-NCP_altitude_home)*111) #In km - this is the average width of a latitudinal degree
##The is the north-south component of the home vector

NS_azimuth<-ifelse(NS_distance==0,NA, as.circular(ifelse(NCP_altitude_displacement>NCP_altitude_home,180,0), type = "angles", units ="degrees", template = "geographics"))
#This is the azimuth of the component - due north (0), or due south (180)

Expected_hour_angle_displacement<-as.circular(Time_from_sunrise_at_displacement*180/Day_length, type = "angles", units ="degrees", template = "geographics")#Expected hor angle calculation

Hour_angle_diff<-as.circular(Observed_hour_angle_displacement - Expected_hour_angle_displacement, type = "angles", units ="degrees", template = "geographics")
#This angle (assuming there is no time shift) represents displacement distance in degrees on the EW axis

EW_distance<-as.numeric(abs(Hour_angle_diff*111.32*as.numeric(cos(NCP_altitude_displacement*pi/180))))#Output in Km, NCP_altitude_displacement converted to radians (pi/180)
#This is the size of the EW component of the home vector

EW_azimuth<-ifelse(EW_distance==0,NA, as.circular(ifelse(Observed_hour_angle_displacement>Expected_hour_angle_displacement,270,90), type = "angles", units ="degrees", template = "geographics"))
#This is the azimuth of this component-due east (90) or due west (270)

#This doesn't need a condition to account for one or both of the compnents being zero becasue the formula hold - it will give the length of the non-zero
#compnent, or zero if both are zeros
Distance_to_home<-as.numeric(sqrt(EW_distance^2+NS_distance^2))
#The straight line distance to home, obtain using Pithagoras as the Hypotenuse of a right angled triangle with EW distancee
#and NS distance as the two other sides

##Convert chord distance (straight line) to arc distance (real distance across a sphere with the earths' average diameter)
##Calculate the arc angle
##I added the condition here because asin can only take [-1,1], and acsc can take anything but that range. The condition basically says if distance to home is larger than the denominator in the brackets of asin, meaning it is larger then 1, (negative values are impossible), use acsc
Sigma<-ifelse(Distance_to_home<12742,2*asin(Distance_to_home/(2*6371)), 2*acsc(Distance_to_home/(2*6371)))##6371 is the earth's reference radius
##Calculate the arc distance into Distance_to_home
Distance_to_home<-Sigma*6371

##Alpha is the angle between the EW component and the homing vector
alpha<-conversion.circular(as.circular(atan(NS_distance/EW_distance), type = "angles", units = "radians", template = "geographics"), type = "angles", units = "degrees", template = "geographics")#Convert to angle separately
#the angle between the NS distance (which is pointing at azimuth 0) and the Hypotenuse 
#To get the Homing azimuth we need a condition set because the calculation varies for the four cases (north+west,north+east,south+west,south+east)

##This conditions outputs the correct azimuth based on alpha and the direction of the North South and East West component. Because azimuth is relative
## To north, to get the right azimuth we need a different caclulation for each of the for component direction combinations
Homing_azimuth<-ifelse (NS_distance == 0 & EW_distance == 0,"Home",
                 ifelse(NS_distance == 0, EW_azimuth , 
                  ifelse (EW_distance == 0 , NS_azimuth, 
                   as.circular(ifelse(EW_azimuth==90&NS_azimuth==0,90-alpha,
                    ifelse(EW_azimuth==90&NS_azimuth==180,90+alpha,
                      ifelse(EW_azimuth==270&NS_azimuth==0,270+alpha,270-alpha))), type = "angles", units = "degrees", template = "geographics"))))
#is the desired flight azimuth to get home on the shortest linear distance "Distance_to_home"


  return(paste("Distance",Distance_to_home,"azimuth",Homing_azimuth)) # If both axis get NA, it means you are home, otherwise, it correctes for the negative output automatically produced by mean.circular by adding 360  
  
}

N_Vector_String_No_Shift_atan2 <- function(NCP_altitude_displacement,NCP_altitude_home,Time_from_sunrise_at_displacement,Day_length,Observed_hour_angle_displacement) {
NS_distance<-as.numeric(NCP_altitude_home-NCP_altitude_displacement)*111 #In km - this is the average width of a latitudinal degree
##The is the north-south component of the home vector

NS_azimuth<-ifelse(NS_distance==0,NA, as.circular(ifelse(NCP_altitude_displacement>NCP_altitude_home,180,0), type = "angles", units ="degrees", template = "geographics"))
#This is the azimuth of the component - due north (0), or due south (180)

Expected_hour_angle_displacement<-as.circular(Time_from_sunrise_at_displacement*180/Day_length, type = "angles", units ="degrees", template = "geographics")#Expected hor angle calculation

Hour_angle_diff<-as.circular(  Expected_hour_angle_displacement-Observed_hour_angle_displacement, type = "angles", units ="degrees", template = "geographics")
#This angle (assuming there is no time shift) represents displacement distance in degrees on the EW axis

EW_distance<-as.numeric(Hour_angle_diff*111.32*as.numeric(cos(NCP_altitude_displacement*pi/180)))#Output in Km, NCP_altitude_displacement converted to radians (pi/180)
#This is the size of the EW component of the home vector

EW_azimuth<-ifelse(EW_distance==0,NA, as.circular(ifelse(Observed_hour_angle_displacement>Expected_hour_angle_displacement,270,90), type = "angles", units ="degrees", template = "geographics"))
#This is the azimuth of this component-due east (90) or due west (270)

#This doesn't need a condition to account for one or both of the compnents being zero becasue the formula hold - it will give the length of the non-zero
#compnent, or zero if both are zeros
Distance_to_home<-as.numeric(sqrt(EW_distance^2+NS_distance^2))
#The straight line distance to home, obtain using Pithagoras as the Hypotenuse of a right angled triangle with EW distancee
#and NS distance as the two other sides

##Convert chord distance (straight line) to arc distance (real distance across a sphere with the earths' average diameter)
##Calculate the arc angle
##I added the condition here because asin can only take [-1,1], and acsc can take anything but that range. The condition basically says if distance to home is larger than the denominator in the brackets of asin, meaning it is larger then 1, (negative values are impossible), use acsc
Sigma<-ifelse(Distance_to_home<12742,2*asin(Distance_to_home/(2*6371)), 2*acsc(Distance_to_home/(2*6371)))##6371 is the earth's reference radius
##Calculate the arc distance into Distance_to_home
Distance_to_home<-Sigma*6371

##Alpha is the angle between the EW component and the homing vector
alpha<-atan2(EW_distance,NS_distance)#Convert to angle separately
#the angle between the NS distance (which is pointing at azimuth 0) and the Hypotenuse 
#To get the Homing azimuth we need a condition set because the calculation varies for the four cases (north+west,north+east,south+west,south+east)

##This conditions outputs the correct azimuth based on alpha and the direction of the North South and East West component. Because azimuth is relative
## To north, to get the right azimuth we need a different caclulation for each of the for component direction combinations
#Homing_azimuth<-(alpha*180/pi)+180
Homing_azimuth<-((alpha*180/pi) %% 360 + 360)%% 360


  return(paste("Distance",Distance_to_home,"azimuth",Homing_azimuth)) # If both axis get NA, it means you are home, otherwise, it correctes for the negative output automatically produced by mean.circular by adding 360  
  
}



#No time shift - Separate function for distance to home and home azimuth for use is analysis 
###########################################################################################

#Straight line distance to home function
#######################################

N_Vector_distance_No_Shift <- function(NCP_altitude_displacement,NCP_altitude_home,Time_from_sunrise_at_displacement,Day_length,Observed_hour_angle_displacement) {
  
NS_distance<-as.numeric(abs(NCP_altitude_displacement-NCP_altitude_home)*111) #In km - this is the average width of a latitudinal degree
##The is the north-south component of the home vector


Expected_hour_angle_displacement<-as.circular(Time_from_sunrise_at_displacement*180/Day_length, type = "angles", units ="degrees", template = "geographics")#Expected hor angle calculation

Hour_angle_diff<-as.circular(Observed_hour_angle_displacement - Expected_hour_angle_displacement, type = "angles", units ="degrees", template = "geographics")
#This angle (assuming there is no time shift) represents displacement distance in degrees on the EW axis

EW_distance<-as.numeric(abs(Hour_angle_diff*111.32*as.numeric(cos(NCP_altitude_displacement*pi/180))))#Output in Km, NCP_altitude_displacement converted to radians (pi/180)
#This is the size of the EW component of the home vector

#This doesn't need a condition to account for one or both of the compnents being zero becasue the formula hold - it will give the length of the non-zero
#compnent, or zero if both are zeros
Distance_to_home<-as.numeric(sqrt(EW_distance^2+NS_distance^2))
#The straight line distance to home, obtain using Pithagoras as the Hypotenuse of a right angled triangle with EW distancee
#and NS distance as the two other sides


##Convert chord distance (straight line) to arc distance (real distance across a sphere with the earths' average diameter)
##Calculate the arc angle
##I added the condition here because asin can only take [-1,1], and acsc can take anything but that range. The condition basically says if distance to home is larger than the denominator in the brackets of asin, meaning it is larger then 1, (negative values are impossible), use acsc
Sigma<-ifelse(Distance_to_home<12742,2*asin(Distance_to_home/(2*6371)), 2*acsc(Distance_to_home/(2*6371)))##6371 is the earth's reference radius
##Calculate the arc distance into Distance_to_home
Distance_to_home<-Sigma*6371


  return(Distance_to_home)# If both axis get NA, it means you are home, otherwise, it correctes for the negative output automatically produced by mean.circular by adding 360  
  
}


#Straight line azimuth home function
###################################

N_Vector_Azimuth_No_Shift <- function(NCP_altitude_displacement,NCP_altitude_home,Time_from_sunrise_at_displacement,Day_length,Observed_hour_angle_displacement) {
  
NS_distance<-as.numeric(abs(NCP_altitude_displacement-NCP_altitude_home)*111) #In km - this is the average width of a latitudinal degree
##The is the north-south component of the home vector

NS_azimuth<-ifelse(NS_distance==0,NA, as.circular(ifelse(NCP_altitude_displacement>NCP_altitude_home,180,0), type = "angles", units ="degrees", template = "geographics"))
#This is the azimuth of the component - due north (0), or due south (180)

Expected_hour_angle_displacement<-as.circular(Time_from_sunrise_at_displacement*180/Day_length, type = "angles", units ="degrees", template = "geographics")#Expected hor angle calculation

Hour_angle_diff<-as.circular(Observed_hour_angle_displacement - Expected_hour_angle_displacement, type = "angles", units ="degrees", template = "geographics")
#This angle (assuming there is no time shift) represents displacement distance in degrees on the EW axis

EW_distance<-as.numeric(abs(Hour_angle_diff*111.32*as.numeric(cos(NCP_altitude_displacement*pi/180))))#Output in Km, NCP_altitude_displacement converted to radians (pi/180)
#This is the size of the EW component of the home vector

EW_azimuth<-ifelse(EW_distance==0,NA, as.circular(ifelse(Observed_hour_angle_displacement>Expected_hour_angle_displacement,270,90), type = "angles", units ="degrees", template = "geographics"))
#This is the azimuth of this component-due east (90) or due west (270)

#This doesn't need a condition to account for one or both of the compnents being zero becasue the formula hold - it will give the length of the non-zero
#compnent, or zero if both are zeros
Distance_to_home<-as.numeric(sqrt(EW_distance^2+NS_distance^2))
#The straight line distance to home, obtain using Pithagoras as the Hypotenuse of a right angled triangle with EW distancee
#and NS distance as the two other sides

##Convert chord distance (straight line) to arc distance (real distance across a sphere with the earths' average diameter)
##Calculate the arc angle
##I added the condition here because asin can only take [-1,1], and acsc can take anything but that range. The condition basically says if distance to home is larger than the denominator in the brackets of asin, meaning it is larger then 1, (negative values are impossible), use acsc
Sigma<-ifelse(Distance_to_home<12742,2*asin(Distance_to_home/(2*6371)), 2*acsc(Distance_to_home/(2*6371)))##6371 is the earth's reference radius
##Calculate the arc distance into Distance_to_home
Distance_to_home<-Sigma*6371

##Alpha is the angle between the EW component and the homing vector
alpha<-conversion.circular(as.circular(atan(NS_distance/EW_distance), type = "angles", units = "radians", template = "geographics"), type = "angles", units = "degrees", template = "geographics")#Convert to angle separately
#the angle between the NS distance (which is pointing at azimuth 0) and the Hypotenuse 
#To get the Homing azimuth we need a condition set because the calculation varies for the four cases (north+west,north+east,south+west,south+east)

##This conditions outputs the correct azimuth based on alpha and the direction of the North South and East West component. Because azimuth is relative
## To north, to get the right azimuth we need a different caclulation for each of the for component direction combinations
Homing_azimuth<-ifelse (NS_distance == 0 & EW_distance == 0,"Home",
                 ifelse(NS_distance == 0, EW_azimuth , 
                  ifelse (EW_distance == 0 , NS_azimuth, 
                   as.circular(ifelse(EW_azimuth==90&NS_azimuth==0,90-alpha,
                    ifelse(EW_azimuth==90&NS_azimuth==180,90+alpha,
                      ifelse(EW_azimuth==270&NS_azimuth==0,270+alpha,270-alpha))), type = "angles", units = "degrees", template = "geographics"))))
#is the desired flight azimuth to get home on the shortest linear distance "Distance_to_home"


  return(Homing_azimuth) # If both axis get NA, it means you are home, otherwise, it correctes for the negative output automatically produced by mean.circular by adding 360  
  
}

N_Vector_Azimuth_No_Shift_atan2 <- function(NCP_altitude_displacement,NCP_altitude_home,Time_from_sunrise_at_displacement,Day_length,Observed_hour_angle_displacement) {
  
NS_distance<-as.numeric(NCP_altitude_home-NCP_altitude_displacement)*111 #In km - this is the average width of a latitudinal degree
##The is the north-south component of the home vector

NS_azimuth<-ifelse(NS_distance==0,NA, as.circular(ifelse(NCP_altitude_displacement>NCP_altitude_home,180,0), type = "angles", units ="degrees", template = "geographics"))
#This is the azimuth of the component - due north (0), or due south (180)

Expected_hour_angle_displacement<-as.circular(Time_from_sunrise_at_displacement*180/Day_length, type = "angles", units ="degrees", template = "geographics")#Expected hor angle calculation

Hour_angle_diff<-as.circular(  Expected_hour_angle_displacement-Observed_hour_angle_displacement, type = "angles", units ="degrees", template = "geographics")
#This angle (assuming there is no time shift) represents displacement distance in degrees on the EW axis

EW_distance<-as.numeric(Hour_angle_diff*111.32*as.numeric(cos(NCP_altitude_displacement*pi/180)))#Output in Km, NCP_altitude_displacement converted to radians (pi/180)
#This is the size of the EW component of the home vector

EW_azimuth<-ifelse(EW_distance==0,NA, as.circular(ifelse(Observed_hour_angle_displacement>Expected_hour_angle_displacement,270,90), type = "angles", units ="degrees", template = "geographics"))
#This is the azimuth of this component-due east (90) or due west (270)

#This doesn't need a condition to account for one or both of the compnents being zero becasue the formula hold - it will give the length of the non-zero
#compnent, or zero if both are zeros
Distance_to_home<-as.numeric(sqrt(EW_distance^2+NS_distance^2))
#The straight line distance to home, obtain using Pithagoras as the Hypotenuse of a right angled triangle with EW distancee
#and NS distance as the two other sides

##Convert chord distance (straight line) to arc distance (real distance across a sphere with the earths' average diameter)
##Calculate the arc angle
##I added the condition here because asin can only take [-1,1], and acsc can take anything but that range. The condition basically says if distance to home is larger than the denominator in the brackets of asin, meaning it is larger then 1, (negative values are impossible), use acsc
Sigma<-ifelse(Distance_to_home<12742,2*asin(Distance_to_home/(2*6371)), 2*acsc(Distance_to_home/(2*6371)))##6371 is the earth's reference radius
##Calculate the arc distance into Distance_to_home
Distance_to_home<-Sigma*6371

##Alpha is the angle between the EW component and the homing vector
alpha<-atan2(EW_distance,NS_distance)#Convert to angle separately
#the angle between the NS distance (which is pointing at azimuth 0) and the Hypotenuse 
#To get the Homing azimuth we need a condition set because the calculation varies for the four cases (north+west,north+east,south+west,south+east)

##This conditions outputs the correct azimuth based on alpha and the direction of the North South and East West component. Because azimuth is relative
## To north, to get the right azimuth we need a different caclulation for each of the for component direction combinations
#Homing_azimuth<-(alpha*180/pi)+180
Homing_azimuth<-((alpha*180/pi) %% 360 + 360)%% 360


#is the desired flight azimuth to get home on the shortest linear distance "Distance_to_home"


  return(Homing_azimuth) # If both axis get NA, it means you are home, otherwise, it correctes for the negative output automatically produced by mean.circular by adding 360  
  
}

#With time shift 
#############

#String function
###############

#Returns string: "Distance ------, Azimuth-------" For use in testing and isolated cases


N_Vector_String_Shifted <- function(NCP_altitude_displacement,NCP_altitude_home,Time_from_sunrise_at_displacement,Day_length,Observed_hour_angle_displacement_no_shift, Time_shift) {
  
NS_distance<-as.numeric(abs(NCP_altitude_displacement-NCP_altitude_home)*111) #In km - this is the average width of a latitudinal degree
##The is the north-south component of the home vector

NS_azimuth<-ifelse(NS_distance==0,NA, as.circular(ifelse(NCP_altitude_displacement>NCP_altitude_home,180,0), type = "angles", units ="degrees", template = "geographics"))
#This is the azimuth of the component - due north (0), or due south (180)

Expected_hour_angle_displacement<-as.circular(Time_from_sunrise_at_displacement*180/Day_length, type = "angles", units ="degrees", template = "geographics")#Expected hor angle calculation

Observed_hour_angle_displacement<-Observed_hour_angle_displacement_no_shift+Time_shift*180/Day_length#There might be a better way to do this becasue the observed that the function inputs is meant to contain the time shift
    
Hour_angle_diff<-as.circular(Observed_hour_angle_displacement - Expected_hour_angle_displacement, type = "angles", units ="degrees", template = "geographics")
#This angle (assuming there is no time shift) represents displacement distance in degrees on the EW axis

EW_distance<-as.numeric(abs(Hour_angle_diff*111.32*as.numeric(cos(NCP_altitude_displacement*pi/180))))#Output in Km, NCP_altitude_displacement converted to radians (pi/180)
#This is the size of the EW component of the home vector

EW_azimuth<-ifelse(EW_distance==0,NA, as.circular(ifelse(Observed_hour_angle_displacement>Expected_hour_angle_displacement,270,90), type = "angles", units ="degrees", template = "geographics"))
#This is the azimuth of this component-due east (90) or due west (270)

#This doesn't need a condition to account for one or both of the compnents being zero becasue the formula hold - it will give the length of the non-zero
#compnent, or zero if both are zeros
Distance_to_home<-as.numeric(sqrt(EW_distance^2+NS_distance^2))
#The straight line distance to home, obtain using Pithagoras as the Hypotenuse of a right angled triangle with EW distancee
#and NS distance as the two other sides

##Convert chord distance (straight line) to arc distance (real distance across a sphere with the earths' average diameter)
##Calculate the arc angle
##I added the condition here because asin can only take [-1,1], and acsc can take anything but that range. The condition basically says if distance to home is larger than the denominator in the brackets of asin, meaning it is larger then 1, (negative values are impossible), use acsc
Sigma<-ifelse(Distance_to_home<12742,2*asin(Distance_to_home/(2*6371)), 2*acsc(Distance_to_home/(2*6371)))##6371 is the earth's reference radius
##Calculate the arc distance into Distance_to_home
Distance_to_home<-Sigma*6371

##Alpha is the angle between the EW component and the homing vector
alpha<-conversion.circular(as.circular(atan(NS_distance/EW_distance), type = "angles", units = "radians", template = "geographics"), type = "angles", units = "degrees", template = "geographics")#Convert to angle separately
#the angle between the NS distance (which is pointing at azimuth 0) and the Hypotenuse 
#To get the Homing azimuth we need a condition set because the calculation varies for the four cases (north+west,north+east,south+west,south+east)

##This conditions outputs the correct azimuth based on alpha and the direction of the North South and East West component. Because azimuth is relative
## To north, to get the right azimuth we need a different calculation for each of the for component direction combinations
Homing_azimuth<-ifelse (NS_distance == 0 & EW_distance == 0,"Home",
                 ifelse(NS_distance == 0, EW_azimuth , 
                  ifelse (EW_distance == 0 , NS_azimuth, 
                   as.circular(ifelse(EW_azimuth==90&NS_azimuth==0,90-alpha,
                    ifelse(EW_azimuth==90&NS_azimuth==180,90+alpha,
                      ifelse(EW_azimuth==270&NS_azimuth==0,270+alpha,270-alpha))), type = "angles", units = "degrees", template = "geographics"))))
#is the desired flight azimuth to get home on the shortest linear distance "Distance_to_home"


  return(paste("Distance",Distance_to_home,"azimuth",Homing_azimuth)) # If both axis get NA, it means you are home, otherwise, it correctes for the negative output automatically produced by mean.circular by adding 360  
  
}


N_Vector_String_Shifted_atan2 <- function(NCP_altitude_displacement,NCP_altitude_home,Time_from_sunrise_at_displacement,Day_length,Observed_hour_angle_displacement_no_shift, Time_shift) {
  
NS_distance<-as.numeric(NCP_altitude_displacement-NCP_altitude_home)*111 #In km - this is the average width of a latitudinal degree
##The is the north-south component of the home vector


Expected_hour_angle_displacement<-as.circular(Time_from_sunrise_at_displacement*180/Day_length, type = "angles", units ="degrees", template = "geographics")#Expected hor angle calculation

Observed_hour_angle_displacement<-Observed_hour_angle_displacement_no_shift+Time_shift*180/Day_length#There might be a better way to do this becasue the observed that the function inputs is meant to contain the time shift

Hour_angle_diff<-as.circular(  Expected_hour_angle_displacement-Observed_hour_angle_displacement, type = "angles", units ="degrees", template = "geographics")
#This angle (assuming there is no time shift) represents displacement distance in degrees on the EW axis

EW_distance<-as.numeric(Hour_angle_diff*111.32*as.numeric(cos(NCP_altitude_displacement*pi/180)))#Output in Km, NCP_altitude_displacement converted to radians (pi/180)
#This is the size of the EW component of the home vector

EW_azimuth<-ifelse(EW_distance==0,NA, as.circular(ifelse(Observed_hour_angle_displacement>Expected_hour_angle_displacement,270,90), type = "angles", units ="degrees", template = "geographics"))
#This is the azimuth of this component-due east (90) or due west (270)

#This doesn't need a condition to account for one or both of the compnents being zero becasue the formula hold - it will give the length of the non-zero
#compnent, or zero if both are zeros
Distance_to_home<-as.numeric(sqrt(EW_distance^2+NS_distance^2))
#The straight line distance to home, obtain using Pithagoras as the Hypotenuse of a right angled triangle with EW distancee
#and NS distance as the two other sides

##Convert chord distance (straight line) to arc distance (real distance across a sphere with the earths' average diameter)
##Calculate the arc angle
##I added the condition here because asin can only take [-1,1], and acsc can take anything but that range. The condition basically says if distance to home is larger than the denominator in the brackets of asin, meaning it is larger then 1, (negative values are impossible), use acsc
Sigma<-ifelse(Distance_to_home<12742,2*asin(Distance_to_home/(2*6371)), 2*acsc(Distance_to_home/(2*6371)))##6371 is the earth's reference radius
##Calculate the arc distance into Distance_to_home
Distance_to_home<-Sigma*6371

##Alpha is the angle between the EW component and the homing vector
alpha<-atan2(EW_distance,NS_distance)#Convert to angle separately
#the angle between the NS distance (which is pointing at azimuth 0) and the Hypotenuse 
#To get the Homing azimuth we need a condition set because the calculation varies for the four cases (north+west,north+east,south+west,south+east)

##This conditions outputs the correct azimuth based on alpha and the direction of the North South and East West component. Because azimuth is relative
## To north, to get the right azimuth we need a different caclulation for each of the for component direction combinations
#Homing_azimuth<-(alpha*180/pi)+180
Homing_azimuth<-((alpha*180/pi) %% 360 + 360)%% 360


  return(paste("Distance",Distance_to_home,"azimuth",Homing_azimuth)) # If both axis get NA, it means you are home, otherwise, it correctes for the negative output automatically produced by mean.circular by adding 360  
  
}





#With time shift - Separate function for distance to home and home azimuth for use is analysis 
################################################################################################

#Straight line distance to home function
########################################

N_Vector_Distance_Shifted <- function(NCP_altitude_displacement,NCP_altitude_home,Time_from_sunrise_at_displacement,Day_length,Observed_hour_angle_displacement_no_shift, Time_shift) {
  
NS_distance<-as.numeric(abs(NCP_altitude_displacement-NCP_altitude_home)*111) #In km - this is the average width of a latitudinal degree
##The is the north-south component of the home vector

Observed_hour_angle_displacement<-Observed_hour_angle_displacement_no_shift+Time_shift*180/Day_length#There might be a better way to do this becasue the observed that the function inputs is meant to contain the time shift
    

Expected_hour_angle_displacement<-as.circular(Time_from_sunrise_at_displacement*180/Day_length, type = "angles", units ="degrees", template = "geographics")#Expected hor angle calculation

Hour_angle_diff<-as.circular(Observed_hour_angle_displacement - Expected_hour_angle_displacement, type = "angles", units ="degrees", template = "geographics")
#This angle (assuming there is no time shift) represents displacement distance in degrees on the EW axis

EW_distance<-as.numeric(abs(Hour_angle_diff*111.32*as.numeric(cos(NCP_altitude_displacement*pi/180))))#Output in Km, NCP_altitude_displacement converted to radians (pi/180)
#This is the size of the EW component of the home vector

#This doesn't need a condition to account for one or both of the compnents being zero becasue the formula hold - it will give the length of the non-zero
#compnent, or zero if both are zeros
Distance_to_home<-as.numeric(sqrt(EW_distance^2+NS_distance^2))
#The straight line distance to home, obtain using Pithagoras as the Hypotenuse of a right angled triangle with EW distancee
#and NS distance as the two other sides

##Convert chord distance (straight line) to arc distance (real distance across a sphere with the earths' average diameter)
##Calculate the arc angle
##I added the condition here because asin can only take [-1,1], and acsc can take anything but that range. The condition basically says if distance to home is larger than the denominator in the brackets of asin, meaning it is larger then 1, (negative values are impossible), use acsc
Sigma<-ifelse(Distance_to_home<12742,2*asin(Distance_to_home/(2*6371)), 2*acsc(Distance_to_home/(2*6371)))##6371 is the earth's reference radius
##Calculate the arc distance into Distance_to_home
Distance_to_home<-Sigma*6371


  return(Distance_to_home)# If both axis get NA, it means you are home, otherwise, it correctes for the negative output automatically produced by mean.circular by adding 360  
  
}

#Straight line azimuth home function
####################################

N_Vector_Azimuth_Shifted <- function(NCP_altitude_displacement,NCP_altitude_home,Time_from_sunrise_at_displacement,Day_length,Observed_hour_angle_displacement_no_shift, Time_shift) {
  
NS_distance<-as.numeric(abs(NCP_altitude_displacement-NCP_altitude_home)*111) #In km - this is the average width of a latitudinal degree
##The is the north-south component of the home vector

NS_azimuth<-ifelse(NS_distance==0,NA, as.circular(ifelse(NCP_altitude_displacement>NCP_altitude_home,180,0), type = "angles", units ="degrees", template = "geographics"))
#This is the azimuth of the component - due north (0), or due south (180)

Observed_hour_angle_displacement<-Observed_hour_angle_displacement_no_shift+Time_shift*180/Day_length#There might be a better way to do this becasue the observed that the function inputs is meant to contain the time shift
    

Expected_hour_angle_displacement<-as.circular(Time_from_sunrise_at_displacement*180/Day_length, type = "angles", units ="degrees", template = "geographics")#Expected hor angle calculation

Hour_angle_diff<-as.circular(Observed_hour_angle_displacement - Expected_hour_angle_displacement, type = "angles", units ="degrees", template = "geographics")
#This angle (assuming there is no time shift) represents displacement distance in degrees on the EW axis

EW_distance<-as.numeric(abs(Hour_angle_diff*111.32*as.numeric(cos(NCP_altitude_displacement*pi/180))))#Output in Km, NCP_altitude_displacement converted to radians (pi/180)
#This is the size of the EW component of the home vector

EW_azimuth<-ifelse(EW_distance==0,NA, as.circular(ifelse(Observed_hour_angle_displacement>Expected_hour_angle_displacement,270,90), type = "angles", units ="degrees", template = "geographics"))
#This is the azimuth of this component-due east (90) or due west (270)

#This doesn't need a condition to account for one or both of the compnents being zero becasue the formula hold - it will give the length of the non-zero
#compnent, or zero if both are zeros
Distance_to_home<-as.numeric(sqrt(EW_distance^2+NS_distance^2))
#The straight line distance to home, obtain using Pithagoras as the Hypotenuse of a right angled triangle with EW distancee
#and NS distance as the two other sides

##Convert chord distance (straight line) to arc distance (real distance across a sphere with the earths' average diameter)
##Calculate the arc angle
##I added the condition here because asin can only take [-1,1], and acsc can take anything but that range. The condition basically says if distance to home is larger than the denominator in the brackets of asin, meaning it is larger then 1, (negative values are impossible), use acsc
Sigma<-ifelse(Distance_to_home<12742,2*asin(Distance_to_home/(2*6371)), 2*acsc(Distance_to_home/(2*6371)))##6371 is the earth's reference radius
##Calculate the arc distance into Distance_to_home
Distance_to_home<-Sigma*6371

##Alpha is the angle between the EW component and the homing vector
alpha<-conversion.circular(as.circular(atan(NS_distance/EW_distance), type = "angles", units = "radians", template = "geographics"), type = "angles", units = "degrees", template = "geographics")#Convert to angle separately
#the angle between the NS distance (which is pointing at azimuth 0) and the Hypotenuse 
#To get the Homing azimuth we need a condition set because the calculation varies for the four cases (north+west,north+east,south+west,south+east)
#צריך לשנות את הפונקציה שבמקום התנאי הזה יהיה חיבור רכיבים כמו ביו וי של הרוח
##This conditions outputs the correct azimuth based on alpha and the direction of the North South and East West component. Because azimuth is relative
## To north, to get the right azimuth we need a different caclulation for each of the for component direction combinations
Homing_azimuth<-ifelse (NS_distance == 0 & EW_distance == 0,"Home",
                 ifelse(NS_distance == 0, EW_azimuth , 
                  ifelse (EW_distance == 0 , NS_azimuth, 
                   as.circular(ifelse(EW_azimuth==90&NS_azimuth==0,90-alpha,
                    ifelse(EW_azimuth==90&NS_azimuth==180,90+alpha,
                      ifelse(EW_azimuth==270&NS_azimuth==0,270+alpha,270-alpha))), type = "angles", units = "degrees", template = "geographics"))))
#is the desired flight azimuth to get home on the shortest linear distance "Distance_to_home"


  return(Homing_azimuth) # If both axis get NA, it means you are home, otherwise, it correctes for the negative output automatically produced by mean.circular by adding 360  
  
}


N_Vector_Azimuth_Shifted_atan2 <- function(NCP_altitude_displacement,NCP_altitude_home,Time_from_sunrise_at_displacement,Day_length,Observed_hour_angle_displacement_no_shift, Time_shift) {
  
NS_distance<-as.numeric(NCP_altitude_displacement-NCP_altitude_home)*111 #In km - this is the average width of a latitudinal degree
##The is the north-south component of the home vector


Observed_hour_angle_displacement<-Observed_hour_angle_displacement_no_shift+Time_shift*180/Day_length#There might be a better way to do this becasue the observed that the function inputs is meant to contain the time shift
    

Expected_hour_angle_displacement<-as.circular(Time_from_sunrise_at_displacement*180/Day_length, type = "angles", units ="degrees", template = "geographics")#Expected hor angle calculation

Hour_angle_diff<-as.circular(  Expected_hour_angle_displacement-Observed_hour_angle_displacement, type = "angles", units ="degrees", template = "geographics")
#This angle (assuming there is no time shift) represents displacement distance in degrees on the EW axis

EW_distance<-as.numeric(Hour_angle_diff*111.32*as.numeric(cos(NCP_altitude_displacement*pi/180)))#Output in Km, NCP_altitude_displacement converted to radians (pi/180)
#This is the size of the EW component of the home vector

EW_azimuth<-ifelse(EW_distance==0,NA, as.circular(ifelse(Observed_hour_angle_displacement>Expected_hour_angle_displacement,270,90), type = "angles", units ="degrees", template = "geographics"))
#This is the azimuth of this component-due east (90) or due west (270)

#This doesn't need a condition to account for one or both of the compnents being zero becasue the formula hold - it will give the length of the non-zero
#compnent, or zero if both are zeros
Distance_to_home<-as.numeric(sqrt(EW_distance^2+NS_distance^2))
#The straight line distance to home, obtain using Pithagoras as the Hypotenuse of a right angled triangle with EW distancee
#and NS distance as the two other sides

##Convert chord distance (straight line) to arc distance (real distance across a sphere with the earths' average diameter)
##Calculate the arc angle
##I added the condition here because asin can only take [-1,1], and acsc can take anything but that range. The condition basically says if distance to home is larger than the denominator in the brackets of asin, meaning it is larger then 1, (negative values are impossible), use acsc
Sigma<-ifelse(Distance_to_home<12742,2*asin(Distance_to_home/(2*6371)), 2*acsc(Distance_to_home/(2*6371)))##6371 is the earth's reference radius
##Calculate the arc distance into Distance_to_home
Distance_to_home<-Sigma*6371

##Alpha is the angle between the EW component and the homing vector
alpha<-atan2(EW_distance,NS_distance)#Convert to angle separately
#the angle between the NS distance (which is pointing at azimuth 0) and the Hypotenuse 
#To get the Homing azimuth we need a condition set because the calculation varies for the four cases (north+west,north+east,south+west,south+east)

##This conditions outputs the correct azimuth based on alpha and the direction of the North South and East West component. Because azimuth is relative
## To north, to get the right azimuth we need a different caclulation for each of the for component direction combinations
#Homing_azimuth<-(alpha*180/pi)+180
Homing_azimuth<-((alpha*180/pi) %% 360 + 360)%% 360

  return(Homing_azimuth) # If both axis get NA, it means you are home, otherwise, it correctes for the negative output automatically produced by mean.circular by adding 360  
  
}



#Southern Hemisphere
#Full Binary function - No time shift:
########################################


#Function:

S_Binary_No_Shift<- function(SCP_altitude_displacement,SCP_altitude_home,Time_from_sunrise_at_displacement,Day_length,Observed_hour_angle_displacement) {
  
    NS_correction<-as.circular(ifelse(SCP_altitude_displacement>SCP_altitude_home,0,ifelse(SCP_altitude_displacement<SCP_altitude_home,180,NA)) , type = "angles", units = "degrees", template = "geographics")#NS correction, NA marks NCP equivalence, meaning no correction needee ***(Switched 0 & 180)
    
    Expected_hour_angle_displacement<-as.circular(Time_from_sunrise_at_displacement*180/Day_length, type = "angles", units ="degrees", template = "geographics")#Expected hor angle calculation
    
    EW_correction<-as.circular(ifelse(Observed_hour_angle_displacement<Expected_hour_angle_displacement,90,ifelse(Observed_hour_angle_displacement>Expected_hour_angle_displacement,270,NA)), type = "angles", units = "degrees", template = "geographics")# EW correction, NA marks Obsereved-Expected equivalence, meaning no correction needed
    
    t<-mean.circular(as.circular(na.omit(c(NS_correction,EW_correction)),  type = "angles", units = "degrees", template = "geographics"))#na.omit remove NAs from the mean calculation - because they signify ne correction needed on the relevant axis

  return(ifelse(all(is.na(t)),"Home",ifelse( t>0,t,360+t))) # If both axis get NA, it means you are home, otherwise, it correctes for the negative output automatically produced by mean.circular by adding 360  
  
}

#Full Binary function - with time shift
########################################
#Time shift is integrated in the observed hour angle

S_Binary_Shifted <- function(SCP_altitude_displacement,SCP_altitude_home,Time_from_sunrise_at_displacement,Day_length,Observed_hour_angle_displacement_no_shift, Time_shift) {
  
    NS_correction<-as.circular(ifelse(SCP_altitude_displacement>SCP_altitude_home,0,ifelse(SCP_altitude_displacement<SCP_altitude_home,180,NA)) , type = "angles", units = "degrees", template = "geographics")#NS correction, NA marks NCP equivalence, meaning no correction needee ***(Switched 0 & 180)
    
    Expected_hour_angle_displacement<-as.circular(Time_from_sunrise_at_displacement*180/Day_length, type = "angles", units ="degrees", template = "geographics")#Expected hor angle calculation
    
    Observed_hour_angle_displacement<-Observed_hour_angle_displacement_no_shift+Time_shift*180/Day_length#There might be a better way to do this becasue the observed that the function inputs is meant to contain the time shift
    
    EW_correction<-as.circular(ifelse(Observed_hour_angle_displacement<Expected_hour_angle_displacement,90,ifelse(Observed_hour_angle_displacement>Expected_hour_angle_displacement,270,NA)), type = "angles", units = "degrees", template = "geographics")# EW correction, NA marks Obsereved-Expected equivalence, meaning no correction needed
    
    t<-mean.circular(as.circular(na.omit(c(NS_correction,EW_correction)),  type = "angles", units = "degrees", template = "geographics"))#na.omit remove NAs from the mean calculation - because they signify ne correction needed on the relevant axis

  return(ifelse(all(is.na(t)),"Home",ifelse( t>0,t,360+t))) # If both axis get NA, it means you are home, otherwise, it correctes for the negative output automatically produced by mean.circular by adding 360  
  
}

#Separate functions for North-South and East-West components
############################################################
#North-South component
######################
#Time shift is irrelevant for this compoonent

S_Binary_NS <- function(SCP_altitude_displacement,SCP_altitude_home) {
  
  NS_correction<-as.circular(ifelse(SCP_altitude_displacement>SCP_altitude_home,0,ifelse(SCP_altitude_displacement<SCP_altitude_home,180,NA)) , type = "angles", units = "degrees", template = "geographics")#NS correction, NA marks NCP equivalence, meaning no correction needee ***(Switched 0 & 180)
  
  return(NS_correction)
  
}

#East-West component
####################

#No time shift
###############

S_Binary_EW_No_Shift <- function(Time_from_sunrise_at_displacement,Day_length,Observed_hour_angle_displacement) {
  
    Expected_hour_angle_displacement<-as.circular(Time_from_sunrise_at_displacement*180/Day_length)

        EW_correction<-as.circular(ifelse(Observed_hour_angle_displacement<Expected_hour_angle_displacement,90,ifelse(Observed_hour_angle_displacement>Expected_hour_angle_displacement,270,NA)), type = "angles", units = "degrees", template = "geographics")# EW correction, NA marks Obsereved-Expected equivalence, meaning no correction needed
    
  return(EW_correction)
  
}

#With time shift
###############


S_Binary_EW_Shifted <- function(Time_from_sunrise_at_displacement,Day_length,Observed_hour_angle_displacement_no_shift, Time_shift) {
  
    Expected_hour_angle_displacement<-as.circular(Time_from_sunrise_at_displacement*180/Day_length)
    
        Observed_hour_angle_displacement<-Observed_hour_angle_displacement_no_shift+Time_shift*180/Day_length#There might be a better way to do this becasue the observed that the function inputs is meant to contain the time shift
    

        EW_correction<-as.circular(ifelse(Observed_hour_angle_displacement<Expected_hour_angle_displacement,90,ifelse(Observed_hour_angle_displacement>Expected_hour_angle_displacement,270,NA)), type = "angles", units = "degrees", template = "geographics")# EW correction, NA marks Obsereved-Expected equivalence, meaning no correction needed
    
  return(EW_correction)
  
}

###############################################################
###############################################################
#VECTOR HOMING FUNCTIONS
###############################################################
###############################################################

#No time shift 
#############

#String function
###############

#Returns string: "Distance ------, Azimuth-------" For use in testing and isolated cases


S_Vector_String_No_Shift <- function(SCP_altitude_displacement,SCP_altitude_home,Time_from_sunrise_at_displacement,Day_length,Observed_hour_angle_displacement) {
  
NS_distance<-as.numeric(abs(SCP_altitude_displacement-SCP_altitude_home)*111) #In km - this is the average width of a latitudinal degree
##The is the north-south component of the home vector

NS_azimuth<-ifelse(NS_distance==0,NA, as.circular(ifelse(SCP_altitude_displacement>SCP_altitude_home,0,180), type = "angles", units ="degrees", template = "geographics"))
#This is the azimuth of the component - due north (0), or due south (180)

Expected_hour_angle_displacement<-as.circular(Time_from_sunrise_at_displacement*180/Day_length, type = "angles", units ="degrees", template = "geographics")#Expected hor angle calculation

Hour_angle_diff<-as.circular(Observed_hour_angle_displacement - Expected_hour_angle_displacement, type = "angles", units ="degrees", template = "geographics")
#This angle (assuming there is no time shift) represents displacement distance in degrees on the EW axis

EW_distance<-as.numeric(abs(Hour_angle_diff*111.32*as.numeric(cos(SCP_altitude_displacement*pi/180))))#Output in Km, NCP_altitude_displacement converted to radians (pi/180)
#This is the size of the EW component of the home vector

EW_azimuth<-ifelse(EW_distance==0,NA, as.circular(ifelse(Observed_hour_angle_displacement>Expected_hour_angle_displacement,270,90), type = "angles", units ="degrees", template = "geographics"))
#This is the azimuth of this component-due east (90) or due west (270)

#This doesn't need a condition to account for one or both of the compnents being zero becasue the formula hold - it will give the length of the non-zero
#compnent, or zero if both are zeros
Distance_to_home<-as.numeric(sqrt(EW_distance^2+NS_distance^2))
#The straight line distance to home, obtain using Pithagoras as the Hypotenuse of a right angled triangle with EW distancee
#and NS distance as the two other sides

##Convert chord distance (straight line) to arc distance (real distance across a sphere with the earths' average diameter)
##Calculate the arc angle
##I added the condition here because asin can only take [-1,1], and acsc can take anything but that range. The condition basically says if distance to home is larger than the denominator in the brackets of asin, meaning it is larger then 1, (negative values are impossible), use acsc
Sigma<-ifelse(Distance_to_home<12742,2*asin(Distance_to_home/(2*6371)), 2*acsc(Distance_to_home/(2*6371)))##6371 is the earth's reference radius
##Calculate the arc distance into Distance_to_home
Distance_to_home<-Sigma*6371

##Alpha is the angle between the EW component and the homing vector
alpha<-conversion.circular(as.circular(atan(NS_distance/EW_distance), type = "angles", units = "radians", template = "geographics"), type = "angles", units = "degrees", template = "geographics")#Convert to angle separately
#the angle between the NS distance (which is pointing at azimuth 0) and the Hypotenuse 
#To get the Homing azimuth we need a condition set because the calculation varies for the four cases (north+west,north+east,south+west,south+east)

##This conditions outputs the correct azimuth based on alpha and the direction of the North South and East West component. Because azimuth is relative
## To north, to get the right azimuth we need a different caclulation for each of the for component direction combinations
Homing_azimuth<-ifelse (NS_distance == 0 & EW_distance == 0,"Home",
                 ifelse(NS_distance == 0, EW_azimuth , 
                  ifelse (EW_distance == 0 , NS_azimuth, 
                   as.circular(ifelse(EW_azimuth==90&NS_azimuth==0,90-alpha,
                    ifelse(EW_azimuth==90&NS_azimuth==180,90+alpha,
                      ifelse(EW_azimuth==270&NS_azimuth==0,270+alpha,270-alpha))), type = "angles", units = "degrees", template = "geographics"))))
#is the desired flight azimuth to get home on the shortest linear distance "Distance_to_home"


  return(paste("Distance",Distance_to_home,"azimuth",Homing_azimuth)) # If both axis get NA, it means you are home, otherwise, it correctes for the negative output automatically produced by mean.circular by adding 360  
  
}




S_Vector_String_No_Shift_atan2 <- function(SCP_altitude_displacement,SCP_altitude_home,Time_from_sunrise_at_displacement,Day_length,Observed_hour_angle_displacement) {
  
NS_distance<-as.numeric(SCP_altitude_displacement-SCP_altitude_home)*111 #In km - this is the average width of a latitudinal degree
##The is the north-south component of the home vector

NS_azimuth<-ifelse(NS_distance==0,NA, as.circular(ifelse(SCP_altitude_displacement>SCP_altitude_home,0,180), type = "angles", units ="degrees", template = "geographics"))
#This is the azimuth of the component - due north (0), or due south (180)

Expected_hour_angle_displacement<-as.circular(Time_from_sunrise_at_displacement*180/Day_length, type = "angles", units ="degrees", template = "geographics")#Expected hor angle calculation

Hour_angle_diff<-as.circular(  Expected_hour_angle_displacement-Observed_hour_angle_displacement, type = "angles", units ="degrees", template = "geographics")
#This angle (assuming there is no time shift) represents displacement distance in degrees on the EW axis

EW_distance<-as.numeric(Hour_angle_diff*111.32*as.numeric(cos(SCP_altitude_displacement*pi/180)))#Output in Km, NCP_altitude_displacement converted to radians (pi/180)
#This is the size of the EW component of the home vector

EW_azimuth<-ifelse(EW_distance==0,NA, as.circular(ifelse(Observed_hour_angle_displacement>Expected_hour_angle_displacement,270,90), type = "angles", units ="degrees", template = "geographics"))
#This is the azimuth of this component-due east (90) or due west (270)

#This doesn't need a condition to account for one or both of the compnents being zero becasue the formula hold - it will give the length of the non-zero
#compnent, or zero if both are zeros
Distance_to_home<-as.numeric(sqrt(EW_distance^2+NS_distance^2))
#The straight line distance to home, obtain using Pithagoras as the Hypotenuse of a right angled triangle with EW distancee
#and NS distance as the two other sides

##Convert chord distance (straight line) to arc distance (real distance across a sphere with the earths' average diameter)
##Calculate the arc angle
##I added the condition here because asin can only take [-1,1], and acsc can take anything but that range. The condition basically says if distance to home is larger than the denominator in the brackets of asin, meaning it is larger then 1, (negative values are impossible), use acsc
Sigma<-ifelse(Distance_to_home<12742,2*asin(Distance_to_home/(2*6371)), 2*acsc(Distance_to_home/(2*6371)))##6371 is the earth's reference radius
##Calculate the arc distance into Distance_to_home
Distance_to_home<-Sigma*6371

##Alpha is the angle between the EW component and the homing vector
alpha<-atan2(EW_distance,NS_distance)#Convert to angle separately
#the angle between the NS distance (which is pointing at azimuth 0) and the Hypotenuse 
#To get the Homing azimuth we need a condition set because the calculation varies for the four cases (north+west,north+east,south+west,south+east)

##This conditions outputs the correct azimuth based on alpha and the direction of the North South and East West component. Because azimuth is relative
## To north, to get the right azimuth we need a different caclulation for each of the for component direction combinations
#Homing_azimuth<-(alpha*180/pi)+180
Homing_azimuth<-((alpha*180/pi) %% 360 + 360)%% 360



  return(paste("Distance",Distance_to_home,"azimuth",Homing_azimuth)) # If both axis get NA, it means you are home, otherwise, it correctes for the negative output automatically produced by mean.circular by adding 360  
  
}


#No time shift - Separate function for distance to home and home azimuth for use is analysis 
###########################################################################################

#Straight line distance to home function
#######################################


S_Vector_distance_No_Shift <- function(SCP_altitude_displacement,SCP_altitude_home,Time_from_sunrise_at_displacement,Day_length,Observed_hour_angle_displacement) {
  
NS_distance<-as.numeric(abs(SCP_altitude_displacement-SCP_altitude_home)*111) #In km - this is the average width of a latitudinal degree
##The is the north-south component of the home vector


Expected_hour_angle_displacement<-as.circular(Time_from_sunrise_at_displacement*180/Day_length, type = "angles", units ="degrees", template = "geographics")#Expected hor angle calculation

Hour_angle_diff<-as.circular(Observed_hour_angle_displacement - Expected_hour_angle_displacement, type = "angles", units ="degrees", template = "geographics")
#This angle (assuming there is no time shift) represents displacement distance in degrees on the EW axis

EW_distance<-as.numeric(abs(Hour_angle_diff*111.32*as.numeric(cos(SCP_altitude_displacement*pi/180))))#Output in Km, NCP_altitude_displacement converted to radians (pi/180)
#This is the size of the EW component of the home vector

#This doesn't need a condition to account for one or both of the compnents being zero becasue the formula hold - it will give the length of the non-zero
#compnent, or zero if both are zeros
Distance_to_home<-as.numeric(sqrt(EW_distance^2+NS_distance^2))
#The straight line distance to home, obtain using Pithagoras as the Hypotenuse of a right angled triangle with EW distancee
#and NS distance as the two other sides

##Convert chord distance (straight line) to arc distance (real distance across a sphere with the earths' average diameter)
##Calculate the arc angle
##I added the condition here because asin can only take [-1,1], and acsc can take anything but that range. The condition basically says if distance to home is larger than the denominator in the brackets of asin, meaning it is larger then 1, (negative values are impossible), use acsc
Sigma<-ifelse(Distance_to_home<12742,2*asin(Distance_to_home/(2*6371)), 2*acsc(Distance_to_home/(2*6371)))##6371 is the earth's reference radius
##Calculate the arc distance into Distance_to_home
Distance_to_home<-Sigma*6371


  return(Distance_to_home)# If both axis get NA, it means you are home, otherwise, it correctes for the negative output automatically produced by mean.circular by adding 360  
  
}




#Straight line azimuth home function
###################################

S_Vector_Azimuth_No_Shift <- function(SCP_altitude_displacement,SCP_altitude_home,Time_from_sunrise_at_displacement,Day_length,Observed_hour_angle_displacement) {
  
NS_distance<-as.numeric(abs(SCP_altitude_displacement-SCP_altitude_home)*111) #In km - this is the average width of a latitudinal degree
##The is the north-south component of the home vector

NS_azimuth<-ifelse(NS_distance==0,NA, as.circular(ifelse(SCP_altitude_displacement>SCP_altitude_home,0,180), type = "angles", units ="degrees", template = "geographics"))
#This is the azimuth of the component - due north (0), or due south (180)

Expected_hour_angle_displacement<-as.circular(Time_from_sunrise_at_displacement*180/Day_length, type = "angles", units ="degrees", template = "geographics")#Expected hor angle calculation

Hour_angle_diff<-as.circular(Observed_hour_angle_displacement - Expected_hour_angle_displacement, type = "angles", units ="degrees", template = "geographics")
#This angle (assuming there is no time shift) represents displacement distance in degrees on the EW axis

EW_distance<-as.numeric(abs(Hour_angle_diff*111.32*as.numeric(cos(SCP_altitude_displacement*pi/180))))#Output in Km, NCP_altitude_displacement converted to radians (pi/180)
#This is the size of the EW component of the home vector

EW_azimuth<-ifelse(EW_distance==0,NA, as.circular(ifelse(Observed_hour_angle_displacement>Expected_hour_angle_displacement,270,90), type = "angles", units ="degrees", template = "geographics"))
#This is the azimuth of this component-due east (90) or due west (270)

#This doesn't need a condition to account for one or both of the compnents being zero becasue the formula hold - it will give the length of the non-zero
#compnent, or zero if both are zeros
Distance_to_home<-as.numeric(sqrt(EW_distance^2+NS_distance^2))
#The straight line distance to home, obtain using Pithagoras as the Hypotenuse of a right angled triangle with EW distancee
#and NS distance as the two other sides

##Convert chord distance (straight line) to arc distance (real distance across a sphere with the earths' average diameter)
##Calculate the arc angle
##I added the condition here because asin can only take [-1,1], and acsc can take anything but that range. The condition basically says if distance to home is larger than the denominator in the brackets of asin, meaning it is larger then 1, (negative values are impossible), use acsc
Sigma<-ifelse(Distance_to_home<12742,2*asin(Distance_to_home/(2*6371)), 2*acsc(Distance_to_home/(2*6371)))##6371 is the earth's reference radius
##Calculate the arc distance into Distance_to_home
Distance_to_home<-Sigma*6371

##Alpha is the angle between the EW component and the homing vector
alpha<-conversion.circular(as.circular(atan(NS_distance/EW_distance), type = "angles", units = "radians", template = "geographics"), type = "angles", units = "degrees", template = "geographics")#Convert to angle separately
#the angle between the NS distance (which is pointing at azimuth 0) and the Hypotenuse 
#To get the Homing azimuth we need a condition set because the calculation varies for the four cases (north+west,north+east,south+west,south+east)

##This conditions outputs the correct azimuth based on alpha and the direction of the North South and East West component. Because azimuth is relative
## To north, to get the right azimuth we need a different caclulation for each of the for component direction combinations
Homing_azimuth<-ifelse (NS_distance == 0 & EW_distance == 0,"Home",
                 ifelse(NS_distance == 0, EW_azimuth , 
                  ifelse (EW_distance == 0 , NS_azimuth, 
                   as.circular(ifelse(EW_azimuth==90&NS_azimuth==0,90-alpha,
                    ifelse(EW_azimuth==90&NS_azimuth==180,90+alpha,
                      ifelse(EW_azimuth==270&NS_azimuth==0,270+alpha,270-alpha))), type = "angles", units = "degrees", template = "geographics"))))
#is the desired flight azimuth to get home on the shortest linear distance "Distance_to_home"


  return(Homing_azimuth) # If both axis get NA, it means you are home, otherwise, it correctes for the negative output automatically produced by mean.circular by adding 360  
  
}


S_Vector_Azimuth_No_Shift_atan2 <- function(SCP_altitude_displacement,SCP_altitude_home,Time_from_sunrise_at_displacement,Day_length,Observed_hour_angle_displacement) {
NS_distance<-as.numeric(SCP_altitude_displacement-SCP_altitude_home)*111 #In km - this is the average width of a latitudinal degree
##The is the north-south component of the home vector

NS_azimuth<-ifelse(NS_distance==0,NA, as.circular(ifelse(SCP_altitude_displacement>SCP_altitude_home,0,180), type = "angles", units ="degrees", template = "geographics"))
#This is the azimuth of the component - due north (0), or due south (180)

Expected_hour_angle_displacement<-as.circular(Time_from_sunrise_at_displacement*180/Day_length, type = "angles", units ="degrees", template = "geographics")#Expected hor angle calculation

Hour_angle_diff<-as.circular(  Expected_hour_angle_displacement-Observed_hour_angle_displacement, type = "angles", units ="degrees", template = "geographics")
#This angle (assuming there is no time shift) represents displacement distance in degrees on the EW axis

EW_distance<-as.numeric(Hour_angle_diff*111.32*as.numeric(cos(SCP_altitude_displacement*pi/180)))#Output in Km, NCP_altitude_displacement converted to radians (pi/180)
#This is the size of the EW component of the home vector

EW_azimuth<-ifelse(EW_distance==0,NA, as.circular(ifelse(Observed_hour_angle_displacement>Expected_hour_angle_displacement,270,90), type = "angles", units ="degrees", template = "geographics"))
#This is the azimuth of this component-due east (90) or due west (270)

#This doesn't need a condition to account for one or both of the compnents being zero becasue the formula hold - it will give the length of the non-zero
#compnent, or zero if both are zeros
Distance_to_home<-as.numeric(sqrt(EW_distance^2+NS_distance^2))
#The straight line distance to home, obtain using Pithagoras as the Hypotenuse of a right angled triangle with EW distancee
#and NS distance as the two other sides

##Convert chord distance (straight line) to arc distance (real distance across a sphere with the earths' average diameter)
##Calculate the arc angle
##I added the condition here because asin can only take [-1,1], and acsc can take anything but that range. The condition basically says if distance to home is larger than the denominator in the brackets of asin, meaning it is larger then 1, (negative values are impossible), use acsc
Sigma<-ifelse(Distance_to_home<12742,2*asin(Distance_to_home/(2*6371)), 2*acsc(Distance_to_home/(2*6371)))##6371 is the earth's reference radius
##Calculate the arc distance into Distance_to_home
Distance_to_home<-Sigma*6371

##Alpha is the angle between the EW component and the homing vector
alpha<-atan2(EW_distance,NS_distance)#Convert to angle separately
#the angle between the NS distance (which is pointing at azimuth 0) and the Hypotenuse 
#To get the Homing azimuth we need a condition set because the calculation varies for the four cases (north+west,north+east,south+west,south+east)

##This conditions outputs the correct azimuth based on alpha and the direction of the North South and East West component. Because azimuth is relative
## To north, to get the right azimuth we need a different caclulation for each of the for component direction combinations
#Homing_azimuth<-(alpha*180/pi)+180
Homing_azimuth<-((alpha*180/pi) %% 360 + 360)%% 360

  return(Homing_azimuth) # If both axis get NA, it means you are home, otherwise, it correctes for the negative output automatically produced by mean.circular by adding 360  
  
}

#With time shift 
#############


#String function
###############

#Returns string: "Distance ------, Azimuth-------" For use in testing and isolated cases

S_Vector_String_Shifted<- function(SCP_altitude_displacement,SCP_altitude_home,Time_from_sunrise_at_displacement,Day_length,Observed_hour_angle_displacement_no_shift, Time_shift) {
  
NS_distance<-as.numeric(abs(SCP_altitude_displacement-SCP_altitude_home)*111) #In km - this is the average width of a latitudinal degree
##The is the north-south component of the home vector

NS_azimuth<-ifelse(NS_distance==0,NA, as.circular(ifelse(SCP_altitude_displacement>SCP_altitude_home,0,180), type = "angles", units ="degrees", template = "geographics"))
#This is the azimuth of the component - due north (0), or due south (180)

Expected_hour_angle_displacement<-as.circular(Time_from_sunrise_at_displacement*180/Day_length, type = "angles", units ="degrees", template = "geographics")#Expected hor angle calculation

Observed_hour_angle_displacement<-Observed_hour_angle_displacement_no_shift+Time_shift*180/Day_length#There might be a better way to do this becasue the observed that the function inputs is meant to contain the time shift
    
Hour_angle_diff<-as.circular(Observed_hour_angle_displacement - Expected_hour_angle_displacement, type = "angles", units ="degrees", template = "geographics")
#This angle (assuming there is no time shift) represents displacement distance in degrees on the EW axis

EW_distance<-as.numeric(abs(Hour_angle_diff*111.32*as.numeric(cos(SCP_altitude_displacement*pi/180))))#Output in Km, NCP_altitude_displacement converted to radians (pi/180)
#This is the size of the EW component of the home vector

EW_azimuth<-ifelse(EW_distance==0,NA, as.circular(ifelse(Observed_hour_angle_displacement>Expected_hour_angle_displacement,270,90), type = "angles", units ="degrees", template = "geographics"))
#This is the azimuth of this component-due east (90) or due west (270)

#This doesn't need a condition to account for one or both of the compnents being zero becasue the formula hold - it will give the length of the non-zero
#compnent, or zero if both are zeros
Distance_to_home<-as.numeric(sqrt(EW_distance^2+NS_distance^2))
#The straight line distance to home, obtain using Pithagoras as the Hypotenuse of a right angled triangle with EW distancee
#and NS distance as the two other sides

##Convert chord distance (straight line) to arc distance (real distance across a sphere with the earths' average diameter)
##Calculate the arc angle
##I added the condition here because asin can only take [-1,1], and acsc can take anything but that range. The condition basically says if distance to home is larger than the denominator in the brackets of asin, meaning it is larger then 1, (negative values are impossible), use acsc
Sigma<-ifelse(Distance_to_home<12742,2*asin(Distance_to_home/(2*6371)), 2*acsc(Distance_to_home/(2*6371)))##6371 is the earth's reference radius
##Calculate the arc distance into Distance_to_home
Distance_to_home<-Sigma*6371

##Alpha is the angle between the EW component and the homing vector
alpha<-conversion.circular(as.circular(atan(NS_distance/EW_distance), type = "angles", units = "radians", template = "geographics"), type = "angles", units = "degrees", template = "geographics")#Convert to angle separately
#the angle between the NS distance (which is pointing at azimuth 0) and the Hypotenuse 
#To get the Homing azimuth we need a condition set because the calculation varies for the four cases (north+west,north+east,south+west,south+east)

##This conditions outputs the correct azimuth based on alpha and the direction of the North South and East West component. Because azimuth is relative
## To north, to get the right azimuth we need a different calculation for each of the for component direction combinations
Homing_azimuth<-ifelse (NS_distance == 0 & EW_distance == 0,"Home",
                 ifelse(NS_distance == 0, EW_azimuth , 
                  ifelse (EW_distance == 0 , NS_azimuth, 
                   as.circular(ifelse(EW_azimuth==90&NS_azimuth==0,90-alpha,
                    ifelse(EW_azimuth==90&NS_azimuth==180,90+alpha,
                      ifelse(EW_azimuth==270&NS_azimuth==0,270+alpha,270-alpha))), type = "angles", units = "degrees", template = "geographics"))))
#is the desired flight azimuth to get home on the shortest linear distance "Distance_to_home"


  return(paste("Distance",Distance_to_home,"azimuth",Homing_azimuth)) # If both axis get NA, it means you are home, otherwise, it correctes for the negative output automatically produced by mean.circular by adding 360  
  
}


S_Vector_String_Shifted_atan2 <- function(SCP_altitude_displacement,SCP_altitude_home,Time_from_sunrise_at_displacement,Day_length,Observed_hour_angle_displacement_no_shift, Time_shift) {
  
NS_distance<-as.numeric(abs(SCP_altitude_displacement-SCP_altitude_home)*111) #In km - this is the average width of a latitudinal degree
##The is the north-south component of the home vector

NS_azimuth<-ifelse(NS_distance==0,NA, as.circular(ifelse(SCP_altitude_displacement>SCP_altitude_home,0,180), type = "angles", units ="degrees", template = "geographics"))
#This is the azimuth of the component - due north (0), or due south (180)

Expected_hour_angle_displacement<-as.circular(Time_from_sunrise_at_displacement*180/Day_length, type = "angles", units ="degrees", template = "geographics")#Expected hor angle calculation

Observed_hour_angle_displacement<-Observed_hour_angle_displacement_no_shift+Time_shift*180/Day_length#There might be a better way to do this becasue the observed that the function inputs is meant to contain the time shift
    

Hour_angle_diff<-as.circular(  Expected_hour_angle_displacement-Observed_hour_angle_displacement, type = "angles", units ="degrees", template = "geographics")
#This angle (assuming there is no time shift) represents displacement distance in degrees on the EW axis

EW_distance<-as.numeric(Hour_angle_diff*111.32*as.numeric(cos(SCP_altitude_displacement*pi/180)))#Output in Km, NCP_altitude_displacement converted to radians (pi/180)
#This is the size of the EW component of the home vector

EW_azimuth<-ifelse(EW_distance==0,NA, as.circular(ifelse(Observed_hour_angle_displacement>Expected_hour_angle_displacement,270,90), type = "angles", units ="degrees", template = "geographics"))
#This is the azimuth of this component-due east (90) or due west (270)

#This doesn't need a condition to account for one or both of the compnents being zero becasue the formula hold - it will give the length of the non-zero
#compnent, or zero if both are zeros
Distance_to_home<-as.numeric(sqrt(EW_distance^2+NS_distance^2))
#The straight line distance to home, obtain using Pithagoras as the Hypotenuse of a right angled triangle with EW distancee
#and NS distance as the two other sides

##Convert chord distance (straight line) to arc distance (real distance across a sphere with the earths' average diameter)
##Calculate the arc angle
##I added the condition here because asin can only take [-1,1], and acsc can take anything but that range. The condition basically says if distance to home is larger than the denominator in the brackets of asin, meaning it is larger then 1, (negative values are impossible), use acsc
Sigma<-ifelse(Distance_to_home<12742,2*asin(Distance_to_home/(2*6371)), 2*acsc(Distance_to_home/(2*6371)))##6371 is the earth's reference radius
##Calculate the arc distance into Distance_to_home
Distance_to_home<-Sigma*6371

##Alpha is the angle between the EW component and the homing vector
alpha<-atan2(EW_distance,NS_distance)#Convert to angle separately
#the angle between the NS distance (which is pointing at azimuth 0) and the Hypotenuse 
#To get the Homing azimuth we need a condition set because the calculation varies for the four cases (north+west,north+east,south+west,south+east)

##This conditions outputs the correct azimuth based on alpha and the direction of the North South and East West component. Because azimuth is relative
## To north, to get the right azimuth we need a different caclulation for each of the for component direction combinations
#Homing_azimuth<-(alpha*180/pi)+180
Homing_azimuth<-((alpha*180/pi) %% 360 + 360)%% 360

  return(paste("Distance",Distance_to_home,"azimuth",Homing_azimuth)) # If both axis get NA, it means you are home, otherwise, it correctes for the negative output automatically produced by mean.circular by adding 360  
  
}


#With time shift - Separate function for distance to home and home azimuth for use is analysis 
################################################################################################

#Straight line distance to home function
########################################

S_Vector_Distance_Shifted <- function(SCP_altitude_displacement,SCP_altitude_home,Time_from_sunrise_at_displacement,Day_length,Observed_hour_angle_displacement_no_shift, Time_shift) {
  
NS_distance<-as.numeric(abs(SCP_altitude_displacement-SCP_altitude_home)*111) #In km - this is the average width of a latitudinal degree
##The is the north-south component of the home vector

Observed_hour_angle_displacement<-Observed_hour_angle_displacement_no_shift+Time_shift*180/Day_length#There might be a better way to do this becasue the observed that the function inputs is meant to contain the time shift
    

Expected_hour_angle_displacement<-as.circular(Time_from_sunrise_at_displacement*180/Day_length, type = "angles", units ="degrees", template = "geographics")#Expected hor angle calculation

Hour_angle_diff<-as.circular(Observed_hour_angle_displacement - Expected_hour_angle_displacement, type = "angles", units ="degrees", template = "geographics")
#This angle (assuming there is no time shift) represents displacement distance in degrees on the EW axis

EW_distance<-as.numeric(abs(Hour_angle_diff*111.32*as.numeric(cos(SCP_altitude_displacement*pi/180))))#Output in Km, NCP_altitude_displacement converted to radians (pi/180)
#This is the size of the EW component of the home vector

#This doesn't need a condition to account for one or both of the compnents being zero becasue the formula hold - it will give the length of the non-zero
#compnent, or zero if both are zeros
Distance_to_home<-as.numeric(sqrt(EW_distance^2+NS_distance^2))
#The straight line distance to home, obtain using Pithagoras as the Hypotenuse of a right angled triangle with EW distancee
#and NS distance as the two other sides

##Convert chord distance (straight line) to arc distance (real distance across a sphere with the earths' average diameter)
##Calculate the arc angle
##I added the condition here because asin can only take [-1,1], and acsc can take anything but that range. The condition basically says if distance to home is larger than the denominator in the brackets of asin, meaning it is larger then 1, (negative values are impossible), use acsc
Sigma<-ifelse(Distance_to_home<12742,2*asin(Distance_to_home/(2*6371)), 2*acsc(Distance_to_home/(2*6371)))##6371 is the earth's reference radius
##Calculate the arc distance into Distance_to_home
Distance_to_home<-Sigma*6371


  return(Distance_to_home)# If both axis get NA, it means you are home, otherwise, it correctes for the negative output automatically produced by mean.circular by adding 360  
  
}



#Straight line azimuth home function
####################################

S_Vector_Azimuth_Shifted <- function(SCP_altitude_displacement,SCP_altitude_home,Time_from_sunrise_at_displacement,Day_length,Observed_hour_angle_displacement_no_shift, Time_shift) {
  
NS_distance<-as.numeric(abs(SCP_altitude_displacement-SCP_altitude_home)*111) #In km - this is the average width of a latitudinal degree
##The is the north-south component of the home vector

NS_azimuth<-ifelse(NS_distance==0,NA, as.circular(ifelse(SCP_altitude_displacement>SCP_altitude_home,0,180), type = "angles", units ="degrees", template = "geographics"))
#This is the azimuth of the component - due north (0), or due south (180)

Observed_hour_angle_displacement<-Observed_hour_angle_displacement_no_shift+Time_shift*180/Day_length#There might be a better way to do this becasue the observed that the function inputs is meant to contain the time shift
    

Expected_hour_angle_displacement<-as.circular(Time_from_sunrise_at_displacement*180/Day_length, type = "angles", units ="degrees", template = "geographics")#Expected hor angle calculation

Hour_angle_diff<-as.circular(Observed_hour_angle_displacement - Expected_hour_angle_displacement, type = "angles", units ="degrees", template = "geographics")
#This angle (assuming there is no time shift) represents displacement distance in degrees on the EW axis

EW_distance<-as.numeric(abs(Hour_angle_diff*111.32*as.numeric(cos(SCP_altitude_displacement*pi/180))))#Output in Km, NCP_altitude_displacement converted to radians (pi/180)
#This is the size of the EW component of the home vector

EW_azimuth<-ifelse(EW_distance==0,NA, as.circular(ifelse(Observed_hour_angle_displacement>Expected_hour_angle_displacement,270,90), type = "angles", units ="degrees", template = "geographics"))
#This is the azimuth of this component-due east (90) or due west (270)

#This doesn't need a condition to account for one or both of the compnents being zero becasue the formula hold - it will give the length of the non-zero
#compnent, or zero if both are zeros
Distance_to_home<-as.numeric(sqrt(EW_distance^2+NS_distance^2))
#The straight line distance to home, obtain using Pithagoras as the Hypotenuse of a right angled triangle with EW distancee
#and NS distance as the two other sides

##Convert chord distance (straight line) to arc distance (real distance across a sphere with the earths' average diameter)
##Calculate the arc angle
##I added the condition here because asin can only take [-1,1], and acsc can take anything but that range. The condition basically says if distance to home is larger than the denominator in the brackets of asin, meaning it is larger then 1, (negative values are impossible), use acsc
Sigma<-ifelse(Distance_to_home<12742,2*asin(Distance_to_home/(2*6371)), 2*acsc(Distance_to_home/(2*6371)))##6371 is the earth's reference radius
##Calculate the arc distance into Distance_to_home
Distance_to_home<-Sigma*6371

##Alpha is the angle between the EW component and the homing vector
alpha<-conversion.circular(as.circular(atan(NS_distance/EW_distance), type = "angles", units = "radians", template = "geographics"), type = "angles", units = "degrees", template = "geographics")#Convert to angle separately
#the angle between the NS distance (which is pointing at azimuth 0) and the Hypotenuse 
#To get the Homing azimuth we need a condition set because the calculation varies for the four cases (north+west,north+east,south+west,south+east)

##This conditions outputs the correct azimuth based on alpha and the direction of the North South and East West component. Because azimuth is relative
## To north, to get the right azimuth we need a different caclulation for each of the for component direction combinations
Homing_azimuth<-ifelse (NS_distance == 0 & EW_distance == 0,"Home",
                 ifelse(NS_distance == 0, EW_azimuth , 
                  ifelse (EW_distance == 0 , NS_azimuth, 
                   as.circular(ifelse(EW_azimuth==90&NS_azimuth==0,90-alpha,
                    ifelse(EW_azimuth==90&NS_azimuth==180,90+alpha,
                      ifelse(EW_azimuth==270&NS_azimuth==0,270+alpha,270-alpha))), type = "angles", units = "degrees", template = "geographics"))))
#is the desired flight azimuth to get home on the shortest linear distance "Distance_to_home"


  return(Homing_azimuth) # If both axis get NA, it means you are home, otherwise, it correctes for the negative output automatically produced by mean.circular by adding 360  
  
}


S_Vector_Azimuth_Shifted_atan2 <- function(SCP_altitude_displacement,SCP_altitude_home,Time_from_sunrise_at_displacement,Day_length,Observed_hour_angle_displacement_no_shift, Time_shift) {
  
NS_distance<-as.numeric(abs(SCP_altitude_displacement-SCP_altitude_home)*111) #In km - this is the average width of a latitudinal degree
##The is the north-south component of the home vector

NS_azimuth<-ifelse(NS_distance==0,NA, as.circular(ifelse(SCP_altitude_displacement>SCP_altitude_home,0,180), type = "angles", units ="degrees", template = "geographics"))
#This is the azimuth of the component - due north (0), or due south (180)

Observed_hour_angle_displacement<-Observed_hour_angle_displacement_no_shift+Time_shift*180/Day_length#There might be a better way to do this becasue the observed that the function inputs is meant to contain the time shift
    

Expected_hour_angle_displacement<-as.circular(Time_from_sunrise_at_displacement*180/Day_length, type = "angles", units ="degrees", template = "geographics")#Expected hor angle calculation


Hour_angle_diff<-as.circular(  Expected_hour_angle_displacement-Observed_hour_angle_displacement, type = "angles", units ="degrees", template = "geographics")
#This angle (assuming there is no time shift) represents displacement distance in degrees on the EW axis

EW_distance<-as.numeric(Hour_angle_diff*111.32*as.numeric(cos(SCP_altitude_displacement*pi/180)))#Output in Km, NCP_altitude_displacement converted to radians (pi/180)
#This is the size of the EW component of the home vector

EW_azimuth<-ifelse(EW_distance==0,NA, as.circular(ifelse(Observed_hour_angle_displacement>Expected_hour_angle_displacement,270,90), type = "angles", units ="degrees", template = "geographics"))
#This is the azimuth of this component-due east (90) or due west (270)

#This doesn't need a condition to account for one or both of the compnents being zero becasue the formula hold - it will give the length of the non-zero
#compnent, or zero if both are zeros
Distance_to_home<-as.numeric(sqrt(EW_distance^2+NS_distance^2))
#The straight line distance to home, obtain using Pithagoras as the Hypotenuse of a right angled triangle with EW distancee
#and NS distance as the two other sides

##Convert chord distance (straight line) to arc distance (real distance across a sphere with the earths' average diameter)
##Calculate the arc angle
##I added the condition here because asin can only take [-1,1], and acsc can take anything but that range. The condition basically says if distance to home is larger than the denominator in the brackets of asin, meaning it is larger then 1, (negative values are impossible), use acsc
Sigma<-ifelse(Distance_to_home<12742,2*asin(Distance_to_home/(2*6371)), 2*acsc(Distance_to_home/(2*6371)))##6371 is the earth's reference radius
##Calculate the arc distance into Distance_to_home
Distance_to_home<-Sigma*6371

##Alpha is the angle between the EW component and the homing vector
alpha<-atan2(EW_distance,NS_distance)#Convert to angle separately
#the angle between the NS distance (which is pointing at azimuth 0) and the Hypotenuse 
#To get the Homing azimuth we need a condition set because the calculation varies for the four cases (north+west,north+east,south+west,south+east)

##This conditions outputs the correct azimuth based on alpha and the direction of the North South and East West component. Because azimuth is relative
## To north, to get the right azimuth we need a different caclulation for each of the for component direction combinations
#Homing_azimuth<-(alpha*180/pi)+180
Homing_azimuth<-((alpha*180/pi) %% 360 + 360)%% 360
  return(Homing_azimuth) # If both axis get NA, it means you are home, otherwise, it correctes for the negative output automatically produced by mean.circular by adding 360  
  
}

```

Apllying the functions n the same order as they appear in the previous chunck
```{r chunk-name, warning=FALSE, message=FALSE}

###Table100
##Change parameters types

#Table$Latitude_home<-as.circular(Table$Latitude_home, type = "angles", units ="degrees", template = "geographics")
#Table$Latitude_displacement<-as.circular(Table$Latitude_displacement, type = "angles", units ="degrees", template = "geographics")
#Table$Day_length<-as.numeric(Table$Day_length)#in hours
#Table$Time_from_sunrise<-as.numeric(Table$Time_from_sunrise)#in hours
#Table$Observed_hour_angle<-as.circular(Table$Observed_hour_angle, type = "angles", units ="degrees", template = "geographics")
#Table$Observed_hour_angle_no_shift<-as.circular(Table$Observed_hour_angle, type = "angles", units ="degrees", template = "geographics")
#Table$Time_shift<-as.numeric(Table$Time_shift)#in hours

##Sample the data its happenning slow
set.seed(123)
Table2<-Table_100[sample(nrow(Table_100), 100000, replace = FALSE), ]

#Table2<-Table
rm(Table)
##N_Binary_No_Shift####
#######################

N_Binary_No_Shift<-Vectorize(N_Binary_No_Shift)## This steps changes the functions so that they can take vectors (table columns), rather then scalar numbers

Table2$N_Binary_No_Shift<-N_Binary_No_Shift(NCP_altitude_displacement = Table2$Latitude_displacement,
                                          NCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement = Table2$Observed_hour_angle)



##N_Binary_Shifted##
####################
N_Binary_Shifted<-Vectorize(N_Binary_Shifted)

Table2$N_Binary_Shifted<-N_Binary_Shifted(NCP_altitude_displacement = Table2$Latitude_displacement,
                                          NCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement_no_shift = Table2$Observed_hour_angle,
                                          Time_shift = Table2$Time_shift )







#Separate functions for North-South and East-West components
############################################################
#North-South component
######################
#Time shift is irrelevant for this compoonent
N_Binary_NS<-Vectorize(N_Binary_NS )
Table2$N_Binary_NS<-N_Binary_NS(NCP_altitude_displacement = Table2$Latitude_displacement,
                                NCP_altitude_home = Table2$Latitude_home) 

#East-West component
####################
#No time shift
###############
N_Binary_EW_No_Shift<-Vectorize(N_Binary_EW_No_Shift)
Table2$N_Binary_EW_No_Shift<-N_Binary_EW_No_Shift(Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                                  Day_length = Table2$Day_length,
                                                  Observed_hour_angle_displacement = Table2$Observed_hour_angle
                                                  )

#With time shift
###############
N_Binary_EW_Shifted<-Vectorize(N_Binary_EW_Shifted)
Table2$N_Binary_EW_Shifted<-N_Binary_EW_Shifted(Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                                  Day_length = Table2$Day_length,
                                                  Observed_hour_angle_displacement_no_shift =  Table2$Observed_hour_angle,
                                                  Time_shift = Table2$Time_shift
                                                )

###############################################################
###############################################################
#VECTOR HOMING FUNCTIONS
###############################################################
###############################################################

#No time shift 
#############

#String function
###############

#Returns string: "Distance ------, Azimuth-------" For use in testing and isolated cases
N_Vector_String_No_Shift<-Vectorize(N_Vector_String_No_Shift)
Table2$N_Vector_String_No_Shift<-N_Vector_String_No_Shift(NCP_altitude_displacement = Table2$Latitude_displacement,
                                          NCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement = Table2$Observed_hour_angle)

N_Vector_String_No_Shift_atan2<-Vectorize(N_Vector_String_No_Shift_atan2)
Table2$N_Vector_String_No_Shift_atan2<-N_Vector_String_No_Shift_atan2(NCP_altitude_displacement = Table2$Latitude_displacement,
                                          NCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement = Table2$Observed_hour_angle)
#No time shift - Separate function for distance to home and home azimuth for use is analysis 
###########################################################################################

#Straight line distance to home function
#######################################

N_Vector_distance_No_Shift<-Vectorize(N_Vector_distance_No_Shift)
Table2$N_Vector_distance_No_Shift<-N_Vector_distance_No_Shift(NCP_altitude_displacement = Table2$Latitude_displacement,
                                          NCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement = Table2$Observed_hour_angle)
#Straight line azimuth home function
###################################
N_Vector_Azimuth_No_Shift<-Vectorize(N_Vector_Azimuth_No_Shift)
Table2$N_Vector_Azimuth_No_Shift<-N_Vector_Azimuth_No_Shift(NCP_altitude_displacement = Table2$Latitude_displacement,
                                          NCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement = Table2$Observed_hour_angle)


N_Vector_Azimuth_No_Shift_atan2<-Vectorize(N_Vector_Azimuth_No_Shift_atan2)
Table2$N_Vector_Azimuth_No_Shift_atan2<-N_Vector_Azimuth_No_Shift_atan2(NCP_altitude_displacement = Table2$Latitude_displacement,
                                          NCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement = Table2$Observed_hour_angle)


#With time shift 
#############

#String function
###############

#Returns string: "Distance ------, Azimuth-------" For use in testing and isolated cases

N_Vector_String_Shifted<-Vectorize(N_Vector_String_Shifted)
Table2$N_Vector_String_Shifted<-N_Vector_String_Shifted(NCP_altitude_displacement = Table2$Latitude_displacement,
                                          NCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement_no_shift = Table2$Observed_hour_angle,
                                          Time_shift = Table2$Time_shift )

N_Vector_String_Shifted_atan2<-Vectorize(N_Vector_String_Shifted_atan2)
Table2$N_Vector_String_Shifted_atan2<-N_Vector_String_Shifted_atan2(NCP_altitude_displacement = Table2$Latitude_displacement,
                                          NCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement_no_shift = Table2$Observed_hour_angle,
                                          Time_shift = Table2$Time_shift )
#With time shift - Separate function for distance to home and home azimuth for use is analysis 
################################################################################################

#Straight line distance to home function
########################################
N_Vector_Distance_Shifted<-Vectorize(N_Vector_Distance_Shifted)
Table2$N_Vector_Distance_Shifted<-N_Vector_Distance_Shifted(NCP_altitude_displacement = Table2$Latitude_displacement,
                                          NCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement_no_shift = Table2$Observed_hour_angle,
                                          Time_shift = Table2$Time_shift )

#Straight line azimuth home function
####################################
N_Vector_Azimuth_Shifted<-Vectorize(N_Vector_Azimuth_Shifted)
Table2$N_Vector_Azimuth_Shifted<-N_Vector_Azimuth_Shifted(NCP_altitude_displacement = Table2$Latitude_displacement,
                                          NCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement_no_shift = Table2$Observed_hour_angle,
                                          Time_shift = Table2$Time_shift )

N_Vector_Azimuth_Shifted_atan2<-Vectorize(N_Vector_Azimuth_Shifted_atan2)
Table2$N_Vector_Azimuth_Shifted_atan2<-N_Vector_Azimuth_Shifted_atan2(NCP_altitude_displacement = Table2$Latitude_displacement,
                                          NCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement_no_shift = Table2$Observed_hour_angle,
                                          Time_shift = Table2$Time_shift )
####Southern function were not yet adressed - The entire comparison df is for northern based on the geosphere function definition. Run the functions and then
###Figure it out

#Southern Hemisphere
#Full Binary function - No time shift:
########################################

S_Binary_No_Shift<-Vectorize(S_Binary_No_Shift)
Table2$S_Binary_No_Shift<-S_Binary_No_Shift(SCP_altitude_displacement = Table2$Latitude_displacement,
                                          SCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement = Table2$Observed_hour_angle)

#Full Binary function - with time shift
########################################
#Time shift is integrated in the observed hour angle
S_Binary_Shifted<-Vectorize(S_Binary_Shifted)
Table2$S_Binary_Shifted<-S_Binary_Shifted(SCP_altitude_displacement = Table2$Latitude_displacement,
                                          SCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement = Table2$Observed_hour_angle,
                                          Time_shift = Table2$Time_shift)
#Separate functions for North-South and East-West components
############################################################
#North-South component
######################
#Time shift is irrelevant for this compoonent

S_Binary_NS<-Vectorize(S_Binary_NS)
Table2$S_Binary_NS<-S_Binary_NS(SCP_altitude_displacement = Table2$Latitude_displacement,
                                          SCP_altitude_home = Table2$Latitude_home)

#East-West component
####################

#No time shift
###############
S_Binary_EW_No_Shift<-Vectorize(S_Binary_EW_No_Shift)
Table2$S_Binary_EW_No_Shift<-S_Binary_EW_No_Shift(Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement = Table2$Observed_hour_angle)

#With time shift
###############

S_Binary_EW_Shifted<-Vectorize(S_Binary_EW_Shifted)
Table2$S_Binary_EW_Shifted<-S_Binary_EW_Shifted(Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement = Table2$Observed_hour_angle,
                                          Time_shift = Table2$Time_shift)

###############################################################
###############################################################
#VECTOR HOMING FUNCTIONS
###############################################################
###############################################################

#No time shift 
#############

#String function
###############

#Returns string: "Distance ------, Azimuth-------" For use in testing and isolated cases

S_Vector_String_No_Shift<-Vectorize(S_Vector_String_No_Shift)
Table2$S_Vector_String_No_Shift<-S_Vector_String_No_Shift(SCP_altitude_displacement = Table2$Latitude_displacement,
                                          SCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement = Table2$Observed_hour_angle)


S_Vector_String_No_Shift_atan2<-Vectorize(S_Vector_String_No_Shift_atan2)
Table2$S_Vector_String_No_Shift_atan2<-S_Vector_String_No_Shift_atan2(SCP_altitude_displacement = Table2$Latitude_displacement,
                                          SCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement = Table2$Observed_hour_angle)

#No time shift - Separate function for distance to home and home azimuth for use is analysis 
###########################################################################################

#Straight line distance to home function
#######################################
S_Vector_distance_No_Shift<-Vectorize(S_Vector_distance_No_Shift)
Table2$S_Vector_distance_No_Shift<-S_Vector_distance_No_Shift(SCP_altitude_displacement = Table2$Latitude_displacement,
                                          SCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement = Table2$Observed_hour_angle)

#Straight line azimuth home function
###################################
S_Vector_Azimuth_No_Shift<-Vectorize(S_Vector_Azimuth_No_Shift)
Table2$S_Vector_Azimuth_No_Shift<-S_Vector_Azimuth_No_Shift(SCP_altitude_displacement = Table2$Latitude_displacement,
                                          SCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement = Table2$Observed_hour_angle)

S_Vector_Azimuth_No_Shift_atan2<-Vectorize(S_Vector_Azimuth_No_Shift_atan2)
Table2$S_Vector_Azimuth_No_Shift_atan2<-S_Vector_Azimuth_No_Shift_atan2(SCP_altitude_displacement = Table2$Latitude_displacement,
                                          SCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement = Table2$Observed_hour_angle)
#With time shift 
#############

#String function
###############

#Returns string: "Distance ------, Azimuth-------" For use in testing and isolated cases
S_Vector_String_Shifted<-Vectorize(S_Vector_String_Shifted)
Table2$S_Vector_String_Shifted<-S_Vector_String_Shifted(SCP_altitude_displacement = Table2$Latitude_displacement,
                                          SCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement = Table2$Observed_hour_angle,
                                          Time_shift = Table2$Time_shift)

S_Vector_String_Shifted_atan2<-Vectorize(S_Vector_String_Shifted_atan2)
Table2$S_Vector_String_Shifted_atan2<-S_Vector_String_Shifted_atan2(SCP_altitude_displacement = Table2$Latitude_displacement,
                                          SCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement = Table2$Observed_hour_angle,
                                          Time_shift = Table2$Time_shift)

#With time shift - Separate function for distance to home and home azimuth for use is analysis 
################################################################################################

#Straight line distance to home function
########################################
S_Vector_Distance_Shifted<-Vectorize(S_Vector_Distance_Shifted)
Table2$S_Vector_Distance_Shifted<-S_Vector_Distance_Shifted(SCP_altitude_displacement = Table2$Latitude_displacement,
                                          SCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement = Table2$Observed_hour_angle,
                                          Time_shift = Table2$Time_shift)

#Straight line azimuth home function
####################################
S_Vector_Azimuth_Shifted<-Vectorize(S_Vector_Azimuth_Shifted)
Table2$S_Vector_Azimuth_Shifted<-S_Vector_Azimuth_Shifted(SCP_altitude_displacement = Table2$Latitude_displacement,
                                          SCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement = Table2$Observed_hour_angle,
                                          Time_shift = Table2$Time_shift)


S_Vector_Azimuth_Shifted_atan2<-Vectorize(S_Vector_Azimuth_Shifted_atan2)
Table2$S_Vector_Azimuth_Shifted_atan2<-S_Vector_Azimuth_Shifted_atan2(SCP_altitude_displacement = Table2$Latitude_displacement,
                                          SCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement = Table2$Observed_hour_angle,
                                          Time_shift = Table2$Time_shift)

Table_100<-Table2


set.seed(123)

Table2<-Table_1000[sample(nrow(Table_1000), 100000, replace = FALSE), ]


##N_Binary_No_Shift####
#######################

N_Binary_No_Shift<-Vectorize(N_Binary_No_Shift)## This steps changes the functions so that they can take vectors (table columns), rather then scalar numbers

Table2$N_Binary_No_Shift<-N_Binary_No_Shift(NCP_altitude_displacement = Table2$Latitude_displacement,
                                          NCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement = Table2$Observed_hour_angle)



##N_Binary_Shifted##
####################
N_Binary_Shifted<-Vectorize(N_Binary_Shifted)

Table2$N_Binary_Shifted<-N_Binary_Shifted(NCP_altitude_displacement = Table2$Latitude_displacement,
                                          NCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement_no_shift = Table2$Observed_hour_angle,
                                          Time_shift = Table2$Time_shift )







#Separate functions for North-South and East-West components
############################################################
#North-South component
######################
#Time shift is irrelevant for this compoonent
N_Binary_NS<-Vectorize(N_Binary_NS )
Table2$N_Binary_NS<-N_Binary_NS(NCP_altitude_displacement = Table2$Latitude_displacement,
                                NCP_altitude_home = Table2$Latitude_home) 

#East-West component
####################
#No time shift
###############
N_Binary_EW_No_Shift<-Vectorize(N_Binary_EW_No_Shift)
Table2$N_Binary_EW_No_Shift<-N_Binary_EW_No_Shift(Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                                  Day_length = Table2$Day_length,
                                                  Observed_hour_angle_displacement = Table2$Observed_hour_angle
                                                  )

#With time shift
###############
N_Binary_EW_Shifted<-Vectorize(N_Binary_EW_Shifted)
Table2$N_Binary_EW_Shifted<-N_Binary_EW_Shifted(Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                                  Day_length = Table2$Day_length,
                                                  Observed_hour_angle_displacement_no_shift =  Table2$Observed_hour_angle,
                                                  Time_shift = Table2$Time_shift
                                                )

###############################################################
###############################################################
#VECTOR HOMING FUNCTIONS
###############################################################
###############################################################

#No time shift 
#############

#String function
###############

#Returns string: "Distance ------, Azimuth-------" For use in testing and isolated cases
N_Vector_String_No_Shift<-Vectorize(N_Vector_String_No_Shift)
Table2$N_Vector_String_No_Shift<-N_Vector_String_No_Shift(NCP_altitude_displacement = Table2$Latitude_displacement,
                                          NCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement = Table2$Observed_hour_angle)

N_Vector_String_No_Shift_atan2<-Vectorize(N_Vector_String_No_Shift_atan2)
Table2$N_Vector_String_No_Shift_atan2<-N_Vector_String_No_Shift_atan2(NCP_altitude_displacement = Table2$Latitude_displacement,
                                          NCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement = Table2$Observed_hour_angle)
#No time shift - Separate function for distance to home and home azimuth for use is analysis 
###########################################################################################

#Straight line distance to home function
#######################################

N_Vector_distance_No_Shift<-Vectorize(N_Vector_distance_No_Shift)
Table2$N_Vector_distance_No_Shift<-N_Vector_distance_No_Shift(NCP_altitude_displacement = Table2$Latitude_displacement,
                                          NCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement = Table2$Observed_hour_angle)
#Straight line azimuth home function
###################################
N_Vector_Azimuth_No_Shift<-Vectorize(N_Vector_Azimuth_No_Shift)
Table2$N_Vector_Azimuth_No_Shift<-N_Vector_Azimuth_No_Shift(NCP_altitude_displacement = Table2$Latitude_displacement,
                                          NCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement = Table2$Observed_hour_angle)


N_Vector_Azimuth_No_Shift_atan2<-Vectorize(N_Vector_Azimuth_No_Shift_atan2)
Table2$N_Vector_Azimuth_No_Shift_atan2<-N_Vector_Azimuth_No_Shift_atan2(NCP_altitude_displacement = Table2$Latitude_displacement,
                                          NCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement = Table2$Observed_hour_angle)


#With time shift 
#############

#String function
###############

#Returns string: "Distance ------, Azimuth-------" For use in testing and isolated cases

N_Vector_String_Shifted<-Vectorize(N_Vector_String_Shifted)
Table2$N_Vector_String_Shifted<-N_Vector_String_Shifted(NCP_altitude_displacement = Table2$Latitude_displacement,
                                          NCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement_no_shift = Table2$Observed_hour_angle,
                                          Time_shift = Table2$Time_shift )

N_Vector_String_Shifted_atan2<-Vectorize(N_Vector_String_Shifted_atan2)
Table2$N_Vector_String_Shifted_atan2<-N_Vector_String_Shifted_atan2(NCP_altitude_displacement = Table2$Latitude_displacement,
                                          NCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement_no_shift = Table2$Observed_hour_angle,
                                          Time_shift = Table2$Time_shift )
#With time shift - Separate function for distance to home and home azimuth for use is analysis 
################################################################################################

#Straight line distance to home function
########################################
N_Vector_Distance_Shifted<-Vectorize(N_Vector_Distance_Shifted)
Table2$N_Vector_Distance_Shifted<-N_Vector_Distance_Shifted(NCP_altitude_displacement = Table2$Latitude_displacement,
                                          NCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement_no_shift = Table2$Observed_hour_angle,
                                          Time_shift = Table2$Time_shift )

#Straight line azimuth home function
####################################
N_Vector_Azimuth_Shifted<-Vectorize(N_Vector_Azimuth_Shifted)
Table2$N_Vector_Azimuth_Shifted<-N_Vector_Azimuth_Shifted(NCP_altitude_displacement = Table2$Latitude_displacement,
                                          NCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement_no_shift = Table2$Observed_hour_angle,
                                          Time_shift = Table2$Time_shift )

N_Vector_Azimuth_Shifted_atan2<-Vectorize(N_Vector_Azimuth_Shifted_atan2)
Table2$N_Vector_Azimuth_Shifted_atan2<-N_Vector_Azimuth_Shifted_atan2(NCP_altitude_displacement = Table2$Latitude_displacement,
                                          NCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement_no_shift = Table2$Observed_hour_angle,
                                          Time_shift = Table2$Time_shift )
####Southern function were not yet adressed - The entire comparison df is for northern based on the geosphere function definition. Run the functions and then
###Figure it out

#Southern Hemisphere
#Full Binary function - No time shift:
########################################

S_Binary_No_Shift<-Vectorize(S_Binary_No_Shift)
Table2$S_Binary_No_Shift<-S_Binary_No_Shift(SCP_altitude_displacement = Table2$Latitude_displacement,
                                          SCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement = Table2$Observed_hour_angle)

#Full Binary function - with time shift
########################################
#Time shift is integrated in the observed hour angle
S_Binary_Shifted<-Vectorize(S_Binary_Shifted)
Table2$S_Binary_Shifted<-S_Binary_Shifted(SCP_altitude_displacement = Table2$Latitude_displacement,
                                          SCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement = Table2$Observed_hour_angle,
                                          Time_shift = Table2$Time_shift)
#Separate functions for North-South and East-West components
############################################################
#North-South component
######################
#Time shift is irrelevant for this compoonent

S_Binary_NS<-Vectorize(S_Binary_NS)
Table2$S_Binary_NS<-S_Binary_NS(SCP_altitude_displacement = Table2$Latitude_displacement,
                                          SCP_altitude_home = Table2$Latitude_home)

#East-West component
####################

#No time shift
###############
S_Binary_EW_No_Shift<-Vectorize(S_Binary_EW_No_Shift)
Table2$S_Binary_EW_No_Shift<-S_Binary_EW_No_Shift(Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement = Table2$Observed_hour_angle)

#With time shift
###############

S_Binary_EW_Shifted<-Vectorize(S_Binary_EW_Shifted)
Table2$S_Binary_EW_Shifted<-S_Binary_EW_Shifted(Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement = Table2$Observed_hour_angle,
                                          Time_shift = Table2$Time_shift)

###############################################################
###############################################################
#VECTOR HOMING FUNCTIONS
###############################################################
###############################################################

#No time shift 
#############

#String function
###############

#Returns string: "Distance ------, Azimuth-------" For use in testing and isolated cases

S_Vector_String_No_Shift<-Vectorize(S_Vector_String_No_Shift)
Table2$S_Vector_String_No_Shift<-S_Vector_String_No_Shift(SCP_altitude_displacement = Table2$Latitude_displacement,
                                          SCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement = Table2$Observed_hour_angle)

S_Vector_String_No_Shift_atan2<-Vectorize(S_Vector_String_No_Shift_atan2)
Table2$S_Vector_String_No_Shift_atan2<-S_Vector_String_No_Shift_atan2(SCP_altitude_displacement = Table2$Latitude_displacement,
                                          SCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement = Table2$Observed_hour_angle)

#No time shift - Separate function for distance to home and home azimuth for use is analysis 
###########################################################################################

#Straight line distance to home function
#######################################
S_Vector_distance_No_Shift<-Vectorize(S_Vector_distance_No_Shift)
Table2$S_Vector_distance_No_Shift<-S_Vector_distance_No_Shift(SCP_altitude_displacement = Table2$Latitude_displacement,
                                          SCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement = Table2$Observed_hour_angle)

#Straight line azimuth home function
###################################
S_Vector_Azimuth_No_Shift<-Vectorize(S_Vector_Azimuth_No_Shift)
Table2$S_Vector_Azimuth_No_Shift<-S_Vector_Azimuth_No_Shift(SCP_altitude_displacement = Table2$Latitude_displacement,
                                          SCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement = Table2$Observed_hour_angle)

S_Vector_Azimuth_No_Shift_atan2<-Vectorize(S_Vector_Azimuth_No_Shift_atan2)
Table2$S_Vector_Azimuth_No_Shift_atan2<-S_Vector_Azimuth_No_Shift_atan2(SCP_altitude_displacement = Table2$Latitude_displacement,
                                          SCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement = Table2$Observed_hour_angle)
#With time shift 
#############

#String function
###############

#Returns string: "Distance ------, Azimuth-------" For use in testing and isolated cases
S_Vector_String_Shifted<-Vectorize(S_Vector_String_Shifted)
Table2$S_Vector_String_Shifted<-S_Vector_String_Shifted(SCP_altitude_displacement = Table2$Latitude_displacement,
                                          SCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement = Table2$Observed_hour_angle,
                                          Time_shift = Table2$Time_shift)

S_Vector_String_Shifted_atan2<-Vectorize(S_Vector_String_Shifted_atan2)
Table2$S_Vector_String_Shifted_atan2<-S_Vector_String_Shifted_atan2(SCP_altitude_displacement = Table2$Latitude_displacement,
                                          SCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement = Table2$Observed_hour_angle,
                                          Time_shift = Table2$Time_shift)

#With time shift - Separate function for distance to home and home azimuth for use is analysis 
################################################################################################

#Straight line distance to home function
########################################
S_Vector_Distance_Shifted<-Vectorize(S_Vector_Distance_Shifted)
Table2$S_Vector_Distance_Shifted<-S_Vector_Distance_Shifted(SCP_altitude_displacement = Table2$Latitude_displacement,
                                          SCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement = Table2$Observed_hour_angle,
                                          Time_shift = Table2$Time_shift)

#Straight line azimuth home function
####################################
S_Vector_Azimuth_Shifted<-Vectorize(S_Vector_Azimuth_Shifted)
Table2$S_Vector_Azimuth_Shifted<-S_Vector_Azimuth_Shifted(SCP_altitude_displacement = Table2$Latitude_displacement,
                                          SCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement = Table2$Observed_hour_angle,
                                          Time_shift = Table2$Time_shift)

S_Vector_Azimuth_Shifted_atan2<-Vectorize(S_Vector_Azimuth_Shifted_atan2)
Table2$S_Vector_Azimuth_Shifted_atan2<-S_Vector_Azimuth_Shifted_atan2(SCP_altitude_displacement = Table2$Latitude_displacement,
                                          SCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement = Table2$Observed_hour_angle,
                                          Time_shift = Table2$Time_shift)


Table_1000<-Table2



set.seed(123)

Table2<-Table_20000[sample(nrow(Table_20000), 100000, replace = FALSE), ]


##N_Binary_No_Shift####
#######################

N_Binary_No_Shift<-Vectorize(N_Binary_No_Shift)## This steps changes the functions so that they can take vectors (table columns), rather then scalar numbers

Table2$N_Binary_No_Shift<-N_Binary_No_Shift(NCP_altitude_displacement = Table2$Latitude_displacement,
                                          NCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement = Table2$Observed_hour_angle)



##N_Binary_Shifted##
####################
N_Binary_Shifted<-Vectorize(N_Binary_Shifted)

Table2$N_Binary_Shifted<-N_Binary_Shifted(NCP_altitude_displacement = Table2$Latitude_displacement,
                                          NCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement_no_shift = Table2$Observed_hour_angle,
                                          Time_shift = Table2$Time_shift )







#Separate functions for North-South and East-West components
############################################################
#North-South component
######################
#Time shift is irrelevant for this compoonent
N_Binary_NS<-Vectorize(N_Binary_NS )
Table2$N_Binary_NS<-N_Binary_NS(NCP_altitude_displacement = Table2$Latitude_displacement,
                                NCP_altitude_home = Table2$Latitude_home) 

#East-West component
####################
#No time shift
###############
N_Binary_EW_No_Shift<-Vectorize(N_Binary_EW_No_Shift)
Table2$N_Binary_EW_No_Shift<-N_Binary_EW_No_Shift(Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                                  Day_length = Table2$Day_length,
                                                  Observed_hour_angle_displacement = Table2$Observed_hour_angle
                                                  )

#With time shift
###############
N_Binary_EW_Shifted<-Vectorize(N_Binary_EW_Shifted)
Table2$N_Binary_EW_Shifted<-N_Binary_EW_Shifted(Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                                  Day_length = Table2$Day_length,
                                                  Observed_hour_angle_displacement_no_shift =  Table2$Observed_hour_angle,
                                                  Time_shift = Table2$Time_shift
                                                )

###############################################################
###############################################################
#VECTOR HOMING FUNCTIONS
###############################################################
###############################################################

#No time shift 
#############

#String function
###############

#Returns string: "Distance ------, Azimuth-------" For use in testing and isolated cases
N_Vector_String_No_Shift<-Vectorize(N_Vector_String_No_Shift)
Table2$N_Vector_String_No_Shift<-N_Vector_String_No_Shift(NCP_altitude_displacement = Table2$Latitude_displacement,
                                          NCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement = Table2$Observed_hour_angle)

N_Vector_String_No_Shift_atan2<-Vectorize(N_Vector_String_No_Shift_atan2)
Table2$N_Vector_String_No_Shift_atan2<-N_Vector_String_No_Shift_atan2(NCP_altitude_displacement = Table2$Latitude_displacement,
                                          NCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement = Table2$Observed_hour_angle)
#No time shift - Separate function for distance to home and home azimuth for use is analysis 
###########################################################################################

#Straight line distance to home function
#######################################

N_Vector_distance_No_Shift<-Vectorize(N_Vector_distance_No_Shift)
Table2$N_Vector_distance_No_Shift<-N_Vector_distance_No_Shift(NCP_altitude_displacement = Table2$Latitude_displacement,
                                          NCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement = Table2$Observed_hour_angle)
#Straight line azimuth home function
###################################
N_Vector_Azimuth_No_Shift<-Vectorize(N_Vector_Azimuth_No_Shift)
Table2$N_Vector_Azimuth_No_Shift<-N_Vector_Azimuth_No_Shift(NCP_altitude_displacement = Table2$Latitude_displacement,
                                          NCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement = Table2$Observed_hour_angle)


N_Vector_Azimuth_No_Shift_atan2<-Vectorize(N_Vector_Azimuth_No_Shift_atan2)
Table2$N_Vector_Azimuth_No_Shift_atan2<-N_Vector_Azimuth_No_Shift_atan2(NCP_altitude_displacement = Table2$Latitude_displacement,
                                          NCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement = Table2$Observed_hour_angle)


#With time shift 
#############

#String function
###############

#Returns string: "Distance ------, Azimuth-------" For use in testing and isolated cases

N_Vector_String_Shifted<-Vectorize(N_Vector_String_Shifted)
Table2$N_Vector_String_Shifted<-N_Vector_String_Shifted(NCP_altitude_displacement = Table2$Latitude_displacement,
                                          NCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement_no_shift = Table2$Observed_hour_angle,
                                          Time_shift = Table2$Time_shift )

N_Vector_String_Shifted_atan2<-Vectorize(N_Vector_String_Shifted_atan2)
Table2$N_Vector_String_Shifted_atan2<-N_Vector_String_Shifted_atan2(NCP_altitude_displacement = Table2$Latitude_displacement,
                                          NCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement_no_shift = Table2$Observed_hour_angle,
                                          Time_shift = Table2$Time_shift )
#With time shift - Separate function for distance to home and home azimuth for use is analysis 
################################################################################################

#Straight line distance to home function
########################################
N_Vector_Distance_Shifted<-Vectorize(N_Vector_Distance_Shifted)
Table2$N_Vector_Distance_Shifted<-N_Vector_Distance_Shifted(NCP_altitude_displacement = Table2$Latitude_displacement,
                                          NCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement_no_shift = Table2$Observed_hour_angle,
                                          Time_shift = Table2$Time_shift )

#Straight line azimuth home function
####################################
N_Vector_Azimuth_Shifted<-Vectorize(N_Vector_Azimuth_Shifted)
Table2$N_Vector_Azimuth_Shifted<-N_Vector_Azimuth_Shifted(NCP_altitude_displacement = Table2$Latitude_displacement,
                                          NCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement_no_shift = Table2$Observed_hour_angle,
                                          Time_shift = Table2$Time_shift )

N_Vector_Azimuth_Shifted_atan2<-Vectorize(N_Vector_Azimuth_Shifted_atan2)
Table2$N_Vector_Azimuth_Shifted_atan2<-N_Vector_Azimuth_Shifted_atan2(NCP_altitude_displacement = Table2$Latitude_displacement,
                                          NCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement_no_shift = Table2$Observed_hour_angle,
                                          Time_shift = Table2$Time_shift )
####Southern function were not yet adressed - The entire comparison df is for northern based on the geosphere function definition. Run the functions and then
###Figure it out

#Southern Hemisphere
#Full Binary function - No time shift:
########################################

S_Binary_No_Shift<-Vectorize(S_Binary_No_Shift)
Table2$S_Binary_No_Shift<-S_Binary_No_Shift(SCP_altitude_displacement = Table2$Latitude_displacement,
                                          SCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement = Table2$Observed_hour_angle)

#Full Binary function - with time shift
########################################
#Time shift is integrated in the observed hour angle
S_Binary_Shifted<-Vectorize(S_Binary_Shifted)
Table2$S_Binary_Shifted<-S_Binary_Shifted(SCP_altitude_displacement = Table2$Latitude_displacement,
                                          SCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement = Table2$Observed_hour_angle,
                                          Time_shift = Table2$Time_shift)
#Separate functions for North-South and East-West components
############################################################
#North-South component
######################
#Time shift is irrelevant for this compoonent

S_Binary_NS<-Vectorize(S_Binary_NS)
Table2$S_Binary_NS<-S_Binary_NS(SCP_altitude_displacement = Table2$Latitude_displacement,
                                          SCP_altitude_home = Table2$Latitude_home)

#East-West component
####################

#No time shift
###############
S_Binary_EW_No_Shift<-Vectorize(S_Binary_EW_No_Shift)
Table2$S_Binary_EW_No_Shift<-S_Binary_EW_No_Shift(Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement = Table2$Observed_hour_angle)

#With time shift
###############

S_Binary_EW_Shifted<-Vectorize(S_Binary_EW_Shifted)
Table2$S_Binary_EW_Shifted<-S_Binary_EW_Shifted(Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement = Table2$Observed_hour_angle,
                                          Time_shift = Table2$Time_shift)

###############################################################
###############################################################
#VECTOR HOMING FUNCTIONS
###############################################################
###############################################################

#No time shift 
#############

#String function
###############

#Returns string: "Distance ------, Azimuth-------" For use in testing and isolated cases
S_Vector_String_No_Shift<-Vectorize(S_Vector_String_No_Shift)
Table2$S_Vector_String_No_Shift<-S_Vector_String_No_Shift(SCP_altitude_displacement = Table2$Latitude_displacement,
                                          SCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement = Table2$Observed_hour_angle)

S_Vector_String_No_Shift_atan2<-Vectorize(S_Vector_String_No_Shift_atan2)
Table2$S_Vector_String_No_Shift_atan2<-S_Vector_String_No_Shift_atan2(SCP_altitude_displacement = Table2$Latitude_displacement,
                                          SCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement = Table2$Observed_hour_angle)

#No time shift - Separate function for distance to home and home azimuth for use is analysis 
###########################################################################################

#Straight line distance to home function
#######################################
S_Vector_distance_No_Shift<-Vectorize(S_Vector_distance_No_Shift)
Table2$S_Vector_distance_No_Shift<-S_Vector_distance_No_Shift(SCP_altitude_displacement = Table2$Latitude_displacement,
                                          SCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement = Table2$Observed_hour_angle)

#Straight line azimuth home function
###################################
S_Vector_Azimuth_No_Shift<-Vectorize(S_Vector_Azimuth_No_Shift)
Table2$S_Vector_Azimuth_No_Shift<-S_Vector_Azimuth_No_Shift(SCP_altitude_displacement = Table2$Latitude_displacement,
                                          SCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement = Table2$Observed_hour_angle)

S_Vector_Azimuth_No_Shift_atan2<-Vectorize(S_Vector_Azimuth_No_Shift_atan2)
Table2$S_Vector_Azimuth_No_Shift_atan2<-S_Vector_Azimuth_No_Shift_atan2(SCP_altitude_displacement = Table2$Latitude_displacement,
                                          SCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement = Table2$Observed_hour_angle)
#With time shift 
#############

#String function
###############

#Returns string: "Distance ------, Azimuth-------" For use in testing and isolated cases
S_Vector_String_Shifted<-Vectorize(S_Vector_String_Shifted)
Table2$S_Vector_String_Shifted<-S_Vector_String_Shifted(SCP_altitude_displacement = Table2$Latitude_displacement,
                                          SCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement = Table2$Observed_hour_angle,
                                          Time_shift = Table2$Time_shift)

S_Vector_String_Shifted_atan2<-Vectorize(S_Vector_String_Shifted_atan2)
Table2$S_Vector_String_Shifted_atan2<-S_Vector_String_Shifted_atan2(SCP_altitude_displacement = Table2$Latitude_displacement,
                                          SCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement = Table2$Observed_hour_angle,
                                          Time_shift = Table2$Time_shift)

#With time shift - Separate function for distance to home and home azimuth for use is analysis 
################################################################################################

#Straight line distance to home function
########################################
S_Vector_Distance_Shifted<-Vectorize(S_Vector_Distance_Shifted)
Table2$S_Vector_Distance_Shifted<-S_Vector_Distance_Shifted(SCP_altitude_displacement = Table2$Latitude_displacement,
                                          SCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement = Table2$Observed_hour_angle,
                                          Time_shift = Table2$Time_shift)

#Straight line azimuth home function
####################################
S_Vector_Azimuth_Shifted<-Vectorize(S_Vector_Azimuth_Shifted)
Table2$S_Vector_Azimuth_Shifted<-S_Vector_Azimuth_Shifted(SCP_altitude_displacement = Table2$Latitude_displacement,
                                          SCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement = Table2$Observed_hour_angle,
                                          Time_shift = Table2$Time_shift)

S_Vector_Azimuth_Shifted_atan2<-Vectorize(S_Vector_Azimuth_Shifted_atan2)
Table2$S_Vector_Azimuth_Shifted_atan2<-S_Vector_Azimuth_Shifted_atan2(SCP_altitude_displacement = Table2$Latitude_displacement,
                                          SCP_altitude_home = Table2$Latitude_home,
                                          Time_from_sunrise_at_displacement = Table2$Time_from_sunrise,
                                          Day_length = Table2$Day_length,
                                          Observed_hour_angle_displacement = Table2$Observed_hour_angle,
                                          Time_shift = Table2$Time_shift)

Table_20000<-Table2

#save(Table2,file="Test_data_V1_lat0_90_lon0_180_shift_0_6")
Table2<-rbind(Table_100,Table_1000,Table_20000)
save(Table2,file="Simulation_data_V9_shift_min6_6_by_0.5")
```

Analisys and presentation

Figure 3


```{r}
load("Supplamentary_data")

#Northern Hemisphere

Table2$N_Binary_No_Shift<-as.circular(Table2$N_Binary_No_Shift, type = "angles", units ="degrees", template = "geographics")
Table2$Azimuth_rhumb_geo<-as.circular(Table2$Azimuth_rhumb_geo, type = "angles", units ="degrees", template = "geographics")
Table2$angular_diff<-Table2$N_Binary_No_Shift-Table2$Azimuth_rhumb_geo
Table2$angular_diff<-as.numeric(Table2$angular_diff)
Table2$abs_angular_diff<-abs(Table2$angular_diff)
Table2$N_Binary_No_Shift<-as.numeric(Table2$N_Binary_No_Shift)
Table2$Azimuth_rhumb_geo<-as.numeric(Table2$Azimuth_rhumb_geo)
Table2$abs_angular_diff<-as.numeric(Table2$abs_angular_diff)
Table2$angular_diff<-as.numeric(Table2$angular_diff)
Table2$abs_angular_diff<-as.numeric(abs(Table2$angular_diff))
Table2$angular_diff<-ifelse(Table2$angular_diff>359,0,Table2$angular_diff)#These mean 0 difference
Table2$angular_diff<-ifelse(Table2$angular_diff>179,0,Table2$angular_diff)#single one, milimeters apart
t<-Table2%>%filter(angular_diff>170)
#Mechanism output vs real distance
Hundred<-Table2%>%filter(range==100)
range(Hundred$angular_diff)
sum(Hundred$angular_diff>170)
#100km

Hundreds_error_distribution_azimuth_N<-ggplot(Hundred, aes(x = angular_diff)) +
 geom_histogram(aes(y = after_stat(count / sum(count))),fill="#7aa0c0",colour = "black",size = 0.01) +
  scale_y_continuous(labels = scales::percent) +theme_minimal()+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1),axis.title.y = element_blank())+labs(title =  "> 0 km\n< 100 km")+ theme(plot.title = element_text(size = 10,hjust = 0.5),axis.text.y = element_blank(),axis.title.x = element_blank())
range(Hundred$angular_diff)
sum(Hundred$angular_diff>50)
##Use thousand circle only - 100, 20,000 are identical, will be reported verbally
#1000km
Thousand<-Table2%>%filter(range==1000)##Maybe increase sample size

Thousand_error_distribution_azimuth_N<-ggplot(Thousand, aes(x = angular_diff)) +
 geom_histogram(aes(y = after_stat(count / sum(count))),fill="#7aa0c0",colour = "black") +
  scale_y_continuous(labels = scales::percent) +theme_minimal()+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1),axis.title.y = element_blank())+labs(title =  "> 100 km\n< 1,000 km")+ theme(plot.title = element_text(size = 10,hjust = 0.5),axis.text.y = element_blank(),axis.title.x = element_blank())

Thousand_plot_azimuth_N<-ggplot(Thousand,aes(x = Azimuth_rhumb_geo, y = N_Binary_No_Shift)) +
  geom_point(aes(color=abs_angular_diff),size=1) +
  geom_smooth(method = "lm", col = "blue") + # Adds linear regression line
  labs(x = "Angle = Rhumb azimuth (deg)",
       y = "Radial distance = Mechanism output (deg)")+
  theme_minimal()+coord_polar()+scale_y_continuous(expand = (0),breaks = c(0, 45, 90,135,180,224,270,315,360),limits = c(0, 360))+theme(legend.position = "bottom")+
  scale_color_continuous(
    name = "Error\nsize (°)", # Set the legend title
    guide = guide_colorbar(
  barwidth = unit(0.2, "cm"),
  barheight = unit(2, "cm")
))+
  theme(
    legend.position = "right"     # Aligns the legend to the right edge
  )+scale_x_continuous(breaks = c(0, 45, 90,135,180,224,270,315,360)) + 
  expand_limits(x = 0)


#10000km
Ten_thousand<-Table2%>%filter(range==20000)

Ten_thousand_error_distribution_azimuth_N<-ggplot(Ten_thousand, aes(x = angular_diff)) +
 geom_histogram(aes(y = after_stat(count / sum(count))),fill="#7aa0c0",colour = "black") +
  scale_y_continuous(labels = scales::percent)  +theme_minimal()+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1),axis.title.y = element_blank())+labs(title =  "> 1,000 km\n< 20,000 km")+ theme(plot.title = element_text(size = 10,hjust = 0.5),axis.text.y = element_blank(),axis.title.x = element_blank())

model_1000<-lm(Thousand$N_Binary_No_Shift~Thousand$Azimuth_rhumb_geo)
summary(model_1000)

model_100<-lm(Hundred$N_Binary_No_Shift~Hundred$Azimuth_rhumb_geo)
summary(model_100)

model_10000<-lm(Ten_thousand$N_Binary_No_Shift~Ten_thousand$Azimuth_rhumb_geo)
summary(model_10000)



###Southern Hemisphere

Table2$S_Binary_No_Shift<-as.circular(Table2$S_Binary_No_Shift, type = "angles", units ="degrees", template = "geographics")
Table2$Azimuth_rhumb_geo<-as.circular(Table2$Azimuth_rhumb_geo, type = "angles", units ="degrees", template = "geographics")
Table2$angular_diff<-Table2$S_Binary_No_Shift-Table2$Azimuth_rhumb_geo_S
Table2$angular_diff<-as.numeric(Table2$angular_diff)
Table2$abs_angular_diff<-abs(Table2$angular_diff)
Table2$S_Binary_No_Shift<-as.numeric(Table2$S_Binary_No_Shift)
Table2$Azimuth_rhumb_geo<-as.numeric(Table2$Azimuth_rhumb_geo)
Table2$abs_angular_diff<-as.numeric(Table2$abs_angular_diff)
Table2$angular_diff<-as.numeric(Table2$angular_diff)

Table2$angular_diff<-ifelse(Table2$angular_diff>359,0,Table2$angular_diff)#These mean 0 difference
Table2$angular_diff<-ifelse(Table2$angular_diff>179,0,Table2$angular_diff)#single one, milimeters apart
Table2$abs_angular_diff<-as.numeric(abs(Table2$angular_diff))
range(Table2$angular_diff)

#Mechanism output vs real distance
Hundred<-Table2%>%filter(range==100)
#100km

Hundreds_error_distribution_azimuth_S<-ggplot(Hundred, aes(x = angular_diff)) +
 geom_histogram(aes(y = after_stat(count / sum(count))),fill="#7aa0c0",colour = "black",size = 0.01) +
  scale_y_continuous(labels = scales::percent) +theme_minimal()+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1),axis.title.y = element_blank())+labs(title =  "> 0 km\n< 100 km")+ theme(plot.title = element_text(size = 10,hjust = 0.5),axis.text.y = element_blank(),axis.title.x = element_blank())


##Use thousand circle only - 100, 20,000 are identical, will be reported verbally
#1000km
Thousand<-Table2%>%filter(range==1000)##Maybe increase sample size

Thousand_error_distribution_azimuth_S<-ggplot(Thousand, aes(x = angular_diff)) +
 geom_histogram(aes(y = after_stat(count / sum(count))),fill="#7aa0c0",colour = "black") +
  scale_y_continuous(labels = scales::percent) +theme_minimal()+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1),axis.title.y = element_blank())+labs(title =  "> 100 km\n< 1,000 km")+ theme(plot.title = element_text(size = 10,hjust = 0.5),axis.text.y = element_blank(),axis.title.x = element_blank())

Thousand_plot_azimuth_S<-ggplot(Thousand,aes(x = Azimuth_rhumb_geo_S, y = S_Binary_No_Shift)) +
  geom_point(aes(color=abs_angular_diff),size=1) +
  geom_smooth(method = "lm", col = "blue") + # Adds linear regression line
  labs(x = "Angle = Rhumb azimuth (deg)",
       y = "Radial distance = Mechanism output (deg)")+scale_y_continuous(expand = (0),breaks = c(0, 45, 90,135,180,224,270,315,360),limits = c(0, 360))+
  theme_minimal()+coord_polar()+theme(legend.position = "bottom")+
  scale_color_continuous(
    name = "Error\nsize (°)", # Set the legend title
    guide = guide_colorbar(
  barwidth = unit(0.2, "cm"),
  barheight = unit(2, "cm")
))+
  theme(
    legend.position = "right"     # Aligns the legend to the right edge
  )+scale_x_continuous(breaks = c(0, 45, 90,135,180,224,270,315,360))



#10000km
Ten_thousand<-Table2%>%filter(range==20000)

Ten_thousand_error_distribution_azimuth_S<-ggplot(Ten_thousand, aes(x = angular_diff)) +
 geom_histogram(aes(y = after_stat(count / sum(count))),fill="#7aa0c0",colour = "black") +
  scale_y_continuous(labels = scales::percent)  +theme_minimal()+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1),axis.title.y = element_blank())+labs(title =  "> 1,000 km\n< 20,000 km")+ theme(plot.title = element_text(size = 10,hjust = 0.5),axis.text.y = element_blank(),axis.title.x = element_blank())

model_100<-lm(Hundred$S_Binary_No_Shift~Hundred$Azimuth_rhumb_geo_S)
summary(model_100)

model_1000<-lm(Thousand$S_Binary_No_Shift~Thousand$Azimuth_rhumb_geo_S)
summary(model_1000)



model_10000<-lm(Ten_thousand$S_Binary_No_Shift~Ten_thousand$Azimuth_rhumb_geo_S)
summary(model_10000)




Binary_mechanism_summary<-ggdraw() +
  theme(plot.background = element_rect(fill="white", color = NA))+
   draw_label("Northern Hemisphere", x = 0.5, y = 0.98, size = 18, fontface = "bold") +
    draw_plot(Thousand_plot_azimuth_N, x = 0., y = 0.45, width =0.6 , height = .6) +
   draw_plot(Hundreds_error_distribution_azimuth_N , x = 0.55, y = 0.5, width =0.15 , height = .2) +
   draw_plot(Thousand_error_distribution_azimuth_N, x = 0.70, y = 0.5, width =0.15 , height = .2) +
   draw_plot(Ten_thousand_error_distribution_azimuth_N, x = 0.85, y = 0.5, width =0.15 , height = .2) +
  draw_line(x = c(0, 1),  y = c(.5, .5), color = "black",size = 0.5)+
   draw_label("Regression equations: \nMechanism output ~ Homing azimuth", x = 0.78, y = 0.94, size = 12, fontface = "bold")+
  draw_line(x = c(0.59, .97),  y = c(.915, .915), color = "black",size = 0.5)+
    draw_label("Homing Distance < 100km ", x = 0.78, y = 0.9, size = 12,fontface = "bold")+
   draw_label(expression(paste("y = 11.23 + 0.93x , p<0.0001, R"^2,"=0.94 ")), x = 0.79, y = 0.878, size = 12)+
    draw_label("100km < Homing Distance < 1,000km ", x = 0.79, y = 0.845, size = 12,fontface = "bold")+
draw_label(expression(paste("y = 11.15 + 0.93x , p<0.0001, R"^2,"=0.94 ")), x = 0.79, y = 0.82, size = 12)+
    draw_label("1,000km < Homing Distance < 20,000km ", x = 0.79, y = 0.79, size = 12,fontface = "bold")+
draw_label(expression(paste("y = 9.95 + 0.95x , p<0.0001, R"^2,"=0.89 ")), x = 0.79, y = 0.765, size = 12)+
  draw_label("Azimuth Error Distributions (deg)", x = 0.78, y = 0.72, size = 12, fontface = "bold") +
  draw_line(x = c(0.61, .95),  y = c(.705, .705), color = "black",size = 0.5)+
    draw_label("Southern Hemisphere", x = 0.5, y = 0.48, size = 18, fontface = "bold") +
   draw_plot(Thousand_plot_azimuth_S, x = 0., y = -0.05, width =0.6 , height = .6) +
   draw_plot(Hundreds_error_distribution_azimuth_S , x = 0.55, y = 0., width =0.15 , height = .2) +
   draw_plot(Thousand_error_distribution_azimuth_S, x = 0.70, y = 0., width =0.15 , height = .2) +
   draw_plot(Ten_thousand_error_distribution_azimuth_S, x = 0.85, y = 0, width =0.15 , height = .2) +
   draw_label("Regression equations: \nMechanism output ~ Homing azimuth", x = 0.78, y = 0.44, size = 12, fontface = "bold")+
  draw_line(x = c(0.59, .97),  y = c(.415, .415), color = "black",size = 0.5)+
    draw_label("Homing Distance < 100km ", x = 0.78, y = 0.4, size = 12,fontface = "bold")+
   draw_label(expression(paste("y = 11.5 + 0.93x , p<0.0001, R"^2,"=0.93 ")), x = 0.79, y = 0.378, size = 12)+
    draw_label("100km < Homing Distance < 1,000km ", x = 0.79, y = 0.345, size = 12,fontface = "bold")+
draw_label(expression(paste("y = 11.92 + 0.98x , p<0.0001, R"^2,"=0.94 ")), x = 0.79, y = 0.32, size = 12)+
    draw_label("1,000km < Homing Distance < 20,000km ", x = 0.79, y = 0.29, size = 12,fontface = "bold")+
draw_label(expression(paste("y = 16.3 + 0.91x , p<0.0001, R"^2,"=0.96 ")), x = 0.79, y = 0.265, size = 12)+
  draw_label("Azimuth Error Distributions (deg)", x = 0.78, y = 0.22, size = 12, fontface = "bold") +
  draw_line(x = c(0.61, .95),  y = c(.205, .205), color = "black",size = 0.5)+
  draw_label(" A.", x = 0.01, y = 0.98, size = 12, fontface = "bold")+
  draw_label(" B.", x = 0.01, y = 0.48, size = 12, fontface = "bold")




ggsave("Binary_no_shift_summary.3.png", width =8 , height =9,dpi=300)

ggplot(Hundred, aes(x = Longitude_home, y = angular_diff)) +
  geom_point(size=0.1,aes()) +
  geom_smooth(method = "lm", col = "blue") + # Adds linear regression line
  labs( )+scale_y_continuous(expand = (0))+
  stat_regline_equation(label.y = 100, aes(label =  paste(..eq.label..,  sep = "~~~~"))) +
   stat_regline_equation(label.y = 95, aes(label =  paste(..rr.label.., sep = "~~~~"))) +
  theme_minimal()+
  scale_color_continuous(
    name = "Error\nsize (%)", # Set the legend title
    labels = scales::percent,# Format labels as percentages
    guide = guide_colorbar(
  barwidth = unit(0.2, "cm"),
  barheight = unit(2, "cm")))

  
  
```
Figure 4
Vector test

```{r}
load("Supplamentary_data")
#load("Simulation_data_V9.1_Time_shift_with_decimal_corrected")
#plots:
#Mechanism output vs real distance


Table2$Distance_diff<-Table2$N_Vector_distance_No_Shift-Table2$distance_km_geo
hist(Table2$Distance_diff)

#Looks like our function estiamte has more larger then shorter distances compared - maybe because its the geosphere outputs shortest difference and we scalculate fixed azimuth. Test geosphere's rhumb distance as well

##Calculate the relative size of the error - relative to the geosphere output of shortest distance
Table2$Distance_diff_rel<-Table2$Distance_diff/Table2$distance_km_geo
Table2$Distance_diff_rel_abs<-abs(Table2$Distance_diff/Table2$distance_km_geo)

Table2$range<-factor(Table2$range)

#100km
Hundred<-Table2%>%filter(range==100)

Hundreds_error_distribution_distance<-ggplot(Hundred, aes(x = Distance_diff_rel)) +
  geom_histogram(aes(y = after_stat(count / sum(count)))) +
  scale_y_continuous(labels = scales::percent)  +
   scale_x_continuous(labels = scales::percent) +theme_minimal()+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1),axis.title.y = element_blank())+labs(title =  "Error\ndistribution")+ theme(plot.title = element_text(size = 10,hjust = 0.5),axis.text.y = element_blank(),axis.title.x = element_blank())


Hundred_plot_distance<-ggplot(Hundred, aes(x = distance_km_rhumb_geo, y = N_Vector_distance_No_Shift)) +
  geom_point(size=0.1,aes(color=Distance_diff_rel_abs)) +
  geom_smooth(method = "lm", col = "blue") + # Adds linear regression line
  labs( x = "Rhumb distance (km)",
       y = "Mechanism output (km)")+scale_y_continuous(expand = (0))+
  stat_regline_equation(label.y = 100, aes(label =  paste(..eq.label..,  sep = "~~~~"))) +
   stat_regline_equation(label.y = 95, aes(label =  paste(..rr.label.., sep = "~~~~"))) +
  theme_minimal()+
  scale_color_continuous(
    name = "Error\nsize (%)", # Set the legend title
    labels = scales::percent,# Format labels as percentages
    guide = guide_colorbar(
  barwidth = unit(0.2, "cm"),
  barheight = unit(2, "cm")
),
  breaks = seq(0, 0.25, by = 0.05))+
  theme(
    legend.position = "right",             # Places the legend above the plot
    legend.justification = "top",
    legend.title = element_text(hjust = 0.5)# Aligns the legend to the right edge
  )



Combined_hundered_distance<-ggdraw() +
  theme(plot.background = element_rect(fill="white", color = NA))+
  draw_plot(Hundred_plot_distance, x = 0, y = 0.0, width =1 , height = 1) +
  draw_plot(Hundreds_error_distribution_distance, x = 0.8, y = 0.08, width =.2 , height = 0.5)



#1000km
Thousand<-Table2%>%filter(range==1000)##Maybe increase sample size

Thousands_error_distribution_distance<-ggplot(Thousand, aes(x = Distance_diff_rel)) +
 geom_histogram(aes(y = after_stat(count / sum(count)))) +
  scale_x_continuous(labels = scales::percent)  +theme_minimal()+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1),axis.title.y = element_blank())+labs(title =  "Error\ndistribution")+ theme(plot.title = element_text(size = 10,hjust = 0.5),axis.text.y = element_blank(),axis.title.x = element_blank())


Thousand_plot_distance<-ggplot(Thousand, aes(x = distance_km_rhumb_geo, y = N_Vector_distance_No_Shift)) +
  geom_point(size=0.1,aes(color=Distance_diff_rel)) +
  geom_smooth(method = "lm", col = "blue") + # Adds linear regression line
  labs(y = "Mechanism output (km)",
       x = "Rhumb distance (km)")+scale_y_continuous(expand = (0))+
  stat_regline_equation(label.y = 2000, aes(label =  paste(..eq.label..,  sep = "~~~~"))) +
   stat_regline_equation(label.y = 1900, aes(label =  paste(..rr.label.., sep = "~~~~"))) +
  theme_minimal()+
  scale_color_continuous(
    name = "Error\nsize (%)", # Set the legend title
    labels = scales::percent, # Format labels as percentages
    guide = guide_colorbar(
  barwidth = unit(0.2, "cm"),
  barheight = unit(2, "cm")),
  breaks = seq(0, 1, by = 0.2)
)+
  theme(
    legend.position = "right",             # Places the legend above the plot
    legend.justification = "top",
    legend.title = element_text(hjust = 0.5)# Aligns the legend to the right edge
  )


Combined_Thousand_distance<-ggdraw() +
  theme(plot.background = element_rect(fill="white", color = NA))+
  draw_plot(Thousand_plot_distance, x = 0, y = 0.0, width =1 , height = 1) +
  draw_plot(Thousands_error_distribution_distance, x = 0.8, y = 0.08, width =.2 , height = 0.5)


#10000km
Ten_thousand<-Table2%>%filter(range==20000)

Ten_thousands_error_distribution_distance<-ggplot(Ten_thousand, aes(x = Distance_diff_rel)) +
geom_histogram(aes(y = after_stat(count / sum(count)))) +
  scale_x_continuous(labels = scales::percent)  +theme_minimal()+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1),axis.title.y = element_blank())+labs(title =  "Error\ndistribution")+ theme(plot.title = element_text(size = 10,hjust = 0.5),axis.text.y = element_blank(),axis.title.x = element_blank())

hist(Ten_thousand$Distance_diff_rel)
Ten_thousand_plot_distance<-ggplot(Ten_thousand, aes(x = distance_km_rhumb_geo, y = N_Vector_distance_No_Shift)) +
  geom_point(size=0.1,aes(color=Distance_diff_rel)) +
  geom_smooth(method = "lm", col = "blue") + # Adds linear regression line
  labs( x = "Rhumb distance (km)",
       y = "Mechanism output (km)")+scale_y_continuous(expand = (0))+
  stat_regline_equation(label.y = 20000,label.x = 100, aes(label =  paste(..eq.label..,  sep = "~~~~"))) +
   stat_regline_equation(label.y = 19000,label.x = 100, aes(label =  paste(..rr.label.., sep = "~~~~"))) +
  theme_minimal()+
  scale_color_continuous(
    name = "Error\nsize (%)", # Set the legend title
    labels = scales::percent, # Format labels as percentages
    guide = guide_colorbar(
  barwidth = unit(0.2, "cm"),
  barheight = unit(2, "cm")
))+
  theme(
    legend.position = "right",             # Places the legend above the plot
    legend.justification = "top" ,
    legend.title = element_text(hjust = 0.5)# Aligns the legend to the right edge
  )

Combined_Ten_thousand_distance<-ggdraw() +
  theme(plot.background = element_rect(fill="white", color = NA))+
  draw_plot(Ten_thousand_plot_distance, x = 0, y = 0.0, width =1 , height = 1) +
  draw_plot(Ten_thousands_error_distribution_distance, x = 0.8, y = 0.08, width =.2 , height = 0.5)





##Now same fpr azimuth

Table2$N_Vector_Azimuth_No_Shift_atan2<-as.circular(Table2$N_Vector_Azimuth_No_Shift, type = "angles", units ="degrees", template = "geographics")
Table2$Azimuth_rhumb_geo<-as.circular(Table2$Azimuth_rhumb_geo, type = "angles", units ="degrees", template = "geographics")
Table2$angular_diff<-Table2$N_Vector_Azimuth_No_Shift_atan2-Table2$Azimuth_rhumb_geo
Table2$angular_diff<-as.numeric(Table2$angular_diff)
Table2$abs_angular_diff<-abs(Table2$angular_diff)
Table2$N_Vector_Azimuth_No_Shift_atan2<-as.numeric(Table2$N_Vector_Azimuth_No_Shift_atan2)
Table2$Azimuth_rhumb_geo<-as.numeric(Table2$Azimuth_rhumb_geo)
Table2$abs_angular_diff<-as.numeric(Table2$abs_angular_diff)
Table2$angular_diff<-as.numeric(Table2$angular_diff)
Table2$angular_diff<-ifelse(Table2$angular_diff>359,0,Table2$angular_diff)#These mean 0 difference
Table2$angular_diff<-ifelse(Table2$angular_diff>179,0,Table2$angular_diff)#single one, milimeters apart
Table2$abs_angular_diff<-as.numeric(abs(Table2$angular_diff))
#Mechanism output vs real distance
Hundred<-Table2%>%filter(range==100)
#100km

Hundreds_error_distribution_azimuth<-ggplot(Hundred, aes(x = angular_diff)) +
 geom_histogram(aes(y = after_stat(count / sum(count)))) +
  scale_y_continuous(labels = scales::percent) +theme_minimal()+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1),axis.title.y = element_blank())+labs(title =  "Error\ndistribution")+ theme(plot.title = element_text(size = 10,hjust = 0.5),axis.text.y = element_blank(),axis.title.x = element_blank())

Hundred_plot_azimuth<-ggplot(Hundred, aes(x = Azimuth_rhumb_geo, y = N_Vector_Azimuth_No_Shift)) +
  geom_point(aes(color=abs_angular_diff),size=0.1) +
  geom_smooth(method = "lm",color = "blue") + # Adds linear regression line
  labs(x = "Angle = Rhumb azimuth (deg)",
       y = "Radial distance = Mechanism output (deg)")+scale_y_continuous(expand = (0),breaks = c(0, 45, 90,135,180,224,270,315,360),limits = c(0, 360))+
  stat_regline_equation(label.y = 360,label.x = 230, aes(label =  paste(..eq.label.., ..rr.label.., sep = "~~~~"))) +
  theme_minimal()+coord_polar()+
  scale_color_continuous(
    name = "Error\nsize (°)", # Set the legend title
    guide = guide_colorbar(
  barwidth = unit(0.2, "cm"),
  barheight = unit(2, "cm")
))+
  theme(
    legend.position = "right",             # Places the legend above the plot
    legend.justification = "top"       # Aligns the legend to the right edge
  )+scale_x_continuous(breaks = c(0, 45, 90,135,180,224,270,315,360)) + 
  expand_limits(x = 0)



Combined_hundred_azimuth<-ggdraw() +
  theme(plot.background = element_rect(fill="white", color = NA))+
  draw_plot(Hundred_plot_azimuth, x = 0, y = 0.0, width =1 , height = 1) +
  draw_plot(Hundreds_error_distribution_azimuth, x = 0.8, y = 0.08, width =.2 , height = 0.45)



#1000km
Thousand<-Table2%>%filter(range==1000)##Maybe increase sample size

Thousand_error_distribution_azimuth<-ggplot(Thousand, aes(x = angular_diff)) +
 geom_histogram(aes(y = after_stat(count / sum(count)))) +
  scale_y_continuous(labels = scales::percent) +theme_minimal()+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1),axis.title.y = element_blank())+labs(title =  "Error\ndistribution")+ theme(plot.title = element_text(size = 10,hjust = 0.5),axis.text.y = element_blank(),axis.title.x = element_blank())

Thousand_plot_azimuth<-ggplot(Thousand,aes(x = Azimuth_rhumb_geo, y = N_Vector_Azimuth_No_Shift)) +
  geom_point(aes(color=abs_angular_diff),size=0.1) +
  geom_smooth(method = "lm", col = "blue") + # Adds linear regression line
 labs(x = "Angle = Rhumb azimuth (deg)",
       y = "Radial distance = Mechanism output (deg)")+scale_y_continuous(expand = (0),,breaks = c(0, 45, 90,135,180,224,270,315,360),limits = c(0, 360))+
  stat_regline_equation(label.y = 360,label.x = 230, aes(label =  paste(..eq.label.., ..rr.label.., sep = "~~~~"))) +
  theme_minimal()+coord_polar()+
  scale_color_continuous(
    name = "Error\nsize (°)", # Set the legend title
    guide = guide_colorbar(
  barwidth = unit(0.2, "cm"),
  barheight = unit(2, "cm")
))+
  theme(
    legend.position = "right",             # Places the legend above the plot
    legend.justification = "top"       # Aligns the legend to the right edge
  )+scale_x_continuous(breaks = c(0, 45, 90,135,180,224,270,315,360)) + 
  expand_limits(x = 0)


Combined_thousand_azimuth<-ggdraw() +
  theme(plot.background = element_rect(fill="white", color = NA))+
  draw_plot(Thousand_plot_azimuth, x = 0, y = 0.0, width =1 , height = 1) +
  draw_plot(Thousand_error_distribution_azimuth, x = 0.8, y = 0.08, width =.2 , height = 0.45)


#10000km
Ten_thousand<-Table2%>%filter(range==20000)

Ten_thousand_error_distribution_azimuth<-ggplot(Ten_thousand, aes(x = angular_diff)) +
 geom_histogram(aes(y = after_stat(count / sum(count)))) +
  scale_y_continuous(labels = scales::percent)  +theme_minimal()+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1),axis.title.y = element_blank())+labs(title =  "Error\ndistribution")+ theme(plot.title = element_text(size = 10,hjust = 0.5),axis.text.y = element_blank(),axis.title.x = element_blank())

Ten_thousand_plot_azimuth<-ggplot(Ten_thousand,aes(x = Azimuth_rhumb_geo, y = N_Vector_Azimuth_No_Shift)) +
  geom_point(aes(color=abs_angular_diff),size=0.1) +
  geom_smooth(method = "lm", col = "blue") + # Adds linear regression line
 labs(x = "Angle = Rhumb azimuth (deg)",
       y = "Radial distance = Mechanism output (deg)")+scale_y_continuous(expand = (0),,breaks = c(0, 45, 90,135,180,224,270,315,360),limits = c(0, 360))+
  stat_regline_equation(label.y = 360,label.x = 230, aes(label =  paste(..eq.label.., ..rr.label.., sep = "~~~~"))) +
  theme_minimal()+coord_polar()+
  scale_color_continuous(
    name = "Error\nsize (°)", # Set the legend title
    guide = guide_colorbar(
  barwidth = unit(0.2, "cm"),
  barheight = unit(2, "cm")
))+
  theme(
    legend.position = "right",             # Places the legend above the plot
    legend.justification = "top"       # Aligns the legend to the right edge
  )+scale_x_continuous(breaks = c(0, 45, 90,135,180,224,270,315,360)) + 
  expand_limits(x = 0)

Combined_ten_thousand_azimuth<-ggdraw() +
  theme(plot.background = element_rect(fill="white", color = NA))+
  draw_plot(Ten_thousand_plot_azimuth, x = 0, y = 0.0, width =1 , height = 1) +
  draw_plot(Ten_thousand_error_distribution_azimuth, x = 0.8, y = 0.08, width =.2 , height = 0.45)


##How many are accuarte here?
hist(Ten_thousand$angular_diff)
t<-Ten_thousand%>%filter(angular_diff<1)
t<-t%>%filter(angular_diff>-1)

filtered_df <- Ten_thousand %>% 
  filter(angular_diff >= -5 & angular_diff <= 5)

filtered_df <- Ten_thousand %>% 
  filter(Distance_diff_rel >= -5 & angular_diff <= 5)


#make the plot
vector_mechanism_summary<-ggdraw() +
  theme(plot.background = element_rect(fill="white", color = NA))+
  draw_plot(Combined_hundered_distance, x = 0, y = 0.66, width =0.5 , height = .33) +
    draw_plot(Combined_hundred_azimuth, x = 0.5, y = 0.66, width =0.5 , height = .33) +
     draw_plot(Combined_Thousand_distance, x = 0., y = 0.33, width =0.5 , height = .33) +
     draw_plot(Combined_thousand_azimuth, x = 0.5, y = 0.33, width =0.5 , height = .33) +
       draw_plot(Combined_Ten_thousand_distance, x = 0., y = 0., width =0.5 , height = .33) +
     draw_plot(Combined_ten_thousand_azimuth, x = 0.5, y = 0., width =0.5 , height = .33) +
  draw_line(x = c(0, 1),  y = c(.66, .66), color = "black",size = 0.5)+
  draw_line(x = c(0, 1),  y = c(.33, .33), color = "black",size = 0.5)+
  draw_label("100 KM,", x = 0.5, y = 0.675, size = 18, fontface = "bold") +
   draw_label("1,000 KM,", x = 0.5, y = 0.345, size = 18, fontface = "bold") +
   draw_label("20,000 KM,", x = 0.5, y = 0.015, size = 18, fontface = "bold")+
    draw_label(" N=100,000", x = 0.61, y = 0.672, size = 12, fontface = "bold") +
    draw_label(" N=100,000", x = 0.625, y = 0.342, size = 12, fontface = "bold") +
    draw_label(" N=100,000", x = 0.63, y = 0.012, size = 12, fontface = "bold") +
  draw_label(" A.", x = 0.01, y = 0.98, size = 12, fontface = "bold")+ 
   draw_label(" i.", x = 0.025, y = 0.978, size = 8, fontface = "bold")+
    draw_label(" ii.", x = 0.5, y = 0.978, size = 8, fontface = "bold")+
  draw_label(" B.", x = 0.01, y = 0.65, size = 12, fontface = "bold")+ 
   draw_label(" i.", x = 0.025, y = 0.648, size = 8, fontface = "bold")+
    draw_label(" ii.", x = 0.5, y = 0.648, size = 8, fontface = "bold")+
  draw_label(" C.", x = 0.01, y = 0.32, size = 12, fontface = "bold")+
   draw_label(" i.", x = 0.025, y = 0.318, size = 8, fontface = "bold")+ 
  draw_label(" ii.", x = 0.5, y = 0.318, size = 8, fontface = "bold")
ggsave("vector_no_shift_summary_3.5.png", width =8 , height =10,dpi=300)

#Composit plot of other parameters vs relative error
#Date
#Day length
#Time from sunrise
#Lats
#Lons

```

Figure 5
Time Shift
```{r}
load("Supplamentary_data")
unique(Table2$Time_shift)
##אני חושב שמה שקורה זה שב1 ו2 אנחנו רואים בערך את האפקט של השיפט + ממוצע הסטייה בלי שיפט. כך שעולים בשיפט
##Add time shift driven change in hour angle 
Table2$Time_shift<-as.numeric(Table2$Time_shift)
unique(Table2$Time_shift)

Table2$Shift_effect<-Table2$Time_shift*180/Table2$Day_length##Added change in ohur angle - Allways positive


##Calculate the observed hour angle based on the Expected formula (which uses time) and the difference in longitude between displacement and home.
Table2$Observed_hour_angle_shift<-Table2$Time_from_sunrise*180/Table2$Day_length+(Table2$Longitude_displacement-Table2$Longitude_home)+Table2$Shift_effect

Table2$Expected_hour_angle_displacement<-Table2$Time_from_sunrise*180/Table2$Day_length
##Filter only home-displacement combinations that produce obsered hour angles between 0-180, below which and above which the sun is not in the sky
##at dispalcement. Expected hour angle doesnt need filtering becuase the forula restricts it to 0-180
Table2<-Table2%>%filter(as.numeric(Observed_hour_angle)>0)
Table2<-Table2%>%filter(as.numeric(Observed_hour_angle)<180)

###Scenarios with shift effect larger than the remainig hour angle range are unrealistic and should be dropped from teh dataset
##If the shift effect is larger than the daylight time remaining -expreseed as hour angle range left to move
##through the day given time from sunrise and day length, then the sun will set and its not in the game.
###############################################################################################################
Table2$Shift_effect_real<-ifelse(Table2$Observed_hour_angle_shift<180,"Yes","No")

##Real
Table2<-Table2%>%filter(Shift_effect_real=="Yes")
Table2$angular_diff_vec<-Table2$N_Vector_Azimuth_Shifted-Table2$Azimuth_rhumb_geo
range(Table2$angular_diff_vec)
#Table2$angular_diff_vec<-ifelse(Table2$angular_diff_vec<0,Table2$angular_diff_vec+360,Table2$angular_diff_vec)
#Table2$angular_diff_vec<-ifelse(Table2$angular_diff_vec>180,360-Table2$angular_diff_vec,Table2$angular_diff_vec)


Table2$angular_diff_bin<-Table2$N_Binary_Shifted-Table2$Azimuth_rhumb_geo
range(Table2$angular_diff_bin)

#Table2$angular_diff_bin<-ifelse(Table2$angular_diff_bin<0,Table2$angular_diff_bin+360,Table2$angular_diff_bin)
#Table2$angular_diff_bin<-ifelse(Table2$angular_diff_bin>180,360-Table2$angular_diff_bin,Table2$angular_diff_bin)
#range(Table2$angular_diff2)

##בוקס שמראה את האאוטפוט של השיפט מול האאופוט של הנון שיפטד
##Calculate the relative size of the error - relative to the geosphere output of shortest distance
Table2<-Table2%>%filter(Latitude_home>=0)
Table2<-Table2%>%filter(Longitude_home>=0)

test_results<-Table2%>%
  group_by(range,Time_shift)%>%
  dplyr::summarise(sample_size=n(),
                   Shift_effect=mean(Shift_effect),
                   median_angular_diff_vec=median(as.numeric(angular_diff_vec)),
                   median_angular_diff_bin=median(as.numeric(angular_diff_bin))
      )

#Table3<-Table2


t<-test_results
t$vec_bin_diff<-t$median_angular_diff_vec
t$mechanism<-rep("vec",nrow(t)) 

t2<-test_results
t2$vec_bin_diff<-t2$median_angular_diff_bin
t2$mechanism<-rep("bin",nrow(t2)) 
test_results<-rbind(t,t2)

test_results$range<-factor(test_results$range)
test_results$mechanism<-ifelse(test_results$mechanism=="bin","Binary Mechanism", "Vector Mechanism")
test_results$range<-ifelse(test_results$range==100,"100",ifelse(test_results$range==1000,"1,000","20,000"))
test_results$range<-factor(test_results$range,levels=c("100","1,000","20,000"))
Mean_angular_diff_plot<-ggplot(test_results, aes(x = Time_shift, y =vec_bin_diff,fill = range )) +
  geom_point(shape=21,color="black",size = 1) +geom_smooth(aes(color=range),se=F)+scale_fill_manual(values =  c("lightblue1","lightblue3","lightblue4"))+scale_color_manual(values =  c("lightblue1","lightblue3","lightblue4"))+theme_minimal()+theme(legend.position = "bottom")+facet_wrap(.~mechanism,ncol=2)+theme(strip.text = element_text(size = 14,face = "bold"))+scale_y_continuous(breaks = seq(-110, 160, by = 10),expand = c(0,0))+scale_x_continuous(breaks=seq(-6,6,by=0.5))+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))+labs(fill= "Displacement range (km)",color="Displacement range (km)")+ylab("[Shited function azimuth output - Real homing azimuth] (degrees)")+xlab("Clock shift (hours)")



Table2$range<-factor(Table2$range)
 Table2$Time_shift_round<-round(Table2$Time_shift,digits = 0)
 Table2$Time_shift_round<-factor(Table2$Time_shift_round,levels = c(-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6))
Distance_vec<-ggplot(data = Table2, aes(x = distance_km_rhumb_geo, y = angular_diff_vec, color = Time_shift_round)) +
  geom_smooth(method = "lm",se=F)  +
  theme_minimal()+scale_color_manual(values=c("#191970","#2b418a","#4f71a4","#7aa0c0","#b6d6e6","#bfefff","black","#FF8DA1","#FF708C","#FF3E60","#DC143C","#A52A2A","#8B0000"))+theme(legend.position = "right")+labs(color = "Clock\nshift\n(hours)")


Distance_bin<-ggplot(data = Table2, aes(x = distance_km_rhumb_geo, y = angular_diff_bin, color = Time_shift_round)) +
  geom_smooth(method = "lm")  + guides(color = guide_legend(override.aes = list(fill = NA), reverse = TRUE))+
  theme_minimal()+scale_color_manual(values=c("#191970","#2b418a","#4f71a4","#7aa0c0","#b6d6e6","#bfefff","black","#FF8DA1","#FF708C","#FF3E60","#DC143C","#A52A2A","#8B0000"))+theme(legend.position = "right")+labs(color = "Clock\nshift\n(hrs)")+xlab("Displacement distace \n(km)")+ylab("[Shited function azimuth output - Real homing azimuth] (degrees)")+scale_x_continuous(breaks=seq(0,20000,by=2000) ,labels = label_number(scale_cut = cut_short_scale()))+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))+scale_y_continuous(breaks = seq(-125, 130, by = 10),expand = c(0,0))+coord_cartesian(ylim = c(-125, 125))+ 
  theme(legend.key.width = unit(0.2, "cm"))

Table2$shift<-ifelse(Table2$Shift_effect==0,"None",ifelse(Table2$Shift_effect<0,"Slow","Fast"))
Table2$shift<-factor(Table2$shift,levels = c("Slow","None","Fast"))
Shift_effect_latitude<-ggplot(data = Table2, aes(x = Latitude_displacement, y = angular_diff_bin, color = shift)) +
  geom_smooth()  +
  theme_minimal()+scale_color_manual(values=c("#7aa0c0","black","#DC143C"))+theme(legend.position = "none")+labs(color = "Shift type")+xlab("Displacement Latitude \n(deg)")+scale_y_continuous(breaks = seq(-125, 130, by = 10),expand = c(0,0))+scale_x_continuous(breaks = seq(0, 90, by = 10),expand = c(0,0))+ylab("[Shited function azimuth output - Real homing azimuth] (degrees)")+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))

##Remove instances where the shift effect goes over sunset or befire sunrise
Table3<-Table2
Table3$t<-ifelse(Table3$Shift_effect+Table3$Observed_hour_angle_shift<0,T,F)
Table3<-Table3%>%filter(t!=T)
Table3$t<-ifelse(Table3$Shift_effect+Table3$Observed_hour_angle_shift>180,T,F)
Table3<-Table3%>%filter(t!=T)


  Shift_effect_plot<-ggplot(Table3, aes(x = factor(Time_shift_round), y = Shift_effect)) +
  geom_boxplot(outlier.size = 0.5,fill="lightblue1") +theme_minimal()+theme(legend.position = "none")+ylim(-180,180)+xlab("Clock shift \n(hours)")+ylab("Shift effect (hour angle degrees)")+scale_y_continuous(breaks = seq(-180, 180, by = 10),expand = c(0,0))
  
  
 plot<-ggdraw() +
  draw_plot(Shift_effect_plot, x = 0, y = 0.5, width =0.28, height = 0.5) +
     draw_plot(Shift_effect_latitude, x = 0.285, y = 0.5, width =0.28, height = 0.5) +
  draw_plot(Distance_bin, x = 0.57, y = 0.5, width =0.44, height = .5) +
     draw_plot( Mean_angular_diff_plot, x = 0., y = 0., width =1, height = .5) +  theme(plot.background = element_rect(fill="white", color = NA))+
   draw_label(" A.", x = 0.01, y = 0.99, size = 12, fontface = "bold")+ 
   draw_label(" B.", x = 0.28, y = 0.99, size = 12, fontface = "bold")+
      draw_label(" C.", x = 0.56, y = 0.99, size = 12, fontface = "bold")+
      draw_label(" D.", x = 0.04, y = 0.48, size = 12, fontface = "bold") +
   draw_label(" Fast shifts (+)", x = 0.45, y = 0.98, size = 10, fontface = "bold",color ="#8B0000" ) +
   draw_label(" No shift", x = 0.45, y = 0.79, size = 10, fontface = "bold") +
   draw_label(" Slow shifts (-)", x = 0.45, y = 0.58, size = 10, fontface = "bold",color="#191970") 
 
  
ggsave("Shift_plot2a.png", width = 7, height = 10,dpi=300)  

 


```

Table 1

```{r}


load("Supplamentary_data")
#load("Simulation_data_V9.1_Time_shift_with_decimal_corrected")
#plots:
#Mechanism output vs real distance


##North

Table2$Distance_diff<-Table2$N_Vector_distance_No_Shift-Table2$distance_km_geo

#Looks like our function estiamte has more larger then shorter distances compared - maybe because its the geosphere outputs shortest difference and we scalculate fixed azimuth. Test geosphere's rhumb distance as well

##Calculate the relative size of the error - relative to the geosphere output of shortest distance
Table2$Distance_diff_rel<-Table2$Distance_diff/Table2$distance_km_geo
Table2$Distance_diff_rel_abs<-abs(Table2$Distance_diff/Table2$distance_km_geo)



##Now same fpr azimuth

Table2$N_Vector_Azimuth_No_Shift_atan2<-as.circular(Table2$N_Vector_Azimuth_No_Shift, type = "angles", units ="degrees", template = "geographics")
Table2$Azimuth_rhumb_geo<-as.circular(Table2$Azimuth_rhumb_geo, type = "angles", units ="degrees", template = "geographics")
Table2$angular_diff<-Table2$N_Vector_Azimuth_No_Shift_atan2-Table2$Azimuth_rhumb_geo
Table2$angular_diff<-as.numeric(Table2$angular_diff)
Table2$abs_angular_diff<-abs(Table2$angular_diff)
Table2$N_Vector_Azimuth_No_Shift_atan2<-as.numeric(Table2$N_Vector_Azimuth_No_Shift_atan2)
Table2$Azimuth_rhumb_geo<-as.numeric(Table2$Azimuth_rhumb_geo)
Table2$abs_angular_diff<-as.numeric(Table2$abs_angular_diff)
Table2$angular_diff<-as.numeric(Table2$angular_diff)
Table2$angular_diff<-ifelse(Table2$angular_diff>359,0,Table2$angular_diff)#These mean 0 difference
Table2$angular_diff<-ifelse(Table2$angular_diff>179,0,Table2$angular_diff)#single one, milimeters apart
Table2$abs_angular_diff<-as.numeric(abs(Table2$angular_diff))




Table2$N_Binary_No_Shift<-as.circular(Table2$N_Binary_No_Shift, type = "angles", units ="degrees", template = "geographics")
Table2$Azimuth_rhumb_geo<-as.circular(Table2$Azimuth_rhumb_geo, type = "angles", units ="degrees", template = "geographics")
Table2$angular_diff_binary<-Table2$N_Binary_No_Shift-Table2$Azimuth_rhumb_geo
Table2$angular_diff_binary<-as.numeric(Table2$angular_diff_binary)
Table2$abs_angular_diff_binary<-abs(Table2$angular_diff_binary)
Table2$N_Binary_No_Shift<-as.numeric(Table2$N_Binary_No_Shift)
Table2$Azimuth_rhumb_geo<-as.numeric(Table2$Azimuth_rhumb_geo)
Table2$abs_angular_diff_binary<-as.numeric(Table2$abs_angular_diff_binary)
Table2$angular_diff_binary<-as.numeric(Table2$angular_diff_binary)
Table2$abs_angular_diff_binary<-as.numeric(abs(Table2$angular_diff_binary))
Table2$angular_diff_binary<-ifelse(Table2$angular_diff_binary>359,0,Table2$angular_diff_binary)#These mean 0 difference
Table2$angular_diff_binary<-ifelse(Table2$angular_diff_binary>179,0,Table2$angular_diff_binary)#single one, milimeters apart

Table_North<-Table2 %>%
  group_by(range) %>%
  summarize(n=n(),
  Mean_distance_rel=mean(Distance_diff_rel_abs),
   SD_distance_rel=sd(Distance_diff_rel_abs),
  max_distance_rel=max(Distance_diff_rel_abs),
   Percentile_95_distance_rel = quantile(Distance_diff_rel_abs, 0.95, na.rm = TRUE),
  Mean_distance=mean(abs(Distance_diff)),
   SD_distance=sd(abs(Distance_diff)),
  max_distance=max(abs(Distance_diff)),
   Percentile_95_distance = quantile(abs(Distance_diff), 0.95, na.rm = TRUE),
  Mean_azimuth_diff=mean(abs_angular_diff),
  SD_azimuth_diff=sd(abs_angular_diff),
  max_azimuth_diff=max(abs_angular_diff),
     Percentile_95_azimuth = quantile(abs_angular_diff, 0.95, na.rm = TRUE),
   Mean_azimuth_diff_binary=mean(abs_angular_diff_binary),
  SD_azimuth_diff_binary=sd(abs_angular_diff_binary),
  max_azimuth_diff_binary=max(abs_angular_diff_binary),
     Percentile_95_azimuth_binary = quantile(abs_angular_diff_binary, 0.95, na.rm = TRUE),)



##South


Table2$Distance_diff<-Table2$S_Vector_distance_No_Shift-Table2$distance_km_geo

#Looks like our function estiamte has more larger then shorter distances compared - maybe because its the geosphere outputs shortest difference and we scalculate fixed azimuth. Test geosphere's rhumb distance as well

##Calculate the relative size of the error - relative to the geosphere output of shortest distance
Table2$Distance_diff_rel<-Table2$Distance_diff/Table2$distance_km_geo
Table2$Distance_diff_rel_abs<-abs(Table2$Distance_diff/Table2$distance_km_geo)



##Now same fpr azimuth

Table2$S_Vector_Azimuth_No_Shift<-as.circular(Table2$S_Vector_Azimuth_No_Shift, type = "angles", units ="degrees", template = "geographics")
Table2$Azimuth_rhumb_geo<-as.circular(Table2$Azimuth_rhumb_geo_S, type = "angles", units ="degrees", template = "geographics")
Table2$angular_diff<-Table2$S_Vector_Azimuth_No_Shift-Table2$Azimuth_rhumb_geo_S
Table2$angular_diff<-as.numeric(Table2$angular_diff)
Table2$abs_angular_diff<-abs(Table2$angular_diff)
Table2$S_Vector_Azimuth_No_Shift<-as.numeric(Table2$S_Vector_Azimuth_No_Shift)
Table2$Azimuth_rhumb_geo<-as.numeric(Table2$Azimuth_rhumb_geo_S)
Table2$abs_angular_diff<-as.numeric(Table2$abs_angular_diff)
Table2$angular_diff<-as.numeric(Table2$angular_diff)
Table2$angular_diff<-ifelse(Table2$angular_diff>359,0,Table2$angular_diff)#These mean 0 difference
Table2$angular_diff<-ifelse(Table2$angular_diff>179,0,Table2$angular_diff)#single one, milimeters apart
Table2$abs_angular_diff<-as.numeric(abs(Table2$angular_diff))




Table2$S_Binary_No_Shift<-as.circular(Table2$S_Binary_No_Shift, type = "angles", units ="degrees", template = "geographics")
Table2$Azimuth_rhumb_geo<-as.circular(Table2$Azimuth_rhumb_geo, type = "angles", units ="degrees", template = "geographics")
Table2$angular_diff_binary<-Table2$S_Binary_No_Shift-Table2$Azimuth_rhumb_geo_S
Table2$angular_diff_binary<-as.numeric(Table2$angular_diff_binary)
Table2$abs_angular_diff_binary<-abs(Table2$angular_diff_binary)
Table2$S_Binary_No_Shift<-as.numeric(Table2$S_Binary_No_Shift)
Table2$Azimuth_rhumb_geo<-as.numeric(Table2$Azimuth_rhumb_geo)
Table2$abs_angular_diff_binary<-as.numeric(Table2$abs_angular_diff_binary)
Table2$angular_diff_binary<-as.numeric(Table2$angular_diff_binary)

Table2$angular_diff_binary<-ifelse(Table2$angular_diff>359,0,Table2$angular_diff_binary)#These mean 0 difference
Table2$angular_diff_binary<-ifelse(Table2$angular_diff>179,0,Table2$angular_diff_binary)#single one, milimeters apart
Table2$abs_angular_diff_binary<-as.numeric(abs(Table2$angular_diff_binary))
range(Table2$angular_diff_binary)
Table_South<-Table2 %>%
  group_by(range) %>%
  summarize(n=n(),
  Mean_distance_rel=mean(Distance_diff_rel_abs),
   SD_distance_rel=sd(Distance_diff_rel_abs),
  max_distance_rel=max(Distance_diff_rel_abs),
   Percentile_95_distance_rel = quantile(Distance_diff_rel_abs, 0.95, na.rm = TRUE),
  Mean_distance=mean(abs(Distance_diff)),
   SD_distance=sd(abs(Distance_diff)),
  max_distance=max(abs(Distance_diff)),
   Percentile_95_distance = quantile(abs(Distance_diff), 0.95, na.rm = TRUE),
  Mean_azimuth_diff=mean(abs_angular_diff),
  SD_azimuth_diff=sd(abs_angular_diff),
  max_azimuth_diff=max(abs_angular_diff),
     Percentile_95_azimuth = quantile(abs_angular_diff, 0.95, na.rm = TRUE),
   Mean_azimuth_diff_binary=mean(abs_angular_diff_binary),
  SD_azimuth_diff_binary=sd(abs_angular_diff_binary),
  max_azimuth_diff_binary=max(abs_angular_diff_binary),
     Percentile_95_azimuth_binary = quantile(abs_angular_diff_binary, 0.95, na.rm = TRUE),)



```