aacgmidl.pro (Documentation for /home/davin/idl/socware/)

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

@default.pro

end