Using Kirchhoff’s and the high-frequency assumptions, we assume the local points on the planetary surface to be a tangent plane approximation to the actual surface and we assume the roughness is smaller than the wavelength thus the scattering of the local surface is equivalent to the specular reflection of a smooth surface. In this case, the direction of the reflected wave will obey Huygens’ principle such that the incident and reflected (in this case scattered) angle and points are equal and the incident ray, reflected ray, and surface normal all lie in the same plane as seen in
Figure 1. The magnitude of the reflected ray can be determined using Fresnel’s reflection coefficients. For this model, we will assume all the Fresnel reflection coefficients are unity. Thus, to determine the reflected ray, the incident direction, the incident surface location, and the normal to the surface at the interception point are necessary. As shown in
Figure 2, the position of the satellite at an epoch time
, is designated as
, the position of the receiver antenna at epoch time is
, and
is the centre position of the planetary surface. The position of the central ray intercepting the surface at epoch time is
, which can be obtained by referring to
,
and
or else by using the SPICE software functions.
While illuminating the surface with a large number of radio rays would lead to a higher resolution it would also lead to slower ray-tracing computation due to the increased need for ray interception checking. To limit the ray tracing computational time, the model initially sets the number of rays equal to the number of surface elements, , with N being equal to the number of divisions along each side of the half-power beamwidth, which generates a square grid of elements.
3.1. The Simulation of Reflected Signals
To trace the ray through the model, we assume that the satellite transmits a bundle of rays at a certain half-beam angle, α. For the sake of the algorithm, the beam footprint can be assumed to consist of a central curve along with a left and right curve. This algorithm processes the ray interception points along the central curve first and then the lateral points going left to right along the footprint. We use the Law of Reflection which states that the incident ray, reflected ray, and surface normal all exist in the same plane and that the incident angle equals the reflected angle. In short, a point on the central curve exists within a plane defined by that point and the positions of the satellite and planet and this plane has a normal which goes through the point and the center planet position as seen in
Figure 2. Thus, the reflected ray position and angle can be found rotating the point around the normal. While in reality the reflected ray position equals the surface interception ray position, for the sake of this algorithm it is defined as another point. This new point, similar to the initial point, is defined as existing within a plane defined by itself and the positions of the satellite and planet with a normal going through itself and the planet’s position. At the local ray interception position, on the central curve,
, the reflected ray position and direction can be obtained by the rotation of
around the normal to the plane on the central curve determined by the surface interception point,
along with
and
, denoted as
. The normalized vector, denoted as
, that
rotates around, is the normal to this plane,
, and goes through point
. The range of the rotation angle,
, is between
, and the
rotated by
is denoted as
, shown in
Figure 2; this point represents the reflected/scattered point. This point, similar to the initial ray interception point,
, is also defined as existing within a plane defined by the reflected point,
along with
and
. The normal to
at
is denoted as
. The vector from
to
, denoted as
, is perpendicular to the normal at
,
.
To determine a position laterally offset from the central curve, denoted as
, the offset length between the corresponding point on the central curve,
, and
can be determined using the geometrical relationships between the length of the vector from
to
, denoted as
, which is perpendicular to the normal at
,
and the angle,
, formed between the vectors
and
as shown in
Figure 2. The distance between
and
, denoted as
, can be obtained by
The position of the offset point is thus:
where
,
is normalized
.
Using a ray-tracing algorithm we can determine the position and direction of scattered/reflected rays. However, these reflection points are generated by assuming a smooth surface. In practice, the surface is rough, and each illuminated surface element can be described as having a height profile with a particular root mean square height distribution, denoted here as
, and a particular correlation length which determines the spatial frequencies over the surface in the x and y directions, denoted herein as
and
. In this model, we consider the local illuminated surface to be a square area with a local surface coordinate system
. The origin of coordinate system is denoted as
. We use the spacecraft as a reference point and maintain that the directions from north to south and from left to right are positive. The characteristics of the
system are shown in
Figure 3.
Overlaid upon this coordinate system is an auxiliary square grid, with Gaussian elevation associated to the nodes of the grid, used to simulate the rough surface; thus, the number of auxiliary elements in the or direction is , resulting in total elements. The number must conform to the Nyquist theorem in the spatial domain such that the number of auxiliary elements per unit length is inversely related to the correlation length for its lower limit. This relationship also provides constraints for the side length of the square facet (N/SL > 2/CoX).
The element side length visible to the transmitter at a certain time epoch is approximately:
To simplify the description, the generated auxiliary elements overlaid on the coordinate system are numbered as where and are row number and column number respectively. The elevation data of is initialized by .
Figure 4a shows an example of an
grid of auxiliary elements. To generate a random matrix with a Gaussian distribution, as seen in
Figure 4b, we use MATLAB’s random number generator function
.
To produce a spatially correlated surface, this surface is convolved with a Gaussian filter using MATLAB’s discrete Fast Fourier Transform (FFT) function. When
, the surface is considered isotropic, meaning the roughness is the same in the
and
directions and thus scatters light evenly in all directions. To achieve a spatially correlated isotropic surface, the Gaussian filter with a low pass frequency response is
Where
and
are coordinates of
on
and
respectively. Using this Gaussian filter in addition to the inverse
, to change from the spectral domain to the spatial one, we can obtain the rmsH for the facet,
where * is the dot product between two vectors. The spatially correlated surface seen in
Figure 4c, is the result of the convolution of H and the Gaussian filter for an isotropic surface. On the other hand, When
and
are different, the surface is considered anisotropic, meaning the roughness varies in the
and
directions thus scattering light in different directions with different intensities. To achieve a spatially correlated surface for an anisotropic, the Gaussian filter with a low pass frequency response is
With this auxiliary grid we can correct the initial reflection point by adding the rms height (
)
Now that the reflection points on the rough surface elements can be calculated there is the possibility that an overlapping or shadowed area possibly exists between two adjacent illuminated surface elements when the moving step of the spacecraft during the simulation is short. Considering that the elevation data of each illuminated surface is independently generated, it is necessary to correct the elevation data of the latter illuminated surface located in the overlapping part. Supposing the moving step is denoted as , the point belonging to the latter illuminated surface and overlapping part can be denoted as . After projecting the overlapping part on the system, it is easy to find which element belonging to the former illuminated surface that will be located in. Supposing this element is , can be corrected and replaced by the intersection of vector and element .
Once the reflection positions are determined, the direction vector of the reflected ray must be determined. The rough surface is meshed using a collection of non-overlapping triangular facets approximating the surface shape. In this way, the intersection of a ray with each triangular facet is calculated. The points selected to create the surface element are
,
and
, where
is selected as the reflection point of this surface element. An arbitrary facet denoted as
is seen in
Figure 5 projected above the local surface coordinate system.
To determine the normal to the triangular facet plane,
, the cross product of two edges is calculated. To find the vector of the reflected ray
, the component of the incident ray is reversed parallel to the triangle facet normal such that
where * is the dot product,
is the incident ray.