From: turner@bigbang.astro.indiana.edu (George Wm Turner)
Subject: Re: Formulae for projections of sky maps?
Reply-To: turner@bigbang.astro.indiana.edu (George Wm Turner)
Organization: Indiana University
Date: Sun, 16 Aug 1992 19:04:41 GMT


>In article <16ka47INN5u6@agate.berkeley.edu>, ted@physics1 (Emory F. Bunn) writes:
>Can anybody point me towards a reference that lists formulae for
>the various projections of full-sky maps to two dimensions?
>I don't even know what the different projections are called, but
>what I'm looking for are the two or three most common ones.
>
>Thanks a lot,
>
>-Ted
>(ted@physics.berkeley.edu)

any geography book on maps will have the information you seek.

there are usually three mapping problems you come across. mapping
small areas, larger areas up to half spheres and full 4-pi steradian
projections of the full sphere.  the three i use are the gnomonic for
small areas on a sphere (such as a ccd images of the sky), the
oblique-azimuthal equalarea for areas up to a hemisphere and the 
hammer-aitoff projection of the full sphere. 

in the code below, : its fortran, vax fortran at that. 
		     ra/ha = long  and  dec = lat
		     no guarantees, you get what you pay for

the gnomonic projection is just a flat plane tangent to the
sphere; you give it the ra,dec (long,lat in geographical coordinates)
of the tangent point (the center of your image/map) and the second coordinate
you want to map and you get out an x,y coordinate.

the oblique azimuthal equal area projection is good for hemispheres & less
e.g. the night sky. you pass in the ra,dec (long,lat) of the center of 
the helisphere and the point you want to project. the x,y coordinate is 
returned.  this function is valid for a full 4*pi steradian projection. 
the "inner" hemisphere (the 2*pi steradians closest to the center point) 
is mapped into a circle of  radius sqrt(2.).  the "outer" hemisphere 
(2*steradians farthest from the center point) is projected onto a 
ring surrounding the inner circle.  the inner radius being sqrt(2.) 
and the outer radius equal to 2.  its neat to project bright stars and
 see if you can identify constellations in this outer ring.

the hammer-aitoff code maps the full 4 pi steradians onto an ellipse. the
ellipse parameter are adjustable, see the code for further questions. this
is the common projection you see for ellipitical world maps. the hammer-aitoff
projection is an equal area projection.

   george wm turner         turner@bigbang.astro.indiana.edu

	note!!!!   ra/ha = long  and  dec = lat


******************************************************************************
* gnomonic projection good for small spherical areas (<< 2*pi steradians)
	program test gnomonic
        implicit none
        real s2,s22,pi,twopi,pibytwo,t,r,g,b,lon,lat,lat0,lon0
        character*80 str
* setup constants
        s2=sqrt(2.)
        s22=2.*s2
        pi=acos(-1.)
        twopi=2*pi
        pibytwo=pi*.5
* draw meridians                        ! go
        do lon=-pi*.1+lon0,pi*.1+lon0,pi*.1/8.
          do lat=-pibytwo*.1+lat0,pibytwo*.1+lat0,pibytwo*.1/8.
            call gnomonic(lat0,lon0,lat,lon,x,y)
            write(7,*)x,y
          end do
        end do
* draw parallels                        ! go
        do lat=-pibytwo*.1+lat0,pibytwo*.1+lat0,pibytwo*.1/8.
          do lon=-pi*.1+lon0,pi*.1+lon0,pi*.1/8.
            call gnomonic(lat0,lon0,lat,lon,x,y)
            write(7,*)x,y
          end do
        end do
* all done
        end
        subroutine gnomonic(lat0,lon0,lat,lon,x,y)
*         oblique gnomonic projection
*               lat0 : spherical latitude for central latitude
*               lat  : spherical latitude    -pi/2. : pi/2. (dec or geo-latitude)
*               lon  : spherical longatitude  -pi : pi     (ha or geo-longitude)
*               x,y  : planar cartesian output coordinates
*       george turner   1-may-1989      indiana university astronomy department
*       internet : turner@bigbang.astro.indiana.edu
        implicit none
        real lat0,lon0,lat,lon,x,y,a
        real r
        parameter (r=1.)
        a=r/(sin(lat0)*sin(lat)+cos(lat0)*cos(lat)*cos(lon-lon0))
        x=a*cos(lat)*sin(lon-lon0)
        y=a*(cos(lat0)*sin(lat)-sin(lat0)*cos(lat)*cos(lon-lon0))
        return
        end

******************************************************************************
* oblique azimuthal equal area projection  good for hemispheres & less.
	program test ob az equalarea
        implicit none
        real lat,lon,x,y,pi,lon0,lat0
        pi=acos(-1.)
        lat=60.*pi/180.
        lon=60.*pi/180.
        lat0=0.0
        lon0=40.*pi/180.
        do lon=-pi,pi,.1
        do lat=-pi/2.,pi/2.,.1
          call ob_azimuthal_equalarea(lat0,lon0,lat,lon,x,y)
        end do
        end do
        end
        subroutine ob_azimuthal_equalarea(lat0,lon0,lat,lon,x,y)
*       oblique azimuthal equal area projection  good for hemispheres & less
*               lat0,lon0 : central lat & long of projection
*               lat : spherical latitude
*               lon : spherical longatitude
*               x,y : planar cartesian output coordinates
*               x,y  : planar cartesian output coordinates
*       george turner   30-apr-1989     indiana university astronomy department
*       internet : turner@bigbang.astro.indiana.edu
        implicit none
        real x,y,lat,lon,lat0,lon0
        real r
        real z,alpha,dlambda
        parameter (r=318.55)
        dlambda=lon0-lon
        z=acos(sin(lat0)*sin(lat)+cos(lat0)*cos(lat)*cos(dlambda))
        alpha=asin(cos(lat0)*sin(dlambda)*(1./sin(z)))
        x=2.*r*sin(z*.5)*sin(alpha)
        y=2.*r*sin(z*.5)*cos(alpha)
        return
        end
*****************************************************************************
* hammer aitoff projection good for full sphere mapping 
        program test hammer aitoff projection
        implicit none
        real lat,lon,x,y,pi
        pi=acos(-1.)
        do lon=-pi,pi,.1
        do lat=-pi/2.,pi/2.,.1
          call hammer_aitoff(lat,lon,x,y)
        end do
        end do
        end
        subroutine hammer_aitoff(lat,lon,x,y)
*         modified hammer-aitoff projection
*               lat : spherical latitude    -pi/2. : pi/2. (dec or geo-latitude)
*               lon : spherical longatitude  -pi : pi      (ha or geo-longitude)
*               x,y : planar cartesian output coordinates

*       notes:
*       m=sqrt(a/(2*b)) where a & b are ellipse parameters of the
*           projection ellipse.
*       if m=1 then the projection is the original hammer mapping
*       if m=1 then remove m from the equations to avoid an unnessary multiply
*       george turner   1-may-1989      indiana university astronomy department
*       internet : turner@bigbang.astro.indiana.edu
        implicit none
        real x,y,lat,lon,t
        real m
*       parameter (m=0.7071067811865475)  a = b = 1  (map into a circle)
        parameter (m=1.0)
        t=sqrt(1.+cos(lat)*cos(lon*.5))
        x=2.*sqrt(2.)*m*cos(lat)*sin(lon*.5)/t
        y=sqrt(2.)*sin(lat)/(m*t)
        return
        end
Newsgroups: sci.astro