[paparazzi-commits] [4195] EGM96 model from 15'x15' file

[Pascal Brisset]

Revision: 4195
          http://svn.sv.gnu.org/viewvc/?view=rev&root=paparazzi&revision=4195
Author:   hecto
Date:     2009-09-21 19:22:28 +0000 (Mon, 21 Sep 2009)
Log Message:
-----------
 EGM96 model from 15'x15' file

Modified Paths:
--------------
    paparazzi3/trunk/sw/lib/ocaml/latlong.ml
    paparazzi3/trunk/sw/lib/ocaml/latlong.mli

Added Paths:
-----------
    paparazzi3/trunk/data/srtm/WW15MGH.DAC.bz2

Added: paparazzi3/trunk/data/srtm/WW15MGH.DAC.bz2
===================================================================
(Binary files differ)


Property changes on: paparazzi3/trunk/data/srtm/WW15MGH.DAC.bz2
___________________________________________________________________
Added: svn:mime-type
   + application/octet-stream

Modified: paparazzi3/trunk/sw/lib/ocaml/latlong.ml
===================================================================
--- paparazzi3/trunk/sw/lib/ocaml/latlong.ml    2009-09-20 16:20:07 UTC (rev 
4194)
+++ paparazzi3/trunk/sw/lib/ocaml/latlong.ml    2009-09-21 19:22:28 UTC (rev 
4195)
@@ -24,6 +24,8 @@
  *
  *)
 
+let (//) = Filename.concat
+
 module C = struct
   include Complex
   let make x y = {re = x; im = y}
@@ -712,3 +714,32 @@
     (float geoid_data.(ilat1).(ilon2))
     (float geoid_data.(ilat2).(ilon1)) 
     (float geoid_data.(ilat2).(ilon2))
+
+
+
+let egm96_ncols = 1440
+let egm96_nrows = 721
+let egm96_data =
+  lazy (
+    let path = [Env.paparazzi_home // "data" // "srtm"] in
+    let f = Ocaml_tools.open_compress (Ocaml_tools.find_file path 
"WW15MGH.DAC") in
+    let n = egm96_ncols * egm96_nrows * 2 in
+    let buf = String.create n in
+    really_input f buf 0 n;
+    buf)
+
+
+(* 
http://earth-info.nima.mil/GandG/wgs84/gravitymod/egm96/binary/binarygeoid.html 
*)
+let egm96 = fun geo ->
+  let egm96_data = Lazy.force egm96_data in
+
+  let lat = truncate ((Rad>>Deg) (norm_angle geo.posn_lat))
+  and lon = truncate ((Rad>>Deg) (norm_angle geo.posn_long)) in
+  let ilat = (90-lat)*4      (* 15' == 4 entries per degree *)
+  and ilon = (lon+if lon < 0 then 360 else 0)*4 in
+
+  let i = (2*(ilat*egm96_ncols+ilon)) in
+
+  let x = Char.code egm96_data.[i] lsl 8 + Char.code egm96_data.[i+1] in
+
+  float ((x lsl 16) asr 16) /. 100.

Modified: paparazzi3/trunk/sw/lib/ocaml/latlong.mli
===================================================================
--- paparazzi3/trunk/sw/lib/ocaml/latlong.mli   2009-09-20 16:20:07 UTC (rev 
4194)
+++ paparazzi3/trunk/sw/lib/ocaml/latlong.mli   2009-09-21 19:22:28 UTC (rev 
4195)
@@ -199,3 +199,6 @@
 (** [geo_of_ecef ellipsoid geo ecef] Returns llh *)
 
 val wgs84_hmsl : geographic -> fmeter
+
+val egm96 :  geographic -> fmeter
+(** Return geoid height from 15' precomputed file *)