;+
; NAME:
;	HammerAitoff
; PURPOSE:
;	Transforms latitude and longitude to X and Y coordinates in the
;	Hammer-Aitoff map projection.
; CALLING SEQUENCE:
	FUNCTION HammerAitoff, Pos	, $
		degrees 	= degrees	, $
		dabg		= dabg		, $
		zero_phase	= zero_phase
; INPUTS:
;	Pos 		array[2,n]; type: float
;					n locations in the sky
;					Pos[0,n]	longitudes
;					Pos[1,n]	latitudes in [-90.,90.]
; OPTIONAL INPUT PARAMETERS:
;	/degrees	if set then all angles are in degrees (default: radians)
;
;	dabg		array[3]; type: float; default: none
;					Euler angles for a rotation to be applied to
;					'Pos' prior to the Hammer-Aitoff projection.
;					(see href=FishEye=)
;	zero_phase	scalar; type: float; default: none
;					Additional offset applied to phase angle
;					(see href=FishEye=)
; OUTPUTS:
;	Pos			array[2,n]; type: float
;					X,Y	Cartesian coordinates (same structure as input array)
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
; CALLS:
;	IsType, AngleRange, ToRadians, EulerRotate, SyncDims
; RESTRICTIONS:
; PROCEDURE:
;	The Hammer-Aitoff projection is the format used at NOAA SEL; for
;	details, see the book Introduction to Map Projections by Porter W.
;	McDonnell Jr.
;
;	Pos[0,*] > 0 has X > 0	Pos[1,*] > 0 has Y > 0
;	Pos[0,*] < 0 has X < 0	Pos[1,*] < 0 has Y < 0
;
;	The returned coordinates will be in the range -sqrt(2),+sqrt(2) in the vertical
;	(latitudinal) direction and -2*sqrt(2),+2*sqrt(2) in the horizontal (longitudinal
;	direction).
; MODIFICATION HISTORY:
;	NOV-1992, Tom Davidson (UCSD)
;	MAR-1999, Paul Hick (UCSD/CASS); added /degrees keyword)
;	JAN-2002, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;		Added 'dabg' keyword
;-

rpm = ToRadians(degrees=Degrees)

Out = Pos

sz  = size(Out)
Out = reform(Out,sz[1],sz[sz[0]+2]/sz[1], /overwrite)

IF IsType(dabg, /defined) THEN	$
	IF IsType(a, /undefined) THEN a = dabg ELSE a += dabg

IF IsType(zero_phase, /defined) THEN	$
	IF IsType(a, /undefined) THEN a = [zero_phase,0,0] ELSE a += [zero_phase,0,0]

IF IsType(a, /defined) THEN BEGIN
	Out = EulerRotate(a, Out, degrees=degrees)
	Out[0,*] = AngleRange(reform(Out[0,*]), /pi, degrees=degrees)
ENDIF

Out *= rpm

Out[0,*] /= 2.0

K = sqrt(2.0/(1.0+cos(Out[0,*])*cos(Out[1,*])))
Out = [ 2*K*cos(Out[1,*])*sin(Out[0,*]) , K*sin(Out[1,*]) ]

SyncDims, Out, size=sz

RETURN, Out  &  END