;+ ; function hsi_calib_ev2vis, det_index,calib_event_list,$ ; VERB=verb,PLOT=plot,CYCLE_FACTOR=cycle_factor,CONJ_SIGN=conj_sign, $ ; FILTER=lopass_filter,$ ; WEIGHTS=weights,$ ; RUNNING_MEAN=running_mean ; ; NAME: hsi_calib_ev2vis ; ; PURPOSE: ; To convert HESSI count profiles into visibilities using ; calibrated phases at map center ; ; METHOD: ; ; a. Gets the ang_pitch from det_index using hessi_grid_parameters() ; b. Finds the cycle periods from phase_map_ctr(roll_angle) ; ; Explanation: ; If aspect offsets were all zero, the phase_map_ctr would be ; equal to 2*!pi*r_map*cos(roll_angle + const), and on a ; range of roll angles spanning 180 deg, the ; phase_map_ctr function is purely sinusoidal with ; one absolute maximum and a minimum magnitude of 0. ; In general, the slowest modulation occurs when ; abs(phase_map_ctr) = 0, and fastest modulation occurs ; at roll angles 90 degrees earlier or later. ; But when aspect offsets are non zero, phase_map_ctr is ; only approximately sinusoidal, so a fit to the phase values ; is made in the form A*sin(roll_angle)+B*cos(roll_angle), ; and this gets us the cycle lengths. ; ; c. Loops through angles, fitting sinusoids of phase_map_ctr to the ; modulation profile values on each modulation cycle ; d. Uses amplitude and phase to create visibility, taking ; appropriate means of the least-square parameters. ; ; INPUTS: ; det_index = 0-27 index of HESSI detector ; calib_eventlist ; whose tags contain: ; calib_eventlist.roll_angle=array of roll angles spanning 360 ; calib_eventlist.count= array of counts, 1 for each ; angle, same array size as roll ; calib_eventlist.phase_map_ctr = array of map-center phases ; ; OPTIONAL INPUTS: ; CYCLE_FACTOR=cycle_factor (default=1.0) Small values will use ; fewer points in a cycle to do the least-square sinusoidal fit. ; CONJ_SIGN=conj_sign (default=1) If negative, the sign of ; sqrt(-1) will be taken to be negative, corresponding to ; mirror source. ; RUNNING_MEAN If set, vis will be constructed from a running ; mean along the profile, if not, no running mean. ; OUTPUTS: ; visibilties = amp*exp(i*(pha - localmean(phase_map_ctr))) ; ; ; EXAMPLE: ; N=512 ; ; det_index=4 ; roll_angle=findgen(N)*2*!pi/N - !pi/4 ; 360 deg of rotation ; params=hsi_grid_parameters() ; grid_angle=(params.orient)(det_index) ; pitch=(params.pitch)(det_index) ;; Simulate some aspect coords (can be anything): ; x_asp=2.*randomn(seed,N)-100.*cos(roll_angle) ; ; y_asp=2.*randomn(seed,N)-100.*sin(roll_angle) ; ;; Select an inertial map center: ; X_map=600. ; arc sec ; Y_map=200. ;; Compute the spacecraft coords of map center: ; x_map_sc = x_asp + X_map*cos(roll_angle) - Y_map*sin(roll_angle) ; y_map_sc = y_asp + Y_map*cos(roll_angle) + X_map*sin(roll_angle) ; r_map_sc = sqrt(x_map_sc^2+y_map_sc^2) ; ang_sc=atan(y_map_sc,x_map_sc) ; azimuth of map center ; r_map=mean(r_map_sc) ; approx dist from spin axis ;; Choose a value for the peak_resp_offset (normally from hessi_grm): ; peak_resp_offset=!pi/6 ; h=1 ; fundamental ; Harmonics 2,3,... not tested ;; Simulate the array of map-center phases: ; phase_map_ctr=2*!pi*h*(x_map_sc*cos(grid_angle) $ ; +y_map_sc*sin(grid_angle)+peak_resp_offset)/pitch ; phase_map_ctr=2*!pi*h*(r_map_sc*cos(ang_sc-grid_angle) $ ; +peak_resp_offset)/pitch ;; Simulate a modulation profile for a pair of point sources: ; count=100.*cos(1.05*phase_map_ctr+!PI/6)+80*cos(phase_map_ctr)+180 ;; Make a calibrated event list structure: ; ce={dx:x_asp,dy:y_asp,count:count,roll_angle:roll_angle, $ ; phase_map_ctr:phase_map_ctr} ; ;; Convert the calib_event_list to visibilities: ; vis=hsi_calib_ev2vis(det_index,ce) ; plot,vis,xsty=1 & oplot,imaginary(vis),line=1 ; ;; To test the fit, compare vis to the true visibilities: ; i=complex(0,1) ; vis_true=100.*exp(i*1.05*phase_map_ctr+i*!PI/6)+80*exp(i*phase_map_ctr) ; CALLS: ; hsi_grid_parameters ; hsi_ls_sin_cos ; hsi_ls_amp_pha_ctr ; ; RESTRICTIONS: ; The sign of imaginary(visibility) is indeterminate, but can be selected ; ; VERSION HISTORY: ; VERS. 1 Derived from test_lscoscos.pro -- EJS May 26, 1999. ; VERS. 2 Added /verb option, upgraded variable names ; -- EJS Jun 15, 1999. ; VERS. 3 Corrected several bugs -- EJS June 18, 1999 ; VERS. 4 Eliminated r_map, uses phase_map_ctr only ; -- EJS June 21, 1999 ; 4.1 Reduced storage space requirements immensely by ; wrapping intermediate least-sq fits into a less sparse array. ; Generalized averaging interval by adding optional cycle_factor ; Added option to reverse sign of imaginary part of vis. ; -- EJS June 25, 1999. ; 4.2 Forced the sine and cosine parts of the visibility ; to be derived from the same exponential. ; Subtract out mean of count profile ; -- EJS Oct 15, 1999 ; 4.3 Incorporated a shift of counts to force correct ; azimuth of simulated source. -- EJS Oct 25, 1999 ; VERS 5.0 Changed windowing scheme to wrap around the ends of ; the arrays, assuming they are periodic, ; Revised to include aspect phase shifts BEFORE ; sinusoids are fit. Changed name from calibevevents2vis ; EJS Nov 4, 1999. ; 5.1 Corrected bug in which roll_angle shift was applied ; before the phase shift in polar_mapper. Now the ; roll_angle shift is done there. EJS Nov 23, 1999 ; 5.2 Added keyword WEIGHTS to return the COEFFICIENTS ; multiplying the count profile in each window. ; ; NOTES: ; July 20, 1999: Found that the sign of the imaginary part of (vis) ; is opposite to what was expected. This is not a problem since the ; sign is arbitrary, but an explanation is needed. ; ; Harmonics h=2,3,4,... need to be tested ; tested with simulated score data ; Should add optional user input of (rough?) map azimuth to determine ; sign of imaginary part of visibility. ;- function hsi_calib_ev2vis, det_index,calib_event_list,$ VERB=verb,$ PLOT=plot,$ CYCLE_FACTOR=cycle_factor,$ CONJ_SIGN=conj_sign, $ FILTER=lopass_filter, $ WEIGHTS=weights,$ RUNNING_MEAN=running_mean modamp=calib_event_list.modamp gridtran=calib_event_list.gridtran livetime=calib_event_list.livetime count=calib_event_list.count ccount=count/(gridtran*livetime) ; Check this out! if keyword_set(lopass_filter) then begin ccount=hsi_lo_pass_filter(ccount) message,'Applying lo-pass filter to count profile',/inform endif roll=calib_event_list.roll_angle nangles=n_elements(roll) phase_map_ctr=calib_event_list.phase_map_ctr pitch=((hsi_grid_parameters()).pitch)((det_index mod 9)) grid_angle=((hsi_grid_parameters()).orient)((det_index mod 9)) if keyword_set(verb) then message,'Computing amplitudes and phases...',/inform if keyword_set(running_mean) AND keyword_set(weights) then $ message,'WARNING: Weights will be only approximate with running_mean visibilities',/inform amp=roll*0 & A_vector=amp & B_vector=amp & pha=amp & sig=amp zcount=ccount-mean(ccount) ; zero-mean ; plot,roll,count,tit='Count rate',xsty=1 ; plot,roll,smooth(zcount,11,/edge),tit='Smoothed Count rate',xsty=1 imin_cycle_period=hsi_timebin_evaluator(calib_event_list) message,'imin_cycle_period='+string(imin_cycle_period),/info if keyword_set(cycle_factor) then begin message,'Cycle factor='+string(cycle_factor),/info if imin_cycle_period*cycle_factor LT 2.0 then $ message,'WARNING: With current cycle_factor, you will undersample cycles!',/info endif else begin cycle_factor=4.0/float(imin_cycle_period) message,'Setting cycle_factor='+string( cycle_factor),/info endelse ; Find the cycle periods: di=hsi_cycle_bins(roll,phase_map_ctr,cycle_factor) if (keyword_set(plot)) then plot,roll,zcount,xsty=1,tit='MODULATION PROFILE', $ xtit='ROLL ANGLE' i=complex(0,1) if (keyword_set(conj_sign)) then if (conj_sign lt 0) then i=-i dmodulo=3*fix(max(di)) ; could also be anything from max(di)+1 to nangles, ; but the bigger it is, the more wasted storage. if keyword_set(running_mean) then $ vv=complexarr(nangles,nangles)*!values.f_nan ; array of NaN ;weight_arrayRE=fltarr(nangles,nangles) ; arrays to store weights ;weight_arrayIM=fltarr(nangles,nangles) ; for later statistics A_wt=fltarr(nangles,nangles) B_wt=fltarr(nangles,nangles) t_0=systime(1) AiB=complexarr(nangles) AiB_wt=complexarr(nangles,nangles) weights=complexarr(nangles,nangles) ; Loop thru angles fitting sinusoids within window di: for k=0,nangles-1 do begin ; Create array of window indices, wrapping around ends if necessary w=(k+findgen(2*di[k]+1)-di[k]+nangles) mod nangles ph_map_ctr=phase_map_ctr(w) badfit=0 wts=[[zcount(w)],[zcount(w)]]; where ls_amp_pha_ctr stores the AB weights sigma=1. AB=hsi_ls_amp_pha_ctr(roll(w),zcount(w),ph_map_ctr,yfit,sigma,flag=badfit, WTS=wts) ; AB[s]=total(wts[*,s]*zcount(w)), s=0,1 AiB[k]=(AB[0]-i*AB[1])/modamp[k] ; This must be checked! ; plot,AiB & oplot,AiB_wt#zcount,thick=2 ; these should be identical if (badfit eq 1) then message,'bad fit at k='+string(k),/continue if keyword_set(running_mean) then $ vv[k,w]=(AB[0]-i*AB[1])*exp(i*ph_map_ctr) ; one column of vv array AiB_wt[k,w]=wts[*,0]-i*wts[*,1] weights[k,w]=AiB_wt[k,w]*exp(i*phase_map_ctr[k]) endfor ;k angle loop ; Create the visibility (unless ROLLING_MEAN is set) vis0=AiB*exp(i*phase_map_ctr) ; from single-point in window ; plot,vis0 & oplot,ab_wt#zcount,thick=2 ; these should be almost identical ;message,'Done with vv, doing weights...',/info ;d=fltarr(nangles,nangles) ;d(where(finite(vv)))=1. ;e=d ;for j=0,nangles-1 do e[*,j]=d[*,j]/total(d[*,j]) ;message,'Done with e, doing matrix multiplation...' ;W=e#( A_wt-i*B_wt)--impractical--a million dot products of 1024-elt vectors ;message,'Done with first part of W, doing 2nd part...',/info ;for n=0,nangles-1 do W[*,n]= exp(i*phase_map_ctr)*W[*,n] ;w[j,n] = exp(i*phi[j]) * sum_k { e[j,k]* (Awt[k,n] -i* Bwt[k,n]) } if keyword_set(RUNNING_MEAN) then begin tvscl,vv print,'Time to construct vv=',systime(1)-t_0,' seconds.' if (keyword_set(plot)) then begin window,/free,title="hsi_calib_ev2vis.pro" tvscl,float(vv),0 & tvscl,imaginary(vv),2 endif vis=total(vv,1,/NaN)/total(finite(vv),1) ; average unflagged row elements sigma=fltarr(nangles) for k=0,nangles-1 do begin sig=moment(abs(vv(*,k)-vis[k]),/NaN,sdev=sigma_k) sigma[k]=sigma_k endfor endif else begin vis=vis0 endelse !x.style=1 & !p.multi=0 ; wtIM_sparse=sprsin(weight_arrayIM) !p.multi=0 & !x.style=1 window,0,xsize=1200,ysize=256 plot,smooth(zcount(*),3,/edge),tit='Smoothed zcount and vis',psym=10 if keyword_set(running_mean) then oplot,vis,thick=2 else plot,vis0,thick=2 print,'3-pt smooth corr=',correlate(smooth(zcount(*),3,/edge),float(vis)) ;stop wlarge=where(abs(vis) gt 10*mean(abs(vis)),nlarge) if nlarge gt 0 then begin plot,abs(vis),tit='ABS(vis)' message,strcompress(nlarge)+' visibilities exceeding 10*mean',/continue message,'nlarge='+string(nlarge) endif z=check_math() ; clear the NaN error flag return,vis end ; A[k] = sum_n { Awt[k,n] * Count[n] ; B[k] = sum_n { Bwt[k,n] * Count[n] ; ; VV[k,j] = (A[k]-i*B[k]) * exp(i*phi[j]) ; = exp(i*phi[j]) * sum_n { ( Awt[k,n] -i* Bwt[k,n] ) * Count[n] } ; = sum_n {WW[j,k,n] * Count[n] } ; ; WW[k,j,n] = exp(i*phi[j])*( Awt[k,n] -i* Bwt[k,n] ) ; ; V[j] = sum_k { e[k,j] * vv[k,j] } ; = sum_k { sum_n {WW[k,j,n]* Count[n]} } ; = sum_n { sum_k {d[k,j] *WW[k,j,n] } * Count[n] } ; = sum_n { w[j,n] * Count[n] } ; ; d[k,j] = 1 if vv[k,j] NE NaN else d[k,j] = 0 ; e[k,j]= d[k,j]/sum_k d[k,j] ; ; w[j,n] = exp(i*phi[j]) * sum_k { e[j,k]* (Awt[k,n] -i* Bwt[k,n]) } ; ; ; ; ; ; ; ; ; ; ; ; ; ; ; ;