; 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." ; ; ; ; ; ; ; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;+ ; 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 ; ; common aacgm_com,coef,cint,height_old,first_coeff_old common mlt_com,sol_dec_old,told,mslon1,mslon2 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 @default.pro end