R/asaf.R

In AlexCast/apsfs: Monte Carlo code for simulation of atmospheric Point Spread Functions and spatially-resolved Spherical Albedo Functions

#' Monte Carlo spatially-resolved spherical albedo
#'
#' Backward Monte Carlo Code to calculate the elastic scattering 
#' spatially-resolved spherical albedo function (integrated in bin area) of a 
#' vertically stratified but horizontally homogeneous atmosphere at monocromatic
#' regime. The spatially-resolved spherical albedo is the angular integral of 
#' the bi-directional scattering-surface reflectance distribution function 
#' (BSSRDF). The spherical albedo is defined as the normalized irradiance 
#' backscattered by the atmosphere when the input irradiance at the bottom is 
#' isotropic (Vermote et al., 1997).
#'
#' Vermote, E. F., Tanre, D., Deuze, J. L., Herman, M., Morcette, J.-J. (1997). 
#' Second Simulation of the Satellite Signal in the Solar Spectrum, 6S: an 
#' overview. IEEE Transactions on Geoscience and Remote Sensing, 35(3), 675–686. 
#'
#' @param atm     Atmosphere model, see Details. 
#' @param cdf_aer Aerosol scattering CDF LUT (unitless).
#' @param geom    One of "annular", "sectorial" or "grid", see Details.
#' @param res     Resolution (km), [0,10]. See Details.
#' @param ext     Maximum spatial extent at geom resolution (km), see Details. 
#' @param np      Number of photons to trace, [0,2^64 - 1].
#' @param mnw     Minimal weight to keep tracing photon package (unitless), 
#'                [0,1].
#' @param press   Atmospheric pressure at surface (mbar). Can be omitted if 
#'                present on the attributes of atm. See Details.
#'
#' @details
#' This backward Monte Carlo radiative transfer code calculates the area 
#' integrals of the spatially-resolved spherical albedo of the atmosphere due to
#' elastic scattering in a layered atmosphere. Turbulence effects and inelastic 
#' scattering are not included. Additional simplifications of the current 
#  version include scalar solution only (polarization not included), Lambertian 
#' surface reflectance (Bidirectional Reflectance Distribution Function effects 
#' not included). The surface is considered totally absorbing, such that 
#' multiple interaction between surface and atmosphere are not included.
#'
#' The sensor position is at the origin, at the bottom of the atmosphere, and 
#' the zenith view angle is 180 degrees (looking at the zenith). The FOV is 
#' hemispherical.
#'
#' The atmosphere profile atm should be a named data frame as the one created by 
#' function \code{get_opt_atm_prfl}. If the profile is create manually or with 
#' user functions, it must contain at least the same variables and column names.
#' It may also include pressure in the attributes (attr(atm, "press") <- value),
#' and in this case the parameter press can be omitted.
#'
#' The initial direction of the photon package is sampled assuming equal 
#' sensibility of the sensor to all directions within its FOV (hemispherical 
#' in the case of irradiance).
#'
#' At each iteration, a free optical path is sampled and the photon package is 
#' moved in metric scale. A check is made to verify if photon reached the 
#' surface or if it escaped the system. If it reached the surface, its intensity
#' ("photon package weight") is added to the appropriate accumulator bin. The 
#' accumulator geometry can be: (1) 'annular', which only keep track of the 
#' annular bin in which the photon reached the surface; (2) 'sectorial', wich 
#' divide the annuli in sectors for azimuthal assymetry; or (3) 'grid', which 
#' will keep track of the {x,y} position. The 'annular' accumulator is 
#' appropriate for symmetrical conditions, as nadir view over a Lambertian 
#' surface. Note that regardless of extent of the accumulator, a final bin will 
#' be added between ext and Inf. The resolution of the accumulator will have an
#' impact on the performance, as too high resolution, specially in 'grid' 
#' geometry will require more photons to be traced to reduce statistical 
#' sampling variation. Note that when using 'sectorial' geometry, angular 
#' resolution in azimuth is fixed to 1 degree and break points returned by the 
#' function refer only to the radius break points (angular break points will be
#' 0:360 always).
#'
#' The reduced performance with high resolution grid and the need for high 
#' resolution close to the target pixel while low resolution is sufficient for 
#' larger distances suggests two separate runs. Those grids can be combined
#' later with the function \code{fit_grid_dpsf} to produce a Zernike polynomial
#' fit to the density PSF. This model can then be used to generate a discrete 
#' PSF kernel at any sensor resolution.
#'
#' To reduce the random sampling fluctuation in the sectorial and grid 
#' geometries, the grid is 'folded' on its axis of symmetry (sensor azimuth 
#' line), averaged and then 'unfolded'.
#'
#' Note that the returned values are fractional contribution per accumulation 
#' bin (integrals), and so are not normalized per area. Conversion to per area
#' contribution can be made with function \code{an_psf}.
#'
#' @return
#' A list with elements 'bin_phtw' containing the fractional contribution of 
#' downward reflected photons per accumulation bin (depends on geom), 'bin_brks'
#' containing the break points of the accumulation bins, 'bin_mid' with the mid 
#' point of the bins and 'metadata' containing the input parameters.
#'
#' @examples
#' # Using the C version of the code:
#' data(us76)
#' atm <- rayleigh_od(us76, lambda = 550) %>%
#'        get_opt_atm_prfl(atm = us76, tau_aer = 0.5, H_aer = 2, w0_aer = 0.89, 
#'          tau_ray_z = ., a_mol_z = rep(0, nrow(us76)))
#' cdf_aer <- cdf_lut(continental_ph_6sv[, c(1, 8)])
#' psf_int <- mc_saf(atm     = atm, 
#'                   geom    = 'annular', 
#'                   res     = 0.03, 
#'                   ext     = 10, 
#'                   np      = 1E5, 
#'                   mnw     = 1E-6,
#'                   cdf_aer = cdf_aer)
#' par(mar = c(5, 6, 3, 2))
#' x <- psf_int$bin_brks
#' plot(x, cum_psf(psf_int), xaxs = "i",  yaxs = "i", xlab = "Radius (km)", 
#' ylab = expression(2*pi*integral("r'"*f("r'")*d*"r'", 0, r)), xlim = c(0, 10))
#'
#' @useDynLib apsfs C_mc_saf
#' @export

mc_saf <- function(
    atm, 
    geom, 
    res, 
    ext, 
    np, 
    mnw,
    cdf_aer, 
    cdf_ray = NULL, 
    press = NULL) {

    ain <- atm
    atm <- data.frame(
        km     = ain$km,
        tau    = ain$tau,
        tau_u  = c(ain$tau[-length(ain$tau)], NA),
        tau_d  = c(ain$tau[-1], NA),
        dkm    = c(abs(diff(ain$km)), NA),
        dtau   = c(abs(diff(ain$tau)), NA),
        w0_tot = ain$w0_tot,
        b_ray  = ain$b_ray,
        b_aer  = ain$b_aer,
        b_tot  = ain$b_ray + ain$b_aer,
        c_tot  = ain$c_tot
    )

    snspos <- c(0, 0, 0) 
    snsznt <- pi

    if(is.null(cdf_ray)) cdf_ray <- rayleigh_cdf

    # Check pressure information. It is not necessary on the MC code per se, but 
    # it has to be included on the metadata of the PSF object for further 
    # processing.
    if(is.null(press)) {
        press <- attr(ain, "press")
        if(is.null(press)) {
            "Pressure mst be defined in press is not present as an attribute of atm" %>%
            stop(call. = FALSE)
        }
    }

    if(geom == "annular") {
        geom <- 1
    } else if(geom == "sectorial") {
        geom <- 2
    } else if(geom == "grid") {
        geom <- 3
    } else {
        stop("geom must be 'annular', 'sectorial' or 'grid', see Details.", 
              call. = FALSE)
    }

    if(np > (2^64 - 1)) {
        warning("Maximum number of photons to simulate is 2^64 - 1 (18,446,744,073,709,551,615). np set to that limit.", call. = FALSE)
        np <- (2^64 - 1)
    }

    saf <- .Call(
        "C_mc_saf", 
        atm, 
        geom, 
        res, 
        ext, 
        snspos, 
        snsznt, 
        np, 
        mnw, 
        cdf_aer, 
        cdf_ray
    )

    names(saf) <- c("bin_phtw", "bin_brks", "bin_mid")

    # Averaging along axis of symmetry to reduce effect of statistical 
    # fluctuations:
    if(geom == 2) {
        saf$bin_phtw <- matrix(saf$bin_phtw, nrow = length(saf$bin_mid))

        nd  <- ncol(saf$bin_phtw)
        ndh <- nd / 2 

        saf$bin_phtw[, 1:ndh] <- (saf$bin_phtw[, 1:ndh] + 
                                  saf$bin_phtw[, nd:(ndh+1)]) / 2
        saf$bin_phtw[, nd:(ndh+1)] <- saf$bin_phtw[, 1:ndh]
    }

    if(geom == 3) {
        saf$bin_phtw <- matrix(saf$bin_phtw, ncol = length(saf$bin_mid))

        # Averaging along axis of symmetry to reduce effect of statistical 
        # fluctuations:
        nd  <- nrow(saf$bin_phtw)
        ndh <- (nd - 1) / 2 

        saf$bin_phtw[1:ndh,] <- (saf$bin_phtw[1:ndh,] + 
                                 saf$bin_phtw[nd:(ndh+2),]) / 2
        saf$bin_phtw[nd:(ndh+2),] <- saf$bin_phtw[1:ndh,]
    }

    geom <- switch(geom, "annular", "sectorial", "grid")

    saf$metadata <- list(
        press   = press,
        atm     = ain,
        cdf_aer = cdf_aer, 
        snspos  = snspos, 
        snsfov  = pi,
        snsznt  = snsznt, 
        geom    = geom,
        res     = res, 
        ext     = ext, 
        np      = np, 
        mnw     = mnw
    )

    class(saf) <- c("saf", class(saf))

    saf
}