This document compares the results of the the DWBA implementation in ZooScatR to an anlytical scattering model of a sphere [@anderson_sound_1950], following routines described in Jech et al. [-@jech_comparisons_2015].
This vignette requires the ZooScatR package for the implementation of the DWBA model, ggplot2 for plotting, viridis for nicer colors.
library(ZooScatR) library(ggplot2) library(viridis)
For the purpose of this example we will build the general anlytical model for a scattered wave by a sphere [@anderson_sound_1950]:
$$p_{scat}=p_0\sum^{\infty}_{n=0}A_n(2n+1)P_n (cos\theta)h^1_n (k_w r)$$
with
For backscatter $\theta=\pi$; therefore $P_n$ becomes $P_n(-1)=(-1)^n$.
In the farfield
$$
\begin{aligned}
k_wr \gg 1, \
h^1_n(k_wr) \xrightarrow[k_wr\gg1]{} \frac{[(-i)^{n+1}]}{(k_wr)}e^{ik_wr}
\end{aligned}
$$
the acoustic backscatter pressure $p_{bs}$ can be expressed as a function of backscattering amplitude ($f_{bs}$, incident pressure, phase and range:
$$p_{bs}=p_0 \frac{e^{ik_w r}}{r}f_{bs}$$
with the simplified $f_bs$ becoming:
$$f_{bs}=-\frac{i}{k_w}\sum^{\infty}_{n=0}(-1)^n(2n+1)A_n $$
following Jech et al. [-@jech_comparisons_2015], $n$ was limited to the nearest integer of $||k_wa||+20$, providing a precision of 0.1 dB.
From $f_bs$, target strength $TS$ in $dB \text{ re }m^2$ can be computed as [@anderson_sound_1950]:
$$TS = 10log_{10}(|f_bs|^2)$$
For a fluid filled sphere, the coefficient, the coefficient $A_n$ can be computed as:
$$A_n=\frac{-1}{1+iC_n}$$
and
$$C_n = \frac{\frac{[j'_n(k_1a)y_n(ka)]}{[j_n(k_1a)j'_n(ka)]} - gh\frac{y'_n(ka)}{j'_n(ka)}}{\frac{[j'_n(k_1a)j_n(ka)]}{[j_n(k_1a)j'_n(ka)]}-gh}$$
with:
Letting $gh \rightarrow \infty$ results in a rigid sphere and:
$$A_n=-\frac{j'_n(k_wa)}{h'_n(k_wa)}$$
Similarly, letting $gh \rightarrow 0$ results in a pressure release sphere and:
$$A_n=-\frac{j_n(k_wa)}{h_n(k_wa)}$$
These functions have been implemented within ZooScatR and can be called using:
library(ZooScatR) r <- 10 #range a <- 0.01 # radius c <- 1477.4 #soundspeed surrounding fluid rho <- 1026.8 #density surrounding fluid g <- 1028.9/rho #density contrast h <- 1480.3/c #soundspeed contrast f <- 200 * 1000 #Frequency in Hz TS.sphere2(f=f,r=r,a=a,c=c,h=h,g=g,rho=rho)
Note: Results involving Bessel functions have to be used with a bit of caution. Given the diffrent ways available to compute the Bessel functions slight variations occur for example between results obtained in Matlab and R. Generally the fluctuations are relatively small (within a fraction of 1%), but they are noticeable.
The following settings were used for a weakly scattering sphere, as defined in [@jech_comparisons_2015]:
r <- 10 #range a <- 0.01 # radius c <- 1477.4 #soundspeed surrounding fluid rho <- 1026.8 #density surrounding fluid g <- 1028.9/rho #density contrast h <- 1480.3/c #soundspeed contrast #Frequency range fmin=12 fmax=400 freqs <- seq(fmin,fmax, by=1)
Set the settings inside ZooScatR:
fname <- paste0(system.file(package="ZooScatR"),"/extdata/configs/config_0.dat") #Loacation of the parameters file para = read_para(fname) #Read parameters file #Create list with soundspeed info misc <- list(cw=c) #Set soundspeed in the surrounding fluid #Set the material properties para$phy$g0 <- g #set the density contrast para$phy$h0 <- h #set the soundspeed contrast #set the simulation parameters para$simu$var0<-fmin #Set the minimum frequency para$simu$var1<-fmax #Set the maximum frequency para$simu$n=length(freqs) #Set the number of output points para$simu$ni=2000 #any high number
#Set the shape parameters para$shape$L <- 2*1000*a # 2 * the radius in m to get length in mm para$shape$order <- 2 #Set the tapering order para$shape$L_a <- 2 #Set L/a para$shape$rho_L <- 2000 #any high number to get a high curvature #Build the sphere and plot the result sp <- buildpos(para) #build the shape sp$plot #show the shape plot
Run the analytical model:
fs <- as.list(freqs*1000) #Create a list of frequencies TS.MJ <- sapply(fs,TS.sphere2,r=r,a=a,c=c,h=h,g=g,rho=rho) #Apply frequency list to the analytical model
Run the ZooScatR DWBA:
#Run DWBA based on config file DWBA <- bscat(para=para, misc=misc) #Show the model outcome DWBA$rplot
Build a function to compute and plot the error:
The error is computed following Jech et al. [-@jech_comparisons_2015] as:
$$|\overline{\Delta TS}|=\frac{1}{N}\sum^N_{i=1}|TS_{i. predicted}-TS_{i-benchmark}| $$
with $N$ the number of sample points (frequency points).
#Compare results TS.comp<-function(TS1,TS2,freqs){ require(ggplot2) require(viridis) TScomp <- as.data.frame(cbind(TS = c(TS1,TS2), Freq = rep(freqs,2), Model=rep(c("Analytical","DWBA"),each=length(freqs)))) names(TScomp)<- c("TS","Frequency","Model") TScomp$TS <- as.numeric(as.character(TScomp$TS)) TScomp$Frequency <- as.numeric(as.character(TScomp$Frequency)) err <- sum(abs(TS1-TS2))/length(freqs) TScomp$diff <- rep(abs(TS1-TS2),2) message(paste("Overall error:", round(mean(err),2),"dB")) #ggplot remove backgorund rmb <- theme(axis.line = element_line(colour = "black"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.border = element_blank(), panel.background = element_blank()) print(ggplot(data=TScomp,aes(x=Frequency,y=TS, group=Model,col=Model))+ geom_line(lwd=1)+ scale_color_brewer(palette="Set1")+ theme(legend.position = "top")+ xlab("Frequency [kHz]")+ rmb) print(ggplot(data=TScomp,aes(x=Frequency,y=diff, col=diff))+ geom_line(lwd=1.2)+ scale_color_viridis(name="dB difference")+ theme(legend.position="top")+xlab("Frequency [kHz]")+rmb) return(TScomp) }
Run the error function for the anylitcal and modelled solutions:
TS.compare <- TS.comp(TS.MJ,DWBA$y,freqs)
For a weak scattering sphere, the error is mostly below 1%. Only within the null locations, the error increases, while th elocations of the nulls are good.
Running the model for material properties of krill as reported in Calie et al. [-@calise_sensitivity_2011] in this example we are looking at a hypothetical perfectly round krill with a radius of 3 mm:
a <- 0.01 g <- 1.0357 #density contrast h <- 1.0279 #soundspeed contrast c <- 1456 #Set the material properties in ZooScatR para$shape$L <- 2*1000*a # 2 * the radius in m to get length in mm para$phy$g0 <- g #set the density contrast para$phy$h0 <- h #set the soundspeed contrast misc$cw <- c #Run the analytical model TS.MJ <- sapply(fs,TS.sphere2,r=r,a=a,c=c,h=h,g=g,rho=rho) #Apply frequency list to the analytical model #Run the ZooScatR DWBA: DWBA <- bscat(para=para, misc=misc)
Compare the results:
TS.compare <- TS.comp(TS.MJ,DWBA$y,freqs)
Compare for TS > -80:
th <- - 80 tmp1 <- TS.MJ tmp1[which(tmp1 < th)] <- th tmp2 <- DWBA$y tmp2[which(tmp2 < th)] <- th Ts.compare2 <- TS.comp(tmp1,tmp2,freqs)
The ZooScatR DWBA model performs well for the material preperties of krill, when compared to the analytical model.
Running the model for very strong scatterer like a 38.1 mm tungsten-carbide sphere:
a <- 0.0381 g <- 14900/rho #density contrast h <- 6853/c #soundspeed contrast #Set the material properties in ZooScatR para$shape$L <- 2*1000*a # 2 * the radius in m to get length in mm para$phy$g0 <- g #set the density contrast para$phy$h0 <- h #set the soundspeed contrast #Run the analytical model TS.MJ <- sapply(fs,TS.sphere2,r=r,a=a,c=c,h=h,g=g,rho=rho) #Apply frequency list to the analytical model #Run the ZooScatR DWBA: DWBA <- bscat(para=para, misc=misc)
Compare the results:
TS.compare <- TS.comp(TS.MJ,DWBA$y,freqs)
The ZooScatR DWBA model holds no validity for such strong scatterers.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.