R/TSsphere.R

In AustralianAntarcticDivision/ZooScatR: ZooScatR - Acoustic Scattering of zooplankton using DWBA

#'Calculate the acoustic wave length lambda
#' @param c Soundspeed in m/s
#' @param f frequency in Hz (s^-1)
#' @examples
#' lambda(c=1500,f=200000)
#' @export
lambda <- function(c,f) {c/f}

#'Calculate the acoustic wavenumber k
#' @param c Soundspeed in m/s
#' @param f frequency in Hz (s^-1)
#' @examples
#' k(c=1500,f=200000)
#' @export
k <- function(c,f){2 * pi / lambda(c,f)}

#' Spherical bessel function first kind
#' @param x numeric data vector
#' @param order order of the bessel funciton
#' @examples
#' Sbessj(x=1:10,order=1)
#' @export
Sbessj <- function(x,order){
  sqrt(pi/2) * 1/sqrt(x) * besselJ(x,(order + (1/2)))
}

#' Spherical bessel function second kind (Neumann)
#' @param x numeric data vector
#' @param order order of the bessel funciton
#' @examples
#' Sneum(x=1:10,order=1)
#' @export
Sneum <- function(x,order){
  sqrt(pi/2) * 1/sqrt(x) * besselY(x,(order + 1/2))
}

#' Derivative of the Spherical bessel function second kind
#' @param x numeric data vector
#' @param order order of the bessel funciton
#' @examples
#' Sneum_p(x=1:10,order=1)
#' Sneum_p(x=1:10,order=0)
#' @export
Sneum_p <- function(x,order){
  sn <- Sneum(x,order)

  if(order==0){
    sp1 <- sqrt(pi/(2*x))*besselY(x, order+3/2)
    return((order/x) * sn - sp1)
  }else{
    sm1 <- sqrt(pi/(2*x))*besselY(x, order-1/2)
    return(sm1 - ((order+1)/x) * sn)
  }
}

#' Derivative of the Spherical bessel function first kind
#' @param x numeric data vector
#' @param order order of the bessel funciton
#' @examples
#' Sbessj_p(x=1:10,order=1)
#' Sbessj_p(x=1:10,order=0)
#' @export
Sbessj_p <- function(x,order){
  bn <- Sbessj(x,order)
  bn <- (sqrt(pi/(2*x)))*besselJ(x, order+0.5)

  if(order==0){
    bp1 <- (sqrt(pi/(2*x)))*besselJ(x, order+3/2)
    return((order/x) * bn - bp1)
  }else{
    bm1 <- sqrt(pi/(2*x))*besselJ(x, order-1/2)
    return(bm1 - ((order+1)/x) * bn)
  }
}

#' Hankel function
#' @param k kind of Hankel function must be 1 or 2
#' @param n order of the Hankel funciton
#' @param z numeric data vector
#' @examples
#' besselH(k=1,n=1,z=1:10)
#' @export
besselH <- function(k,n,z){
  if(!(k %in% c(1,2))){
    message("The kind of the Hankel function must be 1 or 2")
    return()
  }
  if(k==1){bh<-besselJ(z,n) + 1i*besselY(z,n)}
  if(k==2){bh<-besselJ(z,n) - 1i*besselY(z,n)}
  return(bh)
}



#' Spherical Hankel function
#' @param z numeric data vector
#' @param order order of the bessel funciton
#' @examples
#' SHankel(x=1:10,order=1)
#' SHankel(x=1:10,order=0)
#' @export
SHankel<- function(z,order){
  (sqrt(pi/(2*x))) * besselH(k=1,z=z,rep(order+1/2,length(x)))
}

#' Derivative of the Spherical Hankel function
#' @param x numeric data vector
#' @param order order of the bessel funciton
#' @examples
#' SHankel_p(x=1:10,order=1)
#' SHankel_p(x=1:10,order=0)
#' @export
SHankel_p <- function(x,order){
  if(order == 0){
    Hankel_n = (sqrt(pi/(2*x)))* besselJ(x,order+1/2) +
      1i*(sqrt(pi/(2*x)))* besselY(x,order+1/2)
    Hankel_nplus1 = (sqrt(pi/(2*x)))* besselJ(x,order+3/2) +
      1i*(sqrt(pi/(2*x)))* besselY(x,order+3/2)
    SphHankel_deriv = (order/x)*Hankel_n - Hankel_nplus1
  }else if(order > 0){
    #Spherical Bessel functions of other order are directly calculated
    #here to try and reduce the number of function calls
    Hankel_n = (sqrt(pi/(2*x)))* besselJ(x,order+1/2)+
      1i*(sqrt(pi/(2*x)))* besselY(x,order+1/2);
    Hankel_ntake1 = (sqrt(pi/(2*x)))* besselJ(x,order-1/2) +
      1i*(sqrt(pi/(2*x)))* besselY(x,order-1/2);

    SphHankel_deriv <- Hankel_ntake1 - ((order+1)/x)*Hankel_n;
  }
  return(SphHankel_deriv)
}

#' Calculate Legendre Polynomial
#' @description Function to calculate the Legendre Polynomial,
#' based on the recurrence relation in Numerical Recipes in C (Eq. 4.6.10) with j = j-1
#' Therefore:
#' P_j = (2j-1) * P_(j-1) - (j-1) * P_(j-2)
#' P_0 = 1; P_1 = x
#'
#' @param degree The degree of the desired polynomial
#' @param coeffs the values at which the p[olynomial will be evaluated
#' @return Legendre polynomials
#' @examples
#' coeffs=2;degree=3
#' LegPoly(degree,coeffs)
#' @export
LegPoly <- function(degree, coeffs){
  lp <- matrix(1,dim(matrix(coeffs))[2], dim(matrix(coeffs))[1])
  lp1 <- coeffs * lp

  if(degree==1){
    return(lp1)
  }

  if (degree > 1){
    for(cntr in 2:degree){
      lp2 <- ((coeffs * (2*cntr-1) * lp1 - (cntr - 1) * lp)/cntr)
      lp <- lp1
      lp1 <- lp2
    }
    return(lp2)
  }else{
    return(matrix(1,dim(matrix(coeffs))[2], dim(matrix(coeffs))[1]))
    }
}




#' Calculate the analytical solution for a weakly scattering sphere
#' @description Scattered pressure from an incident plane wave upon a fluid sphere. Based on:
#' Anderson, V. C. (1950). "Sound scattering from a fluid sphere."
#' The Journal of the Acoustical Society of America, 22(4), 426-431.
#' @param f Frequency in Hz (s^-1)
#' @param r Range in m from center of sphere
#' @param a Radius of sphere
#' @param c Soundspeed in surrounding fluid m/s
#' @param h Soundspeed contrast inside/surrounding fluid
#' @param g Density contrast inside/surrounding fluid
#' @param rho Density of surrounding fluid
#' @param theta Scattering angle
#' @examples
#' f <- 200000
#' r <- 10
#' a <- 0.01
#' c <- 1477.4
#' rho <- 1026.8
#' g <- 1028.9/rho
#' h <- 1480.3/c
#' theta = pi #For backscatter
#' ts<- TS.sphere(f=200000,r,a,c,h,g, rho)
#' #For a range of frequencies
#' ts<- TS.sphere(f=seq(10,400)*1000,r,a,c,h,g, rho, theta=pi)
#' plot(ts~seq(10,400),xlab="Frequency [kHz]", ylab="TS [dB re m2]", type="l",ylim=c(-160,-80))
#' @export
TS.sphere <- function(f,r,a,c,h,g, rho, theta = pi){

  ###Basic calculations

  #External fluid
  k0 <- k(c,f) #Wavenumber surrounding fluid
  kr <- k0 * r #Wavenumber surrounding fluid * range
  ka <- k0 * a #Wavenumber surrounding fluid * radius
  rho_c <- rho * c #rho external fluid * soundspeed

  #Internal fluid
  k1 <- k(h*c,f) #Wavenumber internal fluid
  k1a <- k1 * a  #Wavenumber internal fluid * radius
  #rho1 <- g *rho #density internal fluid

  cosTheta <- cos(theta)

  #Limits
  ConvLim <- 1e-10 #Convergence limit
  MaxCount <- 200 #Max integration points

  m <- 0
  count=1
  #Loop
  done=0
  while(done==0){
    #Bessel functions
    Jm_kr <- Sbessj(kr,m)
    Jm_ka <- Sbessj(ka,m)
    Jm_k1a <- Sbessj(k1a,m)
    Nm_kr <- Sneum(kr,m)
    Nm_ka <- Sneum(ka,m)
    Alpham_k1a <- (2*m+1)*Sbessj_p(k1a,m)
    Alpham_ka <- (2*m+1)*Sbessj_p(ka,m)
    Betam_ka <- (2*m+1)*Sneum_p(ka,m)
    Alpham_kr <- (2*m+1)*Sbessj_p(kr,m)
    Betam_kr <- (2*m+1)*Sneum_p(kr,m)

    #Legendre polynomials - depdendent on Theta
    Pm <- LegPoly(m, cosTheta)

    #Compute Cm with current m - independent of theta
    Cm <- ((Alpham_k1a/Alpham_ka)*(Nm_ka/Jm_k1a) - (Betam_ka/Alpham_ka)*g*h) /
      ((Alpham_k1a/Alpham_ka)*(Jm_ka/Jm_k1a)-g)

    #Use Cm to compute Am
    Am <- -((-1i)^m) * (2*m+1)/(1+1i*Cm)

    #Use Cm to calculate current 'm' value scattered pressure
    p.sm <- -(((-1i)^m)*(2*m+1)/(1+1i*Cm)) * Pm * (Jm_kr + 1i*Nm_kr)
    p.im <- ((-1i)^m) * (2*m+1) * Pm * Jm_kr
    TS.pm <- ((-1)^m) * (2*m+1) / (1+1i*Cm)

    if(m==0){
      TS = TS.pm
    }else{
      TS = TS + TS.pm
    }

    m <- m+1
    count = count + 1

    if(max(abs(p.sm)) < ConvLim){
      done = 1
      message("Converged")
    }
     if(count > MaxCount){
       done = 1
       message("Convergence failed")
     }
  }
  TS <- (2/ka)^2*abs(TS)^2*(pi*a*a)
  TS <- 10*log10(TS/(4*pi))
  return(TS)
}

#' Calculate the analytical solution for a weakly scattering sphere
#' @description Scattered pressure from an incident plane wave upon a fluid sphere. Based on:
#' Jech, J. M., Horne, J. K., Chu, D., Demer, D. A., Francis, D. T., Gorska, N., ... & Reeder, D. B. (2015).
#' "Comparisons among ten models of acoustic backscattering used in aquatic ecosystem research."
#' The Journal of the Acoustical Society of America, 138(6), 3742-3764.
#' modified from:
#' Anderson, V. C. (1950). "Sound scattering from a fluid sphere."
#' The Journal of the Acoustical Society of America, 22(4), 426-431.
#' @param f Frequency in Hz (s^-1)
#' @param r Range in m from center of sphere
#' @param a Radius of sphere
#' @param c Soundspeed in surrounding fluid m/s
#' @param h Soundspeed contrast inside/surrounding fluid
#' @param g Density contrast inside/surrounding fluid
#' @param rho Density of surrounding fluid
#' @examples
#' fs <- as.list(seq(10,400, by=1)*1000) #Frequencies
#' 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
#' TS <- sapply(fs,TS.sphere2,r=r,a=a,c=c,h=h,g=g,rho=rho)
#' plot(fs,TS, type="l", xlab="Frequency [Hz]",ylab="TS [dB re m2]")
#' @export
TS.sphere2 <- function(f,r,a,c,h,g, rho){

  ###Basic calculations

  #External fluid
  k0 <- k(c,f) #Wavenumber surrounding fluid
  kr <- k0 * r #Wavenumber surrounding fluid * range
  ka <- k0 * a #Wavenumber surrounding fluid * radius
  rho_c <- rho * c #rho external fluid * soundspeed

  #Internal fluid
  k1 <- k(h*c,f) #Wavenumber internal fluid
  k1a <- k1 * a  #Wavenumber internal fluid * radius
  #rho1 <- g *rho #density internal fluid


    Cn_fun <- function(n){
      ((Sbessj_p(k1a,n)*Sneum(ka,n))/(Sbessj(k1a,n)*Sbessj_p(ka,n)) -
         g*h*(Sneum_p(ka,n)/Sbessj_p(ka,n))) /
        ((Sbessj_p(k1a,n)*Sbessj(ka,n))/(Sbessj(k1a,n)*Sbessj_p(ka,n)) - g*h)
    }
    n <- seq(0,round(ka+20,0))
    Cn <- apply(matrix(n),1,Cn_fun)

    A <- -1/(1+1i*Cn)
    fbs <- -1i/k0*sum((-1)^n * (2*n +1) * A)
    TS <- 10*log10(abs(fbs)^2)
    return(TS)
}