;+ ; Contains the dfd_main procedure ; ; :Author: ; Baptiste Cecconi ; ; :History: ; 2006/03/08: Created ; ; 2006/03/08: Last Edit ;- ; ;+ ; dfd_main is a procedure that <behavior desc here> ; ; :Uses: ; read_antenna_set ; ; :Params: ; df_result: in, required, type=sometype ; A parameter named df_result ; data: in, required, type=sometype ; A parameter named data ; 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 dfd_main,df_result,data,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 ; -------------------------------------------------------------------- alZtmp = [antenna_LF.Z.al,antenna_HF.Z.al] alZ = alZtmp(data.f ge f_antenna_lim) beZtmp = [antenna_LF.Z.be,antenna_HF.Z.be] beZ = beZtmp(data.f ge f_antenna_lim) hZtmp = [antenna_LF.Z.h,antenna_HF.Z.h] hZ = hZtmp(data.f ge f_antenna_lim) alXtmp = [[0.,antenna_LF.Xp.al,antenna_LF.Xm.al,antenna_LF.Dip.al],[0.,antenna_HF.Xp.al,antenna_HF.Xm.al,antenna_HF.Dip.al]] alX = alXtmp(data.ant mod 10,data.f ge f_antenna_lim) beXtmp = [[0.,antenna_LF.Xp.be,antenna_LF.Xm.be,antenna_LF.Dip.be],[0.,antenna_HF.Xp.be,antenna_HF.Xm.be,antenna_HF.Dip.be]] beX = beXtmp(data.ant mod 10,data.f ge f_antenna_lim) hXtmp = [[0.,antenna_LF.Xp.h,antenna_LF.Xm.h,antenna_LF.Dip.h],[0.,antenna_HF.Xp.h,antenna_HF.Xm.h,antenna_HF.Dip.h]] hX = hXtmp(data.ant mod 10,data.f ge f_antenna_lim) ; -------------------------------------------------------------------- ; reference frame rotation : S/C to Antenna ; -------------------------------------------------------------------- vhZ = make_vect_sph(1.,alZ,beZ) vhX = make_vect_sph(1.,alX,beX) e3 = [0.,0.,1.] vr1 = crossp(vhZ,e3) vr1 = vr1/sqrt(total(vr1^2.)) qr1 = q_make(alZ,vr1) qhZ = q_vmake(vhZ) qhX = q_vmake(vhX) qhZ1 = q_rot(qr1,qhZ) qhX1 = q_rot(qr1,qhX) alX1 = acos(qhX1(3,*)<1.>(-1.)) ; en radian beX1 = atan(qhX1(2,*),qhX1(1,*)) ; en radian ; --------------------------------------------------------------------------- ; Computing ephemeris in antenna frame ; --------------------------------------------------------------------------- th = ephem.th ph = ephem.ph veph = make_vect_sph(1.,th,ph) qeph = q_vmake(veph) qeph1 = q_rot(qr1,qeph) ph1 = (atan(qeph1(2,*),qeph1(1,*)) + 2*!pi) mod (2.*!pi) th1 = acos(qeph1(3,*)<1.>(-1.)) ; --------------------------------------------------------------------------- ; Preparing Useful Quantities... ; --------------------------------------------------------------------------- sx = hX*(cos(alX1)*sin(th1) - sin(alX1)*cos(th1)*cos(ph1-beX1)) sy = -hX*sin(alX1)*sin(ph1-beX1) sz = hZ*sin(th1) S = fltarr(ndim) V = S & U = S & Q = S dataValid = bytarr(ndim) ; --------------------------------------------------------------------------- ; Determination of Flux S (in antenna frame) ; --------------------------------------------------------------------------- S = (data.autoX*sz^2. - 2.*data.crossR*sx*sz + $ data.autoZ*(sx^2.+sy^2.)) / (sy^2.*sz^2.) ws = where(s ne 0.) if keyword_set(verbose) then message,'% S=0 for '+string(ndim-n_elements(ws))+' data points ('+string(format='(F5.1)',100.*(ndim-n_elements(ws))/ndim)+'%).',/info if ws(0) ne -1 then dataValid(ws) = 1b ; --------------------------------------------------------------------------- ; Determination of Stokes Parameter V (in antenna frame) ; --------------------------------------------------------------------------- ww = where(dataValid) V(ww) = -2.*data(ww).crossI/(S(ww)*sz(ww)*sy(ww)) ; signe modifie / Baptiste !!! cf. equations DF ; --------------------------------------------------------------------------- ; Determination of Stokes Parameters Q,U (linear polarization) in SRC frame ; --------------------------------------------------------------------------- ; ------------------------------------- ; reference frame rotation : S/C to SRC ; ------------------------------------- vk = -veph qk = q_vmake(vk) vr1 = make_vect_sph(1.,!pi/2.,!pi/2.+ph) qr1 = q_make(!pi-th,vr1) vZp = make_vect_sph(1.,ephem.thZp,ephem.phZp) qZp = q_vmake(vZp) qZp1 = q_rot(qr1,qZp) phZp1=(atan(qZp1(2,*),qZp1(1,*)) + 2*!pi) mod (2.*!pi) qr2 = q_make(reform(phZp1),e3) qhX2 = q_rot(qr2,q_rot(qr1,qhX)) qhZ2 = q_rot(qr2,q_rot(qr1,qhZ)) alX2 = acos(qhX2(3,*)<1.>(-1.)) beX2 = (atan(qhX2(2,*),qhX2(1,*)) + 2*!pi) mod (2.*!pi) alZ2 = acos(qhZ2(3,*)<1.>(-1.)) beZ2 = (atan(qhZ2(2,*),qhZ2(1,*)) + 2*!pi) mod (2.*!pi) ; +---------------------+ ; | NB : th=-!pi, ph=0. | ; +---------------------+ ; | -> sin(th) = 0. | ; | cos(th) = -1. | ; +---------------------+ sxz = hZ*sin(alZ2)*cos(beZ2) sxx = hX*sin(alX2)*cos(beX2) syz = hZ*sin(alZ2)*sin(beZ2) syx = hX*sin(alX2)*sin(beX2) U(ws) = (data(ws).autoZ*(sxx(ws)^2.-syx(ws)^2.) - data(ws).autoX*(sxz(ws)^2.-syz(ws)^2.))/ $ (S(ws)*(syz(ws)*sxx(ws)-syx(ws)*sxz(ws))*(sxz(ws)*sxx(ws)+syz(ws)*syx(ws))) $ - (syz(ws)*sxx(ws)+syx(ws)*sxz(ws))/ $ (sxz(ws)*sxx(ws)+syz(ws)*syx(ws)) Q(ws) = 2.*(data(ws).autoX*sxz(ws)*syz(ws) - data(ws).autoZ*sxx(ws)*syx(ws))/ $ (S(ws)*(syz(ws)*syx(ws)+sxz(ws)*sxx(ws))*(syz(ws)*sxx(ws)-syx(ws)*sxz(ws))) + $ (syz(ws)*syx(ws)-sxz(ws)*sxx(ws))/(syz(ws)*syx(ws)+sxz(ws)*sxx(ws)) ;stop ; --------------------------------------------------------------------------- ; Filling data_n3d structure with DFd results ; --------------------------------------------------------------------------- df_result.ydh = data.ydh df_result.num = data.num df_result.s = reform(s,ndim) df_result.q = reform(q,ndim) df_result.u = reform(u,ndim) df_result.v = reform(v,ndim) df_result.th = reform(th,ndim) df_result.ph = reform(ph,ndim) df_result.sn = sn return end