#####################################################
# Marco FW Gauger
# 2021/02/23
#
# This is the script to obtain fig 3 of the manuscript submitted to PeerJ Gauger et al 2021 Diel influences on bottlenose dolphin acoustic detection in a coastal lagoon in the southwestern Gulf of California
#
# personel computer Intel Core I7-7500 U, @2.6 GHz, RAM 8 GB, Windows 10 
# start 09:20
# end 09:24
#####################################################

ifelse("NPCirc" %in% rownames(installed.packages()) == FALSE, 
       install.packages("NPCirc", dependencies = T),
       ("NPCirc is already installed"))
ifelse("chron" %in% rownames(installed.packages()) == FALSE, 
       install.packages("chron", dependencies = T),
       ("chron is already installed"))
ifelse("circular" %in% rownames(installed.packages()) == FALSE, 
       install.packages("circular", dependencies = T),
       ("circular is already installed"))
ifelse("CircMLE" %in% rownames(installed.packages()) == FALSE, 
       install.packages("CircMLE", dependencies = T),
       ("CircMLE is already installed"))
ifelse("lunar" %in% rownames(installed.packages()) == FALSE, 
       install.packages("lunar", dependencies = T),
       ("lunar is already installed"))

require(circular)
require(NPCirc)
require(CircMLE)


analysed_data_hour <- read.csv("Supplemental Data S1.csv", sep=";")
str(analysed_data_hour)

analysed_data_hour$fcampaign <- as.factor(analysed_data_hour$month_campaign)

analysed_data_hour$h2= chron::hours(as.POSIXlt(paste0(analysed_data_hour$datehour,":00"), format="%d.%m.%Y %H:%M:%S", tz="MST"))

analysed_data_hour$eventhour <- 
  circular(analysed_data_hour$h2%%24,
           # convert to 24 hrs
           units="hours", template="clock24")
is.circular(analysed_data_hour$eventhour)

analysed_data_hour$eventmoon <- 
  circular(lunar::lunar.phase(as.Date(analysed_data_hour$date, format="%d.%m.%Y"), shift= analysed_data_hour$h), template="none", rotation = "clock")
is.circular(analysed_data_hour$eventmoon)


write.csv(analysed_data_hour[analysed_data_hour$dp10m_ELAP1 > 0,], file="Suppl_fig4.csv", sep = ",")
str(analysed_data_hour)

analysed_data_hour$eventhour <- 
  circular(analysed_data_hour$h2%%24,
           # convert to 24 hrs
           units="hours", template="clock24")
is.circular(analysed_data_hour$eventhour)
#CircMLE::HR_test(fig4_data$eventhour, iter = 999)
# hour is in utc
# h is in local time MST

analysed_data_hour$eventmoon <- 
  circular(lunar::lunar.phase(as.Date(analysed_data_hour$date, format="%d.%m.%Y"), shift= analysed_data_hour$h2)/2%%pi, 
           template="geographics")
is.circular(analysed_data_hour$eventmoon)

CircMLE::HR_test(analysed_data_hour$eventmoon[analysed_data_hour$dp10m_ELAP1 > 0], iter = 999)
CircMLE::HR_test(analysed_data_hour$eventhour[analysed_data_hour$dp10m_ELAP1 > 0], iter = 999)





analysed_data_hour$eventmoon <- 
  circular(lunar::lunar.phase(as.Date(analysed_data_hour$date, format="%d.%m.%Y"), shift= analysed_data_hour$h), template="none", rotation = "clock")

sizer_hour_dp10m <- circsizer.density(analysed_data_hour$eventhour[analysed_data_hour$dp10m_ELAP1 > 0], bws=seq(0,50,length=12))
sizer_moon_dp10m <- circsizer.density(analysed_data_hour$eventmoon[analysed_data_hour$dp10m_ELAP1 > 0], bws=seq(0,50,length=12))

dev.off()
png("Figure 3.png", width = 150, height = 80, units = "mm", pointsize = 12, res=300)
par(mfrow=c(1,2), oma=c(0,0,0,0), mar=c(0,0,0,0))
#### ######
circsizer.map(sizer_hour_dp10m,type=2,zero=pi/2,clockwise=TRUE,raw.data=F, labels=rep("",8), label.pos=seq(0,7/4*pi,pi/4))
mtext(side=3, line=-1, adj=-.15, "A")
axis.circular(at=circular(seq(pi/2,2*pi,pi/4)), 
              labels=paste(c(0,21,18,15,12,9,6),"h"), tcl.text =-.65)
###
circsizer.map(sizer_moon_dp10m, type=3, zero=pi, clockwise=T, raw.data=F, labels=c("","","",""), label.pos=seq(0,7/4*pi,pi/2))
#circsizer.map(sizer_moon_dp10m, type=3, zero=pi, clockwise=TRUE,raw.data=F)
par(xpd=T)
#lines(circular(seq(0,pi*2,.1)), rep(.7,length(seq(0,pi*2,.1))), lwd=25, offset=1.1, col="white")
axis.circular(at=circular(seq(pi/2,2*pi,pi/2)), 
              labels=c("NM","Wan","FM","Wax"), tcl.text =-.75)
par(xpd=F)
mtext(side=3, line=-1, adj=-.15, "B")
#### ######
dev.off()