mc_psf: Monte Carlo Point Spread Function

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

Monte Carlo Point Spread Function

Description

Backward Monte Carlo Code to calculate the elastic scattering Point Spread Function (integrated in bin area) of a vertically stratified but horizontally homogeneous atmosphere at monocromatic regime.

Usage

mc_psf(
  atm,
  geom,
  res,
  ext,
  snspos,
  snsfov,
  snsznt,
  np,
  mnw,
  cdf_aer,
  cdf_ray = NULL,
  press = NULL,
  Rmc = FALSE
)

Arguments

atm	Atmosphere model, see Details.
geom	One of "annular", "sectorial" or "grid", see Details.
res	Resolution (km), [0,10]. See Details.
ext	Maximum spatial extent at geom resolution (km), see Details.
snspos	Sensor position (km; x, y, z), [0,Inf).
snsfov	Sensor field of view (radians), [0,2*pi].
snsznt	Sensor view zenith angle (radians), [0,pi].
np	Number of photons to trace, [0,2^64 - 1].
mnw	Minimal weight to keep tracing photon package (unitless), [0,1].
cdf_aer	Aerosol scattering CDF LUT (unitless).
press	Atmospheric pressure at surface (mbar). Can be omitted if present on the attributes of atm. See Details.
Rmc	Logical. Should the R version of the MC code be used? Defaults to FALSE. See Details.

Details

This backward Monte Carlo radiative transfer code calculates the area integrals of the Point Spread Function of the atmosphere due to elastic scattering in a layered atmosphere. Turbulence effects and inelastic scattering are not included. Additional simplifications of the current 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 atmosphere profile atm should be a named data frame as the one created by function 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. As an example, the FOV of the Operational Land Imager (OLI/Landsat 8) is ~ 8.6E-5(radians). To simulate a sensor with infinitesimal FOV, set snsfov = 0. The FOV is centered in sensor zenith angle. The photon package is initiated at the upper boundary of the layer where the sensor is located. The sensor azimuth is set to 90 degrees, such that it is to the right of the center of the PSF grid.

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 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 an_psf.

An R version of the code is also provided. It is present mainly for validation reasons but can be useful in a system without a compiler to compile the C code. It will run many times slower though.

Value

A list with elements 'dirtw' containing the fractional contribution of the upward directly transmitted photons, 'bin_phtw' containing the fractional contribution of upward diffusely transmitted 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_psf(atm     = atm, 
                  geom    = 'annular', 
                  res     = 0.03, 
                  ext     = 10, 
                  snspos  = c(0, 0, 800), 
                  snsfov  = 0, 
                  snsznt  = 0, 
                  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))

AlexCast/apsfs documentation built on Feb. 1, 2024, 9:48 p.m.