dfb_main_full.pro (Kronos IDL library)

;+
; Contains the dfb_main_full procedure
;
; :Author:
;    Baptiste Cecconi
;
; :History:
;    2009/06/03: Created
;
;    2009/11/24: Last Edit
;-
;
;+
; dfb_main_full is a procedure that <behavior desc here>
;
; :Uses:
;    crossp1, read_antenna_set, ti_t97
;
; :Params:
;    df_result:    in, required, type=sometype
;       A parameter named df_result
;    data_p:       in, required, type=sometype
;       A parameter named data_p
;    data_m:       in, required, type=sometype
;       A parameter named data_m
;    ephem:        in, required, type=sometype
;       A parameter named ephem
;    antenna_file: in, required, type=sometype
;       A parameter named antenna_file
;    sn:           in, required, type=sometype
;       A parameter named sn
;
; :Keywords:
;    VERBOSE:      in, optional, type=sometype
;       A keyword named VERBOSE
;    TEST:         in, optional, type=sometype
;       A keyword named TEST
;-
PRO dfb_main_full,df_result,data_p,data_m,ephem,antenna_file,sn, $
             VERBOSE=VERBOSE,TEST=TEST
; ====================================================================
; all angles are in radian
; ====================================================================
; --------------------------------------------------------------------

if antenna_file eq 'calDec04' then begin 
  read_antenna_set,antenna_HF,'calDec04_H12',/rad
  read_antenna_set,antenna_LF,'calDec04_ABC',/rad
endif else begin
  read_antenna_set,antenna_HF,antenna_file,/rad
  antenna_LF = antenna_HF
endelse

f_antenna_lim = 320.

ndim = n_elements(df_result)

; --------------------------------------------------------------------
;                    antenna parameters initialization
; --------------------------------------------------------------------

al_tmp  = [antenna_LF.Z.al,antenna_HF.Z.al]
be_tmp  = [antenna_LF.Z.be,antenna_HF.Z.be]
h__tmp  = [antenna_LF.Z.h, antenna_HF.Z.h ]
alZ = al_tmp( data_p.f ge f_antenna_lim )
beZ = be_tmp( data_p.f ge f_antenna_lim )
hZ  = h__tmp( data_p.f ge f_antenna_lim )

al_tmp = [antenna_LF.Xp.al,antenna_HF.Xp.al]
be_tmp = [antenna_LF.Xp.be,antenna_HF.Xp.be]
h__tmp = [antenna_LF.Xp.h, antenna_HF.Xp.h ]
alXp = al_tmp( data_p.f ge f_antenna_lim )
beXp = be_tmp( data_p.f ge f_antenna_lim )
hXp  = h__tmp( data_p.f ge f_antenna_lim )

al_tmp = [antenna_LF.Xm.al,antenna_HF.Xm.al]
be_tmp = [antenna_LF.Xm.be,antenna_HF.Xm.be]
h__tmp = [antenna_LF.Xm.h, antenna_HF.Xm.h ]
alXm = al_tmp( data_p.f ge f_antenna_lim )
beXm = be_tmp( data_p.f ge f_antenna_lim )
hXm  = h__tmp( data_p.f ge f_antenna_lim )

; --------------------------------------------------------------------
; reference frame rotation : S/C to Antenna
; --------------------------------------------------------------------

vhZ  = make_vect_sph(1.,alZ,beZ)
vhXp = make_vect_sph(1.,alXp,beXp)
vhXm = make_vect_sph(1.,alXm,beXm)

e3 = [0.,0.,1.] 

vr1 = crossp1(vhZ,e3,/norm)
qr1 = q_make(alZ,vr1)

qhZ = q_vmake(vhZ)
qhXp = q_vmake(vhXp)
qhXm = q_vmake(vhXm)

qhZ1 = q_rot(qr1,qhZ)
qhXp1 = q_rot(qr1,qhXp)
qhXm1 = q_rot(qr1,qhXm)

beXp1 = reform(atan(qhXp1(2,*),qhXp1(1,*)))     ; in radian
beXm1 = reform(atan(qhXm1(2,*),qhXm1(1,*)))     ; in radian

beta = ( !pi - (beXm1 - beXp1) )/2. ; in radian

qr2 = q_make(-beXp1+beta,e3)

qhZ2 = q_rot(qr2,qhZ1)
qhXp2 = q_rot(qr2,qhXp1)
qhXm2 = q_rot(qr2,qhXm1)

alXp2 = reform(acos(qhXp2(3,*)))                ; in radian
alXm2 = reform(acos(qhXm2(3,*)))                ; in radian

beXp2 = reform(atan(qhXp2(2,*),qhXp1(1,*)))     ; in radian
beXm2 = reform(atan(qhXm2(2,*),qhXm1(1,*)))     ; in radian

;if keyword_set(verbose) then begin
;  print,'Antenna Frame : '
;  print,'Alpha (X+) = ',alXp2(0)*!radeg
;  print,'Beta (X+)  = ',beXp2(0)*!radeg
;  print,'Alpha (X-) = ',alXm2(0)*!radeg
;  print,'Beta (X-)  = ',beXm2(0)*!radeg
;endif

; Computing ephemeris in antenna frame
; ------------------------------------

veph = make_vect_sph(1.,ephem.th,ephem.ph)
qeph = q_vmake(veph)

qeph1 = q_rot(qr1,qeph)
qeph2 = q_rot(qr2,qeph1)

ph_eph2 = (atan(qeph2(2,*),qeph2(1,*)) + 2*!pi) mod (2.*!pi)
th_eph2 = acos(qeph2(3,*))

;if keyword_set(test) then begin
  ph_eph2i = ph_eph2
  th_eph2i = th_eph2
;endif else begin
;  ph_eph2i = interpol(ph_eph2,ephem.time,ti_t97(data.ti))
;  th_eph2i = interpol(th_eph2,ephem.time,ti_t97(data.ti))
;endelse

; --------------------------------------------------------------------------
; Determination of Phi in the Antenna frame
; --------------------------------------------------------------------------

ph2 = (atan((hXp*sin(alXp2)*data_m.CrossI - $
             hXm*sin(alXm2)*data_p.CrossI)*tan(beta), $
            (hXp*sin(alXp2)*data_m.CrossI + $
             hXm*sin(alXm2)*data_p.CrossI) ) + 2.*!pi) mod (2.*!pi)

ph2_a = reform(ph2,ndim) 
ph2_b = reform(ph2,ndim) + !pi

; ---------------------------------------------------------------------------
; Determination of Theta in the antenna frame 
; ---------------------------------------------------------------------------

vvz = (data_p.AutoZ+data_m.AutoZ)/2.

Th2_a = atan(vvz*hXp*hXm*sin(alXp2)*sin(alXm2)*sin(2.*beta), $
           ((hXp*vvz*cos(alXp2)-hZ*data_p.crossR)*hXm*sin(alXm2)* $
            sin(ph2_a+beta)+(hXm*vvz*cos(alXm2)-hZ*data_m.crossR)* $
            hXp*sin(alXp2)*sin(ph2_a-beta)))

Th2_a =( Th2_a + !pi ) mod !pi

Th2_b = atan(vvz*hXp*hXm*sin(alXp2)*sin(alXm2)*sin(2.*beta), $
           ((hXp*vvz*cos(alXp2)-hZ*data_p.crossR)*hXm*sin(alXm2)* $
            sin(ph2_b+beta)+(hXm*vvz*cos(alXm2)-hZ*data_m.crossR)* $
            hXp*sin(alXp2)*sin(ph2_b-beta)))

Th2_b =( Th2_b + !pi ) mod !pi

; ---------------------------------------------------------------------------
; Preparing Useful Quantities...
; ---------------------------------------------------------------------------

sxp_a = hXp*(cos(alXp2)*sin(th2_a) - sin(alXp2)*cos(th2_a)*cos(ph2_a-beta))
sxm_a = hXm*(cos(alXm2)*sin(th2_a) + sin(alXm2)*cos(th2_a)*cos(ph2_a+beta))

sxp_b = hXp*(cos(alXp2)*sin(th2_b) - sin(alXp2)*cos(th2_b)*cos(ph2_b-beta))
sxm_b = hXm*(cos(alXm2)*sin(th2_b) + sin(alXm2)*cos(th2_b)*cos(ph2_b+beta))

syp_a = -hXp*sin(alXp2)*sin(ph2_a-beta)
sym_a = +hXm*sin(alXm2)*sin(ph2_a+beta)

syp_b = -hXp*sin(alXp2)*sin(ph2_b-beta)
sym_b = +hXm*sin(alXm2)*sin(ph2_b+beta)

sz_a  = hZ*sin(th2_a)

sz_b  = hZ*sin(th2_b)

; ---------------------------------------------------------------------------
; Determination of Flux (S)
; ---------------------------------------------------------------------------

Sp_a = (data_p.autoX*sz_a^2. - 2.*data_p.crossR*sxp_a*sz_a + $
        data_p.autoZ*(sxp_a^2.+syp_a^2.)) / (syp_a^2.*sz_a^2.)
Sm_a = (data_m.autoX*sz_a^2. - 2.*data_m.crossR*sxm_a*sz_a + $
        data_m.autoZ*(sxm_a^2.+sym_a^2.)) / (sym_a^2.*sz_a^2.)

Sp_b = (data_p.autoX*sz_b^2. - 2.*data_p.crossR*sxp_b*sz_b + $
        data_p.autoZ*(sxp_b^2.+syp_b^2.)) / (syp_b^2.*sz_b^2.)
Sm_b = (data_m.autoX*sz_b^2. - 2.*data_m.crossR*sxm_b*sz_b + $
        data_m.autoZ*(sxm_b^2.+sym_b^2.)) / (sym_b^2.*sz_b^2.)

; ---------------------------------------------------------------------------
; Determination of Stokes Parameter V (in antenna frame)
; ---------------------------------------------------------------------------

Vp_a = -2.*data_p.crossI/(Sp_a*sz_a*syp_a)
Vm_a = -2.*data_m.crossI/(Sm_a*sz_a*sym_a)

Vp_b = -2.*data_p.crossI/(Sp_b*sz_b*syp_b)
Vm_b = -2.*data_m.crossI/(Sm_b*sz_b*sym_b)

; ---------------------------------------------------------------------------
; Theta and Phi in S/C frame
; ---------------------------------------------------------------------------

; solution (a)
; ---------------------------------------------------------------------------

vsrc2 = make_vect_sph(1.,th2_a,ph2_a)
qsrc2 = Q_vmake(vsrc2)

qsrc1 = Q_rot(Q_conj(qr2),qsrc2)
qsrc = Q_rot(Q_conj(qr1),qsrc1)

th_a = reform(acos(qsrc(3,*)),ndim)
ph_a = reform((atan(qsrc(2,*),qsrc(1,*)) + 2*!pi) mod (2.*!pi),ndim)

ith = where(abs(th_a - !pi/2.) - !pi/2. ge -0.01, cnt)
if cnt ne 0 then ph_a(ith) = 0.

; solution (b)
; ---------------------------------------------------------------------------

vsrc2 = make_vect_sph(1.,th2_b,ph2_b)
qsrc2 = Q_vmake(vsrc2)

qsrc1 = Q_rot(Q_conj(qr2),qsrc2)
qsrc = Q_rot(Q_conj(qr1),qsrc1)

th_b = reform(acos(qsrc(3,*)),ndim)
ph_b = reform((atan(qsrc(2,*),qsrc(1,*)) + 2*!pi) mod (2.*!pi),ndim)

ith = where(abs(th_b - !pi/2.) - !pi/2. ge -0.01, cnt)
if cnt ne 0 then ph_b(ith) = 0.

; ---------------------------------------------------------------------------
; Determination of Stokes Parameters Q,U (linear polarization) in S/C frame
; ---------------------------------------------------------------------------

; solution (a)
; ---------------------------------------------------------------------------

sxz = hZ*(cos(alZ)*sin(th_a)-sin(alZ)*cos(th_a)*cos(ph_a-beZ))
sxp = hXp*(cos(alXp)*sin(th_a)-sin(alXp)*cos(th_a)*cos(ph_a-beXp))
sxm = hXm*(cos(alXm)*sin(th_a)-sin(alXm)*cos(th_a)*cos(ph_a-beXm))

syz = -hZ*sin(alZ)*sin(ph_a-beZ)
syp = -hXp*sin(alXp)*sin(ph_a-beXp)
sym = -hXm*sin(alXm)*sin(ph_a-beXm)

Up_a = (data_p.autoZ*(sxp^2.-syp^2.) - data_p.autoX*(sxz^2.-syz^2.))/ $
       (Sp_a*(syz*sxp-syp*sxz)*(sxz*sxp+syz*syp)) - (syz*sxp+syp*sxz)/ $
       (sxz*sxp+syz*syp)
Um_a = (data_m.autoZ*(sxm^2.-sym^2.) - data_m.autoX*(sxz^2.-syz^2.))/ $
       (Sm_a*(syz*sxm-sym*sxz)*(sxz*sxm+syz*sym)) - (syz*sxm+sym*sxz)/ $
       (sxz*sxm+syz*sym)


Qp_a = 2.*(data_p.autoX*sxz*syz - data_p.autoZ*sxp*syp)/ $
       (Sp_a*(syz*syp+sxz*sxp)*(syz*sxp-syp*sxz)) + $
       (syz*syp-sxz*sxp)/(syz*syp+sxz*sxp)
Qm_a = 2.*(data_m.autoX*sxz*syz - data_m.autoZ*sxm*sym)/ $
       (Sm_a*(syz*sym+sxz*sxm)*(syz*sxm-sym*sxz)) + $
       (syz*sym-sxz*sxm)/(syz*sym+sxz*sxm)

; solution (b)
; ---------------------------------------------------------------------------

sxz = hZ*(cos(alZ)*sin(th_b)-sin(alZ)*cos(th_b)*cos(ph_b-beZ))
sxp = hXp*(cos(alXp)*sin(th_b)-sin(alXp)*cos(th_b)*cos(ph_b-beXp))
sxm = hXm*(cos(alXm)*sin(th_b)-sin(alXm)*cos(th_b)*cos(ph_b-beXm))

syz = -hZ*sin(alZ)*sin(ph_b-beZ)
syp = -hXp*sin(alXp)*sin(ph_b-beXp)
sym = -hXm*sin(alXm)*sin(ph_b-beXm)

Up_b = (data_p.autoZ*(sxp^2.-syp^2.) - data_p.autoX*(sxz^2.-syz^2.))/ $
       (Sp_b*(syz*sxp-syp*sxz)*(sxz*sxp+syz*syp)) - (syz*sxp+syp*sxz)/ $
       (sxz*sxp+syz*syp)
Um_b = (data_m.autoZ*(sxm^2.-sym^2.) - data_m.autoX*(sxz^2.-syz^2.))/ $
       (Sm_b*(syz*sxm-sym*sxz)*(sxz*sxm+syz*sym)) - (syz*sxm+sym*sxz)/ $
       (sxz*sxm+syz*sym)


Qp_b = 2.*(data_p.autoX*sxz*syz - data_p.autoZ*sxp*syp)/ $
       (Sp_b*(syz*syp+sxz*sxp)*(syz*sxp-syp*sxz)) + $
       (syz*syp-sxz*sxp)/(syz*syp+sxz*sxp)
Qm_b = 2.*(data_m.autoX*sxz*syz - data_m.autoZ*sxm*sym)/ $
       (Sm_b*(syz*sym+sxz*sxm)*(syz*sxm-sym*sxz)) + $
       (syz*sym-sxz*sxm)/(syz*sym+sxz*sxm)

; ---------------------------------------------------------------------------
; Filling output structure
; ---------------------------------------------------------------------------


s = reform([transpose([[sp_a],[sm_a]]),transpose([[sp_b],[sm_b]])],2,2,ndim)  
q = reform([transpose([[qp_a],[qm_a]]),transpose([[qp_b],[qm_b]])],2,2,ndim)  
u = reform([transpose([[up_a],[um_a]]),transpose([[up_b],[um_b]])],2,2,ndim)  
v = reform([transpose([[vp_a],[vm_a]]),transpose([[vp_b],[vm_b]])],2,2,ndim)  

df_result.ydh = data_p.ydh
df_result.num = transpose(reform([data_p.num,data_m.num],ndim,2))
df_result.s  = s
df_result.q  = q
df_result.u  = u
df_result.v  = v
df_result.th = reform(transpose([[th_a],[th_b]]),2,ndim)
df_result.ph = reform(transpose([[ph_a],[ph_b]]),2,ndim)
df_result.zr = reform((data_p.autoZ-data_m.autoZ)/(data_p.autoZ+data_m.autoZ),ndim)
df_result.sn = sn

return
end