Fortress example of WMM2010
Hello Fortress developers,
I have written up a Fortress implementation of the World Magnetic Model
for 2010-2015. (http://www.ngdc.noaa.gov/geomag/WMM). It calculates the
field strength of the earth's magnetic field and is wildly used in
navigation system. See attached file.
This model was chosen because they provide a very nice numerical example
which were used to check the results. Most of the formula was taken
directly from the technical report, without looking at the C
implementation. Only the PlmSchmidt_d1() function was extracted from C code.
Fortify tool can process most the code. The urls in the comments seems
to cause problems because it tries to create escape characters which
later cause LaTeX to break.
Anyway, thought you might like it. Feedback is welcome.
Cheers
Francois
component geomag
export Executable
(*
* This is an example of the World Magnetic Model for 2010-2015.
* Reference: http://www.ngdc.noaa.gov/geomag/WMM/DoDWMM.shtml
*)
run()
= do
println "High-precision numerical example of the WMM2010 model ..."
(*
* Reference:
* [1] "Section 1.5: MODEL EQUATIONS NUMERICAL EXAMPLE",
* The US/UK World Magnetic Model for 2010-2015, Technical Report,
* http://www.ngdc.noaa.gov/geomag/WMM/data/WMM2010/WMM2010_Report.pdf
*)
(* High-precision numerical example:
* ||-------------------------------------------------||
* || Time || 2012.5000 0000 || yr ||
* || Height-above-Ellipsoid || 100.0000 0000 || km ||
* || Latitude || -80.0000 0000 || deg ||
* || Longitude || 240.0000 0000 || deg ||
* ||-------------------------------------------------||
*)
date = 2012.5 (*) years
lat = -80 (*) degrees
lon = 240 (*) degrees
alt = 100 (*) km
geomagnetic(date,lat,lon,alt)
(*
* ||----------------------------------------------||
* || t || 2012.50000 00000 || yr ||
* || lambda || 4.18879 02048 || rad ||
* || phi || -1.39626 34016 || rad ||
* || h || 1 00000.00000 00000 || m ||
* || phi-prime || -1.39512 89589 || rad ||
* || r || 64 57402.34844 73705 || m ||
* || g(1,0,t) || -29467.60000 00000 || nT ||
* || g(1,1,t) || -1545.05000 00000 || nT ||
* || g(2,0,t) || -2426.85000 00000 || nT ||
* || g(2,1,t) || 3015.10000 00000 || nT ||
* || g(2,2,t) || 1673.35000 00000 || nT ||
* || h(1,0,t) || 0.00000 00000 || nT ||
* || h(1,1,t) || 4879.65000 00000 || nT ||
* || h(2,0,t) || 0.00000 00000 || nT ||
* || h(2,1,t) || -2763.95000 00000 || nT ||
* || h(2,2,t) || -605.60000 00000 || nT ||
* || Xprime || 5478.08914 74225 || nT ||
* || Yprime || 14765.37032 43050 || nT ||
* || Zprime || -50632.17770 56324 || nT ||
* || Xprime-dot || 20.58517 51801 || nT/yr ||
* || Yprime-dot || 1.02725 92716 || nT/yr ||
* || Zprime-dot || 83.50809 72670 || nT/yr ||
* || X || 5535.52491 48687 || nT ||
* || Y || 14765.37032 43050 || nT ||
* || Z || -50625.93054 78794 || nT ||
* || Xdot || 20.49042 68023 || nT/yr ||
* || Ydot || 1.02725 92716 || nT/yr ||
* || Zdot || 83.53139 62281 || nT/yr ||
* || F || 53024.92848 40226 || nT ||
* || H || 15768.89967 29956 || nT ||
* || D || 1.21211 40681 || rad ||
* || I || -1.26884 21543 || rad ||
* || Fdot || -77.32705 44552 || nT/yr ||
* || Hdot || 8.15485 76192 || nT/yr ||
* || Ddot || -0.00119 38570 || rad/yr ||
* || Idot || 0.00061 53148 || rad/yr ||
* ||----------------------------------------------||
*)
println "Compare these values with the examples in the code comments."
end
geomagnetic(date:RR64, lat:RR64, lon:RR64, alt:RR64)
= do
(*
* Reference:
* [1] "Section 1.2: RELEVANT MODEL EQUATIONS",
* The US/UK World Magnetic Model for 2010-2015, Technical Report,
* http://www.ngdc.noaa.gov/geomag/WMM/data/WMM2010/WMM2010_Report.pdf
*)
(* input variables *)
println (date,lat,lon,alt) (*) in years, degrees and kilometers
(* (lat,lon,alt) is geodetic coordinates on WGS84 ellipsoid *)
(* some useful constants *)
pi = acos(-1.0)
(* read the magnetic model coefficients from file *)
(nmax,gnm,hnm,gnm_dot,hnm_dot,a,t0) = readMagneticModel()
(* Set up functions that interpolate coefficients *)
G_nm(n:ZZ32,m:ZZ32,t:RR64):RR64
requires { 0 <= m <= n <= nmax }
= do
k = PlmIndex(n,m)
gnm[k] + (t - t0) gnm_dot[k]
end
H_nm(n:ZZ32,m:ZZ32,t:RR64):RR64
requires { 0 <= m <= n <= nmax }
= do
k = PlmIndex(n,m)
hnm[k] + (t - t0) hnm_dot[k]
end
G_dot_nm(n:ZZ32,m:ZZ32,t:RR64):RR64
requires { 0 <= m <= n <= nmax }
= gnm_dot[PlmIndex(n,m)]
H_dot_nm(n:ZZ32,m:ZZ32,t:RR64):RR64
requires { 0 <= m <= n <= nmax }
= hnm_dot[PlmIndex(n,m)]
(* convert input variables to correct units *)
(t, phi,lamda,h) = (date, lat pi/180, lon pi/180, alt 1000) (*) in years, radians, radians and meters
(* (phi,lamda,h) is geodetic coordinates on WGS84 ellipsoid *)
println "(t,lamda,phi,h) = "(t,lamda,phi,h)
(* Step 1: convert to spherical geocentric coordinates *)
(* semi-major axis (A) of the earth *)
A = 6378137 (*) meters
(* reciprocal flattening (1/f) of the earth *)
f = 1/298.257223563
(* eccentricity squared (e^2) *)
e2 = f(2-f)
(* radius of East-West curvature (R_c) *)
R_c = A/(SQRT(1 - e2 (sin phi)^2))
p = (R_c + h) cos phi
z = (R_c(1 - e2) + h) sin phi
r = SQRT(p^2 + z^2)
phi_prime = asin(z/r)
(* geodetic coordinates: (lamda, phi, h) = (longitude, latitude, altitude) *)
(* spherical geocentric coordinates: (lamda, phi_prime, r) *)
println "(lamda,phi-prime,r) = "(lamda,phi_prime,r)
(* show some coefficients for comparison with example values *)
println "Gnm: "(G_nm(1,0,t), G_nm(1,1,t), G_nm(2,0,t), G_nm(2,1,t), G_nm(2,2,t))
println "Hnm: "(H_nm(1,0,t), H_nm(1,1,t), H_nm(2,0,t), H_nm(2,1,t), H_nm(2,2,t))
(* Step 2: Calculate the Schmidt normalized associated Legendre values *)
(pp,qq) = PlmSchmidt_d1(nmax, sin phi_prime)
P_breve_nm(n:ZZ32,m:ZZ32,mu:RR64)
requires { 0 <= m <= n <= nmax, |mu| <= 1.0 }
= pp[PlmIndex(n,m)]
dP_breve_nm(n:ZZ32,m:ZZ32,mu:RR64)
requires { 0 <= m <= n <= nmax, |mu| <= 1.0 }
= qq[PlmIndex(n,m)]
(* Step 3: Compute the field vector components (equations 10,11,12) *)
N = nmax
X_prime = -(SUM[n<-1:N]( (a/r)^(n+2) (SUM[m<-0:n]( (G_nm(n,m,t) cos(m lamda) + H_nm(n,m,t) sin(m
lamda)) dP_breve_nm(n,m,sin phi_prime) ))))
Y_prime = (1/cos(phi_prime)) (SUM[n<-1:N]( (a/r)^(n+2) (SUM[m<-0:n](m((G_nm(n,m,t) sin(m lamda)) -
(H_nm(n,m,t) cos(m lamda))) P_breve_nm(n,m,sin phi_prime) ))))
Z_prime = -(SUM[n<-1:N]( (n+1) (a/r)^(n+2) (SUM[m<-0:n]( (G_nm(n,m,t) cos(m lamda) + H_nm(n,m,t)
sin(m lamda)) P_breve_nm(n,m,sin phi_prime) ))))
println "(Xprime,Yprime,Zprime) = "(X_prime, Y_prime, Z_prime)
(* Compute the secular variation of the field components (equastions 13,14,15) *)
X_dot_prime = -(SUM[n<-1:N]( (a/r)^(n+2) (SUM[m<-0:n]( (G_dot_nm(n,m,t) cos(m lamda) +
H_dot_nm(n,m,t) sin(m lamda)) dP_breve_nm(n,m,sin phi_prime) ))))
Y_dot_prime = 1/cos(phi_prime) (SUM[n<-1:N]( (a/r)^(n+2) (SUM[m<-0:n]( m (G_dot_nm(n,m,t) sin(m
lamda) - H_dot_nm(n,m,t) cos(m lamda)) P_breve_nm(n,m,sin phi_prime) ))))
Z_dot_prime = -(SUM[n<-1:N]( (n+1) (a/r)^(n+2) (SUM[m<-0:n]( (G_dot_nm(n,m,t) cos(m lamda) +
H_dot_nm(n,m,t) sin(m lamda)) P_breve_nm(n,m,sin phi_prime) ))))
println "(Xprime-dot,Yprime-dot,Zprime-dot) = "(X_dot_prime, Y_dot_prime, Z_dot_prime)
(* Step 4: rotate geocentric magnetic field vector components into ellipsoidal reference frame *)
X = X_prime cos(phi_prime - phi) - Z_prime sin(phi_prime - phi) (*) nT
Y = Y_prime (*) nT
Z = X_prime sin(phi_prime - phi) + Z_prime cos(phi_prime - phi) (*) nT
println "(X,Y,Z) = "(X,Y,Z)
X_dot = X_dot_prime cos(phi_prime - phi) - Z_dot_prime sin(phi_prime - phi) (*) nT per year
Y_dot = Y_dot_prime (*) nT per year
Z_dot = X_dot_prime sin(phi_prime - phi) + Z_dot_prime cos(phi_prime - phi) (*) nT per year
println "(Xdot,Ydot,Zdot) = "(X_dot,Y_dot,Z_dot)
(* Step 5: Compute magnetic elements H,F,I and D from orthogonal componants *)
F = SQRT(X^2 + Y^2 + Z^2) (*) nT
H = SQRT(X^2 + Y^2) (*) nT
I = atan2(Z,H) (*) radians
D = atan2(Y,X) (*) radians
println "(F,H,D,I) = "(F,H,D,I)
F_dot = (X DOT X_dot + Y DOT Y_dot + Z DOT Z_dot)/F (*) nT per year
H_dot = (X DOT X_dot + Y DOT Y_dot)/H (*) nT per year
I_dot = (H DOT Z_dot - Z DOT H_dot)/(F^2) (*) radians per year
D_dot = (X DOT Y_dot - Y DOT X_dot)/(H^2) (*) radians per year
println "(Fdot,Hdot,Ddot,Idot) = " (F_dot,H_dot,D_dot,I_dot)
(X,Y,Z, H,I,D, F)
end
readMagneticModel(): (ZZ32,Array[\RR64,ZZ32\],Array[\RR64,ZZ32\],Array[\RR64,ZZ32\],Array[\RR64,ZZ32\],RR64)
= do
(*
* Reference:
* [1] "Section 1.3: THE WMM2010 COEFFICIENTS",
* The US/UK World Magnetic Model for 2010-2015, Technical Report,
* http://www.ngdc.noaa.gov/geomag/WMM/data/WMM2010/WMM2010_Report.pdf
*)
(*) Coefficients as read from WMM.COF file
(* n m g h g_dot h_dot *)
wmm2010data: RR64[90,6] = [
1 0 -29496.6 0.0 11.6 0.0
1 1 -1586.3 4944.4 16.5 -25.9
2 0 -2396.6 0.0 -12.1 0.0
2 1 3026.1 -2707.7 -4.4 -22.5
2 2 1668.6 -576.1 1.9 -11.8
3 0 1340.1 0.0 0.4 0.0
3 1 -2326.2 -160.2 -4.1 7.3
3 2 1231.9 251.9 -2.9 -3.9
3 3 634.0 -536.6 -7.7 -2.6
4 0 912.6 0.0 -1.8 0.0
4 1 808.9 286.4 2.3 1.1
4 2 166.7 -211.2 -8.7 2.7
4 3 -357.1 164.3 4.6 3.9
4 4 89.4 -309.1 -2.1 -0.8
5 0 -230.9 0.0 -1.0 0.0
5 1 357.2 44.6 0.6 0.4
5 2 200.3 188.9 -1.8 1.8
5 3 -141.1 -118.2 -1.0 1.2
5 4 -163.0 0.0 0.9 4.0
5 5 -7.8 100.9 1.0 -0.6
6 0 72.8 0.0 -0.2 0.0
6 1 68.6 -20.8 -0.2 -0.2
6 2 76.0 44.1 -0.1 -2.1
6 3 -141.4 61.5 2.0 -0.4
6 4 -22.8 -66.3 -1.7 -0.6
6 5 13.2 3.1 -0.3 0.5
6 6 -77.9 55.0 1.7 0.9
7 0 80.5 0.0 0.1 0.0
7 1 -75.1 -57.9 -0.1 0.7
7 2 -4.7 -21.1 -0.6 0.3
7 3 45.3 6.5 1.3 -0.1
7 4 13.9 24.9 0.4 -0.1
7 5 10.4 7.0 0.3 -0.8
7 6 1.7 -27.7 -0.7 -0.3
7 7 4.9 -3.3 0.6 0.3
8 0 24.4 0.0 -0.1 0.0
8 1 8.1 11.0 0.1 -0.1
8 2 -14.5 -20.0 -0.6 0.2
8 3 -5.6 11.9 0.2 0.4
8 4 -19.3 -17.4 -0.2 0.4
8 5 11.5 16.7 0.3 0.1
8 6 10.9 7.0 0.3 -0.1
8 7 -14.1 -10.8 -0.6 0.4
8 8 -3.7 1.7 0.2 0.3
9 0 5.4 0.0 0.0 0.0
9 1 9.4 -20.5 -0.1 0.0
9 2 3.4 11.5 0.0 -0.2
9 3 -5.2 12.8 0.3 0.0
9 4 3.1 -7.2 -0.4 -0.1
9 5 -12.4 -7.4 -0.3 0.1
9 6 -0.7 8.0 0.1 0.0
9 7 8.4 2.1 -0.1 -0.2
9 8 -8.5 -6.1 -0.4 0.3
9 9 -10.1 7.0 -0.2 0.2
10 0 -2.0 0.0 0.0 0.0
10 1 -6.3 2.8 0.0 0.1
10 2 0.9 -0.1 -0.1 -0.1
10 3 -1.1 4.7 0.2 0.0
10 4 -0.2 4.4 0.0 -0.1
10 5 2.5 -7.2 -0.1 -0.1
10 6 -0.3 -1.0 -0.2 0.0
10 7 2.2 -3.9 0.0 -0.1
10 8 3.1 -2.0 -0.1 -0.2
10 9 -1.0 -2.0 -0.2 0.0
10 10 -2.8 -8.3 -0.2 -0.1
11 0 3.0 0.0 0.0 0.0
11 1 -1.5 0.2 0.0 0.0
11 2 -2.1 1.7 0.0 0.1
11 3 1.7 -0.6 0.1 0.0
11 4 -0.5 -1.8 0.0 0.1
11 5 0.5 0.9 0.0 0.0
11 6 -0.8 -0.4 0.0 0.1
11 7 0.4 -2.5 0.0 0.0
11 8 1.8 -1.3 0.0 -0.1
11 9 0.1 -2.1 0.0 -0.1
11 10 0.7 -1.9 -0.1 0.0
11 11 3.8 -1.8 0.0 -0.1
12 0 -2.2 0.0 0.0 0.0
12 1 -0.2 -0.9 0.0 0.0
12 2 0.3 0.3 0.1 0.0
12 3 1.0 2.1 0.1 0.0
12 4 -0.6 -2.5 -0.1 0.0
12 5 0.9 0.5 0.0 0.0
12 6 -0.1 0.6 0.0 0.1
12 7 0.5 0.0 0.0 0.0
12 8 -0.4 0.1 0.0 0.0
12 9 -0.4 0.3 0.0 0.0
12 10 0.2 -0.9 0.0 0.0
12 11 -0.8 -0.2 -0.1 0.0
12 12 0.0 0.9 0.1 0.0
]
(* Note: H and H_dot is undefined for m = 0 *)
nmax = 12
sdim = PlmIndex(nmax+1,0)
G = array[\RR64\](sdim)
H = array[\RR64\](sdim)
G_dot = array[\RR64\](sdim)
H_dot = array[\RR64\](sdim)
(* Not used, but there to fill our packed triangle matrix *)
G[0] := 0.0
G[0] := 0.0
G_dot[0] := 0.0
H_dot[0] := 0.0
(* populate our packed triangle matrix *)
for row <- 0#90 do
k = PlmIndex( wmm2010data[row,0], wmm2010data[row,1] )
G[k] := wmm2010data[row,2]
H[k] := wmm2010data[row,3]
G_dot[k] := wmm2010data[row,4]
H_dot[k] := wmm2010data[row,5]
end
(* geomagnetic reference radius of the model *)
a = 6371200 (*) meters
(* base date of the model *)
t0 = 2010.0 (*) years
(nmax, G, H, G_dot, H_dot, a, t0)
end
PlmIndex(l: ZZ32, m: ZZ32): ZZ32
requires { 0 <= m <= l }
= l(l+1)/2 + m
PlmSchmidt_d1(l_max: ZZ32, z: RR64): (Array[\RR64,ZZ32\], Array[\RR64,ZZ32\])
requires { 0 <= l_max <= 16, |z| <= 1 }
= do
p_size: ZZ32 = PlmIndex(l_max+1, 0)
var p: Array[\RR64,ZZ32\]
var dp: Array[\RR64,ZZ32\]
var schmidtQuasiNorm: Array[\RR64,ZZ32\]
do
mu = SQRT((1.0-z)(1.0+z))
(* Compute the Gauss-normalized associated Legendre functions *)
p := array[\RR64\](p_size)
dp := array[\RR64\](p_size)
p[0] := 1.0
dp[0] := 0.0
for n <- seq(1:l_max), m <- seq(0:n) do
index = PlmIndex(n,m)
if n = m then
index1 = PlmIndex(n-1,m-1)
p[index] := mu p[index1]
dp[index] := mu dp[index1] + z p[index1]
elif n = 1 AND m = 0 then
index1 = PlmIndex(n-1,m)
p[index] := z p[index1]
dp[index] := z dp[index1] - mu p[index1]
elif n > 1 AND n =/= m then
index2 = PlmIndex(n-1,m)
if m > (n-2) then
p[index] := z p[index2]
dp[index] := z dp[index2] - mu p[index2]
else
index1 = PlmIndex(n-2,m)
k = ((n-1)^2 - m^2) / ((2 n - 1)(2 n - 3))
p[index] := z p[index2] - k p[index1]
dp[index] := z dp[index2] - mu p[index2] - k dp[index1]
end
else
fail "should never happen"
end
end
also do
schmidtQuasiNorm := array[\RR64\](p_size)
(*
* Compute the ration between the Gauss-normalized associated Legendre
* functions and the Schmidt quasi-normalized version. This is equivalent
* to SQRT((m==0?1:2)*(n-m)!/(n+m!))*(2n-1)!!/(n-m)!
*)
schmidtQuasiNorm[0] := 1.0
for n <- seq(1:l_max), m <- seq(0:n) do
schmidtQuasiNorm[PlmIndex(n,m)] :=
if m = 0 then
schmidtQuasiNorm[PlmIndex(n-1,0)] (2 n - 1) / n
else
schmidtQuasiNorm[PlmIndex(n,m-1)] (SQRT(((n - m + 1)(if m = 1 then 2 else 1 end)) / (n + m)))
end
end
end
(*
* Converts the Gauss-normalized associated Legendre
* functions to the Schmidt quasi-normalized version using pre-computed
* relation stored in the variable schmidtQuasiNorm
*)
for n <- 1:l_max, m <- 0:n do
index = PlmIndex(n,m)
p[index] := schmidtQuasiNorm[index] p[index]
dp[index] := - schmidtQuasiNorm[index] dp[index]
(* The sign is changed since the new WMM routines
* use derivative with respect to latitude instead
* of co-latitude
*)
end
(p,dp)
end
end (* component geomag *)
_______________________________________________
language mailing list
language@...
http://projectfortress.sun.com/mailman/listinfo/language
RSS Feed