R/bscat.R

Defines functions bscat

Documented in bscat

#' @title Backscattering model
#' @description  Run the DWBA model for a given set of parameters. Visit the vignette for more details (\url{../doc/DwbaCommand.html})
#' @param para [list] A list containing all the model parameters (\url{../doc/parameters.html})
#' @param misc [list] A list containing the soundspeed (cw) of the surrounding fluid (to compute cw visit \code{\link{c_Coppens1981}} or \code{\link{c_Leroy08}} or \code{\link{c_Mackenzie1981}})
#' @param app = FALSE [boolean] function call from shiny interface or command line
#' @param  nang [integer] Number of angular ismulations, this will overwrite increment
#' @param nl [integer] Number of Length simulations, will overwrite increment informaiton
#' @param simOut = TRUE [boolean] If TRUE ysim (results for the simulated orientations) and ysimL (results for the simulated lengths based on the results of the mean simulated angle) as well as ang (simulated orientation angles) and L (simulated lengths) are added to the function output (average is kept as well), which contains all simulated data, if FALSE, only the averaged value is kept
#' @param plotOut = TRUE [boolean] If results plot should be produced or not
#' @author Sven Gastauer
#' @return list with all parameters for DWBA
#' @import ggplot2
#' @export bscat
#' @examples
#' #Get filename of the parameters file
#' fname <- paste0(system.file(package = "ZooScatR"),"/extdata/configs/config_0.dat")
#' #Read in teh parameter
#' para = read_para(fname)
#' #Create list with soundspeed info
#' misc <- list()
#' misc$cw <- 1500
#' #Change some of the settings
#' #Set starting frequency over which to run the model
#' para$simu$var0 <- 38
#' #Set end frequency over which to run te model
#' para$simu$var1 <- 300
#' # run DWBA based on the settings defined in the parameters file
#' res <- bscat(para=para, misc=misc, app=FALSE)
#' #plot the results of the model
#' res$rplot
#'



bscat <- function(para, misc, app=FALSE, nang=NULL, nl=NULL,simOut = TRUE, plotOut=TRUE){
  # if(exists("status")==FALSE){status=list()}
  # status$stop = 0

  # limits of angle variation
  if (para$simu$var_indx == 2){ #If output is vs angles
    if(para$orient$ave_flag == 1){ #if average should be computed
      if(para$orient$PDF == 1){	# uniform PDF
        ang_max=para$orient$ang1 #para$simu$var1+para$orient$PDF_para
        ang_min=para$orient$ang0 #para$simu$var0-para$orient$PDF_para
      }else{								# Gaussian PDF
        ang_max=para$simu$var1+3.1*para$orient$PDF_para
        ang_min=para$simu$var0-3.1*para$orient$PDF_para
      }
      if(ang_min < para$orient$ang0){
        para$orient$ang0 = ang_min
        #h = findobj(gcf,'Tag','EditTextTheta0')
        #set(h,'String',num2str(para$orient$ang0)) ################### To be fixed
      }
      if(ang_max > para$orient$ang1){
        para$orient$ang1=ang_max
        #h=findobj(gcf,'Tag','EditTextTheta1')
        #set(h,'String',num2str(para$orient$ang1))################### To be fixed
      }
      ang = seq(ang_min, ang_max, by = para$orient$dang) # get anular range
    }else{ #if no average for vs angles should be calculated
      ang_min = para$simu$var0
      ang_max = para$simu$var1
      ang=seq(ang_min,ang_max,length=para$simu$n)
    }
  }else{ #if output is not vs angles (but frequency or ka)
      if(para$orient$ave_flag == 1){ #average over angles should be computed
      if(para$orient$PDF == 1){		# uniform PDF
        ang_min = para$orient$ang0
        ang_max = para$orient$ang1
        print(ang_min)
    }else{								# Gaussian PDF
      ang_min = para$orient$angm - 3.1 * para$orient$PDF_para
      ang_max = para$orient$angm + 3.1 * para$orient$PDF_para
    }
        if(is.null(nang)){
          ang = seq(ang_min, ang_max, by = para$orient$dang) # get anular range
        }else{
          ang = seq(ang_min, ang_max, length.out = nang) # get anular range
        }
    }else{ #no average should be computed
      ang_min=para$orient$angm
      ang_max=para$orient$angm
      ang=ang_min
    }
}
# limits of ka variation
a = para$shape$L / para$shape$L_a # mm	# a in mm and freq in kHz
if(para$simu$var_indx == 1){
  ka0 = 2 * pi * para$simu$var0 * a / misc$cw		# a (mm) f(kHz)
  ka1 = 2 * pi * para$simu$var1 * a / misc$cw
}
if(para$simu$var_indx == 2){
  ka0 = 2 * pi * para$simu$freq * a / misc$cw
  ka1 = ka0
}
if(para$simu$var_indx == 3){
  ka0 = para$simu$var0
  ka1 = para$simu$var1
}

if ( para$shape$ave_flag == 1){
  ka_min = ka0 * (1 - 3.1 * para$shape$Lstd)
  ka_max=ka1*(1+3.1*para$shape$Lstd)
  if(is.null(nl)){nl=6*para$shape$Lstd/para$shape$dL}
}else{
  ka_min = ka0
  ka_max = ka1
  nl = 1
}

len_ave_para = c(nl, para$shape$Lstd)

Npts = para$simu$n

if(para$simu$var_indx == 2){
  if(para$shape$ave_flag == 1){
    if(is.null(nl)){nl=6*para$shape$Lstd/para$shape$dL}
    kaL = seq(ka0, ka1, length=1)
    ka = seq(ka_min, ka_max,length=nl)
  }else{
    ka = ka0
    kaL = ka
  }
 }else{
   if(Npts == 1){
     ka = ka_min
     kaL = ka0
    }else{
      ka= seq(ka_min, ka_max, length=Npts)
      kaL = seq(ka0, ka1, length=Npts)
  }
}

misc$ang = ang
misc$ka = ka

# conmpute scattering amplitude/L
if(is.null(para$shape$profile)){para$shape$profile = -1}
if(is.null(para$phy$body_ih)){para$phy$body_ih = FALSE}

#define if in script mode or app mode
if(exists('app')==FALSE){app=FALSE}
if(exists('app')==FALSE){app="script"}
dp = ifelse(plotOut==T,1, 0)
dwba_out=DWBAscat2(para, misc, app,shplot=dp)
ka = dwba_out$ka
ang = dwba_out$ang
f = dwba_out$f
ys = f
 # if(status$stop == 1){
 #   return()
 # }

if(para$simu$var_indx == 2){
  angm = seq(para$simu$var0, para$simu$var1, length=Npts)	# mean incident angle
  len_ave_para = c(nl, para$shape$Lstd)
  for(i in 1:Npts){
      orient_ave_para = c(angm[i], para$orient$PDF_para)
      if(para$orient$ave_flag == 1){
        f1 = orient_ave(ang,f,para$orient$PDF,orient_ave_para)
      }else{
        f1 = f[,i]
      }
  if(exists("f2")==FALSE){f2<-NA}
  f2[i] = length_ave(ka,kaL,f1,2,len_ave_para, app=TRUE)[[1]]
  }
}else{
  orient_ave_para = c(para$orient$angm, para$orient$PDF_para)
  f1 = orient_ave(ang,f,para$orient$PDF,orient_ave_para)
  len_ave_para = c(nl, para$shape$Lstd)
  f2 = length_ave(ka,kaL,f1,2,len_ave_para, app=TRUE)
  fsim_l = f2[[2]]
  Lsim = f2[[3]] * para$shape$L
  f2 = f2[[1]]


}

# convert output to the specified quantity
y = abs(f2)
ylab=expression(paste('Normalised Backscattering Amplitude ', f[bs*L^{-1}]))
if(para$simu$out_indx == 2){
  y = y*y
  ylab= expression(paste("Scattering Cross Section ", sigma["bs"]))
}

if(para$simu$out_indx == 3){
  y=20*log10(y*para$shape$L)-60	# L: mm -> m
  ylab=expression(paste('Target Strength (dB re ', m^{2},')'))
}

if(para$simu$out_indx == 4){
  y=20*log10(y)
  ylab=expression(paste('Reduced Target Strength (dB re ', m^{2},')'))
}
if(para$simu$out_indx == 5){
  y=f2
  ylab=expression(paste('Complex scattering'))
}

xlabs <- c('Frequency (kHz)','Orientation angle','ka')
xlab <- xlabs[para$simu$var_indx]

var = seq(para$simu$var0, para$simu$var1, length=para$simu$n)
if(plotOut ==TRUE){
  p = ggplot2::ggplot()+
    ggplot2::geom_line(ggplot2::aes(x=var,y=y),lwd=1.5)+
    ggplot2::xlab(xlab)+
    ggplot2::ylab(ylab)+
    ggplot2::scale_y_continuous(expand=c(0.2,0.2))+
    ggplot2::theme_bw() +
    ggplot2::theme(panel.border =
                     ggplot2::element_blank(),
                       panel.grid.major =
                     ggplot2::element_blank(),
                       panel.grid.minor =
                     ggplot2::element_blank(),
          axis.line =
            ggplot2::element_line(colour = "black"))
}else{p=NULL}

if(simOut==TRUE){
  if(exists("fsim_l")==FALSE){fsim_l<-f2}
  if(exists("Lsim")==FALSE){Lsim <- para$shape$L}
  ys0=ys
  ys = abs(ys)
  ysl =  abs(fsim_l)
  if(para$simu$out_indx == 2){
    ys = ys*ys
    ysl = ysl*ysl
  }

  if(para$simu$out_indx == 3){
    ys=20*log10(ys*para$shape$L)-60	# L: mm -> m
    ysl=20*log10(ysl*para$shape$L)-60	# L: mm -> m
  }

  if(para$simu$out_indx == 4){
    ys=20*log10(ys)
    ysl=20*log10(ysl)
  }
  if(para$simu$out_indx == 5){
    ys=ys0
    ysl=ys0
  }

}
#xlabel(xlab);
#ylabel(ylab);
#dat.var=var;
#dat.fun=y;
#p=0;
if(simOut==TRUE){
  return(list(var = var, y=y,rplot=p, xlab=xlab, ylab=ylab, shplot = dwba_out$shplot, ysim = ys, ang=ang, ysimL = ysl,L=Lsim, comp=ys0))
}else{
  return(list(var = var, y=y,rplot=p, xlab=xlab, ylab=ylab, shplot = dwba_out$shplot))
}
}
AustralianAntarcticDivision/ZooScatR documentation built on Aug. 13, 2022, 1:21 a.m.