F# for real: a troposphere model
F# is a functional language for .NET with some imperative elements. It has full support in Visual Studio 10 and later (except for the Express editions). F# can be a useful complement to C# but you have to keep them in separate projects.
Instead of making a tutorial we are going to look at a real world example where F# can help us. This is mainly thanks to a feature called Units of Measure which lets us introduce measurement units into the type system.
The example I have written is an implementation of a tropospheric error model for GPS as described in the system specifications for SBAS (an augmentation system for GPS). Section 6.3 on page 32 in this PDF has more information.
namespace GpsErrorModels
open System
// Units and equivalences
[<Measure>] type degree
[<Measure>] type kelvin
[<Measure>] type kilogram
[<Measure>] type meter
[<Measure>] type second
[<Measure>] type joule = kilogram meter^2 / second^2
[<Measure>] type millibar = kilogram / meter / second^2
// Tropospheric model (defined in RTCA/DO-229C) for a specific receiver
// position (latitude and height) and time (date)
type TroposphereModel(latitude:float<degree>, height:float<meter>,
time:DateTime) =
let Interpolate_ζ(latitude:float<degree>, values:float<'u> list):float<'u> =
let interpolate(k, k0, k1, v0, v1):float<'u> =
// Interpolate (equations A-4 and A-5)
v0 + (v1 - v0) * (k - k0) / (k1 - k0);
let latitudes = [15.0<degree>; 30.0<degree>; 45.0<degree>;
60.0<degree>; 75.0<degree>;]
// Interpolation according to table A-2
if latitude < latitudes.[0] then values.[0]
elif latitude < latitudes.[1] then
interpolate(latitude, latitudes.[0], latitudes.[1],
values.[0], values.[1])
elif latitude < latitudes.[2] then
interpolate(latitude, latitudes.[1], latitudes.[2],
values.[1], values.[2])
elif latitude < latitudes.[3] then
interpolate(latitude, latitudes.[2], latitudes.[3],
values.[2], values.[3])
else values.[4];
let Get_ζ(latitude:float<degree>, time:DateTime, averages:float<'u> list,
variations:float<'u> list):float<'u> =
let ζ0 = Interpolate_ζ(latitude, averages)
let Δζ = Interpolate_ζ(latitude, variations)
let D = (float)time.DayOfYear
let Dmin = if latitude > 0.0<degree> then 28.0 else 211.0
// Compute meterological parameter (equation A-3)
ζ0 - Δζ * Math.Cos(2.0 * Math.PI * (D - Dmin) / 365.25);
// Pressure values from table A-2
let Get_P(latitude:float<degree>, time:DateTime) =
let averages = [1013.25<millibar>; 1017.25<millibar>; 1015.25<millibar>;
1011.25<millibar>; 1013.00<millibar>;]
let variations = [ 0.00<millibar>; -3.75<millibar>; -2.25<millibar>;
-1.75<millibar>; -0.50<millibar>;]
Get_ζ(latitude, time, averages, variations);
// Temperature values from table A-2
let Get_T(latitude:float<degree>, time:DateTime) =
let averages = [299.65<kelvin>; 294.15<kelvin>; 283.25<kelvin>;
272.15<kelvin>; 263.65<kelvin>;]
let variations = [ 0.00<kelvin>; 7.00<kelvin>; 11.00<kelvin>;
15.00<kelvin>; 14.50<kelvin>;]
Get_ζ(latitude, time, averages, variations);
// Water vapor pressure values from table A-2
let Get_e(latitude:float<degree>, time:DateTime) =
let averages = [26.31<millibar>; 21.79<millibar>; 11.66<millibar>;
6.78<millibar>; 4.11<millibar>;]
let variations = [0.00<millibar>; 8.85<millibar>; 7.24<millibar>;
5.36<millibar>; 3.39<millibar>;]
Get_ζ(latitude, time, averages, variations);
// Temperature lapse rates from table A-2
let Get_β(latitude:float<degree>, time:DateTime) =
let averages = [6.30E-3<kelvin/meter>; 6.05E-3<kelvin/meter>;
5.58E-3<kelvin/meter>; 5.39E-3<kelvin/meter>;
4.53E-3<kelvin/meter>;]
let variations = [0.00E-3<kelvin/meter>; 0.25E-3<kelvin/meter>;
0.32E-3<kelvin/meter>; 0.81E-3<kelvin/meter>;
0.62E-3<kelvin/meter>;]
Get_ζ(latitude, time, averages, variations);
// Water vapor "lapse rates" from table A-2
let Get_λ(latitude:float<degree>, time:DateTime) =
let averages = [2.77<1>; 3.15<1>; 2.57<1>; 1.81<1>; 1.55<1>;]
let variations = [0.00<1>; 0.33<1>; 0.46<1>; 0.74<1>; 0.30<1>;]
Get_ζ(latitude, time, averages, variations);
let ZenithDelay:float<meter> =
// Meteorological values from table A-2
let P:float<millibar> = Get_P(latitude, time)
let T:float<kelvin> = Get_T(latitude, time)
let e:float<millibar> = Get_e(latitude, time)
let β:float<kelvin/meter> = Get_β(latitude, time)
let λ:float<1> = Get_λ(latitude, time)
// Constants defined on page A-9 and A-10
let k1 = 77.604<kelvin/millibar>
let k2 = 382000.0<kelvin^2/millibar>
let Rd = 287.054<joule/kilogram/kelvin>
let gm = 9.784<meter/second^2>
let g = 9.80665<meter/second^2>
// Zero-altitude zenith dry delay (equation A-6)
let z_hyd:float<meter> = (1E-6 * k1 * Rd * P)
/ gm
// Zero-altitude zenith wet delay (equation A-7)
let z_wet:float<meter> = (1E-6 * k2 * Rd)
/ (gm * (λ + 1.0) - β * Rd)
* e / T
// Height-mapped zenith dry delay (equation A-8)
let d_hyd:float<meter> = (1.0 - β * height / T)
** (g / (Rd * β))
* z_hyd
// Height-mapped zenith wet delay (equation A-9)
let d_wet:float<meter> = (1.0 - β * height / T)
** ((λ + 1.0) * g / (Rd * β) - 1.0)
* z_wet
// Combined zenith dry and wet delay (part of equation A-2)
d_hyd + d_wet;
// Calculate the tropospheric delay for a satellite given its elevation
member this.GetCorrection(elevation:float<degree>):float<meter> =
let Sin(angle:float<degree>):float =
Math.Sin(((float)angle * Math.PI / 180.0))
// Pre-calculated zenith delay
let d = ZenithDelay
// Elevation mapping value (equation A-10)
let m = 1.001 / sqrt (0.002001 + Sin(elevation) ** 2.0)
// Tropospheric delay correction (equation A-2)
-d * m;To show how the unit system is used I have highlighted the relevant lines leading up to the z_wet equation. Because the code compiles I can be reasonably sure the equation is correct. If I had made a mistake and the units didn’t match up it would have shown as an error in exactly the same way as if the types didn’t match.
This system isn’t perfect. You can still make mistakes that aren’t caught. But it gives an extra level of confidence that may well be worth it in some applications.
TroposphereModel troposphere = new TroposphereModel(
59.35264, 39.773, DateTime.Parse("2004-07-08"));
Debug.Assert(Math.Abs(troposphere.GetCorrection(60.4) + 2.791) < 0.05);
Debug.Assert(Math.Abs(troposphere.GetCorrection(37.0) + 4.024) < 0.05);
Debug.Assert(Math.Abs(troposphere.GetCorrection(12.3) + 11.144) < 0.05);Here is a piece of C# code where we call the F# code and check the results against a set of known values taken from a Nordnav-R30 GPS receiver.