POLAR MAPPING (October 22, 1999)

Polar Mapping

HESSI is a Fourier-transform imager and its primary data form consists of samples in the Fourier (u,v) plane of the Fourier transform of a flare. There is a rich heritage of Fourier transform imaging, most of it in the radio domain, but most Fourier imaging in astronomy has been done in rectangular, Cartesian coordinates, where the ``dirty map'' D(x,y) is given by the transform of visibilities:

$\begin{displaymath}D(x,y)= \int_{-\infty}^{\infty} du \int_{-\infty}^{\infty}dv\ V(u,v)\ e^{2\pi i (ux+vy)} \eqno{(1)} \end{displaymath}$

But HESSI is a spinning telescope, and the samples in the U,V plane are on circles $u^2+v^2=k^2={\rm constant}$ , one for each subcollimator and harmonic. Therefore it seems to make sense to go into polar coordinates, and it turns out that this is an efficient mapping strategy. (The method does not assume that the spin axis remains fixed.)

We must convert the Fourier transform from rectangular coordinates (x,y) to polar coordinates ( $r,\theta)$ in map space, and (u,v) to $(k,\phi)$ in visibility space:

$\begin{displaymath}D(r,\theta) = \int_{0}^{\infty} k\ dk \int_{0}^{2\pi} d\phi\ V(k,\phi)\ e^{2\pi i r k \ cos(\theta - \phi)} \eqno{(2)} \end{displaymath}$

But the integral over the radial wavenumber k becomes a sum over the number of subcollimators (and possibly their harmonics), because the visibility is sampled on circles:

$\begin{displaymath}D(r,\theta) = \sum_k k \int_{0}^{2\pi} d\phi\ V_k(\phi)\ e^{2\pi i r k \ cos(\theta - \phi)} \eqno{(3)} \end{displaymath}$

The above integral is a convolution, and therefore intrinsically faster to compute and more compact for computation than equation (1). This in a nutshell explains why Polar Mapping for HESSI is more efficient than rectangular mapping: the integrals are one-dimensional, and if the samples of visibility are equally spaced on the (u,v) circles (impossible in rectangular coordinates), FFT methods may be used in an efficient manner.

SIMULATION EXAMPLES

EXAMPLE 1

EXAMPLE 2

EXAMPLE 3

EXAMPLE 4

PROGRAMS (October 25, 1999)

(WARNING: Bugs exist in these programs, and they are likely to be out of date by the time you read this.)

polar_mapper.pro

;+ 
; NAME:  polar_mapper
;  
; PURPOSE:
;        MAKES A POLAR MAP FROM A CALIBRATED EVENT LIST
;              
; USAGE: 
;         1. Make a model map and call sim_calib_eventlist
;         2. Convert the counts to visibilities using calibevents2vis
;         4. Get the aspect structure for the same time bins
;         5. Shift the visibility phases to sun center
;         6. Make the polar map with vis2pmap 
;        
; INPUTS: 
;         det_index = scalar, 0-27
;         ce        = calibrated event list structure
;         xy_offset = [x,y] vector of inertial coords of map center (arcsec) 
;         /plot         Set this if you want plots and images
;        
; OUTPUTS: 
;         dmap = radius and azimuth map
;         dmap_xy = x,y projection (annulus) of dmap
;
; EXAMPLE:
;;   1. MAKE THE MODEL MAP
;       model=fltarr(64,64)
;       model(31,11)=1 
;
;;   2. CREATE THE CALIBRATED EVENT LIST for the selected det_index
;       det_index=5
;       xy_offset=[650,200]
;       sim_xy_offset=xy_offset+25
;       sim_photons=2000
;       det_mask=intarr(27) & det_mask(det_index)=1
;       bin_size=2L^12 ; microseconds  (this is ok except for det_index=0)
;       num_bins=(2L^22)/bin_size   
;       time_range=[0,4]  ; If this is changed, polar_mapper won't work
;       o = obj_new( 'hsi_calib_eventlist')
;       timebindef = o->get( /time_bin_def )  ; get the default array
;       timebindef(det_index)=bin_size
;       ce = o->getcalibeventlist(det_index, SIM_TIME_range=time_range, $
;             TIME_RANGE=time_range, SIM_PHOT=sim_photons, $
;             A2D_INDEX_MASK=det_mask,TIME_BIN_DEF=timebindef,SIM_MODEL=model, $
;             SIM_XYOFFSET=sim_xy_offset, XYOFFSET=xy_offset)   ;,/saszero 
;;   3. Call the mapper:
;       dmap=polar_mapper(det_index,ce,xy_offset,/plot)
;
;       
; RESTRICTIONS:
;       Works only for time_range=[0,4] and for a single subcollimator.
;       If the map position is several pitches away from the flare centroid,
;       the map will be poor.
;
; CALLS:
;                  hsi_grid_parameters
;                  calibevents2vis, (converts counts to visibilities)
;                     ls_sin_cos
;                     ls_amp_pha_ctr
;                  vis2pmap, which makes a polar map calling:
;                     convl (the key to speed)
;                     pmap2xy (interpolates the polar map onto xy coordinates)
;
; VERSION HISTORY:
;        VERS. 1
;          Revised to match compute_model_ce.pro Sept 17, 1999
;        VERS. 2              Sept. 27, 1999 -- EJS
;          Various bugs squashed
;          Forced visibility to be Hermitian
;        VERS. 3
;          Replaced Bessel-fcn array with FFTs
;          Added CONVL program for fast convolution (extends convolve.pro).
;          Made to work with an externally-produced calib_eventlist
;          Lots of bugs remain.                    Oct 26, 1999 -- EJS
;
; NOTES:
;        The grid_angle is used when correcting for aspect, but not in 
;           mapping, so the azimuth of the source comes out wrong.
;        No general times yet: We extract only the first 1/2 of 1
;           rotation from an entire rotation.
;
;        Maps are in the wrong hemisphere (S instead of N), probably
;           the result of negatively increasing roll angles.
;        There is something suspicious about the first dark ring in array
;           dmap_xy -- the scaling may be wrong
;        
;-

function polar_mapper,det_index,ce,xy_offset,plot=plot
;
;       
; RESTRICTIONS:
;       WORKS ONLY FOR TIMERANGE=[0,4]
;

       pitch=((hsi_grid_parameters()).pitch)(det_index)
       grid_angle=((hsi_grid_parameters()).orient)(det_index)

if keyword_set(plot) then plot,ce.dx,ce.dy,tit='ASPECT VALUES'

;  CONVERT THE COUNTS TO VISIBILITIES: 

   message,'Running calibevents2vis',/cont
   vis0=calibevents2vis(det_index,ce, cycle_factor=2)  ; best cycle_factor 


;  GET THE PHASE FACTOR FOR MAP PHASE CENTER

  alpha=2*!pi*((ce.dx)*cos(grid_angle)+(ce.dy)*sin(grid_angle))/pitch
         ; there may need to be a sign change


;  SHIFT THE VISIBILITY PHASES TO SUN CENTER:

   i=complex(0,1)
   vis1=vis0*exp(-i*(alpha)) ; This may be the wrong place to do this.
;                            ; Probably it should be done in calibevents2vis

;  HERMITIANIZE THE VISIBILITIES -- USE ONLY TIMERANGE=[0,2] SECONDS
   vis=vis1
   num_bins=n_elements(vis)
   vis(num_bins/2:num_bins-1)=conj(vis(0:num_bins/2-1))

if keyword_set(plot) then BEGIN
  !p.multi=[0,1,2]
  window,4,xsize=1000,ysize=400
  plot,ce.roll_angle,ce.count,xsty=1,tit='COUNT RATE',xtit='ROLL ANGLE (radians)
  plot,ce.roll_angle,vis0,tit='Visibility (vis0)',xsty=1,xtit='ROLL ANGLE (radians)
  !p.multi=0
  wset,0
endif

; THEN CREATE AN ANNULAR DIRTY MAP:
        mapwidth=256
        dmap_xy=1
        message,'Calling vis2pmap...',/continue
     if keyword_set(plot) then begin
          message,'Calling vis2pmap with /plot set',/cont
          dmap=vis2pmap(vis,ce,det_index,mapwidth ,dmap_xy,/plot) 
     endif else begin
          message,'Calling vis2pmap with /plot NOT set',/cont
          dmap=vis2pmap(vis,ce,det_index,mapwidth ,dmap_xy)
     endelse



return,dmap

end

vis2pmap.pro

function vis2pmap,vis,ce,det_index,map_width,dmap_xy,plot=plot
;+
; NAME: VIS2PMAP
;       
; PURPOSE:  Converts visibilities into a polar coordinate map
;        
;        
;        
; METHOD:
;      Compute the convolution of visibility and exp(i*k*r*cos(roll_angle))
;        
;
; INPUTS:
;       vis=complex array of N_vis values for phi=0 to 2*!pi
;       Note that the vis must be computed with phase center at the spin axis.
;       ce=calib_eventlist structure
;       det_index=0-27
;       map_width = scalar width (and height) of map to be made
;
; OUTPUTS:
;       dmap(j_r,n_theta)  = "dirty" map
;       Lots of plots if optional keyword plot is set
;       Optional return of rectangular map if dmap_xy exists
;
; EXAMPLE:
;       N=512
;       roll_angle=findgen(N)*2*!pi/N 
;       det_index=4
;       pitch=((hsi_grid_parameters()).pitch)(det_index)
;       i=complex(0,1)
;       vis=100*exp(i*2*!pi*605.*cos(roll_angle)/pitch)
;       count=100+float(vis)
;       phase_map_ctr=2*!pi*600.*cos(roll_angle)/pitch
;       ce={roll_angle:roll_angle,count:count,phase_map_ctr:phase_map_ctr}
;       map_width=200
;       dmap=vis2pmap(vis,ce,det_index,map_width,/plot)
;
;
; FUNCTIONS CALLED:
;      convl.pro
;
; RESTRICTIONS:
;       Makes only annular maps at present
;       The visibilities must have been  phase shifted to the polar center

;
; VERSION HISTORY:
;        VERS. 1 EJS Mar 29, 1999
;        VERS. 2 EJS Sept 28, 1999
;           Doubled width of annulus, increased final map to 360-deg coverage.
;        VERS 3 EJS Oct 15, 1999 
;           Replaced spin radius by a call to phase2radius
;           fixed an error in r2
;        VERS 4 EJS Oct 26, 1999
;           Changed Bessel function calls to convl calls.
; NOTES:
;-
; function vis2pmap,vis,ce,det_index,map_width,dmap_xy,plot=plot

message,'Beginning...',/continue

!p.multi=[0,1,2]
; 1. Look at the visibilities and compute the max spatial freq therein.
print,'Look at the visibilities and compute the max spatial freq therein.'
;

roll_angle=ce.roll_angle
phase_map_ctr=ce.phase_map_ctr
sigma=-1.
radius=phase2radius(ce,det_index,sigma)
message,'Deduced spin radius from phase_map_ctr: ='+string(radius),/cont
print,'sigma=',sigma
r0=fix(radius+0.5)          ; radius of the phase center
r1=(r0-map_width)>1 & r2=r1+2*map_width  ; min and max radii to be mapped
if (r1 LE 0) then begin
  message,'r1 LE 0'
endif


phi1=-!pi & phi2=!pi                 ; min and max azimuths to be mapped
if (keyword_set(plot)) then begin
   plot,float(vis),xsty=1,tit='Real Part of Visibilities'
   plot,imaginary(vis),xsty=1,tit='Imaginary Part of Visibilities'
endif
pitch=((hsi_grid_parameters()).pitch)(det_index)

if (max(abs(vis)) LT .0001) then message,'All visibilities are essentially 0'



Nroll=n_elements(roll_angle)
k=2*!pi/pitch              
r=r1+findgen(r2-r1)
kr=k*r
i=complex(0,1)
dmap=fltarr(Nroll,r2-r1)
message,'Doing'+string(r2-r1)+' radii...',/cont

for j=0,r2-r1-1 do begin
  v=exp(-i*kr[j]*cos(roll_angle))
        ; if ((j mod 10) eq 0) then plot,v,tit='r='+string(j+r1),xsty=1
        ; wait,0.1
  d=convl(v,vis)
  dmap[*,r[j]-r1]=d   ; 1st arg is azimuth, 2nd is radius
endfor


if (keyword_set(plot)) then begin
  erase
  ndx=n_elements(dmap(*,0))  & ndy=n_elements(dmap(0,*)) 
  window,1,xsize=ndx,ysize=ndy,title='Real Part of DMAP'
  print,'Pausing for 1 sec..' & wait,1
  erase
  
  tvscl,dmap

endif

; Modified Oct 26, 99
  message,'Calling pmap2xy.pro...' ,/continue 
;              dmap: 1st index is azimuth (0 to 2*!pi), 2nd is  radius

  dmap_xy=pmap2xy(dmap,r1,r2)  
  message,'r1,r0,r2='+string(r1)+string(r0)+string(r2),/continue

if (keyword_set(plot)) then begin
  if (!d.name NE 'PS') then window,0,xsize=960,ysize=960
  tvscl,dmap_xy
endif

message,'Done. Returning',/cont


return,dmap
end

pmap2xy.pro

function pmap2xy,polmap,r1,r2
; NAME pmap2xy
;
; PURPOSE:
;        Finds the x,y points inside an annulus of radii r1, r2, converts
;        to polar coordinates and uses interpolation to extract values from
;        a polar map polmap of size (
; INPUTS:
;   polmap:
;    1st index is azimuth (0 to 360 degrees) 
;    2nd index of polmap is radius (asec)
;   r1,r2 = integer values of inner and outer radii covered by pol
;
; OUTPUT:
;     a square map of size 960 x 960 with 2 asec pixels
;-


x=findgen(960)#replicate(1,960) 
y=transpose(x)
r=2*sqrt((x-479.5)^2+(y-479.5)^2)
n_az=n_elements(polmap(*,0))
i_radius=r-r1
j_azimuth=(n_az/2)*(1-atan((y-479.5),(x-479.5))/!pi)    
;j_azimuth=(n_az/2)*((2+atan((y-479.5),(x-479.5))/!pi) mod 2)
polmap=float(polmap)

d=interpolate(polmap,j_azimuth,i_radius,missing=0)

return,d


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_coscos, EJS, June 15, 1999.
;  Version 2, Added double and bad-fit flag, EJS Sept 30, 1999
;  NOTE: if the sine vector is nearly constant, a bad fit results.
;   This is because the sine is at an extremum and the cosine is 
;   very small, so it has to cover all the y variability, and its
;   coefficient will be large. (Similarly if cosine = const.)
;-

function ls_amp_pha_ctr,angl,y,phase_map_ctr,yfit,sigma, $
               AMP_PHA=amp_pha,verb=verb,double=double,flag=flag
; message,'stopping'
; Compute the sine and cosine arrays

if keyword_set(double) then begin
 if keyword_set(verb) then message,'Going to double precision',/continue
 phase_map_ctr=double(phase_map_ctr)
endif

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-2*(abs(m00)+abs(m11)) then begin
  message,'Ill-conditioned LS fit.',/continue
  ;message,string(del)+string(m00)+string(m01)+string(m11),/continue
  flag=1
endif
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  ; = amp*cos(phase_map_ctr-pha)
; compute the goodness of fit
if (n_elements(sigma) GT 0) then sigma=sqrt(variance(y-yfit))


if keyword_set(amp_pha) then AB=[sqrt(a^2+b^2),atan(B,A)]

return,AB

end

phase2radius.pro

function phase2radius,ce,det_index,rmsdev
;+
; NAME: phase2radius
;
; PURPOSE:
;       
; Finds the average spin radius from the pitch and the
; calib_eventlist tags .phase_map_center and .roll_angle
; for detector det_index
;
; INPUTS:
;         ce=calib_eventlist structure
;         det_index= scalar (0-27)
;
; OUTPUTS:
;         radius = scalar FLOAT
;
; VERSION HISTORY:
;   Upgraded inputs to structure ce
;                      EJS Oct. 26, 1999
;-

 pitch=((hsi_grid_parameters()).pitch)(det_index mod 9)
 if n_elements(rmsdev) NE 0 then rmsdev=-1.0

 yfit=1.
 z=ls_sin_cos(ce.roll_angle,ce.phase_map_ctr,1,yfit,rmsdev)
 radius=sqrt(z(0)^2+z(1)^2)*pitch/(2*!pi)

return,radius

end

convl.pro

FUNCTION CONVL,A,P
;+
;  PURPOSE:
;   IDL routine to convolve array A with a point spread function P
;
;  CALLING SEQUENCE:   C=CONVL(A,P)
;
;  INPUTS:
;    A = map vector or array (REAL OR COMPLEX)
;    P = PSF vector or array, same size as A
;
;  OUTPUT:
;    Convolved matrix or vector, shifted to make the origin the same
;    as A.  If P is 1 at array center, C=A
;
;  RESTRICTIONS:
;    A and P must be vectors or matrices of the same size, dimension LE 2
; 
;  NOTES:
;    The astron program "convolve" does not allow input arrays of
;    equal size, so this program is complementary to it.
;    With this function, 
;        total(C)=total(A)*total(P)
;        convl(A,P)=convl(P,A)
;        convl(A,P1+P2)=convl(A,P1)+convl(A,P2)
;
;  VERSION HISTORY:
;    VERS. 2.0 OCTOBER 25, 1999, EJS
;
;-

sz_a=size(a) & sz_p=size(p)
if (sz_a(0) NE sz_p(0)) then message,"Input arrays don't have same dimensions"
if (sz_a(0) eq 1) then $
   if (sz_a(1) NE sz_p(1)) then message,'Input vectors have different sizes'
if (sz_a(0) eq 2) then begin
   if ( (sz_a(1) NE sz_p(1)) OR (sz_a(2) NE sz_p(2)) ) then $
                message,'Input arrays have different sizes'
   n1=sz_a(1) & n2=sz_a(2)
endif
if (sz_a(0) gt 2) then message,'Arrays of dimension 3 or higher not supported'


c=float(fft(fft(a,1)*fft(p,1),-1) )

if sz_a(0) eq 1 then c=shift(c,sz_a(1)/2) $
else           c=shift(c,sz_a(1)/2,sz_a(2)/2)


return,c

end

Edward J. Schmahl

Last modified: Nov 8 09:55 EDT 1999