modprofile2vis
;+
; NAME: modprofile2vis
;
; PURPOSE:
; To convert HESSI modulation profiles for one half rotation
; into visibilities using calibrated phases at map center
;
; METHOD:
;
; a. Gets the ang_pitch from det_index using hessi_params()
; b. Finds the cycle period function 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:
; roll = array of roll angles spanning 180 degrees
; modprof = array of modulation amplitudes, 1 for each
; angle, same array size as roll
; det_index = 0-17 index of HESSI detector
; phase_map_ctr = array of map-center phases (same size as roll)
;
; OPTIONAL INPUTS:
; cycle_factor=cycle_factor (default=0.5) 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.
;
; OUTPUTS:
; visibilties = amp*exp(i*(pha - localmean(phase_map_ctr)))
;
;
; EXAMPLE:
; N=512 ;
; det_index=4
; roll_angle=findgen(N)*!pi/N - !pi/4 ; 180 deg of rotation
; params=hessi_params()
; 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:
; mod_profile=100.*cos(1.05*phase_map_ctr+!PI/6)+80*cos(phase_map_ctr)
;; Convert the modprof to visibilities:
; vis=modprofile2vis(roll_angle,mod_profile,det_index,phase_map_ctr)
; plot,imaginary(vis),xsty=1 & oplot,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)
;
; 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.
; 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
; Yet to be tested with simulated score data
; Should add optional user input of (rough?) map azimuth to determine
; sign of imaginary part of visibility.
;-
function modprofile2vis,roll,modprof,det_index,phase_map_ctr,$
verb=verb,plo=plo,cycle_factor=cycle_factor,conj_sign=conj_sign
params=hessi_params()
pitch=(params.pitch)((det_index mod 9))
grid_angle=(params.orient)((det_index mod 9))
;really should update .orient to .grid_angle
angle=roll-grid_angle ;
if keyword_set(verb) then message,'Computing amplitudes and phases...',/cont
nangles=n_elements(roll)
amp=roll*0 & A_vector=amp & B_vector=amp & pha=amp & sig=amp & sigma=0.
plot,angle/!dtor,phase_map_ctr,tit='Phase of Map Center'
; Find the parameters of the modulation frequency
amp_pha=ls_sin_cos(angle,phase_map_ctr,1.0,yfit,sigma,/amp_pha)
print,'amp_pha=',amp_pha
modulation_freq=amp_pha(0)*abs(sin(angle-amp_pha[1]))
oplot,angle/!dtor,modulation_freq
; Set up cycle periods:
rmin_cycle_period=2*!pi/amp_pha[0] ; (radians)
imin_cycle_period=rmin_cycle_period*180 ; (indices)
c=0.75
; The modulation freq vanishes at one angle, so we don't just invert
; it to get the period. Instead we use period=const/(const + frequency)
if (NOT keyword_set(cycle_factor)) then cycle_factor=0.5
di=fix(cycle_factor*imin_cycle_period/ $
(1.-c+c*modulation_freq/max(modulation_freq)))
if (min(di) LT 2) then message,'Warning--sample range LT 4',/cont
if (keyword_set(plo)) then plot,angle,modprof,xsty=1
i=complex(0,1)
if (keyword_set(conj_sign)) then if (conj_sign lt 0) then i=-i
i1=(findgen(nangles)-di)>0
i2=(findgen(nangles)+di)<(nangles-1)
dmodulo=2*fix(max(di)) ; could also be anyhing from max(di)+1 to nangles
; This array vv will contain the successive amp/phase fits;
; The NaN's are to flag elements with no data so proper averages
; can be made.
; A side effect is to cause spurious "Floating illegal operand"
; messages.
vv=complexarr(dmodulo,nangles)*!values.f_nan ; "compressed" array of NaN
for k=0,nangles-1 do begin
rolla=angle[k] ; when angle=0, rolla=grid_angle+phase_ctr_angle
; i1[k]=(k-di[k])>0
; i2[k]=(k+di[k])<(nangles-1)
y=modprof[i1[k]:i2[k]]
alpha=angle[i1[k]:i2[k]]
; Do a least square fit of A*cos(phase_map_ctr[])+B*sin(phase_map_ctr[]):
AB=ls_amp_pha_ctr(alpha,y,phase_map_ctr[i1[k]:i2[k]],yfit,sigma)
;if (keyword_set(plo)) then oplot,alpha,yfit,thick=3
A=AB[0] & B=AB[1]
amp=sqrt(A^2+B^2)
pha=atan(B,A)
; EXPLANATION:
; yfit = A*cos(phase_map_ctr[i1:i2])+B*sin(phase_map_ctr[i1:i2])
; = amp[k]*cos(pha - phase_map_ctr[i1:i2])
; The real part of the visibility = yfit:
; vis_real[k]=yfit(di) ; the central value of yfit (Don't use mean!)
; And the imaginary part of vis is found by replacing cos with sin:
; vis_imag[k]= amp[k]*sin(pha - phase_map_ctr[i1+di]) )
; vis[k]=amp[k]*(cos(pha-phase_map_ctr[i1:i2]) $
; +i * sin(pha-phase_map_ctr[i1:i2]))
; =amp[k]*exp(i*(pha - (phase_map_ctr[i1:i2]))--some kind of average
;
; Keep all the fits in a 2D array for subsequent row-averaging:
;vv[k,i1[k]:i2[k]]=yfit(*) +
;i*amp*sin(pha - phase_map_ctr[i1[k]:i2[k]])
; Accumulate all the fits in array vv to average after the loop.
; wrap around modulo dmodulo to save space:
vv((k mod dmodulo),i1[k]:i2[k])= $
yfit(*) + i*amp*sin(pha - phase_map_ctr[i1[k]:i2[k]])
sig[k]=sigma
endfor ;k angle loop
tvscl,vv
if (keyword_set(plo)) then begin
window,0 & tvscl,float(vv),0 & tvscl,imaginary(vv),2
endif
;vis=total(vv,1,/NaN)/total(finite(vv),1) ; old array (too big)
vis=total(vv,1,/NaN)/total(finite(vv),1) ; average unflagged row elements
print,'i=',i
return,vis
end ; end modprofile2vis
modprof2vis_demo
;+
; NAME: modprof2vis_demo.pro
;
; PURPOSE:
; To demonstrate capability of modprofile2vis
;
; METHOD:
; 1. Simulate a modulation profile for hardwired parameters
; det_index, r_map, phase_ctr_angle, offsets
;
; 2. Calls modprofile2vis (for 180 deg of rotation), which:
; a. Gets the ang_pitch from det_index
; b. Finds the shortest cycle period from map_ctr_phase
; c. Loops through angles, fitting sinusoids of phase to each cycle
; d. Uses amplitude and phase to create visibility
;
; 3. Plots the imaginary and real parts of all the visibilities
;
; INPUTS: None (Hardwired parameters go to modprofile2vis)
;
;
; OUTPUTS: Plots of simulated mod profile, real, imag parts of visibility
;
;
; EXAMPLE:
; compile modprofile2vis, ls_sin_cos, and ls_amp_pha_ctr, then
; .run modprof2vis_demo
;
;
; RESTRICTIONS:
;
;
; VERSION HISTORY:
; VERS. 1 Derived from test_lscoscos.pro -- EJS May 26, 1999.
; VERS. 2 Added /verb option, upgraded variable names--
; EJS Jun 14, 1999.
; NOTES:
;-
pro modprof2vis_demo
wset,0
;1. SIMULATE A MODULATION PROFILE FOR SPECIFIC PARAMETERS:
det_index=4
r_map=700.
phase_ctr_angle=!pi*0.75
sigma=1.
r_src_off1=0.0 ; simulated source positions in polar coords
r_src_off2=0.5
phi_src_off1=0.0
phi_src_off2=0.5
srcti='SIMULATED SOURCE OFFSETS: ' + $
'('+ trim(r_src_off1,'(f3.1)')+','+trim(phi_src_off1,'(f3.1)')+') '+ $
'('+ trim(r_src_off2,'(f3.1)')+','+trim(phi_src_off2,'(f3.1)')+')'
subti='R_map,phase_ctr_angle='+strcompress(r_map)+','+strcompress(phase_ctr_angle/!dtor)
sideti='EJS modprof2vis_demo.pro '+systime(0)
N=512
roll_angle=findgen(N)*!pi/N - !pi/4
params=hessi_params()
grid_angle=(params.orient)(det_index)
pitch=(params.pitch)(det_index)
; Create some aspect coords:
x_asp=15.00*smooth(rebin(randomn(seed,N/8),N,/samp),N/8,/edge)+ $
1.0*randomn(seed,N)
y_asp=15.00*smooth(rebin(randomn(seed,N/8),N,/samp),N/8,/edge)+ $
1.0*randomn(seed,N)
; Select an inertial map center:
X_map=600.
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)
; Choose a value for the peak_resp_offset (normally from hessi_grm):
peak_resp_offset=!pi/6
h=1 ; fundamental
; Make 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
angle=roll_angle-phase_ctr_angle-grid_angle ; =0 at pt of slowest modulation
print,phase_ctr_angle,grid_angle
;stop
modprofile=180.+100.*cos(1.05*phase_map_ctr+!PI/6)+80*cos(phase_map_ctr)
mod_profile=modprofile
; Add Poisson noise:
for k=0,n_elements(angle)-1 do $
mod_profile[k]=randomu(seed,1,poisson=modprofile[k])
mod_profile=mod_profile-mean(mod_profile)
xt='ROLL ANGLE (degrees)
plot,roll_angle/!dtor,mod_profile,xsty=1,tit='SIMULATED MODULATION PROFILE',xtit=xt
; 2. Convert the modprof to visibilities:
vis=modprofile2vis(roll_angle,mod_profile,det_index,phase_map_ctr,/plo, $
cycle_factor=0.125)
; 3. Plot the visibilities
if (!d.name eq "X") then window,2,xsize=750,ysize=500
!p.multi=[0,1,2]
ti='REAL PART OF VISIBILITY'
plot,roll_angle/!dtor,float(vis),tit=ti,xtit=xt,subtit=subti
xyouts,0.5,0.9,srcti,/norm,align=0.5
xyouts,0.99,0.5,sideti,/norm,align=0.5,orient=90
oplot,roll_angle/!dtor,mod_profile,psym=2,symsiz=0.3
ti='IMAG PART OF VISIBILITY'
plot,roll_angle/!dtor,imaginary(vis),tit=ti,xtit=xt,subtit=subti
xyouts,0.5,0.9,srcti,/norm,align=0.5
xyouts,0.99,0.5,sideti,/norm,align=0.5,orient=90
!p.multi=0
end
ls_amp_pha_ctr.pro
;+
; NAME: ls_amp_pha_ctr.pro
;
; PURPOSE: Returns the least-square fit of y(angle) to the function
; A*cos(phase_map_ctr(angle))+B*sin(phase_map_ctr(angle))
; or equivalently, AMP*cos(phase_map_ctr(angle)-PHA).
; The result can be used to fit visibilities to HESSI
; countrates, even in the ranges of slow modulation.
;
; METHOD:
; Theta=phase_map_ctr
; Minimize the total squared deviation H of the fit to
; y[i] = A*cos(Theta[i])+B*sin(Theta[i]))
;
; H=total{ [y - A*cos(Theta) - B*sin(Theta)]^2 }
;
; If the angle[i] are chosen to make the cos and sin terms uncorrelated
; then the regression matrix is diagonal. If, they are well correlated
; or nearly anti-correlated, the matrix is singular, and the method fails.
; The results may be meaningless unless r_p has been
; chosen such that phi[angle] ranges over a large fraction of [0,2*!pi].
;
; INPUTS:
; y=modulation profile (countrates), one for each angle
; angle=vector of angles (radians), possibly irregular
; phase_map_ctr=Array of phases of map center relative to line
; of maximum response, one for each roll angle
;
; OUTPUTS:
; [A,B] cos,sin amplitudes, OR
; [amplitude,phase] if optional amp_pha is set,
; If sigma exists as an argument, the RMS deviation of the fit
; is returned in sigma.
; If y has only 2 elements, sigma is necessarily 0 (exact fit).
; If yfit is a defined argument, fit to y will be returned in it.
;
; EXAMPLE:
; ;Simple simulation of map-center phases
; ;Create an array of angles:
; N=512
; angle=findgen(N)*!pi/N - !pi/4
; ;Create some aspect coords:
; x_asp=1.0*randomn(seed,N)
; y_asp=1.0*randomn(seed,N)
; ;Select an inertial map center:
; X_map=600.
; Y_map=200.
; ;Compute the spacecraft coords of map center:
; x_map_sc = x_asp + X_map*cos(angle) - Y_map*sin(angle)
; y_map_sc = y_asp + Y_map*cos(angle) + X_map*sin(angle)
; ;Choose a value for the peak_resp_offset (normally from hessi_grm):
; peak_resp_offset=!pi/6
; ;Pick a sub_collimator:
; ang_pitch=39.45 & det_index=4
; grid_angle=0. & h=1
; ;Now make 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)/ang_pitch
; ;Select a range of indices:
; inx=180+findgen(20)
; y=cos(phase_map_ctr)+0.25*randomn(seed,N)
; z=ls_amppha_pmc(angle(inx),y(inx),phase_map_ctr(inx),yfit,sig)
; plot,angle,y,psym=10,xsty=1 & oplot,angle(inx),yfit,thick=3
;
; RESTRICTIONS:
; If the data set is enormous, this is slow.
; For a way to speed it up,
; see Press and Rybicki, Ap.J., 338, 277, 1989.
;
; angle must =0 or !pi where the modulation is slowest, and
; angle must= +!pi/2 or -!pi/2 where modulation is fastest
;
; REVISION HISTORY:
; Version 1.0, Derived from ls_cosocs, EJS, June 15, 1999.
;
;-
function ls_amp_pha_ctr,angl,y,phase_map_ctr,yfit,sigma, $
AMP_PHA=amp_pha,verb=verb
; message,'stopping'
; Compute the sine and cosine arrays
cosine=cos(phase_map_ctr)
sine=sin(phase_map_ctr)
if (keyword_set(verb)) then message, $
'Working with '+String(n_elements(phi))+' angles',/cont
; The regression matrix
m00=total(cosine^2) & m01=total(sine*cosine)
m11= total(sine^2)
del=m00*m11-m01^2
if abs(del) lt 1.e-5*m00 then message,'Ill-conditioned LS fit.'
c0=total(y*cosine)
c1=total(y*sine)
A=(c0*m11-c1*m01)/del ; Solution to M # [A,B] = [c0,c1]
B=(c1*m00-c0*m01)/del
AB=[A,B]
yfit=A*cosine+B*sine
; compute the goodness of fit
if (n_elements(sigma) GT 0) then sigma=sqrt(variance(y-yfit))
if keyword_set(amp_pha) then AB=[abs(AB),atan(AB[1],AB[0])]
return,AB
end
ls_coscos.pro
; NAME: ls_coscos.pro
;
; PURPOSE: Returns the least-square fit of y(angle) to the function
; A*cos(2*!pi*r_p*cos(angle))+B*sin(2*!pi*r_p*cos(angle))
; or equivalently, AMP*cos(2*!pi*r_p*cos(angle)-PHA).
; The result can be used to fit visibilities to HESSI
; countrates, even in the ranges of slow modulation.
;
; METHOD:
; Let phi[angle] = 2*!pi*r_p*cos(angle[i])
; Minimize the total squared deviation H of the fit to
; y[i] = A*cos(phi[i)+B*sin(phi[i]))
;
; H=total{ [y - A*cos(phi) - B*sin(phi)]^2 }
;
; If the angle[i] are chosen to make the cos and sin terms uncorrelated
; then the regression matrix is diagonal. If, they are well correlated
; or nearly anti-correlated, the matrix is singular, and the method fails.
; The results may be meaningless unless r_p has been
; chosen such that phi[angle] ranges over a large fraction of [0,2*!pi].
;
; INPUTS:
; y=approximate sinusoid (vector)
; angle=vector of angles (radians), possibly irregular
; angle must =0 or !pi where the modulation is slowest, and
; angle must= +!pi/2 or -!pi/2 where modulation is fastest
; angle=roll_angle-grid_angle-phase_ctr_angle
; r_p = radial position of HESSI map center divided by pitch
;
; OUTPUTS:
; [A,B] cos,sin amplitudes, OR
; [amplitude,phase] if optional amp_pha is set,
; If sigma exists as an argument, the RMS deviation of the fit
; is returned in sigma.
; If y has only 2 elements, sigma is necessarily 0 (exact fit).
; If yfit is a defined argument, fit to y will be returned in it.
;
; EXAMPLE:
; Create an array of angles:
; angle=findgen(512)*!pi/512 - !pi/4
; Define variables:
; r_p=20. & sigma=0.
; Define function values (a simulated modulation profile):
; y=100.*cos(2*!pi*r_p*cos(angle)) + $ ; 80.*cos(2*!pi*(r_p-1.0)*cos(angle+0.01)) & plot,angle,y
; AB=LS_coscos(angle(120:160),y(120:160),r_p,yfit,sigma)
; plot,angle,y,xsty=1 & oplot,angle(120:160),yfit,thick=3
; OR:
; AB=LS_coscos(angle(330:340),y(330:340),r_p,yfit,sigma)
; plot,angle,y,xsty=1 & oplot,angle(330:340),yfit,thick=3
; OR:
; AB=LS_coscos(angle(240:250),y(240:250),r_p,yfit,sigma)
; plot,angle,y,xsty=1 & oplot,angle(240:250),yfit,thick=3
;
; RESTRICTIONS:
; If the data set is enormous, this is slow.
; For a way to speed it up,
; see Press and Rybicki, Ap.J., 338, 277, 1989.
;
; angle must =0 or !pi where the modulation is slowest, and
; angle must= +!pi/2 or -!pi/2 where modulation is fastest
;
; REVISION HISTORY:
; Version 1.0, EJS Derived from LS_sinusoid, May 25, 1999.
;
;-
Ed Schmahl
Last modified: Tue Jul 20 11:10:03 EDT 1999