make_localization_error_3a.pro (Kronos IDL library)

;+
; Contains the make_localization_error_3a procedure
;
; :Author:
;    Baptiste Cecconi, Laurent Lamy
;
; :History:
;    2009/02/20: Created
;
;    2010/05/17: Last Edit
;-
;
;+
; make_localization_error_3a is a procedure that <behavior desc here>
;
; :Uses:
;    add_error_on_gp, iso_fc_intersect, load_data_ephem, load_data_n3bi, load_data_n3edi, magnetic_field, mfl_foot, q_prod, xyz_to_rtp
;
; :Params:
;    aaaajjjd:    in, required, type=sometype
;       A parameter named aaaajjjd
;    hd:          in, required, type=sometype
;       A parameter named hd
;    aaaajjjf:    in, required, type=sometype
;       A parameter named aaaajjjf
;    hf:          in, required, type=sometype
;       A parameter named hf
;
; :Keywords:
;    n_error:     in, optional, type=sometype
;       A keyword named n_error
;    mfl_model:   in, optional, type=sometype
;       A keyword named mfl_model
;    source:      in, optional, type=sometype
;       A keyword named source
;    output_path: in, optional, type=sometype
;       A keyword named output_path
;    verbose:     in, optional, type=sometype
;       A keyword named verbose
;    nocrt:       in, optional, type=sometype
;       A keyword named nocrt
;-
PRO MAKE_LOCALIZATION_ERROR_3A,aaaajjjd,hd,aaaajjjf,hf,n_error=n_error,mfl_model=mfl_model,$
                               source=source,output_path=output_path,verbose=verbose,nocrt=nocrt

if ~keyword_set(n_error) then n_error=8
if ~keyword_set(mfl_model) then mfl_model='SPV'
if ~keyword_set(source) then source='sq'

aj1=aj_t97(aaaajjjd) 
aj2=aj_t97(aaaajjjf) 

	for jour=aj1,aj2 do begin

aj=t97_aj(jour)
h1=0. & h2=23.

if jour eq aj1 then h1=hd<23
if jour eq aj2 then h2=hf<24

		for h=h1,h2-1l do begin

load_data_n3bi,aj, h, aj, h+1, n1,n2,n3b,ind,source=source,data_path=data_path,/verb
;load_data_n3edi,aj,h,aj,h+1,n1,n2,n3d,n3e,ind,source=source,data_path=data_path,/verb
if n2(0).t97 eq 0 then goto,suite
load_data_ephem,aj,h,aj,h+1,veph,qeph,veph='vsat',qeph='qssq',t97=n2.t97

;if n_elements(n2) gt 1 then begin

	read_data_src_list,src_list ;,file='/c/rpws/pro/src_list.txt'
	src_data = (src_list(where(src_list.df_name eq source)))(0)

; ua = 1 Unite astronomique en Rp

ua = 150.e6/src_data.body_rad_value


; Preparation des Ephemerides:
;-----------------------------

; qcenter: quaternion containing the body position (with respect to S/C) in S/C frame
qcenter = q_vmake(veph.r)
; qc0    : quaternion containing the body position (with respect to S/C) in body frame
qc0     = q_rot(qeph.q,qcenter)
; rr0    : distance S/C to body (km)
rr0      = sqrt(total(qc0(1:3,*)^2.,1))
; qr1    : rotational quaternion to put the body direction in (XoZ) plane
;          i.e. rotation around Z of angle - azimuth(body)
qr1     = q_make(reform(-atan(qc0(2,*),qc0(1,*))),[0,0,1])
; qr2    : rotational quaternion to put the body direction in the X direction
;          i.e. rotation around Y of angle pi/2 - colatitude
qr2     = q_make(reform(!pi/2-acos(qc0(3,*)/rr0)),[0,1,0])
; qrot0  : rotational quaternion that merges the qr2, qr1 and qeph.q
;          in the resulting frame:
;             - Z is pointing into North
;             - X is pointing to the body center
;
; ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
;  NB:  in the new frame X pointing to the body center, Z is up along the North direction, 
;      and Y toward left.
; ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
;
qrot1   = q_prod(qr2,qr1) 
qrot0   = q_prod(qrot1,qeph.q)
; qcenter2: quaternion containing the body position (with respect to S/C) in the new frame
qcenter2= q_rot(qrot0,qcenter)
; veph0  : just to keep a link to the original ephemeris data
veph0   = veph
; veph   : updating the body position in the new frame (veph.r should be rr*[1,0,0])
veph.r  = qcenter2(1:3,*)

; Selection des donnees correspondant au jour courant:
;-----------------------------------------------------

; Interpolation des Ephemerides:
;-------------------------------

qrot = qrot0
dth = 2./!radeg
add_error_on_gp,dth,ind.rr,ind.ydf,ind.zdf,qrot,n_error,qdferr,axis='N'


; localisation : 
; --------------
read_data_src_list,src_list
src_data = src_list(where(src_list.df_name eq source))
rp = src_data.BODY_RAD_VALUE

qv_SC = q_vmake(-veph0.r)
qtmp  = q_rot(qeph.q,qv_SC)
v_SC  = qtmp(1:3,*)/rp ;(in Rsat)
l_SC  = sqrt(total(v_SC^2.,1))

ndata = n_elements(n2)

pn2 = ptr_new(n2,/allocate)
n2  = replicate({data_n2},n_error,ndata)
for i=0,n_error-1 do n2(i,*) = reform(*pn2,1,ndata)
n2 = reform(n2,n_error*ndata)
ptr_free,pn2

v_SC = reform(rebin(reform(v_SC,3,1,ndata),3,n_error,ndata),3,n_error*ndata)
l_SC = reform(rebin(reform(l_SC,1,ndata),n_error,ndata),n_error*ndata)

qtmp = fltarr(4,n_error,ndata)
for i=0,n_error-1 do qtmp(*,i,*) = reform(q_rot(qeph.q,reform(qdferr(*,i,*),4,ndata)),4,1,ndata)
qtmp = reform(qtmp,4,n_error*ndata)

v_DF = [qtmp(1:3,*)]
a_SP = angular_distance(v_DF,-v_SC)

ndata = n_elements(n2)

xyz_fc  = fltarr(3,ndata)
foot_fc = fltarr(3,ndata)
beam_fc = fltarr(ndata)
dist_fc = fltarr(ndata)
verbose = 0


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)

if cnt_in ne 0 then begin

	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)
	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,nocrt=nocrt)
	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

	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

	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,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_err = replicate({data_loc},ndata)
data_loc_err.ydh  = n2.ydh
data_loc_err.num  = n2.num
data_loc_err.foot = foot_fc
data_loc_err.xyz  = xyz_fc
data_loc_err.beam = beam_fc
data_loc_err.dist = dist_fc

if keyword_set(output_path) then dataFileOut=output_path else dataFileOut=data_path+'loc/'
if strmid(output_path,0,1,/reverse) ne '/' then output_path=output_path+'/'
dataFileOut = dataFileOut+'loc_err_3A_'+mfl_model+'_'+(keyword_set(nocrt)?'nocrt_':'')+string(format='("n",i2.2,"_",i7,".",i2.2)',n_error,aj,h)

write_data_binary,dataFileOut,data_loc_err

suite:
		endfor

	endfor

; 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