;+ ; 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 ; ;- ;==================================================================================== FUNCTION 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 ;===============================================================================