Sim formulas and bibliography
The goal is to take the solar irradiance received by the body, and compute the resulting reflected irradiance for an observer. We need to :
- Apply the inverse square law (quantity is inversely proportional to the square of the distance)
- Take into account the fact that a body doesn't "mirror" the received energy, but "diffuse" it. At this point, we assume it's a lambertian surface and that all energy is reflected isotropically and hemispherically.
- Take into account the fact that the visible fraction of the body sphere is dependent on the altitude. At an altitude close to 0, the reflected solar irradiance should be near 0.
Note : the stock code uses the same formula, excepted that it doesn't take that last point into account, consequently, it gives plain wrong values at medium/low altitudes.
- with r = body radius,
- with a = altitude,
- with i = solar irradiance at the body position,
- with sr = steradians,
- At the current altitude, the visible portion of the body disc surface area is
π * (r * sin(arccos(r /(a + r))))²
- So total energy received on that disc area is
i * π * (r * sin(arccos(r /(a + r))))² - That energy is isotropically re-emitted hemispherically, that is a solid angle of
2 * π * sr - So the energy emitted in watts per steradian can be expressed as
i * π * (r * sin(arccos(r /(a + r))))² / (2 * π * sr) - The half-sphere enclosing the body at the given altitude has a surface area of
2 * π * (r + a)² - This translate in a surface area / steradian of
2 * π * (r + a)² / (4 * π sr) = (r + a)² / (2 * sr) - So the albedo radiosity at the current altitude is
i * π * (r * sin(arccos(r /(a + r))))² / (2 * π * sr) / (((r + a)² / (2 * sr))) - After simplification, this gives
(i * r² * (a² + 2 * a * r)) / (a + r)^4
Visualization for both formulas, in the case of Kerbin : https://www.desmos.com/calculator/zuyyhiqsnu
- X is the altitude in meters, Y the reflected solar irradiance in W/m²
- Orange line is the solar irradiance (1360 W/m² at Kerbin)
- Green curve is the stock formula
- Purple curve is the altitude-corrected formula
This gives the final method:
public double IsotropicIrradianceAtAltitude(double irradiance, double altitude, double bodyRadius)
{
double radiusSquared = bodyRadius * bodyRadius;
double altitudeSquared = altitude * altitude;
double nominator = irradiance * radiusSquared * (altitudeSquared + 2 * altitude * radius);
return nominator / Math.Pow(altitude + radius, 4.0);
}To take into account the phase angle, we apply a simple dot product between the body-sun direction and the body-vessel direction. Stock does this too.
However, bodies aren't perfect lambertian surfaces : the incoming irradiance is not entirely diffused hemispherically, but (usually) more strongly toward the direction of the source (the sun). The effect is described as the opposition surge and expressed by the geometric albedo.
We assume that the CelestialBody.albedo value is a bond albedo (since the stock code assume a [0 ; 1] range). We assume that the "average" geometric albedo is ~1.25 times the bond albedo on average, based on some quick and dirty data crunching (various sources, but notably "Planetary sciences, Imke De Pater / Jack Jonathan Lissauer) :
| Body | Bond albedo | Geometric albedo | Atmo ? | Ratio |
|---|---|---|---|---|
| Mercury | 0,12 | 0,14 | 0 | 1,16 |
| Moon | 0,11 | 0,12 | 0 | 1,09 |
| Enceladus | 0,80 | 1,40 | 0 | 1,75 |
| Pluto | 0,40 | 0,53 | 0 | 1,31 |
| Venus | 0,75 | 0,84 | 1 | 1,12 |
| Venus (alt) | 0,77 | 0,69 | 1 | 0,90 |
| Earth | 0,31 | 0,43 | 1 | 1,42 |
| Earth (alt) | 0,31 | 0,37 | 1 | 1,20 |
| Mars | 0,16 | 0,15 | 1 | 0,94 |
| Mars (alt) | 0,25 | 0,15 | 1 | 0,60 |
| Jupiter | 0,34 | 0,52 | 1 | 1,52 |
| Saturn | 0,34 | 0,50 | 1 | 1,46 |
| Saturn (alt) | 0,34 | 0,47 | 1 | 1,38 |
| Uranus | 0,30 | 0,49 | 1 | 1,63 |
| Uranus (alt) | 0,29 | 0,51 | 1 | 1,76 |
| Neptune | 0,29 | 0,44 | 1 | 1,52 |
| Neptune (alt) | 0,31 | 0,41 | 1 | 1,32 |
| Titan | 0,20 | 0,21 | 1 | 1,05 |
| Average (total) | 1,28 |
Note that in the context of our irradiance representing the flux at all wavelengths (not only visible light), this is likely wrong, as the above data for geometric albedo is for the human-visible portion of the electromagnetic spectrum, and values for other wavelengths (infrared, etc) can be significantly different (source) :
U = ultraviolet, B = "blue", V = visible light, R = red, I = infrared...
To apply that factor, we take the result of the sun-body-observer dot product, and apply a very basic polynomial, ensuring that the total energy is conserved (the [0 ; 1] integral is equal, see https://www.desmos.com/calculator/jtgsj3wc4v )
The formula gives the following method :
public double GeometricAlbedoFactor(double sunBodyObserverAngleFactor)
{
return 1.25 * Math.Pow(sunBodyObserverAngleFactor, 1.5);
}public double AlbedoIrradiance(Vector3d sunToBodyDirection, Vector3d observerToBodyDirection, double altitude, double starIrradianceAtBody, double bodyRadius, double bodyBondAlbedo)
{
double isotropicIrradiance = IsotropicIrradianceAtAltitude(starIrradianceAtBody, bodyRadius);
double sunBodyObserverAngleFactor = (Vector3d.Dot(sunToBodyDirection, observerToBodyDirection) + 1.0) * 0.5;
double geometricFactor = GeometricAlbedoFactor(sunBodyObserverAngleFactor);
return isotropicIrradiance * geometricFactor * bodyBondAlbedo;
}WIP "all-in-one" desmos graph : https://www.desmos.com/calculator/avtgokyyju https://www.desmos.com/calculator/4e73qklv9n
Earth emissive irradiance (source : ECSS-E-ST-10-04C)
The Earth-emitted thermal radiation has a spectrum of a black body with a characteristic
average temperature of 288 K. The Earth infrared radiation also varies across the globe but
less than the albedo. It also shows a diurnal variation which is small over the ocean but can
amount to 20 % for desert areas. The average infrared radiation emitted by Earth is 230 W m-2. For an orbiting spacecraft, it
can vary from 150 W m-2 to 350 W m-2. The diurnal variations can amount to about 20 %
over desert areas while being smaller over the oceans.