; aacgmidl.pro
;
; By R.J.Barnes
; Original version of AACGM by Kile Baker and Simon Wing
;
; Copyright � 2001 The Johns Hopkins University/Applied Physics Laboratory.
; All rights reserved.
;
; This material may be used, modified, or reproduced by or for the U.S.
; Government pursuant to the license rights granted under the clauses at DFARS
; 252.227-7013/7014.
;
; For any other permissions, please contact the Space Department
; Program Office at JHU/APL.
;
; This Distribution and Disclaimer Statement must be included in all copies of
; RST-ROS (hereinafter "the Program").
;
; The Program was developed at The Johns Hopkins University/Applied Physics
; Laboratory (JHU/APL) which is the author thereof under the "work made for
; hire" provisions of the copyright law.
;
; JHU/APL assumes no obligation to provide support of any kind with regard to
; the Program. This includes no obligation to provide assistance in using the
; Program or to provide updated versions of the Program.
;
; THE PROGRAM AND ITS DOCUMENTATION ARE PROVIDED AS IS AND WITHOUT ANY EXPRESS
; OR IMPLIED WARRANTIES WHATSOEVER. ALL WARRANTIES INCLUDING, BUT NOT LIMITED
; TO, PERFORMANCE, MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE ARE
; HEREBY DISCLAIMED. YOU ASSUME THE ENTIRE RISK AND LIABILITY OF USING THE
; PROGRAM TO INCLUDE USE IN COMPLIANCE WITH ANY THIRD PARTY RIGHTS. YOU ARE
; ADVISED TO TEST THE PROGRAM THOROUGHLY BEFORE RELYING ON IT. IN NO EVENT
; SHALL JHU/APL BE LIABLE FOR ANY DAMAGES WHATSOEVER, INCLUDING, WITHOUT
; LIMITATION, ANY LOST PROFITS, LOST SAVINGS OR OTHER INCIDENTAL OR
; CONSEQUENTIAL DAMAGES, ARISING OUT OF THE USE OR INABILITY TO USE THE
; PROGRAM."
;
; NOTE: Minor modifications made by Patrick Cruce at IGPP/UCLA for use with
; THEMIS TDAS Software in July/August 2008
;
; 5/27/2015:
; modified by egrimes@igpp, added a note to aacgmidl to load the coefficients
; before trying to use these routines. also note that this previously
; loaded the coefficients for year = 2000 by batch executing default.pro
;
;ALSO NOTE:
; Haje Korth's AACGM DLM interface to the C code is significantly faster:
; http://ampere.jhuapl.edu/code/idl_aacgm.html
;
;
;
;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;+
; NAME:
; AACGMIDL
;
; PURPOSE:
; This library is a pure IDL version of the AACGM library
;
;
; ROUTINES:
;
; EQN_OF_TIME equation of time
; SOLAR_LOC find location of sun
; CNV_AACGM convert coordinates to/from magnetic/geographic
; CALC_MLT calculate magnetic local time
; LOAD_COEF load a set of AACGM coefficients
;
;----------------------------------------------------------------------------
;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;+
; NAME:
; EQN_OF_TIME
;
; PURPOSE:
; equation of time for a given longitude and year
;
; Calling sequence:
; eqt = eqn_of_time(mean_lon,yr)
;
;
;
;---------------------------------------------------------------------
;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;+
; NAME:
; SOLAR_LOC
;
; PURPOSE:
; location of the sun for given year and time
;
; Calling sequence:
; solar_loc,yr,t1,mean_lon,dec
; t1 is the seconds from the start of year.
; the mean longitude and declination are returned
; in mean_lon and dec as floats.
;
;
;---------------------------------------------------------------------
;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;+
; NAME:
; CNV_AACGM
;
; PURPOSE:
; convert to and from AACGM and Geographic coordinates
;
; Calling sequence:
; CNV_AACGM,in_lat,in_lon,height,out_lat,out_lon,r,error
; the calculated latitude and longitude for the
; given height are returned in out_lat,out_lon.
;
;
;
;---------------------------------------------------------------------
;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;+
; NAME:
; CALC_MLT
;
; PURPOSE:
; calculate magnetic local time for a given longitude
;
; Calling sequence:
; mlt=CALC_MLT(yr,t0,mlong)
; t1 is the seconds from the start of year and
; mlong is the magnetic longitude of the observing
; point.
;
;---------------------------------------------------------------------
;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;+
; NAME:
; LOAD_COEF
;
; PURPOSE:
; load a set of AACGM coefficients
;
; Calling sequence:
; load_coef,fname
; fname is the filename.
;
;
;---------------------------------------------------------------------
;
;
; $Log: aacgmidl.pro,v $
; Revision 1.2 2002/03/29 16:03:43 barnes
; Fixed credits.
;
; Revision 1.1 2002/03/29 15:50:19 barnes
; Initial revision
;
;
; sind and cosd removed to avoid duplicate routines (other versions in ssl_general/misc/math)
;function sind, theta
; return, sin(theta*!PI/180.0)
;end
;
;function cosd, theta
; return, cos(theta*!PI/180.0)
;end
function asind, theta
return, asin(theta)*180.0/!PI
end
function acosd, theta
return, acos(theta)*180.0/!PI
end
function sgn,a,b
if (a ge 0) then x=a else x=-a
if (b ge 0) then return, x
return, -x
end
function modulus,x,y
quotient=x/y
if (quotient ge 0) then quotient=floor(quotient) $
else quotient=-floor(-quotient)
return, (x-y*quotient)
end
pro rylm,colat,lon,order,ylmval
cos_theta = cos(colat)
sin_theta = sin(colat)
cos_lon = cos(lon)
sin_lon = sin(lon)
d1 = -sin_theta;
z2=complex(cos_lon,sin_lon)
z1=d1*z2
q_fac=z1
ylmval[0]=1;
ylmval[2]=cos_theta;
for l=1,order-1 do begin
la = (l - 1) * l + 1
lb = l * (l + 1) + 1
lc = (l + 1) * (l + 2) + 1;
ca =float(l * 2. + 1.) / (l + 1.)
cb= float(l)/(l + 1.)
ylmval[lc-1] = ca * cos_theta * ylmval[lb-1] - cb * ylmval[la-1]
end
q_val=q_fac;
ylmval[3]=float(q_val)
ylmval[1]=-imaginary(q_val)
for l=2,order do begin
d1 = l*2 - 1.
z2=d1*q_fac
z1=z2*q_val
q_val=z1
la = l*l + (2*l) + 1
lb = l*l + 1;
ylmval[la-1] = float(q_val)
ylmval[lb-1] = -imaginary(q_val)
end
for l=2,order do begin
la = l*l
lb = l*l - 2*(l - 1)
lc = l*l + (2*l);
ld = l*l + 2
fac = l*2 - 1
ylmval[lc-1] = fac * cos_theta * ylmval[la-1]
ylmval[ld-1] = fac * cos_theta * ylmval[lb-1]
end
for m=1,order-2 do begin
la = (m+1)*(m+1)
lb = (m+2)*(m+2)-1
lc = (m+3)*(m+3)-2
ld = la - (2*m)
ldd = lb - (2*m)
lf = lc - (2*m)
for l=m+2,order do begin
ca=float(2.*l-1)/(l-m)
cb=float(l+m-1.)/(l-m)
ylmval[lc-1] = ca * cos_theta *ylmval[lb-1] - cb *ylmval[la-1]
ylmval[lf-1] = ca * cos_theta *ylmval[ldd-1] - cb *ylmval[ld-1]
la = lb
lb = lc
lc = lb + (2*l) + 2
ld = la - (2*m)
ldd = lb - (2*m)
lf = lc - (2*m)
end
end
return
end
pro altitude_to_cgm, r_height_in, r_lat_alt, r_lat_adj
eradius=6371.2
eps=1e-9
unim=0.9999999;
r1 = cosd(r_lat_alt)
ra = r1 * r1
if (ra lt eps) then ra = eps
r0 = (r_height_in/eradius+1) / ra
if (r0 lt unim) then r0 = unim
r1 = acos(sqrt(1/r0));
r_lat_adj= sgn(r1, r_lat_alt)*180./!pi;
return
end
pro cgm_to_altitude, r_height_in,r_lat_in, r_lat_adj, error
eradius=6371.2
unim=1
error=0
r1 = cosd(r_lat_in);
ra = (r_height_in/ eradius+1)*(r1*r1);
if (ra gt unim) then begin
ra = unim;
error=1;
end
r1 = acos(sqrt(ra));
r_lat_adj = sgn(r1,r_lat_in)*180./!pi;
return
end
function eqn_of_time,mean_lon,yr
eqcoef=[ [-105.8,596.2,4.4,-12.7,-429.0,-2.1,19.3], $
[-105.9,596.2,4.4,-12.7,-429.0,-2.1,19.3], $
[-106.1,596.2,4.4,-12.7,-428.9,-2.1,19.3], $
[-106.2,596.2,4.4,-12.7,-428.9,-2.1,19.3], $
[-106.4,596.1,4.4,-12.7,-428.9,-2.1,19.3], $
[-106.5,596.1,4.4,-12.7,-428.8,-2.1,19.3], $
[-106.6,596.1,4.4,-12.7,-428.8,-2.1,19.3], $
[-106.7,596.1,4.4,-12.7,-428.7,-2.1,19.3], $
[-106.8,596.1,4.4,-12.7,-428.7,-2.1,19.3], $
[-107.0,596.1,4.4,-12.7,-428.7,-2.1,19.3], $
[-107.2,596.1,4.4,-12.7,-428.6,-2.1,19.3], $
[-107.3,596.1,4.4,-12.7,-428.6,-2.1,19.3] $
]
index=1
if (yr lt 88) then index=yr+2000-1988
if ((yr ge 88) and (yr lt 100)) then index=yr-88 $
else if ((yr ge 100) and (yr lt 1900)) then index=yr-88 $
else index=yr-1988
if (index lt 1) then index=1
if (index gt 12) then index=12;
return, eqcoef[0,index-1]*sind(mean_lon)+ $
eqcoef[1,index-1]*sind(2.0*mean_lon)+ $
eqcoef[2,index-1]*sind(3.0*mean_lon)+ $
eqcoef[3,index-1]*sind(4.0*mean_lon)+ $
eqcoef[4,index-1]*cosd(mean_lon)+ $
eqcoef[5,index-1]*cosd(2.0*mean_lon)+ $
eqcoef[6,index-1]*cosd(3.0*mean_lon)
end
pro solar_loc,yr, t1, mean_lon, dec
L0=[279.642,279.403,279.165,278.926,279.673,279.434, $
279.196,278.957,279.704,279.465,279.226,278.982]
DL=0.985647
G0=[356.892984,356.637087,356.381191,356.125295, $
356.854999,356.599102,356.343206,356.087308, $
356.817011,356.561113,356.31,356.05]
DG=0.98560028
EPS0=[23.440722,23.440592,23.440462,23.440332, $
23.440202,23.440072,23.439942,23.439811, $
23.439682,23.439552,23.439422,23.439292]
DE=-0.00000036;
d = 0;
if (yr lt 1900) then index = yr - 88 $
else index = yr - 1988
if (index le 0) then delta_yr = index - 1 $
else if (index gt 10) then delta_yr = index - 10 $
else delta_yr = 0;
if (index lt 1) then index = 1
if (index gt 12) then index = 12
yr_step = sgn(1,delta_yr)
delta_yr = abs(delta_yr)
for i=1,delta_yr do begin
if (yr_step gt 0) then yrs=98+i-1 $
else yrs=89-i
if (modulus(yrs,4) eq 0) then d = d + 366*yr_step $
else d = d + 365*yr_step
end
d = d + t1/86400
L = L0[index-1] + DL*d
g = G0[index-1] + DG*d
while (L lt 0) do L = L + 360
while (g lt 0) do g = g + 360
L = modulus(L,360.0)
g = modulus(g,360.0)
lambda = L + 1.915*sind(g) + 0.020*sind(2*g)
eps = EPS0[index-1] + DE*d
dec = asind(sind(eps)*sind(lambda))
mean_lon = L
return
end
pro convert_geo_coord, lat_in,lon_in,height_in,lat_out,lon_out, $
order,error, geo=geo
common aacgm_com,coef,cint,height_old,first_coeff_old
flag=keyword_set(GEO);
if lon_in lt 0 then lon_in=lon_in+360
if (first_coeff_old ne coef[0,0,0,0]) then height_old=[-1.0, -1.0]
first_coeff_old=coef[0,0,0,0]
error=-2
if ((height_in lt 0) or (height_in gt 7200)) then return
error=-8;
if (abs(lat_in) gt 90.) then return
error=-16
if ((lon_in lt 0) or (lon_in gt 360)) then return
if (height_in ne height_old(flag)) then begin
alt_var= height_in/7200.0;
alt_var_sq = alt_var * alt_var;
alt_var_cu = alt_var * alt_var_sq;
alt_var_qu = alt_var * alt_var_cu;
for i=0,2 do begin
for j=0,120 do begin
cint[j,i,flag] =coef[j,i,0,flag]+ $
coef[j,i,1,flag]*alt_var+ $
coef[j,i,2,flag]*alt_var_sq+ $
coef[j,i,3,flag]*alt_var_cu+ $
coef[j,i,4,flag]*alt_var_qu
end
end
height_old[flag] = height_in;
end
x=0.
y=0.
z=0.
lon_input =lon_in*!pi/180.0;
if (flag eq 0) then colat_input = (90.-lat_in)*!pi/180.0 $
else begin
error=-64
cgm_to_altitude,height_in, lat_in,lat_adj,errflg
if (errflg ne 0) then return;
colat_input= (90. - lat_adj)*!pi/180.0;
end
ylmval=fltarr(121);
rylm,colat_input,lon_input,order,ylmval;
for l = 0, order do begin
for m = -l,l do begin
k = l * (l+1) + m+1;
x=x+cint[k-1,0,flag]*ylmval[k-1]
y=y+cint[k-1,1,flag]*ylmval[k-1]
z=z+cint[k-1,2,flag]*ylmval[k-1]
end
end
error=-32
r = sqrt(x * x + y * y + z * z)
if ((r lt 0.9) or (r gt 1.1)) then return
z=z / r
x=x / r
y=y / r
if (z ge 1.) then colat_temp=0 $
else if (z lt -1.) then colat_temp =!pi $
else colat_temp= acos(z)
if ((abs(x) lt 1e-8) and (abs(y) lt 1e-8)) then lon_temp =0 $
else lon_temp = atan(y,x)
lon_output = lon_temp
if (flag eq 0) then begin
lat_alt =90 - colat_temp*180.0/!pi;
altitude_to_cgm,height_in, lat_alt,lat_adj
colat_output = (90. - lat_adj) * !pi/180.0;
end else colat_output = colat_temp
lat_out =90. - colat_output*180.0/!pi
lon_out =lon_output*180.0/!pi
error=0
return
end
pro mlt1,t0,solar_dec,mlon,mlt,mslon
common mlt_com,sol_dec_old,told,mslon1,mslon2
if ((abs(solar_dec-sol_dec_old) gt 0.1) or (sol_dec_old eq 0)) then told=1e12
if (abs(mslon2-mslon1) gt 10) then told=1e12;
if ((t0 ge told) and (t0 lt (told+600))) then $
mslon=mslon1+(t0-told)*(mslon2-mslon1)/600.0 $
else begin
told=t0
sol_dec_old=solar_dec
slon1 = (43200.0-t0)*15.0/3600.0
slon2 = (43200.0-t0-600)*15.0/3600.0
height = 450
convert_geo_coord,solar_dec,slon1,height,mslat1,mslon1,4,err
convert_geo_coord,solar_dec,slon2,height,mslat2,mslon2,4,err
mslon=mslon1
end
mlt = (mlon - mslon) /15.0 + 12.0
if (mlt ge 24) then mlt=mlt-24;
if (mlt lt 0) then mlt=mlt+24;
end
function calc_mlt,yr,t0,mlong
if (yr gt 1900) then yr=yr-1900
mean_lon=0.0
dec=0.0
mlt=0.0
solar_loc,yr, t0,mean_lon,dec
et = eqn_of_time(mean_lon, yr)
dy=floor(t0/(24.*3600.))
ut=float(t0-(dy*24*3600));
apparent_time = ut + et;
mlt1,apparent_time, dec, mlong, mlt,mslong
return, mlt
end
pro cnv_aacgm, in_lat,in_lon,height,out_lat,out_lon,r,error, geo=geo
out_lat=0.
out_lon=0.
geo=keyword_set(geo)
convert_geo_coord,in_lat,in_lon,height,out_lat,out_lon,10,error,geo=geo
r=1.0;
end
pro load_coef,filename
common aacgm_com,coef,cint,height_old,first_coeff_old
openr, unit, filename,/GET_LUN
coef=fltarr(121,3,5,2)
readf, unit,coef
close,unit
free_lun,unit
end
pro aacgmidl
common aacgm_com,coef,cint,height_old,first_coeff_old
common mlt_com,sol_dec_old,told,mslon1,mslon2
coef=fltarr(121,3,5,2)
cint=fltarr(121,3,2)
height_old=[-1.,-1.]
first_coeff_old=-1.
sol_dec_old=0
told=1e12
mslon1=0
mslon2=0
; NOTE: must load the coefficients! they can be found in:
; general/cotrans/aacgm/coef
rt_info = routine_info('aacgmidl',/source)
path = file_dirname(rt_info.path) + '/coef/'
load_coef, path+'aacgm_coeffs2010.asc'
;@default.pro
end