;+

; PROCEDURE:	transform_mag_geo

;

; PURPOSE:	Convert lattitudes and longitudes between MAG and GEO 

;		coordinate systems using the eccentric dipole model.

;		Transformation is from MAG to GEO unless INVERSE is set.

;		Transformation is performed at <Re> + 100km = 6472.1 km

;

; INPUTS:

;

;   LAT		The input lattitude array. (Not altered)

;   LNG		The input longitude array. (Not altered)

;   TLAT	Named variable to accept transformed lattitudes.

;   TLNG	Named variable to accept transformed longitudes.

;

; KEYWORDS:

;

;   YEAR	The 4-digit year.

;   DEGREES	If nonzero, input and ouput in degrees.

;   INVERSE	If nonzero, tranformation is from GEO to MAG.

;

; NOTES:

;

;   MAG --> GEO:

;   When transforming from MAG to GEO, LAT and LNG are assumed to

;   refer to points Re=6372.1 km from the earth's dipole center.

;   (This does not necessarily place them on the surface of the earth.)

;   The points are propagated to <Re> + 100km and then their

;   coordinates are converted to the GEO system through a rotation

;   followed by a translation.

;

;   This procedure is inaccurate for points where field lines are nearly

;   orthogonal to lines through the dipole center, i.e. near the equator.

;

;   GEO --> MAG:

;   When transforming from GEO to MAG, LAT and LNG are assumed to

;   refer to points Re + 100km from the center of the earth.

;   Since this one is simply the inverse of the above transformation, we

;   do the inverse translation followed by the inverse rotation.

;   Right now, the GEO to MAG transformation is less accurate than its

;   inverse; it does not propagate along field lines during the

;   translation to the dipole center.

; 

; 

; WRITTEN BY:	J. Rauchleiba	4/14/97

; ALGORITHM BY: M. Temerin, G. Kaplan, J. Rauchleiba

;

;-

pro transform_mag_geo, lat, lng, tlat, tlng, $

	YEAR=year, $

	DEGREES=deg, $

	INVERSE=inv



earth_rad = 6372.1

mean_rad = 6472.1 ;km  ( = <Re> + 100 km)



; Get location of dipole



dpo = dipole_offset(year)	; 3-element vector in km



; Make the rotation matrices (for transformation either way).



yyvar = year - 1995		; IGRF 1995 revision

g10 = -29682.0 + 17.6*yyvar

g11 = -1789.0 + 13.0*yyvar

h11 = 5318.0 - 18.3*yyvar

h2 = g10*g10 + g11*g11 + h11*h11

alpha = acos( -g10 / sqrt(h2) )

beta = atan(h11/g11)

ca = cos(alpha)

sa = sin(alpha)

cb = cos(beta)

sb = sin(beta)



geo2mag = [	[ca*cb, ca*sb, -sa], $

		[-sb  , cb   , 0  ], $

		[sa*cb, sa*sb, ca ]	]



if NOT keyword_set(inv) then mag2geo = transpose(geo2mag)



if keyword_set(deg) then begin	; I/O in degrees, must convert to radians

	tlat = lat*!dtor	; The transformed array

	tlng = lng*!dtor

endif else begin

	tlat = lat		; Copy input to output arrays

	tlng = lng

endelse



tlat = !pi/2. - tlat		; colattitude in rads for spherical coords



; Create array of column vectors (lats and lons -> cartesian coords.)



pos = [ [sin(tlat)*cos(tlng)], $

	[sin(tlat)*sin(tlng)], $

        [cos(tlat)          ]	] ; unit vectors



; Scale the position vectors to touch 100km alt if doing GEO->MAG (inv)

; or Re if doing MAG->GEO.



if keyword_set(inv) then pos=temporary(pos)*mean_rad $

  else pos=temporary(pos)*earth_rad

	

; PERFORM THE TRANSFORMATION



if NOT keyword_set(inv) then begin ;mag -> geo

        

        tpos = mag2geo ## pos       ; Lengths unchanged, axes now || to GEO



        ;; Now scale the vectors to 100 km alt

        ;; This does not entail constraint to field lines



	square = tpos(*,0)^2 + tpos(*,1)^2 + tpos(*,2)^2

	cross = 2*dpo(0)*tpos(*,0) + 2*dpo(1)*tpos(*,1) + 2*dpo(2)*tpos(*,2)

	const = dpo(0)^2 + dpo(1)^2 + dpo(2)^2 - mean_rad^2

	scale = ( -cross + sqrt(cross^2 - 4*square*const)) / (2*square)

        ;; if where(scale LT 0) NE -1 then message, 'Scaling error in T-form.'

	for i=0,2 do tpos(*,i) = tpos(*,i)*scale

        

        ; Find the difference between original vectors and scaled ones



        new_len = sqrt(tpos(*,0)^2 + tpos(*,1)^2 + tpos(*,2)^2)

        diff = new_len - mean_rad



        ; Find the fix in lattitude necessary to constrain scaled

        ; vectors to field lines.  Then fix and recalculate vectors.



        latfix = atan(diff*tan(tlat)/new_len/2)

        tlat = tlat + latfix

        pos = [ [sin(tlat)*cos(tlng)], $

                [sin(tlat)*sin(tlng)], $

                [cos(tlat)          ]	] ; reset pos to MAG unit vectors

        for i=0,2 do pos(*,i) = pos(*,i)*new_len ; scale to 100km individually

        tpos = mag2geo ## pos

        

        ; Translate the coordinate system to earth center



	for i=0,2 do tpos(*,i) = tpos(*,i) + dpo(i)



endif else begin ;do geo -> mag

        tpos = pos ; Set up the target array

	for i=0,2 do tpos(*,i) = pos(*,i) - dpo(i)

	tpos = geo2mag ## tpos

endelse



; Convert back into lattitude and longitude, degrees if necessary



tlat = atan(tpos(*,2)/sqrt(tpos(*,0)^2 + tpos(*,1)^2))	; assume < pi/2

tlng = atan(tpos(*,1)/tpos(*,0))			; can't assume

flip = where(tpos(*,0) LT 0)

if flip(0) NE -1 then tlng(flip) = tlng(flip) + !pi



if keyword_set(deg) then begin

	rerange, /deg, tlng, tlat

endif else begin

	tlat = tlat MOD (2*!pi)

	tlng = tlng MOD (2*!pi)

endelse



return

end