This program, shadow, will create a shaded-relief cube from a digital elevation model (DEM) and a MATCH cube.

We use the sun's position at the center of the input cube or a user-defined observation time. By default, we factor in shadows cast by features. This program operates much like the ISIS 'shade' program which instead requires azimuth/elevation as input.

However, with the shadow application, using the sun's position allows much higher precision shading and enables the possibility of computing shadowed areas. The algorithm description below is provided to help understand the optimization settings.

User-Requirements
The user must supply an elevation model (DEM) and either an observation time or a cube with raw camera geometry (see 'spiceinit'). If providing an observation time, the PCK and SPK SPICE kernels associated with the DEM's target must have been downloaded and be available on the user's machine in the "base" kernels area, or their file paths must be provided as values to the PCK and SPK user parameters.

Understanding the Algorithm
The 'shade' program's algorithm is called 'hillshade' in this algorithm description.

1. Compute the sun's position at the center of the MATCH cube. The MATCH cube must have the same target as the DEM.
2. For every pixel in the input elevation model, compute the hillshade value for the pixel, and
   
1. if the hillshade result is positive (facing towards the sun), then estimate if the pixel is in shadow;
2. if the hillshade value is positive and the pixel is in shadow, then the result is an LRS;
3. if the hillshade value is positive and the pixel is not in shadow, then the result is the hillshade result and
4. if the hillshade value is negative, the result is low resolution saturation (LRS).

The algorithm to estimate if a pixel is in shadow:

1. Optimization: If SHADOW, and this elevation model pixel is known to be shadowed, consider this pixel to be shadowed and stop (no pixels are initially known to be shadowed).
2. Compute the pixel's body-fixed coordinate (XYZ) position.
3. If SUNEDGE, adjust the sun's position to nearer the highest (on the horizon) edge of the sun.
4. Subtract the sun's position from the elevation model pixel's position, providing a vector to the center (or edge) of the sun.
5. Iterate until a solution is found:
1. Step along the 3D ray the current estimate of PRECISION pixels
2. Project the ray back onto the elevation model to find the equivalent radius
3. Update estimate of how far along the 3D ray is equivalent to the full resolution of the elevation model.
   Optimization: Multiply the step by PRECISION.
   Optimization : If SKIPOVERSHADOW, while the next linearly-extrapolated elevation model position is known to be in shadow, increase the next step size by the estimate up to MAXSKIPOVERSHADOWSTEPS times.
4. Check for a solution
• If the equivalent radius is higher than the ray, then the originating pixel is shadowed.
  
• If the ray's elevation is higher than the highest point on the elevation model, the originating pixel is in light.
  Optimization: If LIGHTCURTAIN, and the ray's elevation is higher than a previous ray that intersected this pixel in the elevation model, consider the originating pixel in light.
5. Optimize
   • If LIGHTCURTAIN, and the pixel was determined to be in light, record the elevations of the ray where it projected onto the elevation model.
     Optimization: If LOWERLIGHTCURTAIN, lower the elevations along the ray by subtracting the minimum difference between the ray and the elevation model while the ray was being walked.
     Optimization: If CACHEINTERPOLATEDVALUES, linearly interpolate the points where the ray would have intersected the elevation model to one pixel accuracy.
   • If SHADOWMAP, and the pixel was determined to be in shadow, record all points where the ray projected onto the elevation model to be known shadowed points, excluding the actual intersection point.
     Optimization: If CACHEINTERPOLATEDVALUES, linearly interpolate the points where the ray would have intersected the elevation model to one pixel accuracy.

The caches inherently cause inaccuracies (approximately 1-2 pixels if CACHEINTERPOLATEDVALUES is off) in the output shadow positions because they record sub-pixel values as if they were the center of the pixel. Periodically, the caches are shrunk to lower memory usage. The light curtain cache is shrunk to around BASELIGHTCACHESIZE entries and the shadow map cache is shrunk down to around BASESHADOWCACHESIZE entries. Each cache entry is approximately 24 bytes for the light cache, and 16 bytes for the shadow cache, but the caches are hash-based causing a large amount of potential overhead. The caches are not limited to their specified sizes, only reduced to them periodically, and larger cache sizes result in this program consuming more memory. The larger the cache, the longer it takes to lookup a single value (which happens often). In general, larger caches mean less CPU time for the cost of memory.

This example will cover running this program in balanced mode. The input file "localdem.cub" was produced with the following commands:

          spiceinit from=ab102401.cub
          fx from=ab102401.cub to=ab102401.radii.cub equation="radius(f1)"
          cam2map from=ab102401.radii.cub to=ab102401.radii.lev2.cub
          demprep from=ab102401.radii.lev2.cub to=localdem.cub
        

This example will cover running this program in balanced mode with a lunar DEM and a time.