# R code to reproduce analysis and plots - Pail et al. 2020 (Herpetozoa)

library(circular) # package for plotting and circular analyis 

### prepare some slightly adapted functions (originallly from circular) - i.e. no crosshair in the center of circular plots and a function to double angles and reduce to modulo 360
source("LLCircularFunctions_herpetozoa.R") # only works like this if the functions are in the working directory, otherwise change path to working directory 

#fonts need only to be added once, if already loaded, no need to run every time 
library(extrafont) # in order to use Arial for figures 
font_import() # run this line if fonts are not yet imported 
loadfonts(device="pdf") 

library(shape) # in order to use different kind of arrows for the plots 

# load data - data files need to be placed in the working directory 
FirstDay<-read.csv("data_Firstday.csv")
SecondDay<-read.csv("data_Secondday.csv")

# calculate directions relative to home and d-axis
relhome_1st <- FirstDay$home-FirstDay$direction
FirstDay$relhome <-ifelse(((relhome_1st+360)>(360)),(relhome_1st),(relhome_1st+360))

relhome_2nd <- SecondDay$home-SecondDay$direction
SecondDay$relhome <-ifelse(((relhome_2nd+360)>(360)),(relhome_2nd),(relhome_2nd+360))

relD_1st <- FirstDay$daxis-FirstDay$direction
FirstDay$reldaxis <-ifelse(((relD_1st+360)>(360)),(relD_1st),(relD_1st+360))

relD_2nd <- SecondDay$daxis-SecondDay$direction
SecondDay$reldaxis <-ifelse(((relD_2nd+360)>(360)),(relD_2nd),(relD_2nd+360))

FirstDay$reldaxis_doubled <-DoublingAnglesDeg(FirstDay$reldaxis)
SecondDay$reldaxis_doubled <-DoublingAnglesDeg(SecondDay$reldaxis)

FirstDay$daxis_doubled <-DoublingAnglesDeg(FirstDay$daxis)
SecondDay$daxis_doubled <-DoublingAnglesDeg(SecondDay$daxis)

#make radians values out of degrees 
FirstDay$relhome <-rad(FirstDay$relhome)
SecondDay$relhome <- rad(SecondDay$relhome)
FirstDay$reldaxis <-rad(FirstDay$reldaxis)
SecondDay$reldaxis <-rad(SecondDay$reldaxis)
FirstDay$direction <-rad(FirstDay$direction)
SecondDay$direction <-rad(SecondDay$direction)
FirstDay$reldaxis_doubled <-rad(FirstDay$reldaxis_doubled)
SecondDay$reldaxis_doubled <-rad(SecondDay$reldaxis_doubled)
FirstDay$daxis<-rad(FirstDay$daxis)
SecondDay$daxis<-rad(SecondDay$daxis)

# add lunar phase for MANOVA approach
#First create posixct object
library(lubridate)
library(suncalc)
FirstDay$date <- as.POSIXct(strptime(as.character(FirstDay$date),"%Y-%m-%d", tz="GMT")) 
SecondDay$date <- as.POSIXct(strptime(as.character(SecondDay$date),"%Y-%m-%d", tz="GMT"))

moon <- getMoonIllumination(FirstDay$date,keep = "angle")
FirstDay$moon<- moon[,2]
moon <- getMoonIllumination(SecondDay$date,keep = "angle")
SecondDay$moon <-moon[,2]

FirstDay$cos_moon<-cos(FirstDay$moon); FirstDay$sin_moon<-sin(FirstDay$moon)
SecondDay$cos_moon<-cos(SecondDay$moon); SecondDay$sin_moon<-sin(SecondDay$moon)

# Split data based on sites 
SecondDay_1 <- SecondDay[SecondDay$location==1,]
FirstDay_1 <- FirstDay[FirstDay$location==1,]

SecondDay_2 <- SecondDay[SecondDay$location==2,]
FirstDay_2 <- FirstDay[FirstDay$location==2,]

##################################################
###########MANOVA##########
########################################
library(car) #loads Manova() 
library(qtlmt)# loads mStep()

#First day just direction - add location as a factor 
sin_1st_Dir <- sin(FirstDay$direction); cos_1st_Dir <- cos(FirstDay$direction)
sincos_1st_Dir <-cbind(cos_1st_Dir,sin_1st_Dir)
LinModel_1st_Dir<-lm(sincos_1st_Dir ~ cloud+temp+cos_moon+sin_moon+as.factor(location),
                     data= FirstDay)
LinModel_1st_Dir_DMan=mStep(LinModel_1st_Dir)
LinModel_1st = Manova(LinModel_1st_Dir_DMan)

LinModel_1st_out <- capture.output(LinModel_1st)
write.csv(LinModel_1st_out,"Results_MANOVA_FirstDay.csv") ## save raw table 

#effect size _ added to the raw table by hand 
library(heplots) # loads etasq()
EffectSize_LinModel_1st <- etasq(Manova(LinModel_1st_Dir_DMan))

#Second day just direction
sin_2nd_Dir <- sin(SecondDay$direction); cos_2nd_Dir <- cos(SecondDay$direction)
sincos_2nd_Dir <-cbind(cos_2nd_Dir,sin_2nd_Dir)
LinModel_2nd_Dir<-lm(sincos_2nd_Dir ~ cloud*temp*cos_moon*sin_moon*as.factor(location),
                       data= SecondDay)
LinModel_2nd_Dir_DMan=mStep(LinModel_2nd_Dir)
LinModel_2nd = Manova(LinModel_2nd_Dir_DMan)

LinModel_2nd_out <- capture.output(LinModel_2nd)
write.csv(LinModel_2nd_out,"Results_MANOVA_SecondDay.csv") ## save raw table

#effect size
EffectSize_LinModel_2nd <- etasq(Manova(LinModel_2nd_Dir_DMan)) # added to table by hand 

#################################################################################
###################circular plots ###############################################
#################################################################################
#predfine colors for all plots 
H ="#FF8000";site = "orange4";pcolor = "darkorchid3";testday1 = "#99004C"
testday2 = "#009900";D = "#0000FF";cexp = 0.8

#########################################################################
################### Figure 2
#########################################################################

pdf("Figure 2.pdf", width=4.98, height=3)
par(mar=c(0.5,0,0,0),mfrow = c(1, 2), family = "Arial")

#define circular plotting parameters
sep=0.12;shrink=1.7;cexD = 0.8;start.sep = 0.08; cexsite=0.9

#First day location 1 - d-axis relative to home (292 is home)
ori=circular(FirstDay_1$daxis); rh=rho.circular(ori)
me=circular(mean.circular(ori)); p=rayleigh.test(ori); p=p$p.value

plot.circular(ori,stack = T,axes=FALSE,template = "geographics", units = "degrees",start.sep = start.sep,
              sep=sep,shrink=shrink, cex = cexD,pch=16, tol =0,col=D )
axis.circularLL(at=circular(rad(c(90))), labels=c("N"),tcl.text=0.19,tick=TRUE,cex = 0.8)

arrows.circular(x=me,y=rh,lwd=2.5,length=0.06,zero = rad(90),code = 2,rotation ="clock") 

xs=sin(rad(292)); ys=cos(rad(292));x_ys=abs(ys/xs)#calculating x and y of homeward direction
Arrows(0,0,x1=(xs+0.1),y1=(ys-(0.1*x_ys)),col=H,arr.type = "triangle",
       arr.width=0.3,arr.length=0.2,arr.lwd=1,segment = FALSE)
text((xs+0.33),(ys-(0.33*x_ys)),"H",cex = 1.2,col=H,font=2)

text(-0.35,-0.5,"Site 1",cex = cexsite,col=site,font=2)

text(-1.2,1.2,"A",font=2,cex=1.5) 

#First day location 2 - d-axis relative to north (45 is home) 
ori=circular(FirstDay_2$daxis); rh=rho.circular(ori)
me=circular(mean.circular(ori)); p=rayleigh.test(ori); p=p$p.value

plot.circular(ori,stack = T,axes=F,template = "geographics", units = "degrees",start.sep = start.sep,
              sep=sep,shrink=shrink, cex = cexD,pch=16, tol =0,col=D)
axis.circularLL(at=circular(rad(c(90))), labels=c("N"),tcl.text=0.19,tick=TRUE,cex = 0.8)

arrows.circular(x=me,y=rh,lwd=2.5,length=0.06,zero = rad(90),code = 2,rotation ="clock")

Arrows(0,0,x1=0.63,y1=0.63,col=H,arr.type = "triangle",
       arr.width=0.3,arr.length=0.2,arr.lwd=1,segment = FALSE)
text(0.45,0.45,"H",cex = 1.2,col=H,font=2)


text(-0.35,-0.5,"Site 2",cex = cexsite,col=site,font=2)
text(-1.2,1.2,"B",font=2,cex=1.5) 

dev.off()

##############################################################################
######################### Figure 3
##############################################################################

pdf("Figure 3.pdf", width=4.98, height=4.5)
par(mar=c(0.5,0,0,0),mfrow = c(2, 2), family = "Arial")

#define circular plotting parameters
sep=0.1;shrink=1.4;cexD = 1;start.sep = 0.08;cexsite=0.9 

#First Day - Site 1 - just Direction
ori=circular(FirstDay_1$direction);rh=rho.circular(ori);me=circular(mean.circular(ori))
p=rayleigh.test(ori,mu = rad(292));p=p$p.value

plot.circular(ori,stack = T,axes=FALSE,template = "geographics", units = "degrees",start.sep = start.sep ,
              sep=sep,shrink=shrink, cex = cexD,pch=16, tol =0,col=testday1)
axis.circularLL(at=circular(rad(c(90))), labels=c("N"),tcl.text=0.19,tick=TRUE,cex = 0.8)
arrows.circular(x=me,y=rh,lwd=3,length=0.07,zero = rad(90),code = 2,rotation ="clock")

text2 = paste("p",ifelse(round(p,2)<0.01,"< 0.01",paste("=",round(p,2))))
text(-0.4,-0.4,text2,cex = cexp,col=pcolor)

text(-1.2,1.2,"A",font=2,cex=1.5) 

xs=sin(rad(292)); ys=cos(rad(292));x_ys=abs(ys/xs)#calculating x and y of homeward direction
Arrows(0,0,x1=(xs+0.1),y1=(ys-(0.1*x_ys)),col=H,arr.type = "triangle",
       arr.width=0.3,arr.length=0.2,arr.lwd=1,segment = FALSE)
text((xs+0.28),(ys-(0.28*x_ys)),"H",cex = 1.2,col=H,font=2)

text(0,.5,"Site 1",cex = cexsite,col=site,font=2)


#First Day - Site 2 - just Direction
ori=circular(FirstDay_2$direction);rh=rho.circular(ori)
me=circular(mean.circular(ori));p=rayleigh.test(ori,mu = rad(45)); p=p$p.value

plot.circular(ori,stack = T,axes=FALSE,template = "geographics", units = "degrees",start.sep = start.sep,
              sep=sep,shrink=shrink, cex = cexD,pch=16, tol =0,col=testday1)
axis.circularLL(at=circular(rad(c(90))), labels=c("N"),tcl.text=0.19,tick=TRUE,cex = 0.8)
arrows.circular(x=me,y=rh,lwd=3,length=0.07,zero = rad(90),code = 2,rotation ="clock") 

text2 = paste("p",ifelse(round(p,2)<0.01,"< 0.01",paste("=",round(p,2))))
text(-0.4,-0.4,text2,cex = cexp,col=pcolor)

text(-1.2,1.2,"B",font=2,cex=1.5) # add more descirbtions 

Arrows(0,0,x1=0.63,y1=0.63,col=H,arr.type = "triangle",
       arr.width=0.3,arr.length=0.2,arr.lwd=1,segment = FALSE)
text(0.48,0.48,"H",cex = 1.2,col=H,font=2)

text(0,0.5,"Site 2",cex = cexsite,col=site,font=2)

#Second Day - Site 1 - just Direction
ori=circular(SecondDay_1$direction); rh=rho.circular(ori)
me=circular(mean.circular(ori)); p=rayleigh.test(ori,mu = rad(292)); p=p$p.value

plot.circular(ori,stack = T,axes=FALSE,template = "geographics", units = "degrees",start.sep = start.sep,
              sep=sep,shrink=shrink, cex = cexD,pch=16, tol =0,col=testday2)
axis.circularLL(at=circular(rad(c(90))), labels=c("N"),tcl.text=0.19,tick=TRUE,cex = 0.8)
arrows.circular(x=me,y=rh,lwd=3,length=0.07,zero = rad(90),code = 2,rotation ="clock")

text2 = paste("p",ifelse(round(p,2)<0.01,"< 0.01",paste("=",round(p,2))))
text(-0.4,-0.4,text2,cex = cexp,col=pcolor)

text(-1.2,1.2,"C",font=2,cex=1.5) 

xs=sin(rad(292)); ys=cos(rad(292));x_ys=abs(ys/xs) #calculating x and y of homeward direction
Arrows(0,0,x1=(xs+0.1),y1=(ys-(0.1*x_ys)),col=H,arr.type = "triangle",
       arr.width=0.3,arr.length=0.2,arr.lwd=1,segment = FALSE)
text((xs+0.28),(ys-(0.28*x_ys)),"H",cex = 1.2,col=H,font=2)

text(0,0.5,"Site 1",cex = cexsite,col=site,font=2)


#Second Day - Site 2 - just Direction
ori=circular(SecondDay_2$direction); rh=rho.circular(ori)
me=circular(mean.circular(ori)); p=rayleigh.test(ori,mu = rad(45)); p=p$p.value

plot.circular(ori,stack = T,axes=FALSE,template = "geographics", units = "degrees",start.sep = start.sep,
              sep=sep,shrink=shrink, cex = cexD,pch=16, tol =0,col=testday2)
axis.circularLL(at=circular(rad(c(90))), labels=c("N"),tcl.text=0.19,tick=TRUE,cex = 0.8)
arrows.circular(x=me,y=rh,lwd=3,length=0.07,zero = rad(90),code = 2,rotation ="clock") 

text2 = paste("p",ifelse(round(p,2)<0.01,"< 0.01",paste("=",round(p,2))))
text(-0.4,-0.4,text2,cex = cexp,col=pcolor)

text(-1.2,1.2,"D",font=2,cex=1.5)

Arrows(0,0,x1=0.63,y1=0.63,col=H,arr.type = "triangle",
       arr.width=0.3,arr.length=0.2,arr.lwd=1,segment = FALSE)
text(0.48,0.48,"H",cex = 1.2,col=H,font=2)

text(0,0.5,"Site 2",cex = cexsite,col=site,font=2)

mtext("Test day 1",side = 3, line = -12,cex = 0.8,col=testday1,outer = TRUE,font=2)
mtext("Test day 2",side = 3, line = -25,cex = 0.8,col=testday2,outer = TRUE,font=2)
dev.off()

######################################################################
##################Figure 4
######################################################################

pdf("Figure 4.pdf", width=4.98, height=4.5)
par(mar=c(0.5,0,0,0),mfrow = c(2, 2), family = "Arial")

#define circular plotting parameters
sep=0.1;shrink=1.4;cexD = 1;start.sep = 0.08; cexsite = 0.9

#First Day - Site 1  - relative to d-axis
ori=circular(FirstDay_1$reldaxis); ori2=circular(FirstDay_1$reldaxis_doubled)
rh=rho.circular(ori); rh2=rho.circular(ori2)
me=circular(mean.circular(ori)); me2=(circular(mean.circular(ori2)))/2; me2=ifelse(((me2+2*pi)>(2*pi)),(me2),(me2+2*pi))
p=rayleigh.test(ori,mu=0); p=p$p.value; p2=rayleigh.test(ori2,mu=0); p2=p2$p.value

plot.circular(ori,stack = T,axes=FALSE,template = "geographics", units = "degrees",start.sep = start.sep,
              sep=sep,shrink=shrink, cex = cexD,pch=16, tol =0,col=testday1)


Arrows(-0.03,0.85,x1=-0.03,y1=0.88,col=D,
       arr.type = "triangle",arr.width=0.3,arr.length=0.2,arr.lwd=1,segment = FALSE)
text(-0.03,0.68,"D",cex = 1.2,col=D,font=2)

arrows.circular(x=me2,y=rh2,lwd=2.5,length=0.06,zero = rad(90),code = 2,rotation ="clock") 
arrows.circular(x=me2-pi,y=rh2,lwd=2.5,length=0.06,zero = rad(90),code = 2,rotation ="clock") 

text2 = paste("p",ifelse(round(p2,2)<0.01,"< 0.01",paste("=",round(p2,2))))
text(-0.6,0,text2,cex = cexp,col = pcolor)

plot_bootstrap_CI_doubled_inRad_outDeg(ori2)

text(-1.2,1.2,"A",font=2,cex=1.5)
text(-0.6,0.4,"Site 1",cex = cexsite,col=site,font=2)


#First Day - Site 2  - relative to d-axis
ori=circular(FirstDay_2$reldaxis); ori2=circular(FirstDay_2$reldaxis_doubled)
rh=rho.circular(ori); rh2=rho.circular(ori2); me=circular(mean.circular(ori))
me2=(circular(mean.circular(ori2)))/2; me2=ifelse(((me2+2*pi)>(2*pi)),(me2),(me2+2*pi))
p=rayleigh.test(ori,mu=0); p=p$p.value; p2=rayleigh.test(ori2,mu=0); p2=p2$p.value

plot.circular(ori,stack = T,axes=FALSE,template = "geographics", units = "degrees",start.sep = start.sep,
              sep=sep,shrink=shrink, cex = cexD,pch=16, tol =0,col=testday1)

Arrows(0.03,0.85,x1=0.03,y1=0.88,col=D,arr.type = "triangle",
       arr.width=0.3,arr.length=0.2,arr.lwd=1,segment = FALSE)
text(0.03,0.68,"D",cex = 1.2,col=D,font=2)

arrows.circular(x=me2,y=rh2,lwd=3,length=0.07,zero = rad(90),code = 2,rotation ="clock") 
arrows.circular(x=(me2-pi),y=rh2,lwd=3,length=0.07,zero = rad(90),code = 2,rotation ="clock")

text2 = paste("p",ifelse(round(p2,2)<0.01,"< 0.01",paste("=",round(p2,2))))
text(-0.6,0,text2,cex = cexp,col=pcolor)

plot_bootstrap_CI_doubled_inRad_outDeg(ori2)
text(-1.2,1.2,"B",font=2,cex=1.5)

text(-0.6,0.4,"Site 2",cex = cexsite,col=site,font=2)


#Second Day - Site 1  - relative to d-axis
ori=circular(SecondDay_1$reldaxis); ori2=circular(SecondDay_1$reldaxis_doubled)
rh=rho.circular(ori); rh2=rho.circular(ori2); me=circular(mean.circular(ori))
me2=(circular(mean.circular(ori2)))/2; me2=ifelse(((me2+2*pi)>(2*pi)),(me2),(me2+2*pi))
p=rayleigh.test(ori,mu=0); p=p$p.value; p2=rayleigh.test(ori2,mu=0); p2=p2$p.value

plot.circular(ori,stack = T,axes=FALSE,template = "geographics", units = "degrees",start.sep = start.sep,
              sep=sep,shrink=shrink, cex = cexD,pch=16, tol =0,col=testday2)

Arrows(-0.03,0.85,x1=-0.03,y1=0.88,col=D,
       arr.type = "triangle",arr.width=0.3,arr.length=0.2,arr.lwd=1,segment = FALSE)
text(-0.03,0.68,"D",cex = 1.2,col=D,font=2)

arrows.circular(x=me,y=rh,lwd=2.5,length=0.06,zero = rad(90),code = 2,rotation ="clock") 

plot_bootstrap_CI_inRad_outDeg(ori)

text2 = paste("p",ifelse(round(p,2)<0.01,"< 0.01",paste("=",round(p,2))))
text(-0.6,0,text2,cex = cexp,col = pcolor)

text(-1.2,1.2,"C",font=2,cex=1.5) 

text(-0.6,0.4,"Site 1",cex = cexsite,col=site,font=2)

#Second Day - Site 2  - relative to d-axis
ori=circular(SecondDay_2$reldaxis); ori2=circular(SecondDay_2$reldaxis_doubled)
rh=rho.circular(ori); rh2=rho.circular(ori2);me=circular(mean.circular(ori))
me2=(circular(mean.circular(ori2)))/2; me2=ifelse(((me2+2*pi)>(2*pi)),(me2),(me2+2*pi))
p=rayleigh.test(ori,mu=0); p=p$p.value;p2=rayleigh.test(ori2,mu=0); p2=p2$p.value

plot.circular(ori,stack = T,axes=FALSE,template = "geographics", units = "degrees",start.sep = start.sep,
              sep=sep,shrink=shrink, cex = cexD,pch=16, tol =0,col=testday2)

Arrows(0.03,0.85,x1=0.03,y1=0.88,col=D,arr.type = "triangle",
       arr.width=0.3,arr.length=0.2,arr.lwd=1,segment = FALSE)
text(0.03,0.68,"D",cex = 1.2,col=D,font=2)

arrows.circular(x=me,y=rh,lwd=3,length=0.07,zero = rad(90),code = 2,rotation ="clock") 

text2 = paste("p",ifelse(round(p2,2)<0.01,"< 0.01",paste("=",round(p,2))))
text(-0.6,0,text2,cex = cexp,col=pcolor)

text(-1.2,1.2,"D",font=2,cex=1.5) # add more descirbtions 
text(-0.6,0.4,"Site 2",cex = cexsite,col=site,font=2)

mtext("Test day 1",side = 3, line = -12,cex = 0.8,col=testday1,outer = TRUE,font=2)
mtext("Test day 2",side = 3, line = -25,cex = 0.8,col=testday2,outer = TRUE,font=2)

dev.off()