From mboxrd@z Thu Jan  1 00:00:00 1970
X-Spam-Checker-Version: SpamAssassin 3.4.4 (2020-01-24) on polar.synack.me
X-Spam-Level: 
X-Spam-Status: No, score=0.1 required=5.0 tests=BAYES_05,INVALID_MSGID
	autolearn=no autolearn_force=no version=3.4.4
X-Google-Language: ENGLISH,ASCII-7-bit
X-Google-Thread: 103376,960d1cb75613cad
X-Google-Attributes: gid103376,public
From: Geoff Bull <gbull@acenet.com.au>
Subject: Re: Earth map transformation in Ada ?
Date: 2000/05/07
Message-ID: <3914B7DE.3CCEF170@acenet.com.au>#1/1
X-Deja-AN: 620125911
Content-Transfer-Encoding: 7bit
References: <8eudi9$ck7$1@news.uit.no>
X-Accept-Language: en
Content-Type: text/plain; charset=us-ascii
X-Complaints-To: abuse@telstra.net
X-Trace: nsw.nnrp.telstra.net 957659289 203.35.118.1 (Sun,
 07 May 2000 10:28:09 EST)
Organization: Customer of Telstra Big Pond Direct
MIME-Version: 1.0
NNTP-Posting-Date: Sun, 07 May 2000 10:28:09 EST
Newsgroups: comp.lang.ada
Date: 2000-05-07T00:00:00+00:00
List-Id: <comp.lang.ada>



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;


...