/usr/share/gnudatalanguage/astrolib/geo2eci.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 | ;+
; NAME:
; GEO2ECI
;
; PURPOSE:
; Convert geographic spherical coordinates to Earth-centered inertial coords
;
; EXPLANATION:
; Converts from geographic spherical coordinates [latitude, longitude,
; altitude] to ECI (Earth-Centered Inertial) [X,Y,Z] rectangular
; coordinates. JD time is also needed.
;
; Geographic coordinates are in degrees/degrees/km
; Geographic coordinates assume the Earth is a perfect sphere, with radius
; equal to its equatorial radius.
; ECI coordinates are in km from Earth center at epoch TOD (True of Date)
;
; CALLING SEQUENCE:
; ECIcoord=geo2eci(gcoord,JDtime)
;
; INPUT:
; gcoord: geographic [latitude,longitude,altitude], or a an array [3,n]
; of n such coordinates
; JDtime: Julian Day time, double precision. Can be a 1-D array of n
; such times.
;
; KEYWORD INPUTS:
; None
;
; OUTPUT:
; a 3-element array of ECI [X,Y,Z] coordinates, or an array [3,n] of
; n such coordinates, double precision. The TOD epoch is the
; supplied JDtime.
;
; COMMON BLOCKS:
; None
;
; PROCEDURES USED:
; CT2LST - Convert Local Civil Time to Local Mean Sidereal Time
;
; EXAMPLES:
;
; IDL> ECIcoord=geo2eci([0,0,0], 2452343.38982663D)
; IDL> print,ECIcoord
; -3902.9606 5044.5548 0.0000000
;
; (The above is the ECI coordinates of the intersection of the equator and
; Greenwich's meridian on 2002/03/09 21:21:21.021)
;
; MODIFICATION HISTORY:
; Written by Pascal Saint-Hilaire (shilaire@astro.phys.ethz.ch)
; on 2002/05/14
; Update documentation to specify epoch is TOD.
; R. Redmon NOAA/NGDC April 2014
;
;-
;====================================================================================
FUNCTION geo2eci,incoord,JDtim
Re=6378.137 ; Earth's equatorial radius, in km
lat = DOUBLE(incoord[0,*])*!DPI/180.
lon = DOUBLE(incoord[1,*])*!DPI/180.
alt = DOUBLE(incoord[2,*])
JDtime= DOUBLE(JDtim)
ct2lst,gst,0,0,JDtime
angle_sid=gst*2.*!DPI/24. ; sidereal angle
theta=lon+angle_sid ; azimuth
r=(alt+Re)*cos(lat)
X=r*cos(theta)
Y=r*sin(theta)
Z=(alt+Re)*sin(lat)
RETURN,[X,Y,Z]
END
;====================================================================================
|