horizon_radio.pro (Kronos IDL library)

;+
; Contains the horizon_radio function
;
; :Author:
;    Baptiste Cecconi
;
; :History:
;    2009/05/19: Created
;
;    2009/05/19: Last Edit
;-
;
;+
; horizon_radio is a function that <behavior desc here>
;
; :Uses:
;    magnetic_field, mfl_foot, xyz_to_rtp
;
; :Returns:
;    <return desc here>
;
; :Params:
;    xyz_cassini:       in, required, type=sometype
;       A parameter named xyz_cassini
;    freq:              in, required, type=sometype
;       A parameter named freq
;
; :Keywords:
;    ntheta:            in, optional, type=sometype
;       A keyword named ntheta
;    nphi:              in, optional, type=sometype
;       A keyword named nphi
;    mfl_model:         in, optional, type=sometype
;       A keyword named mfl_model
;    horizon_foot:      in, optional, type=sometype
;       A keyword named horizon_foot
;    nocrt:             in, optional, type=sometype
;       A keyword named nocrt
;    observer_centered: in, optional, type=sometype
;       A keyword named observer_centered
;    normalized_vector: in, optional, type=sometype
;       A keyword named normalized_vector
;-
FUNCTION horizon_radio,xyz_cassini,freq,ntheta=ntheta,nphi=nphi, $
  mfl_model=mfl_model,horizon_foot=horizon_foot, nocrt=nocrt, $
  observer_centered=observer_centered, normalized_vector=normalized_vector

if keyword_set(ntheta) then ntheta = long(ntheta) else ntheta = 1800l
if keyword_set(nphi)   then nphi   = long(nphi)   else nphi   = 3600l
if ~ keyword_set(mfl_model) then mfl_model = 'SPV'


; v_SP : vecteur Saturn -> Point etudie
; v_SC : vecteur Saturn -> Cassini
; v_PC : vecteur Point etudie -> Cassini

; longueur max du segment SP pour l'amorce de la dichotomie
l_SP_max = 20.

; grille meridienne (ntheta valeurs entre 1 et 179 degre) (meridien: y=0)
; => on ne calcule pas le Bmag au pole... 
v_SP0 = fltarr(3,ntheta)
v_SP1 = transpose([[sin((findgen(ntheta)/(ntheta-1)*178.+1.)/!radeg)], $
                   [fltarr(ntheta)                     ], $
                   [cos((findgen(ntheta)/(ntheta-1)*178.+1.)/!radeg)] ]) * l_SP_max

; nombre de pas de l'algoritme de dichotomie
ndic = 30L

; rapport f/b en kHz
fsb=2.79924835996d0 ; converts Gauss to MHz
fsb=fsb*1.d3        ; converts Gauss to kHz

for idic=0L,ndic-1L do begin
  v_SP2 = (v_SP0+v_SP1)/2
  rtp2  = xyz_to_rtp(v_SP2)
  magnetic_field,rtp2(0,*),rtp2(1,*),rtp2(2,*),bvec2,btot2,mfl_model=mfl_model,nocrt=nocrt
  fc2 = btot2*fsb
  ww0 = where(fc2 ge freq,cnt0,compl=ww1,ncompl=cnt1)
  if cnt0 ne 0 then begin
    v_SP0(*,ww0) = v_SP2(*,ww0)
  endif
  if cnt1 ne 0 then begin 
    v_SP1(*,ww1) = v_SP2(*,ww1)
  endif
endfor

v_SP = (v_SP0+v_SP1)/2.

;plot,v_SP(0,*),v_sp(2,*),psym=-4,/iso,xs=3,ys=3

;calcul de la grille 3D en coordonnees cartesiennes

v_SP_grid = fltarr(3,ntheta,nphi)
for iphi=0,nphi-1 do begin 
  phi = float(iphi)/nphi*2.*!pi
  v_SP_grid(0,*,iphi) = v_SP(0,*)*cos(phi)
  v_SP_grid(1,*,iphi) = v_SP(0,*)*sin(phi)
  v_SP_grid(2,*,iphi) = v_SP(2,*)
endfor

ngrid = ntheta*nphi
v_SP_grid = reform(v_SP_grid,3,ngrid)

v_SC = xyz_cassini

v_PC = v_SP_grid - rebin(reform(v_SC,3,1),3,ngrid)

th_SPC = angular_distance(-v_SP_grid,v_PC)

ww = where(abs(th_spc*!radeg-90) le .1) 

v_SP_horizon = v_SP_grid(*,ww)

n_SP_horizon = n_elements(ww)


ff = fltarr(ngrid)+freq
rtp_horizon = xyz_to_rtp(v_SP_horizon)
horizon_foot = mfl_foot(reform(rtp_horizon(0,*)),             $
                       reform(rtp_horizon(2,*))*!radeg,       $
                       reform(rtp_horizon(1,*))*!radeg,       $
                       ff,mfl_model=mfl_model,nocrt=nocrt)

if keyword_set(observer_centered) then begin 
  v_SP_horizon = v_SP_horizon - rebin(reform(v_SC,3,1),3,n_SP_horizon)
  if keyword_set(normalized_vector) then begin
    v_SP_horizon = v_SP_horizon / rebin(reform(sqrt(total(v_SP_horizon^2,1)),1,n_SP_horizon),3,n_SP_horizon)
  endif
endif

return,v_SP_horizon
end