Re: Earth map transformation in Ada ?

[Geoff Bull <gbull@acenet.com.au>]

Reinert Korsnes wrote:
> 
> Does someone know about available (Earth) map stransformations in Ada ?
> 
> like:
> 
> lat,long  <-> x,y (Polar Stereographic)
> 
> reinert

I have a package for doing wgs84 calculations, which you may find useful
as a starting point. I can't remember wehther it produces 
correct answers. Here is an excerpt, I can mail you the whole
package if you like:

generic
   type Float_Type is digits <>;
package WGS84_Package is

     --  first  numerical eccentricity
      e1sqr : constant := (a ** 2 - b ** 2) / (a ** 2) ;

      --  second numerical eccentricity
      e2sqr : constant := (a ** 2 - b ** 2) / (b ** 2);

...
   subtype Metres is Float_Type;
   subtype Radians is Float_Type;
   subtype Degrees is Float_Type;

   type ECEF is -- earth centreed earth fixed rectangular coordinate
system
      record
         X : Metres;  --  though equator at prime meridian
         Y : Metres;  --  though equator 90 degrees east of prime
meridian
         Z : Metres;  --  through north pole
      end record;

 
   type Polar is
      record
         Latitude : Radians;
         Longitude : Radians;
         Altitude : Metres;
      end record;

   type Azimuth_Elevation is
      record
         Azimuth : Radians;
         Elevation : Radians;
         Distance : Metres;
      end record;

...

package body WGS84_Package is

    --  converts Latitude, Longitude and Altitude (above ellipsoid) to
XYZ
   function To_ECEF (From : in Polar) return ECEF is
      N : Float_Type := a / sqrt (1.0 - e1sqr * sin (From.Latitude) **
2);
   begin -- To_ECEF
      return
        (X => (N + From.Altitude) * cos (From.Latitude) * cos
(From.Longitude),

         Y => (N + From.Altitude) * cos (From.Latitude) * sin
(From.Longitude),

         Z => (N * (1.0 - e1sqr) + From.Altitude) * sin (From.Latitude)
         );
   end To_ECEF;


   --  converts XYZ to Latitude, Longitude, Altitude (above ellipsoid)
   function To_Polar (From : in ECEF) return Polar is
      -- length projected onto equatorial plane;
      P : Metres := sqrt (From.X ** 2 + From.Y ** 2);

      T : Float_Type := arctan ((From.Z * a) / (p * b));


      N : Float_Type;

      Latitude, Longitude : Radians;
      Altitude : Metres;
   begin -- To_Polar
      Latitude := arctan ((From.Z + e2sqr * b * Sin (T) ** 3)
                      / (P - e1sqr * a * Cos (T) ** 3));

      if From.X = 0.0 then
         Longitude := Pi / 2.0;
      else
         Longitude := arctan (From.Y / From.X);
      end if;

      N :=  A / Sqrt (1.0 - e1sqr * sin (Latitude) ** 2);
      Altitude := P / cos (Latitude) - N;

      return (Latitude, Longitude, Altitude);

   end To_Polar;


...