mfl_foot.pro (Kronos IDL library)

;+
; Contains the mfl_foot function 
; and associated functions and procedures
;
; :Author:
;    Baptiste Cecconi
;
; :History:
;    2009/06/14: Created
;
;    V1_BC : added H0 and MFL_MODEL keywords.
;
;    V2_BC : added vector analysis  and verbose keyword
;
;    V21_BC : added NOCRT keyword (update for new magnetic_field routines)
;
;    V30_BC : using pre-computed iso-fc grid.
;
;    V31_BC : fixed memory leak
;
;    V32_BC : fixed path for iso_fc saveset file
;
;    2010/03/16: Last Edit
;-
;
;+
; stop function
;
; :Params:
;    y: in, required, type=fltarr
;       == rtp (t is latitude and p is azimuth, both in degree)
;
; :Returns:
;    result is altitude over saturn's surface.
;-
function G_STOP, y
  eps=0.098
  rep=y(0,*)-1.d0/sqrt(1.d0+eps*(2.d0-eps)/(1.d0-eps)^2*(cos(y(1,*)*!dpi/180.d0))^2)
return,rep
end


;+
; derivative function
;
; computes the tangent (sph coord)) to the local magnetic field line 
; at the RTP point
;-
pro DERIV, rtp, drtp_ds, b, bv,mfl_model=mfl_model, nocrt=nocrt
  nn = n_elements(rtp)/3
  MAGNETIC_FIELD,rtp(0,*),rtp(1,*),rtp(2,*),bv,b,mfl_model=mfl_model,/deg,nocrt=nocrt
  drtp_ds=[bv(0,*)/b, $
           bv(1,*)/rtp(0,*)*!radeg/b, $
           bv(2,*)/rtp(0,*)*!radeg/sin(rtp(1,*)/!radeg)/b]
return
end

;+
; mfl_foot is a function that <behavior desc here>
;
; :Uses:
;    magnetic_field
;
; :Returns:
;    position de la ligne de champ en coord sphériques 
;
; :Params:
;    ray:       in, required, type=sometype
;       r initial en Rs (always>0)
;    azi:       in, required, type=sometype
;       azimuth initial en degre
;    col:       in, required, type=sometype
;       colatitude initiale en degre
;    ff:        in, required, type=sometype
;       A parameter named ff
;
; :Keywords:
;    mfl_model: in, optional, type=string
;       mag field model
;    lun:       in, optional, type=sometype
;       A keyword named lun
;    verbose:   in, optional, type=sometype
;       additionnal info output
;    test:      in, optional, type=sometype
;       A keyword named test
;    nocrt:     in, optional, type=sometype
;       A keyword named nocrt
;-
FUNCTION MFL_FOOT, ray, azi, col, ff, mfl_model=mfl_model, lun=lun, $
                   verbose=verbose, test=test, nocrt=nocrt
; rtp (OUTPUT)=position de la ligne de champ en coord sphériques 
; h0 : computation step (sign depends on hemisphere)
nr = n_elements(ray)
nl = n_elements(azi)
nc = n_elements(col)

nk = [nr,nl,nc]
nk = nk(sort(nk))
if nk(0)*nk(1)*nk(2) eq 1 then begin
  nn = 1 
endif else begin
  if nk(0)*nk(1) eq 1 then begin
    nn = nk(2) 
  endif else begin
    if nk(0) eq 1 then begin
      if nk(1) eq nk(2) then begin
        nn = nk(1)
      endif else message,'dimensions of input arrays do not agree, aborting.'
    endif else begin
      if nk(0) eq nk(1) and nk(1) eq nk(2) then begin
        nn = nk(0)
      endif else message,'dimensions of input arrays do not agree, aborting.'
    endelse
  endelse
endelse


root_rpws = getenv('ROOT_RPWS')
path = root_rpws+'/bin/'
iso_fc_file = 'iso_fc_grid_'+mfl_model
;if keyword_set(nocrt) then iso_fc_file += '_nocrt'
iso_fc_file += '.sav'
restore,path+iso_fc_file,verbose=verbose

rtp = fltarr(3,nn)

uff = uniq_list(ff, rev=rff)
nff = n_elements(uff)
for i=0,nff-1l do begin
  wff = rff(rff(i):rff(i+1)-1)
  w_iso_fc = where(ff_iso_fc eq uff(i),cnt)
  if cnt eq 0 then begin 
    message,/info,"Some input frequencies cannot be processed."
    message,/info,"Frequency should be either one of the HFR frequencies or a multiple of 5kHz."
    message,/info,'INFO: Conflicting frequency value = '+string(uff(i))
    message,"Aborting!"
  endif
  if ptr_valid(iso_fc(w_iso_fc(0))) then begin
    rtp_iso_fc = *iso_fc(w_iso_fc(0))
    rtp_iso_fc_foot = *iso_fc_foot(w_iso_fc(0))
    rtp(0,wff) = interpol(rtp_iso_fc_foot(0,*),rtp_iso_fc(1,*),col(wff))
    rtp(1,wff) = interpol(rtp_iso_fc_foot(1,*),rtp_iso_fc(1,*),col(wff))
    rtp(2,wff) = azi(wff)
  endif
endfor

ptr_free,iso_fc
ptr_free,iso_fc_foot
ptr_free,iso_fc_bvect
ptr_free,iso_fc_bmag
ptr_free,iso_fc_gradb

return,rtp
end