pro hsi_bproj2size,obs_time_interval,map_interval,xyoffset,energy_band, $ SAVE=save,$ PLOT=plot,$ VERB=verb,$ FILE=file,$ SINGLE_LEVEL=single_level,$ CROSSCORR_MIN=crosscorr_min,$ MAX_AXES_RATIO=max_axes_ratio,$ CHK_XYOFFSET=chk_xyoffset,$ CHK_CIRCULARITY=chk_circularity,$ CHK_SINGLE=chk_single,$ CHK_OFFSET=chk_offset ;+ ; PURPOSE: ; gets bproj maps for one flare, 9 dets, one energy_band ; (restricted to 64 x 64 maps with 1 arcsec pixels), and creates ; relative amplitude profile, whose transform is the spatial profile. ; ; METHOD: ; 1. Make 9 (8 if E<20 keV) bproj maps for the given inputs ; 2a Find "best" location of maximum (pixel m) ; re-do 1 if necessary. ; 2b Make PSF for SC #1,2, & 3 --> kmax_mask ; 3. Convolve bproj maps with PSF. Find crosscorr max to get k_max ; 4. Make clean maps using dets 3-9 ; 5. Check that source is single and compact using clean maps, and ; that source is near center. ; 6. Get calib eventlists for all detectors ; 7. Determine spinaxis - map center distance, --> kmin_mask(*) ; 8. Kmask=kmin_mask*Kmax_mask. Set Kmax(1)=1 iff energy_band[1] > 20. ; 9. Compute , total_counts, maximum and get relamp(k) ;10. Gaussian interpolate relamp on (log(R),k^2) space ;11. Bessel transform the interpolated relmps to get size profiles ; ; INPUTS: ; time interval in the form ['2002/03/25 20:05:00','2002/03/25 20:15:00'] ; map_interval in the form ['2002/03/25 20:06:30','2002/03/25 20:07:30'] ; energy_band = FLOAT(2) ; e.g. [12.,25.] ; xyoffset = FLOAT(2) ; e.g. [-282., 247.] ; SINGLE_LEVEL = acceptable contour level (percent) for second source ; CROSSCORR_MIM = minimum acceptable correlation with PSF ; (default=0.5) ; CHK_XYOFFSET If set, checks map max, if not centered, re-runs bproj (default) ; CHK_CIRCULARITY If set, checks circularity ; ; OUTPUTS: ; maparr=FLOAT(64,64,9) ; if SAVE is set, a save array containing the parameter structure and ; the maps is saved with a name based on map_interval{0} & energy ; if PLOT is set, various plots will be made ; if VERB is set, useful parameters will be printed ; if FILE is set equal to a filename, params will be printed to it. ; REQUIRED MODULES: ; crosscorr.pro ; hsi_find_spinaxis.pro ; bessel_trans.pro ; check_multiplicity ; check_circularity.pro ; ; CALLING EXAMPLE: ; energy_band = [12.,25.] ; obs_time_interval=['2002/02/26 10:23:56','2002/02/26 10:33:44'] ; map_interval=['2002/02/26 10:26:00','2002/02/26 10:27:00'] ; xyoffset = [929.,-226.] ; hsi_bproj2size,obs_time_interval,map_interval,xyoffset,energy_band,/save,/plot,/verb,file='Feb26.dat' ; VERSION HISTORY: ; ejs january 27, 2003 (schmahl@hessi.gsfc.nasa.gov) ; EJS April 15, 2003 ; ejs Mar 2004 Updated ;- if keyword_set(FILE) then message,'file='+file,/info ERRMSG='' maparr=fltarr(64,64,9) cpsf=maparr ndets=9 if keyword_set(FILE) then begin close,1 openw,1,file printf,1,'TIME INTERVAL:',obs_time_interval printf,1,'MAP INTERVAL:',map_interval printf,1,'XYOFFSET=',xyoffset printf,1,'ENERGY_BAND=',ENERGY_band endif default,CROSSCORR_MIN,0.25 default,SINGLE_LEVEL,50 default,MAX_AXES_RATIO,0.50 default,CHK_OFFSET,1 default,CHK_CIRCULARITY,1 default,CHK_SINGLE,1 message,'CROSSCORR_MIN='+strcompress(crosscorr_min),/info message,'SINGLE_LEVEL='+strcompress(single_level),/info obj_im = hsi_image() obj_im -> set,use_auto_time_bin=0 obj_im -> set,obs_time_interval=obs_time_interval obj_im -> set,energy_band=energy_band ; obj_im -> set,image_algorithm='back_proj' obj_im -> set,pixel_size=[1,1] obj_im -> set,time_bin_def=[1,1,1,1,1,1,1,1,1] obj_im -> set,image_dim=[64,64] obj_im -> set,image_algorithm='back_proj' istart=anytim(map_interval[0])-anytim(obs_time_interval[0]) iend=anytim(map_interval[1])-anytim(obs_time_interval[0]) if iend-istart GT 120 then message,'Map interval exceeds 2 minutes!',/info obsinterval=anytim(obs_time_interval[1])-anytim(obs_time_interval[0]) if istart LT 0 or iend-istart GT obsinterval then ERRMSG='Map interval outside obs_time_interval' if ERRMSG NE '' then message,errmsg obj_im -> set,time_range=[istart,iend] obj_im -> set,det_index_mask=bytarr(9)+1 p=obj_im->get() ; this gets us the parameter settings bproj_iter=0 message,'STEP 1: Make 9 single-SC bproj maps, PSFs, and cbes',/info STEP1: obj_im -> set,XYOFFSET=xyoffset ; This may be reset below. message,'new xyoffset='+strcompress(xyoffset[0])+strcompress(xyoffset[1]),/info if keyword_set(FILE) then printf,1,'new xyoffset=',xyoffset[0],xyoffset[1] for j=0,8 do begin det_index_mask=bytarr(9) det_index_mask[j]=1b im=obj_im->getdata(det_index_mask=det_index_mask) maparr(*,*,j)=im(0:63,0:63) if n_elements(im(*,0)) NE 64 then message,'WARNING: im dimensions not 64x64',/info if n_elements(im(0,*)) NE 64 then message,'WARNING: im dimensions not 64x64',/info endfor if keyword_set(chk_offset) then begin message,'STEP 2a Check for maxima of maps--must be close to map center or phases will be bad',/info STEP2a: ; Check for maxima of maps--must be close to center or phases will be bad xmax=fltarr(9) & ymax=xmax For j=2,7 do begin xymax=locate_val(maparr(*,*,j),max(maparr(*,*,j))) xmax[j]=xymax[0] & ymax[j]=xymax[1] if keyword_set(VERB) then message,'SC,xy max:'+strjoin([j+1,xmax[j],ymax[j]]),/info if keyword_set(FILE) then printf,1,'SC,xy max:',j+1,xmax[j],ymax[j] endfor;j xbar=median(xmax[3:6]) & ybar=median(ymax[2:6]) dim=obj_im->get(/image_dim) if (xbar-(dim[0]+1.)/2)^2 + (ybar-(dim[1]+1.)/2)^2 GT 1.0 then begin xy_offset=xyoffset+[xbar-(dim[0]+1.)/2,ybar-(dim[1]+1.)/2] message,'Median x,y positions:'+strcompress(xbar)+strcompress(ybar),/info if keyword_set(FILE) then printf,1,'Mean x,y positions:',xbar,ybar xy_off=strcompress(xy_offset[0])+strcompress(xy_offset[1]) xyoffset_ = strcompress(xyoffset[0])+strcompress(xyoffset[1]) MSG='The present xyoffset='+xyoffset_+ ' We will run again with xyoffset='+xy_off if keyword_set(VERB) then message,MSG,/info if keyword_set(FILE) then printf,1,MSG xyoffset=xy_offset bproj_iter=bproj_iter+1 if bproj_iter EQ 1 then begin goto,step1 endif else message,'ABORTING...initial xyoffset looks bad. endif else begin message,'Median x,y positions:'+strcompress(xbar)+strcompress(ybar),/info if keyword_set(VERB) then $ message,'The maximum is within 1 pixel of the map center',/info if keyword_set(FILE) then $ printf,1,'The maximum is within 1 pixel of the map center' endelse endif ; chk_offset message,'STEP 2b j='+strcompress(j),/info STEP2b: dim=obj_im->get(/image_dim) obj_cl = hsi_image(image_algorithm='clean' ) for j=0,8 do begin det_index_mask=bytarr(9) det_index_mask[j]=1b ; psf = obj_cl->GetData(CLASS_NAME = 'HSI_PSF', /SUM_PSF, xy_pixel=xy ) ; OBSOLETE: psf_obj=hsi_psf() ; ejs mar 30, 2004 psf=psf_obj->getdata( /SUM_PSF, xy_pixel=dim/2., $ obs_time_interval=map_interval, $ XYOFFSET=xyoffset, $ energy_band=energy_band, $ use_auto_time_bin=0, $ time_bin_def=[1,1,1,1,1,1,1,1,1], $ det_index_mask=det_index_mask ) rdim=psf_obj->get(/rmap_dim) ;image dimension in polar coordinates cpsf1=reform(psf,rdim(0),rdim(1)) ;transform to xy cpsf1= hsi_annsec2xy(cpsf1,psf_obj) cpsf(*,*,j)=cpsf1(0:63,0:63) if n_elements(cpsf1(*,0)) NE 64 then message,'WARNING: im dimension not 64x64',/info if n_elements(cpsf1(0,*)) NE 64 then message,'WARNING: im dimension not 64x64',/info endfor ;j message,'STEP 3 convolve the bproj maps with the corresponding PSFs',/info STEP3: ; convolve the bproj maps with the corresponding PSFs crosscor9=fltarr(9) for j=0,8 do begin if max(maparr(*,*,j)) GT 0 then $ crosscor9[j]=max(float(crosscorr(maparr(*,*,j),cpsf(*,*,j)))) endfor if keyword_set(VERB) then begin message,'SC Max(cross-corr)',/info for j=0,8 do message,strjoin([j,crosscor9[j]]),/info endif if keyword_set(FILE) then begin printf,1,'SC Max(cross-corr)' for j=0,8 do printf,1,j,crosscor9[j] endif Kmax_mask=1b+bytarr(9) for j=0,3 do Kmax_mask[j] = (crosscor9[j] GT crosscorr_min) message,'Kmax_mask='+strjoin(Kmax_mask),/info if keyword_set(FILE) then printf,1,'Kmax_mask=',Kmax_mask if keyword_set(PLOT) then begin window,/free,xsize=n_elements(maparr(*,0,0))*ndets,ysize=2*n_elements(maparr(0,*,0)) loadct, file=getenv('SSW')+'/hessi/dbase/colors_hessi.tbl', 41 for j=0,8 do tvscl,maparr(*,*,j),j for j=0,8 do tvscl,cpsf(*,*,j),j+9 endif if keyword_set(CHK_SINGLE) OR keyword_set(CHK_CIRCULARITY) then begin message,'STEP 4 MAKE CLEAN MAP for SC 3-8',/info STEP4: ; MAKE CLEAN MAP for SC 3-8 det_index_mask=[0,0,1,1,1,1,1,1,0] obj_im -> set,det_index_mask=det_index_mask obj_im -> set,image_algorithm='clean' obj_im -> set,clean_niter=333 ; Usually big enough clnmap=obj_im->getdata(clean_show_maps=0) obj_im -> plotman if keyword_set(CHK_SINGLE) then begin message,'STEP 5 CHECK THAT SOURCE IS SINGLE',/info STEP5: ; CHECK THAT SOURCE IS SINGLE ncontours=fltarr(9) for j=0,8 do ncontours(j)=check_multiplicity(clnmap,(j+1)*0.1) if keyword_set(VERB) then begin message,'# OF CONTOURS AT 10,20,...,90% levels:',/info message,strjoin(ncontours),/info w=max(where(ncontours GT 1 ,nw)) if nw GT 0 then message,'Source is single down to '+strcompress(20+w*10)+' percent level',/info else message,'Source is not single.',/info endif if keyword_set(FILE) then begin printf,1,'# OF CONTOURS AT 10,20,...,90% levels: printf,1,ncontours w=max(where(ncontours GT 1 ,nw)) if nw GT 0 then printf,1,'Source is single down to ',20+w*10,' percent level' else message,'Source is not single.',/info endif ; Stop if the source isn't single at the single_level contour. if ncontours(single_level/10-1) NE 1 then begin ERRMSG='Source not single at 20% level' message,errmsg endif endif ; chk_single if keyword_set(CHK_CIRCULARITY) then begin ; Check the circularity of the clean map source axes_ratio=check_circularity(clnmap,0.50) if keyword_set(VERB) then message,'Axes Ratio (max-min)/(max+min)='+strjoin(axes_ratio),/info if keyword_set(FILE) then printf,1,'Axes Ratio (max-min)/(max+min)=',axes_ratio if axes_ratio GT max_axes_ratio then begin ERRMSG='Axes ratio too large:'+strcompress(axes_ratio) message,errmsg endif endif ; CHK_CIRCULARITY endif ; CHK_SINGLE OR CHK_CIRCULARITY message,'STEP 6 GET CALIBRATED EVENTLISTS',/info STEP6: ; GET CALIBRATED EVENTLISTS ; Create a pointer to the calibrated eventlists cbe=obj_im->getdata(class_name='hsi_calib_eventlist',$ energy_band=energy_band,det_index=[1,1,1,1,1,1,1,1,1],$ OBS_TIME_INTERVAL=obs_time_interval, $ XYOFFSET=xyoffset, $ use_auto_time_bin=0, $ time_bin_def=[1,1,1,1,1,1,1,1,1], $ /NO_SPATIAL_FREQUENCY_WEIGHT ) message,'Got cbe',/info ;stop message,'STEP 7 FIND & USE DISTANCE TO SPIN AXIS',/info STEP7: ; FIND & USE DISTANCE TO SPIN AXIS ti='PLOT OF SUN-CENTER VECTOR IN SPACECRAFT FRAME if keyword_set(PLOT) then begin window,/free plot,(*cbe[4,0]).dx,(*cbe[4,0]).dy,tit=ti endif message,'Finding spin axis...',/info abxy=hsi_find_spinaxis(cbe,4,/verb) xspin=abxy[2] & yspin=abxy[3] message,'Spin axis at'+strcompress(abxy[2])+strcompress(abxy[3]),/info message,'Map ctr at'+strcompress(xyoffset[0])+strcompress(xyoffset[1]),/info spindist=sqrt((abxy[2]-xyoffset[0])^2+(abxy[3]-xyoffset[1])^2) message,'Distance of map center from spin axis='+strcompress(spindist)+' asec',/info if keyword_set(FILE) then begin printf,1,'Spin axis at'+strcompress(abxy[2])+strcompress(abxy[3]) printf,1,'Map ctr at'+strcompress(xyoffset[0])+strcompress(xyoffset[1]) printf,1,'Distance of map center from spin axis='+strcompress(spindist)+' asec' endif Kmin_mask=1b+bytarr(ndets) ; 1 if spindist > 2*pitch, 0 otherwise pitch=(hsi_grid_parameters()).pitch for j=6,8 do Kmin_mask[j] = ( spindist GT 2*pitch[j] ) ; 1 if true, 0 if false ; The subcollimator should not be used for bproj if kmin_mask[j] = 0 message,'Kmin_mask='+strjoin(fix(Kmin_mask)),/info message,'STEP 8 FINALIZE KMASK (WAVENUMBER MASK)',/info STEP8: ; FINALIZE KMASK (WAVENUMBER MASK)' ; If energy band doesn't extend above 20 keV, set kmask[1]=0 Kmask = Kmin_mask AND Kmax_mask ; If det 2 is on, it can be used if energy > 20 keV: if max(maparr(*,*,1)) GT 0 then Kmask[1]=(energy_band[1] GT 20.) message,'Kmask='+strjoin(fix(Kmask)),/info if keyword_set(FILE) then printf,1,'Kmask=',Kmask message,'STEP 9 Get relative amplitudes',/info STEP9: ; Get relative amplitudes count=fltarr(9) & pbar=fltarr(9) & bmax=fltarr(9) & relamp=fltarr(9) for i=0,8 do begin j=i gtran=(*cbe[j,0]).gridtran phase=(*cbe[j,0]).phase_map_ctr mamp=(*cbe[j,0]).modamp pbar[j]=mean(gtran*(1+mamp*cos(phase))) count[j]=total((*cbe[j,0]).count) bmax[j]=max(maparr(*,*,j)) relamp[j]=bmax[j]*pbar[j]/count[j] endfor if keyword_set(FILE) then begin printf,1,'RELATIVE AMPLITUDES:' FOR J=0,8 DO printf,1,j,relamp[j] endif message,'STEP 10: GAUSSIAN INTERPOLATE RELAMP ON (LOG(R),K^2) SPACE',/info STEP10: ;GAUSSIAN INTERPOLATE RELAMP ON (LOG(R),K^2) SPACE ; Create set of wavenumbers kk: nn=800. ; number of wavenumbers ; Check to see if this goes high enough! kk=findgen(nn)*4./nn ; wavenumbers to span RHESSI (0.017-1.389) k=2*!pi/pitch idet=where(kmax_mask NE 0); Do this using only legal subcollimators: ; Interpolate in (log(relamp), k^2) plane: rel_int=exp(interpol(alog([relamp(idet),1]),[k(idet)^2,0],kk^2)) if keyword_set(PLOT) then plot,kk,rel_int,tit='INTERPOLATED RELATIVE AMPLITUDE',$ XTIT='WAVENUMBER (INVERSE ARCSEC)',YTIT='RELATIVE AMPLITUDE' message,'STEP 11: BESSEL TRANSFORM THE INTERPOLATED RELMPS TO GET SIZE PROFILES',/info STEP11: ; BESSEL TRANSFORM THE INTERPOLATED RELMPS TO GET SIZE PROFILES s=fltarr(nn) r=findgen(nn)*80./nn ; if rel_int=exp(-0.5*(kk*a)^2) then s=0.25*exp(-0.5*(r/a)^2) for i=0,nn-1 do s(i)=bessel_trans(rel_int,kk,r[i]) s=s/s[0] w1=max(where(s GE 0.5)) & w2=min(where(s LE 0.5)) fwhm=r[w1]+r[w2] fwhm=strcompress(strmid(fix((fwhm+.05)*10)/10.,3,6),/rem) message,'fwhm='+strcompress(fwhm)+' asec',/info if keyword_set(FILE) then printf,1,'FWHM=',fwhm,' asec' cumulative=fltarr(nn) for k1=0,nn-1 do cumulative[k1]=total(s(k1:nn-1)*r(k1:nn-1)) cumulative=cumulative/cumulative[0] ti='SPATIAL PROFILE '+map_interval[0]+'-'+strmid(map_interval[1],11)+ $ ' '+strcompress(fix(energy_band(0)),/rem)+'-'+strcompress(fix(energy_band(1)),/rem)+" keV' plot,r,s,tit=ti,xtit='ARCSEC',YTIT='RELATIVE INTENSITY',XRA=[0,20] oplot,r,cumulative,line=1 oplot,[0,r[nn-1]],[0,0] xyouts,r[15],0.80,'FWHM = '+fwhm+' ARCSEC' xyouts,r[15],0.75,'PROFILE' & oplot,[r[25],r[35]],0.76*[1,1] xyouts,r[15],0.70,'INTEGRAL' & oplot,[r[25],r[35]],0.71*[1,1],LINE=1 stop if keyword_set(SAVE) then begin ; fname=str_replace(map_interval[0],'/','-') ;; fname=str_replace(fname,' ','_') ;; fname=str_replace(fname,':','') ;; erange=strcompress(fix(energy_band[0]),/rem)+'-'+strcompress(fix(energy_band[1]),/rem) ;; fname='bproj_'+fname+'E'+erange+'.sav' ;; message,'Writing sav file for '+fname,/info endif ; ; if ERRMSG NE '' then message,'ABNORMAL EXIT. RETURNING.',/info message,errmsg,/info if keyword_set(FILE) then begin if ERRMSG NE '' then printf,1,'ABNORMAL EXIT. printf,1,ERRMSG close,1 endif return end