;+
; NAME:
;       GEO2MAG()
;
; PURPOSE:
;       Convert from geographic to geomagnetic coordinates
; EXPLANATION:
;       Converts from GEOGRAPHIC (latitude,longitude) to GEOMAGNETIC (latitude, 
;       longitude).   (Altitude remains the same)
;
;       Latitudes and longitudes are expressed in degrees.
;
; CALLING SEQUENCE:
;       mcoord=geo2mag(gcoord)
;
; INPUT:
;       gcoord = a 2-element array of geographic [latitude,longitude], or an 
;                array [2,n] of n such coordinates.
;
; KEYWORD INPUTS:
;       None
;
; OUTPUT:
;       a 2-element array of magnetic [latitude,longitude], or an array [2,n] 
;         of n such coordinates                     
;
; COMMON BLOCKS:
;       None
;
; EXAMPLES:
;       geographic coordinates of magnetic south pole
;
;       IDL> mcoord=geo2mag([79.3,288.59])      
;       IDL> print,mcoord
;       89.999992      -173.02325
;
; MODIFICATION HISTORY:
;       Written by Pascal Saint-Hilaire (Saint-Hilaire@astro.phys.ethz.ch), 
;            May 2002
;       Changed the name to avoid conflicts, jmm, 2017-12-13, no other
;       changes from the SSW version               
;-

;====================================================================================
FUNCTION barrel_geo2mag,incoord

        ; SOME 'constants'...
        Dlong=288.59D   ; longitude (in degrees) of Earth's magnetic south pole
                        ;(which is near the geographic north pole!) (1995)
        Dlat=79.30D     ; latitude (in degrees) of same (1995)
        R = 1D          ; distance from planet center (value unimportant -- 
                 ;just need a length for conversion to rectangular coordinates)

        ; convert first to radians
        Dlong=Dlong*!DPI/180.
        Dlat=Dlat*!DPI/180.

        glat=DOUBLE(incoord[0,*])*!DPI/180.
        glon=DOUBLE(incoord[1,*])*!DPI/180.
        galt=glat * 0. + R
        
        coord=[glat,glon,galt]

        ;convert to rectangular coordinates
        ;       X-axis: defined by the vector going from Earth's center towards
        ;            the intersection of the equator and Greenwitch's meridian.
        ;       Z-axis: axis of the geographic poles
        ;       Y-axis: defined by Y=Z^X
        x=coord[2,*]*cos(coord[0,*])*cos(coord[1,*])
        y=coord[2,*]*cos(coord[0,*])*sin(coord[1,*])
        z=coord[2,*]*sin(coord[0,*])

        ;Compute 1st rotation matrix : rotation around plane of the equator,
        ;from the Greenwich meridian to the meridian containing the magnetic
        ;dipole pole.
        geolong2maglong=dblarr(3,3)
        geolong2maglong[0,0]=cos(Dlong)
        geolong2maglong[0,1]=sin(Dlong)
        geolong2maglong[1,0]=-sin(Dlong)
        geolong2maglong[1,1]=cos(Dlong)
        geolong2maglong[2,2]=1.
        out=geolong2maglong # [x,y,z]

        ;Second rotation : in the plane of the current meridian from geographic
        ;                  pole to magnetic dipole pole.
        tomaglat=dblarr(3,3)
        tomaglat[0,0]=cos(!DPI/2-Dlat)
        tomaglat[0,2]=-sin(!DPI/2-Dlat)
        tomaglat[2,0]=sin(!DPI/2-Dlat)
        tomaglat[2,2]=cos(!DPI/2-Dlat)
        tomaglat[1,1]=1.
        out= tomaglat # out

        ;convert back to latitude, longitude and altitude
        mlat=atan(out[2,*],sqrt(out[0,*]^2+out[1,*]^2))
        mlat=mlat*180./!DPI
        mlon=atan(out[1,*],out[0,*])
        mlon=mlon*180./!DPI
        ;malt=sqrt(out[0,*]^2+out[1,*]^2+out[2,*]^2)-R 
;  I don't care about that one...just put it there for completeness' sake

        RETURN,[mlat,mlon]
END
;===============================================================================