Nothing
## ---- include = FALSE---------------------------------------------------------
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
## ----setup, eval=FALSE--------------------------------------------------------
# #Stable - Install package from CRAN
# install.packages("bioSNR")
#
# #Unstable - Install package from Github repository
# devtools::install_github("MattyD797/bioSNR")
#
# #Attach package namespace to active libraries in Rstudio
# library(bioSNR)
## ----setupR, eval=TRUE, include=FALSE-----------------------------------------
library(bioSNR)
## ----wenz, echo=FALSE, fig.align = 'center', message=FALSE--------------------
freqBand <- c(0,0)
shipT <- -1
seaState <- -1
wSpeed <- 0
boolR <- T
# This code re-traces the empirical Wenz curves, that describe Oceanic ambient noise as a function of frequency.
freqMax <- 10000
freq <- pracma::logspace(0,6,freqMax)
#Geophysical contribution for Freq<10 Hz
#N_geophys = 107 - 30log10(f)
intervalG <- freq[1:tail(which(freq<10)+1, 1)]
NL_geophys <- 107 - 30 * log10(intervalG)
NL_geophys2 <- cbind(intervalG, NL_geophys)
colnames(NL_geophys2)[1] <- "freq"
## Contribution of ship traffic (5 <= f <= 500 Hz)
# - 1 - 2 low ship traffic
# - 3-4-5 standard ship traffic
# - 6 - 7 heavy ship traffic
# - 8 - 9 intense ship traffic
# N_ship = 75 -40log10(f/30) + 5(n_shiptraffic - 4)
shipLines <- function(trafficIndex, freqI=freq, freqmax = 10000){
interval <- freqI[head(which(freqI > 5), 1):tail(which(freqI < 500), 1)]
NL_ship<- 76 - (20*(log10(interval/30))^2) + 5*(trafficIndex-4)
return(NL_ship)
}
shipTraffic <- sapply(c(1:9), shipLines)
shipTraffic2 <- cbind(freq[head(which(freq > 5), 1):tail(which(freq < 500), 1)], shipTraffic)
colnames(shipTraffic2) <- c("freq",
"shipTraffic.1",
"shipTraffic.2",
"shipTraffic.3",
"shipTraffic.4",
"shipTraffic.5",
"shipTraffic.6",
"shipTraffic.7",
"shipTraffic.8",
"shipTraffic.9")
## Wind speed and sea state (80 <= f <= 300 kHz)
# Wind speed in knots - related to sea state as Vkn = 5 * sea state
windSpeed <- function(windSpeed, freqI=freq, freqmax = 10000){
#between 80 - 1000 Hz
# < 1000 Hz
interval1 <- head(which(freqI > 80)-1, 1):(tail(which(freqI < 1000)-1, 1))
NL_wind <- 44 + sqrt(21*windSpeed) + 17 * (3 - log10(freqI[interval1])) * (log10(freqI[interval1])-2)
# >= 1000 Hz
interval2 <- head(which(freqI > 1000)-1, 1):tail(which(freqI < 300000)-1, 1)
NL_wind <- append(NL_wind, 95 + sqrt(21*windSpeed) - 17*log10(freqI[interval2]))
return(NL_wind)
}
if(seaState != -1){
wSpeed = seaState * 5
}
wind <- sapply(c(0, 10, 15, 20, 25, 30, 40, 45), windSpeed)
intervalW <- freq[head(which(freq > 80)-1, 1):tail(which(freq < 300000)-1, 1)]
wind2 <- cbind(intervalW, wind)
colnames(wind2) <- c("freq", "seaState.0", "seaState.10", "seaState.15", "seaState.20", "seaState.25", "seaState.30", "seaState.40", "seaState.45")
## Thermal noise (f > 50 kHz)
# NLthermal = -75 + 20log10(f)
intervalT <- freq[head(which(freq > 50000), 1):freqMax]
NL_thermal <- -75 + 20 * log10(intervalT)
NL_thermal2 <- cbind(intervalT, NL_thermal)
colnames(NL_thermal2)[1] <- "freq"
freq <- as.data.frame(freq)
colnames(freq) <- "freq"
#Solve for given bandwidth
bandW <- freq[freq>=freqBand[1] & freq <= freqBand[2]]
gpInputF <- vector()
stInputF <- vector()
wInputF <- vector()
tnInputF <- vector()
gpInput <- vector()
stInput <- vector()
wInput <- vector()
tnInput <- vector()
for(f in bandW){
#Only geophysical SL
if(f < 10){
gpInput<- c(gpInput, 107 - 30 * log10(f))
gpInputF <- c(gpInputF, f)
}
if (f >= 5 & f <= 500){
if(shipT < 0 ){
if(f <= 10){
stop("ERROR: shipT is not specified")
}
} else {
stInput <- c(stInput, 76 - (20*(log10(f/30))^2) + 5*(shipT-4))
stInputF <- c(stInputF, f)
}
}
if(f>=80 & f < 1000){
wInput <- c(wInput, 44 + sqrt(21*wSpeed) + 17 * (3 - log10(f)) * (log10(f)-2))
wInputF <- c(wInputF, f)
}
if(f>=1000 & f <= 300000){
wInput <- c(wInput, 95 + sqrt(21*wSpeed) - 17*log10(f))
wInputF <- c(wInputF, f)
}
if(f>50000){
tnInput <- c(tnInput,-75 + 20 * log10(f))
tnInputF <- c(tnInputF, f)
}
}
gpInput2 <- cbind(gpInputF, gpInput)
colnames(gpInput2)[1] <- "freq"
stInput2 <- cbind(stInputF, stInput)
colnames(stInput2)[1] <- "freq"
wInput2 <- cbind(wInputF, wInput)
colnames(wInput2)[1] <- "freq"
tnInput2 <- cbind(tnInputF, tnInput)
colnames(tnInput2)[1] <- "freq"
df_list <- list(freq, NL_geophys2, shipTraffic2, wind2, NL_thermal2, gpInput2, stInput2, wInput2, tnInput2)
dt <- Reduce(function(x, y) merge(x, y, by = "freq", all = TRUE), df_list)
dt <- dplyr::mutate(dt, gpInput = as.numeric(gpInput),
stInput = as.numeric(stInput),
wInput = as.numeric(wInput),
tnInput = as.numeric(tnInput))
suppressWarnings(print(ggplot2::ggplot(dt, ggplot2::aes(freq)) +
ggplot2::geom_line(ggplot2::aes(y=shipTraffic.1))+
ggplot2::geom_line(ggplot2::aes(y=shipTraffic.2))+
ggplot2::geom_line(ggplot2::aes(y=shipTraffic.3))+
ggplot2::geom_line(ggplot2::aes(y=shipTraffic.4))+
ggplot2::geom_line(ggplot2::aes(y=shipTraffic.5))+
ggplot2::geom_line(ggplot2::aes(y=shipTraffic.6))+
ggplot2::geom_line(ggplot2::aes(y=shipTraffic.7))+
ggplot2::geom_line(ggplot2::aes(y=shipTraffic.8))+
ggplot2::geom_line(ggplot2::aes(y=shipTraffic.9))+
ggplot2::scale_x_continuous(trans='log10', ,
breaks = scales::trans_breaks("log10", function(x) 10^x),
labels = scales::trans_format("log10", scales::math_format(10^.x)),
n.breaks = 25)+
ggplot2::geom_line(ggplot2::aes(y=NL_geophys))+
ggplot2::geom_line(ggplot2::aes(y=seaState.0))+
ggplot2::geom_line(ggplot2::aes(y=seaState.10))+
ggplot2::geom_line(ggplot2::aes(y=seaState.15))+
ggplot2::geom_line(ggplot2::aes(y=seaState.20))+
ggplot2::geom_line(ggplot2::aes(y=seaState.25))+
ggplot2::geom_line(ggplot2::aes(y=seaState.30))+
ggplot2::geom_line(ggplot2::aes(y=seaState.40))+
ggplot2::geom_line(ggplot2::aes(y=seaState.45))+
ggplot2::geom_line(ggplot2::aes(y=NL_thermal))+
ggplot2::ylab(as.expression(bquote("Spectrum Level (dB ref 1 \u00B5" ~ Pa^2 ~ ')')))+
ggplot2::xlab("Frequency (Hz)") +
ggplot2::geom_point(ggplot2::aes(y=gpInput))+
ggplot2::geom_point(ggplot2::aes(y=stInput))+
ggplot2::geom_point(ggplot2::aes(y=wInput))+
ggplot2::geom_point(ggplot2::aes(y=tnInput))+
ggplot2::geom_vline(xintercept = freqBand[1])+
ggplot2::geom_vline(xintercept = freqBand[2])+
ggplot2::annotation_logticks(sides = "b")))
## ----ex1, fig.align = 'center', warning=FALSE---------------------------------
specLvlGraph(c(28,33), ship=4,seaState = 1, wSpeed = 10, boolR = T)[[2]]
## ----ex2, fig.align = 'center'------------------------------------------------
#Source Level
SL <- 195
#Noise Level - from above
NL <- 82.9897
#Detection Threshold
DT <- 10
#Transition Range <- depth/2
TR <- 2500
rmax(sl=SL,nl=NL,dt=DT,d=5000, xaxis=10000000)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.