;+ ; Contains the dfc_main function ; ; :Author: ; Baptiste Cecconi ; ; :History: ; 2008/01/28: Created ; ; 2008/01/28: Last Edit ;- ; ;+ ; dfc_main is a procedure that <behavior desc here> ; ; :Uses: ; crossp1, read_antenna_set ; ; :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 ; nan_exception: in, required, type=sometype ; A parameter named nan_exception ; ; :Keywords: ; VERBOSE: in, optional, type=sometype ; A keyword named VERBOSE ; TEST: in, optional, type=sometype ; A keyword named TEST ;- PRO dfc_main,df_result,data_p,data_m,ephem,antenna_file,sn,nan_exception, $ 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) nan_exception = bytarr(ndim) ;t97 = reform(rebin(reform(data.t97,2,ndim),1,ndim),ndim) ;wp = where(data.ant eq 11b) ;wm = where(data.ant eq 12b) ;vv = transpose(reform([data(wp).autoX,data(wp).autoZ,data(wm).autoX,data(wm).autoZ, $ ; data(wp).crossR,data(wp).crossI,data(wm).crossR,data(wm).crossI],ndim,8)) ; ------------------------------------------------------------------------------ ; 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,*)) ; -------------------------------------------------------------------------- ; Determination of Phi in the Antenna frame ; -------------------------------------------------------------------------- wzp0 = where(data_p.autoZ ne 0.) wzm0 = where(data_m.autoZ ne 0.) ; -------------------------------------------------------------------------- ; B+ and B- variables ; -------------------------------------------------------------------------- ; B+*h+x^2*sin^2(alXp) = A+xx - Cr+xz^2/Azz ; B-*h-x^2*sin^2(alXm) = A-xx - Cr-xz^2/Azz ; -------------------------------------------------------------------------- Bp = fltarr(ndim) Bm = fltarr(ndim) Bp(wzp0) = 2.*(data_p(wzp0).autoX - $ data_p(wzp0).crossR/data_p(wzp0).autoZ*data_p(wzp0).crossR)/$ (hXp(wzp0)*sin(alXp2(wzp0)))^2. Bm(wzm0) = 2.*(data_m(wzm0).autoX - $ data_m(wzm0).crossR/data_m(wzm0).autoZ*data_m(wzm0).crossR)/$ (hXm(wzm0)*sin(alXm2(wzm0)))^2. TwoTheta = atan((Bp+Bm)*sin(2.*beXp2),(Bp-Bm)*cos(2.*beXp2)) w = where(abs(cos(TwoTheta)/cos(2.*beXp2)) gt 1.) if w(0) ne -1 then begin message,'warning : [Ph2] '+string(n_elements(w))+' points may be wrong see nan_exception report...',/info nan_exception(w) = nan_exception(w) + 1 endif ph2p = 0.5*(acos(cos(TwoTheta)/cos(2.*beXp2)<1.>(-1.))- TwoTheta) ph2m = !pi - 0.5*(acos(cos(TwoTheta)/cos(2.*beXp2)<1.>(-1.))+ TwoTheta) ;stop delta_ph = abs(reform([ [ph2p-ph_eph2-2.*!pi], $ [ph2p-ph_eph2-!pi], $ [ph2p-ph_eph2], $ [ph2p-ph_eph2+!pi], $ [ph2p-ph_eph2+2.*!pi], $ [ph2m-ph_eph2-2.*!pi], $ [ph2m-ph_eph2-!pi], $ [ph2m-ph_eph2], $ [ph2m-ph_eph2+!pi], $ [ph2m-ph_eph2+2.*!pi] $ ],ndim,10)) idelta_ph = bytarr(ndim) for i=0l,ndim-1 do idelta_ph(i)=(where(delta_ph(i,*) eq min(delta_ph(i,*))))(0) ph2 = (idelta_ph le 4)*ph2p + (idelta_ph gt 5)*ph2m + ((idelta_ph mod 5) - 2)*!pi ; -------------------------------------------------------------------------- ; Determination of Flux S ; -------------------------------------------------------------------------- S = Bp+Bm - (Bp-Bm)/tan(2.*ph2)/tan(2.*beXp2) ; Q = 0. ; U = 0. ; -------------------------------------------------------------------------- ; Determination of Theta in the Antenna frame ; -------------------------------------------------------------------------- wws0 = where(s ne 0.) wwp = where(sqrt(2.*data_p(wws0).autoZ/(S(wws0)*hZ(wws0)^2.)) gt 1.) wwm = where(sqrt(2.*data_m(wws0).autoZ/(S(wws0)*hZ(wws0)^2.)) gt 1.) if wwp(0) ne -1 then begin message,'warning : [Th2] '+string(n_elements(wwp))+' points may be wrong see nan_exception report...',/info nan_exception(wwp) = nan_exception(wwp)+1 endif if wwm(0) ne -1 then begin message,'warning : [Th2] '+string(n_elements(wwm))+' points may be wrong see nan_exception report...',/info nan_exception(wwm) = nan_exception(wwm)+1 endif th2p = fltarr(ndim) th2m = fltarr(ndim) th2p(wws0) = asin(sqrt(2.*data_p(wws0).autoZ/(S(wws0)*hZ(wws0)^2.))<1.) th2m(wws0) = asin(sqrt(2.*data_m(wws0).autoZ/(S(wws0)*hZ(wws0)^2.))<1.) delta_thp = abs(reform([[th2p-th_eph2],[!pi-th2p-th_eph2]],ndim,2)) delta_thm = abs(reform([[th2m-th_eph2],[!pi-th2m-th_eph2]],ndim,2)) idelta_thp = bytarr(ndim) idelta_thm = bytarr(ndim) for i=0l,ndim-1 do idelta_thp(i)=(where(delta_thp(i,*) eq min(delta_thp(i,*))))(0) for i=0l,ndim-1 do idelta_thm(i)=(where(delta_thm(i,*) eq min(delta_thm(i,*))))(0) th2p = th2p + (idelta_thp eq 1)*(!pi - 2.*th2p) th2m = th2m + (idelta_thm eq 1)*(!pi - 2.*th2m) th2p = reform(th2p) th2m = reform(th2m) ; -------------------------------------------------------------------------- ; Determination of Circular Polarization V [added in v1.0] ; -------------------------------------------------------------------------- svp = ( (sin(th2p)*sin(alXp2)*sin(ph2-beXp2)) ge 0. )*2. -1. svm = ( (sin(th2m)*sin(alXm2)*sin(ph2-beXm2)) ge 0. )*2. -1. Vp = fltarr(ndim) Vm = fltarr(ndim) wwp0 = where(data_p.autoX ne 0 and data_p.autoZ ne 0) wwm0 = where(data_m.autoX ne 0 and data_m.autoZ ne 0) Vp(wwp0) = float(data_p(wwp0).crossI/sqrt(data_p(wwp0).autoX*data_p(wwp0).autoZ-double(data_p(wwp0).crossR)^2.)) * svp(wwp0) Vm(wwm0) = float(data_m(wwm0).crossI/sqrt(data_m(wwm0).autoX*data_m(wwm0).autoZ-double(data_m(wwm0).crossR)^2.)) * svm(wwm0) ; --------------------------------------------------------------------------- ; Theta and Phi in S/C frame ; --------------------------------------------------------------------------- vsrc2p = make_vect_sph(1.,th2p,ph2) vsrc2m = make_vect_sph(1.,th2m,ph2) qsrc2p = Q_vmake(vsrc2p) qsrc2m = Q_vmake(vsrc2m) qsrc1p = Q_rot(Q_conj(qr2),qsrc2p) qsrcp = Q_rot(Q_conj(qr1),qsrc1p) qsrc1m = Q_rot(Q_conj(qr2),qsrc2m) qsrcm = Q_rot(Q_conj(qr1),qsrc1m) thp = reform(acos(qsrcp(3,*))) php = reform((atan(qsrcp(2,*),qsrcp(1,*)) + 2*!pi) mod (2.*!pi)) thm = reform(acos(qsrcm(3,*))) phm = reform((atan(qsrcm(2,*),qsrcm(1,*)) + 2*!pi) mod (2.*!pi)) ith = where(abs(thp - !pi/2.) - !pi/2. ge -0.01) if ith(0) ne -1 then php(ith) = 0. ith = where(abs(thm - !pi/2.) - !pi/2. ge -0.01) if ith(0) ne -1 then phm(ith) = 0. df_result.ydh = data_p.ydh df_result.num = transpose(reform([data_p.num,data_m.num],ndim,2)) df_result.s = reform(S,ndim) df_result.q = 0. df_result.u = 0. df_result.v = transpose(reform([Vp,Vm],ndim,2)) df_result.th = transpose(reform([thp,thm],ndim,2)) df_result.ph = transpose(reform([php,phm],ndim,2)) df_result.zr = reform((data_p.autoZ-data_m.autoZ)/(data_p.autoZ+data_m.autoZ),ndim) df_result.sn = sn return end