;+ ; Contains the dfe_main procedure ; ; :Author: ; Baptiste Cecconi ; ; :History: ; 2009/06/02: Created ; ; 2009/06/02: Last Edit ;- ; ;+ ; dfe_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: 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 ; 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 dfe_main,df_result,data,ephem,antenna_file,sn,nan_exception, $ VERBOSE=VERBOSE,TEST=TEST ; ----------------------------------------------------------------------------------- ; all angles are in radian ; v1.0 [BC, 2005-jul-06] : first version ; v1.1 [BC, 2009-jun-02] : best solution selection modified ; ----------------------------------------------------------------------------------- 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) ; -------------------------------------------------------------------- ; 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) ; --------------------------------------------------------------------------- ; Computing ephemeris in 2-antenna system frame ; NB: ; 2-antenna system frame = hX and hZ in xOy plane, y bissecting antenna directions ; --------------------------------------------------------------------------- vhZ = make_vect_sph(1.,alZ,beZ) vhX = make_vect_sph(1.,alX,beX) qhZ = q_vmake(vhZ) qhX = q_vmake(vhX) qsrc = q_vmake(make_vect_sph(1,ephem.th,ephem.ph)) e3 = [0.,0.,1.] vz1 = crossp1(vhz,vhx,/norm) vr1 = crossp1(vz1,e3,/norm) be1 = angular_distance(vz1,e3) qr1 = q_make(be1,vr1) qZ1 = q_rot(qr1,qhZ) qX1 = q_rot(qr1,qhX) phx1 = atan(qx1(2,*),qx1(1,*)) phZ1 = atan(qz1(2,*),qz1(1,*)) be2 = reform((!pi - (phZ1-phx1))/2.-phX1) qr2 = q_make(be2,e3) qz2 = q_rot(qr2,qz1) qx2 = q_rot(qr2,qx1) alX2 = acos(qX2(3,*)) alZ2 = acos(qZ2(3,*)) beX2 = atan(qX2(2,*),qX2(1,*)) beZ2 = atan(qZ2(2,*),qZ2(1,*)) qsrc2 = q_rot(qr2,q_rot(qr1,qsrc)) theph2 = acos(qsrc2(3,*)) pheph2 = atan(qsrc2(2,*),qsrc2(1,*)) ;!p.multi=[0,2,2] ;plot,alz2,psym=1,ytit='alz' ;plot,alx2,psym=1,ytit='alx' ;plot,bez2,psym=1,ytit='bez' ;plot,bex2,psym=1,ytit='bex' ;stop ; --------------------------------------------------------------------------- ; Preparing Useful Quantities... ; --------------------------------------------------------------------------- P = fltarr(ndim) & Q=P & R=P ws = where(data.autoZ ne 0. or data.autoX ne 0.) P(ws) = 2.*hX(ws)*hZ(ws)*data(ws).crossR/ $ (hx(ws)^2.*data(ws).autoZ+hZ(ws)^2.*data(ws).autoX) Q(ws) = (hx(ws)^2.*data(ws).autoZ-hZ(ws)^2.*data(ws).autoX)/ $ (hx(ws)^2.*data(ws).autoZ+hZ(ws)^2.*data(ws).autoX) R(ws) = 2.*hX(ws)*hZ(ws)*data(ws).crossI/ $ (hx(ws)^2.*data(ws).autoZ+hZ(ws)^2.*data(ws).autoX) c = sin(alX2)*sin(alZ2)*cos(beX2-beZ2)+cos(alX2)*cos(alZ2) Pc = (P-c)/(1.-P*c) Qc = Q*sqrt(1.-c^2.)/(1.-P*c) Rc = R*sqrt(1.-c^2.)/(1.-P*c) ; --------------------------------------------------------------------------- ; Determination of azimuth Ph ; --------------------------------------------------------------------------- ph2 = fltarr(ndim) ph2(ws) = 0.5*atan(-Qc(ws),Pc(ws)) ph2a = ph2 ph2b = ph2 ph2c = ph2+!pi ph2d = ph2+!pi ; --------------------------------------------------------------------------- ; Determination of colatitude Th ; --------------------------------------------------------------------------- th2 = acos(sqrt((2./(1.+sqrt(Pc^2.+Qc^2.))-1.)>0.)<1.>(-1.)) th2a = th2 th2b = !pi-th2 th2c = th2 th2d = !pi-th2 ; --------------------------------------------------------------------------- ; Selecting source with input guessed ephemeris ; --------------------------------------------------------------------------- dtha = angular_distance(make_vect_sph(1,th2a,ph2a),qsrc2(1:3,*)) dthb = angular_distance(make_vect_sph(1,th2b,ph2b),qsrc2(1:3,*)) dthc = angular_distance(make_vect_sph(1,th2c,ph2c),qsrc2(1:3,*)) dthd = angular_distance(make_vect_sph(1,th2d,ph2d),qsrc2(1:3,*)) dth = reform([dtha,dthb,dthc,dthd],ndim,4) dthmin = min(dth,idthmin,dimension=2) th2 = ([th2a,th2b,th2c,th2d])(idthmin) ph2 = ([ph2a,ph2b,ph2c,ph2d])(idthmin) ; --------------------------------------------------------------------------- ; Determination of flux S ; --------------------------------------------------------------------------- S = (data.autoX/hX^2.+data.autoZ/hZ^2.)*(1.+sqrt(Pc^2.+Qc^2.))/(1.+Pc*c) ; --------------------------------------------------------------------------- ; Determination of Circular polarization degree V ; --------------------------------------------------------------------------- V = Rc/(1.+sqrt(Pc^2.+Qc^2.))/cos(th2) ; --------------------------------------------------------------------------- ; Back into S/C frame from angles ; --------------------------------------------------------------------------- src2 = make_vect_sph(1.,th2,ph2) qsrc2 = q_vmake(src2) qsrc = q_rot(q_conj(qr1),q_rot(q_conj(qr2),qsrc2)) th = reform(acos(qsrc(3,*)<1.>(-1.))) ph = reform(atan(qsrc(2,*),qsrc(1,*))) mod (2.*!pi) ; --------------------------------------------------------------------------- ; Filling data_n3e structure with DFe results ; --------------------------------------------------------------------------- df_result.ydh = data.ydh df_result.num = data.num df_result.s = reform(s,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