;+
; NAME:
;	FishEye
; PURPOSE:
;	Converts spherical angles to coordinates for a 'fish-eye' plot
; CATEGORY:
;	Math: projections
; CALLING SEQUENCE:
	FUNCTION FishEye, Pos		, $
		degrees 	= degrees	, $
		maxelo		= maxelo	, $
		polar		= polar 	, $
		dabg		= dabg		, $
		zero_phase	= zero_phase, $
		inverse 	= inverse
; INPUTS:
;	Pos				array[2,*]; type: float
;						pos[0,*]: phase angle, longitudes
;						pos[1,*]: latitude
; OPTIONAL INPUT PARAMETERS:
;	/degrees		if set, all angles are in degrees (default: radians)
;	maxelo=maxelo	scalar; type: float; default: !pi
;						elongation 'maxelo' is scaled to 2*sqrt(2)
;	dabg=dabg		array[3]; type: float; default: [0,0,0]
;						Euler angles to be added to the rotation
;						discussed in PROCEDURE.
;	zero_phase		scalar; type: float; default: none
;						Additional offset applied to phase angle
;						i.e. adds Euler triplet [zero_phase,0,0] to
;						the rotation discussed in PROCEDURE.
;	/inverse		if set, the inverse transformation (from (x,y) in the
;					fisheye to phase angle/lat) is done.
; OUTPUTS:
;	Result			array[2,*]; type: float
;						x and y coordinates in fish-eye plot
; OPTIONAL OUTPUT PARAMETERS:
;	polar=polar		array[2,*]; type: float
;						polar coordinates in fish-eye plot
; INCLUDE:
	@compile_opt.pro			; On error, return to caller
; CALLS:
;	IsType, ToRadians, EulerRotate, SyncDims
; PROCEDURE:
; >	The input angles are usually celestial spherical coordinates,
;	e.g. ecliptic longitude and latitude, or right ascension and
;	declination.
; >	The fish-eye plot uses phase and polar angles relative to a
;	coordinate system with the z-axis along longitude=0,latitude=0,
;	and the x-axis along longitude=90, latitude=0. The transformation
;	to this coordinate system takes Euler angles E=[0,90,90] degrees.
; >	The Euler angle triplet specified in keyword dabg will be
;	added to E, i.e. the total rotation applied to the input
;	coordinates is E+dabg.
;	dabg is used to change the direction for the center of the
;	fish-eye plot. E.g. if the input angles are ecliptic coordinates
;	and the ecliptic coordinates of the Sun are [lng,dec], then
;	dabg = [lng,-dec,0] would put the Sun in the center of the
;	fish-eye (the total rotation would be [lng,90-dec,90]).
;	Note that in this example the 3rd angle in dabg (set to zero)
;	rotates the plot around the center of the fisheye.
; >	In the fish-eye coordinate system the following transformation is used:
;		x = polar*cos( phase )*(2*sqrt(2)/maxelo)
;		y = polar*sin( phase )*(2*sqrt(2)/maxelo)
;	i.e. the maximum polar angle 'maxelo' is scaled to 2*sqrt(2)
; >	 scale -sqrt(2),sqrt(2)
;		X-axis: east (left)
;		Y-axis: north (up)
; MODIFICATION HISTORY:
;	SEP-1999, Paul Hick (UCSD/CASS)
;	JUL-2008, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;		Added keyword /inverse
;-

InitVar, inverse, /key

rpm = ToRadians(degrees=Degrees)
pi  = !dpi/rpm
pi2 = pi/2.0d0

IF IsType(maxelo, /undefined) THEN edge = pi ELSE edge = maxelo
edge /= 2.0*sqrt(2.0) 			; Scale edge to 2*sqrt(2)

a = [0.0d0,pi2,pi2]
IF IsType(dabg      , /defined) THEN a    += dabg
IF IsType(zero_phase, /defined) THEN a[0] += zero_phase

sz  = size(Pos)

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

CASE inverse OF

0: BEGIN

	Out = EulerRotate(a, Out, degrees=degrees)

	Out[1,*] = pi2-Out[1,*]			; Convert to colatitude

	IF arg_present(polar) THEN polar = Out

	Out[0,*] *= rpm					; Phase angle to radians
	Out  = [Out[1,*]*cos(Out[0,*]),Out[1,*]*sin(Out[0,*])]
	Out /= edge 			; Scale to 2*sqrt(2)

END

1: BEGIN

	Out = [ atan(Out[1,*],Out[0,*])/rpm, reform(sqrt(total(Out*Out,1)),size(Out[1,*],/dim))*edge ]

	IF arg_present(polar) THEN polar = Out

	Out[1,*] = pi2-Out[1,*]			; Convert to latitude

	Out = EulerRotate(-reverse(a), Out, degrees=degrees)
	Out[0,*] = AngleRange(Out[0,*],degrees=degrees)
END

ENDCASE

SyncDims, Out, size=sz
IF IsType(polar, /defined) THEN SyncDims, polar, size=sz

RETURN, Out  &  END