/usr/share/gnudatalanguage/astrolib/geo2mag.pro is in gdl-astrolib 2018.02.16+dfsg-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 | ;+
; 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
;===============================================================================
|