;+ ; Contains the localization procedure ; ; :Author: ; Baptiste Cecconi, Laurent Lamy ; ; :History: ; 2009/02/20: Created ; ; 2010/08/26: Last Edit ;- ; ;+ ; Instead of computing a grid of magnetic field lines in advance, we use a ; dichotomic algorithm to obtain the intersection of each direction of arrival ; with the correponding iso-fc surface, on its visible side as seen from ; Cassini. If there is no intersection, we determine the line which is ; perpendicular to the direction of arrival and which passes through the center ; of Saturn. The "intersection point" is then the intersection between the ; iso-fc surface and that latter line. We also store the distance between the ; direction of arrival and the iso-fc surface, measured on that line. ; Once the "intersection point" is determined, we follow he magnetic field line ; down to the surface of Saturn (top of the ionosphere) to find the foot print ; position of that field line, using the MFL_FOOT(V30) routine. ; ; we use iso-fc and not iso-fx as the locus of SKR emission regions. ; ; :Uses: ; iso_fc_intersect, load_data_ephem, load_data_n3edi, magnetic_field, mfl_foot, xyz_to_rtp ; ; :Params: ; aj1: in, required, type=sometype ; A parameter named aj1 ; h1: in, required, type=sometype ; A parameter named h1 ; ; :Keywords: ; mfl_model: in, optional, type=sometype ; A keyword named mfl_model ; test: in, optional, type=sometype ; A keyword named test ; VERBOSE: in, optional, type=sometype ; A keyword named VERBOSE ; harmonic: in, optional, type=sometype ; HARMONIC keyword to compute localization supposing we observe the ; harmonic. ; output_path: in, optional, type=sometype ; A keyword named output_path ; nocrt: in, optional, type=sometype ; A keyword named nocrt ; source: in, optional, type=sometype ; A keyword named source ; debug: in, optional, type=sometype ; A keyword named debug ;- PRO LOCALIZATION, aj1, h1, mfl_model=mfl_model,test=test,$ VERBOSE=verbose, harmonic=harmonic, $ output_path=output_path, nocrt=nocrt, source=source, $ debug=debug ; +++++++++++++++++++++ ; + localization_V30 + ; +++++++++++++++++++++ ; localization,2006016,0,mfl_model='Z3',/verb t0 = systime(/sec) if keyword_set(verbose) then verbose=1b else verbose=0b if keyword_set(debug) then debug=1b else debug=0b if keyword_set(mfl_model) then mfl_model=mfl_model else mfl_model='SPV' if keyword_set(test) then test_enabled=1b else test_enabled=0b if ~keyword_set(source) then source='sq' if aj1 eq 0 then begin read_data_binary,'test/R0000000.00',n1,level='n1' read_data_binary,'test/P0000000.00',n2,level='n2' read_data_binary,'test/N3d_xxx0000000.00',n3d,level='n3d' read_data_binary,'test/N3e_xxx0000000.00',n3e,level='n3e' read_data_binary,'test/INDEX0000000.00',index,level='index' read_data_binary,'test/0000000.vxxx',veph,level='veph' read_data_binary,'test/0000000.qxxx',qeph,level='qeph' ndata = n_elements(n2) endif else begin if verbose then message,/info,'Loading RPWS data.' load_data_n3edi,aj1, h1, aj1, h1+1, n1,n2,n3d,n3e,index,source='sq',verbose=verbose if debug then help,/heap ndata = n_elements(n2) if ndata eq 0 || (ndata eq 1 and n2(0).t97 eq 0) then return if verbose then message,/info,'Loading EPHEM data.' load_data_ephem,aj1, h1, aj1, h1+1,veph,qeph,veph='vsat',qeph='qssq',t97_ramp = n2.t97 if debug then help,/heap endelse if keyword_set(harmonic) then n2.f = n2.f/2. read_data_src_list,src_list src_data = src_list(where(src_list.df_name eq source)) rp = src_data.BODY_RAD_VALUE if verbose then message,/info,'Preparing GP Vectors.' qv_SC = q_vmake(-veph.r) qtmp = q_rot(qeph.q,qv_SC) v_SC = qtmp(1:3,*)/rp ;(in Rsat) l_SC = sqrt(total(v_SC^2.,1)) qv_DF = q_vmake(make_vect_sph(1.,n3e.th,n3e.ph)) qtmp = q_rot(qeph.q,qv_DF) v_DF = qtmp(1:3,*) a_SP = angular_distance(v_DF,-v_SC) fsb=2.79924835996d0 ; converts Gauss to MHz fsb=fsb*1.d3 ; converts Gauss to kHz ; looking for maximum Fc between V_SP0 (ie Cassini) and V_SP1 v_SP0 = v_SC v_SP1 = v_SC + 2.*v_DF*rebin(reform(-(v_SC(0,*)^2.+v_SC(1,*)^2.)/(v_DF(0,*)*v_SC(0,*)+v_DF(1,*)*v_SC(1,*)),1,ndata),3,ndata) ndic = 30L for idic=0L,ndic-1L do begin v_SP2 = (v_SP0+v_SP1)/2 v_SP2a = v_SP2 + v_DF/10./float(idic+1) v_SP2b = v_SP2 - v_DF/10./float(idic+1) rtp2a = xyz_to_rtp(v_SP2a) magnetic_field,rtp2a(0,*),rtp2a(1,*),rtp2a(2,*),bvec2a,btot2a,mfl_model=mfl_model,nocrt=nocrt fc2a = btot2a*fsb rtp2b = xyz_to_rtp(v_SP2b) magnetic_field,rtp2b(0,*),rtp2b(1,*),rtp2b(2,*),bvec2b,btot2b,mfl_model=mfl_model,nocrt=nocrt fc2b = btot2b*fsb dfc2 = fc2a-fc2b wwa = where(dfc2 ge 0, cnta, compl=wwb, ncompl=cntb) if cnta ne 0 then begin v_SP0(*,wwa) = v_SP2(*,wwa) endif if cntb ne 0 then begin v_SP1(*,wwb) = v_SP2(*,wwb) endif endfor v_SP = (v_SP0+v_SP1)/2. ; fixing v_SP0 to where Fc is maximum on DF direction : v_SP0 = V_SP rtp0 = xyz_to_rtp(v_SP0) magnetic_field,rtp0(0,*),rtp0(1,*),rtp0(2,*),bvec0,btot0,mfl_model=mfl_model,nocrt=nocrt fc0 = btot0*fsb wdata_in = where(fc0 ge n2.f,cnt_in,compl=wdata_out,ncompl=cnt_out) xyz_fc = fltarr(3,ndata) foot_fc = fltarr(3,ndata) beam_fc = fltarr(ndata) dist_fc = fltarr(ndata) if verbose then message,/info,'Analyzing GP results:' if cnt_in ne 0 then begin if verbose then message,/info,'+ when GP in iso_fc ('+strtrim(string(cnt_in),2)+' data points).' tc=systime(/sec) f_in = n2(wdata_in).f iso_fc_intersect,v_SC(*,wdata_in),v_SP(*,wdata_in),f_in, $ beam_fc_in,dist_fc_in,xyz_fc_in,30,/in,mfl_model=mfl_model,nocrt=nocrt xyz_fc(*,wdata_in) = xyz_fc_in beam_fc(wdata_in) = beam_fc_in dist_fc(wdata_in) = dist_fc_in rtp_fc_in = xyz_to_rtp(xyz_fc_in) if verbose then message,/info,' determining mfl foot print.' rtp_mfl_in = mfl_foot(reform(rtp_fc_in(0,*)), $ reform(rtp_fc_in(2,*))*!radeg, $ reform(rtp_fc_in(1,*))*!radeg, $ f_in,mfl_model=mfl_model,test=test_enabled,nocrt=nocrt) if debug then help,/heap foot_fc(*,wdata_in) = rtp_mfl_in if verbose then message,/info,' done in '+string(format='(F5.2," sec.")',systime(/sec)-tc) endif if cnt_out ne 0 then begin if verbose then message,/info,'+ when GP out of iso_fc ('+strtrim(string(cnt_out),2)+' data points).' tc=systime(/sec) f_out = n2(wdata_out).f iso_fc_intersect,v_SC(*,wdata_out),v_SP(*,wdata_out),f_out, $ beam_fc_out,dist_fc_out,xyz_fc_out,25,/out,mfl_model=mfl_model,nocrt=nocrt if debug then help,/heap xyz_fc(*,wdata_out) = xyz_fc_out beam_fc(wdata_out) = beam_fc_out dist_fc(wdata_out) = dist_fc_out if verbose then message,/info,' determining mfl foot print.' rtp_fc_out = xyz_to_rtp(xyz_fc_out) rtp_mfl_out = mfl_foot(reform(rtp_fc_out(0,*)), $ reform(rtp_fc_out(2,*))*!radeg, $ reform(rtp_fc_out(1,*))*!radeg, $ f_out,mfl_model=mfl_model,test=test_enabled,nocrt=nocrt) foot_fc(*,wdata_out) = rtp_mfl_out if verbose then message,/info,' done in '+string(format='(F5.2," sec.")',systime(/sec)-tc) endif data_loc = replicate({data_loc},n_elements(n2)) data_loc.ydh = n2.ydh data_loc.num = n2.num data_loc.foot = foot_fc data_loc.xyz = xyz_fc data_loc.beam = beam_fc data_loc.dist = dist_fc if ~keyword_set(output_path) then output_path='./' else $ if strmid(output_path,0,1,/reverse) ne '/' then output_path=output_path+'/' file = output_path+'loc_'+mfl_model+ $ (keyword_set(nocrt)?'_nocrt':'')+ $ (keyword_set(harmonic)?'_h':'')+ $ string(format='("_",I7.7,".",I2.2)',aj1,h1) if verbose then message,/info,"Outputing data in "+file write_data_binary,file,data_loc if verbose then message,/info,string(format='("Everything done in ",F6.2," sec.")',systime(/sec)-t0) ; Update the Kronos Database with the new files spawn, [getenv('ROOT_RPWS') + '/pro/kronosdb/upsert.sh', '-l', 'loc', string(format='(I7)', aaaajjjd), string(format='(I7)', aaaajjjf)], /NOSHELL END