;+ ; PROJECT: SPEX ; NAME: drm_correct_albedo ; ; ; PURPOSE: Function to include photospheric albedo into DRM for original spex or object ; spex (OSPEX). In SPEX call through drm_albedo. In OSPEX, is called in getdata for ; spex_drm object if spex_albedo_correct is set. ; ; ; CATEGORY: SPEX, spectral analysis ; ; ; CALLING SEQUENCE: ; ; drm_correct_albedo,theta=theta,anisotropy=anisotropy, drm=drm, ph_edges=ph_edges ; ; ; INPUTS KEYWORDS: ; theta (degrees) - heliocentric angle of the flare (default is theta =45 degrees) ; anisotropy - a coeficient showing the ratio of the flux in observer direction ; to the flux downwards ; if anisotropy=1 (default) the source is isotropic ; drm - uncorrected drm ; ph_edges - energy edges in photon space in keV (2xN) ; verbose - if set, print some informational messages ; ; USES precomputed green functions from files green_compton_mu***.dat, where *** is cos(theta) ; These files are stored in $SSWDB_XRAY/albedo. ; ; OUTPUTS: ; drm - detector response matrix per input photon/cm2 in cnts/cm2/keV with albedo ; correction ; ; SIDE EFFECTS: ; none ; ; RESTRICTIONS: ; works only for integer binning of edges and E_in ; ; PROCEDURE: ; none ; ; MODIFICATION HISTORY: ; eduard@astro.gla.ac.uk, 5-July-2004 ; 6-Aug-2004, Kim Tolbert, ; - This was called drm_albedo. Now drm_albedo calls this. ; - Changed angle and anisotropy input to keywords ; - Added input keywords drm and ph_edges so that this ; could be called for spex or ospex. ; - Default location for files is in $SSWDB_XRAY/albedo. Also tries current directory. ; If tries to read a file that's not found, catch will trap the error. ; - Added catch for errors so won't crash. Returns -1 if error. ; - Abort if ph_edges aren't integers ; eduard@astro.gla.ac.uk, 10-Nov-2004 ; - integerpolation of Green's matrix added ; - integer energy edges check removed ;- function drm_correct_albedo,theta=angle,anisotropy=anisotropy, drm=drm_in, ph_edges=ph_edges, verbose=verbose ;This will catch any errors (e.g. can't find the data file) print error message, ; and return -1 in drm_out. catch,error if error ne 0 then begin catch,/cancel print, !error_state.msg return, -1 endif ; Currently program only works if ph_edges are integers. ;diff = ph_edges - round(ph_edges) ;if min(diff) gt .001 or max(diff) gt .001 then $ ; message,'Error. Photon energy edges must be integers for albedo correction (currently).' verbose = keyword_set(verbose) drm = drm_in angle=fcheck(angle,45.) anisotropy=fcheck(anisotropy,1.) ; the value that mimics anisotropy of the source mu=cos(angle*!PI/180.) if verbose then message,/info,$ 'Albedo correction for a source at ' + trim(acos(mu)*180./!PI) + $ ' degrees cos(theta) =' + trim(mu) IF ((mu GT 0.05) AND (mu LT 0.95)) THEN BEGIN if verbose then message,/info,'reading data from files ..............................................' file1='green_compton_mu'+string(format='(I3.3)',5*FLOOR(mu*20.),/print)+'.dat' temp = loc_file( path='$SSWDB_XRAY/albedo', file1, count=count) if count gt 0 then file1 = temp file2='green_compton_mu'+string(format='(I3.3)',5*CEIL(mu*20.),/print)+'.dat' temp = loc_file( path='$SSWDB_XRAY/albedo', file2, count=count,err=err) if count gt 0 then file2 = temp restore,file1 p1=p restore,file2 p2=p alb=p1.albedo+(p2.albedo-p1.albedo)*(mu - float(floor(mu*20))/20.) END ELSE BEGIN IF (mu LT 0.05) THEN BEGIN message,/info,'Warning! Assuming heliocentric angle =' + trim(acos(0.05)*180./!PI) file='green_compton_mu'+string(format='(I3.3)',5*FLOOR(0.05*20.),/print)+'.dat' temp = loc_file( path='$SSWDB_XRAY/albedo', file, count=count) if count gt 0 then file = temp restore,file alb=p.albedo ENDIF IF (mu GT 0.95) THEN BEGIN message,/info,'Warning! Assuming heliocentric angle =' + trim(acos(0.95)*180./!PI) file='green_compton_mu'+string(format='(I3.3)',5*FLOOR(0.95*20.),/print)+'.dat' temp = loc_file( path='$SSWDB_XRAY/albedo', file, count=count) if count gt 0 then file = temp restore,file alb=p.albedo ENDIF ENDELSE alb=transpose(alb) Nc =N_elements(drm(*,0)) Nph=N_elements(drm(0,*)) ;Anew=fltarr(Nph,Nph) ; rebinning of green matrix e_ph_m =(ph_edges(0,*)+ph_edges(1,*))/2. de_ph =(ph_edges(1,*)-ph_edges(0,*)) IX=e_ph_m-min(e_ph_m) Anew=INTERPOLATE(Alb, IX, IX,/GRID) for i=0, N_elements(e_ph_m)-1 do begin Anew(*,i)=Anew(*,i)*(Anew(*,i) GT 1e-10)*de_ph/anisotropy endfor ;Anew=0.+Anew*(Anew GT 1e-8)/anisotropy ;for i=0, Nph-1 do begin ; irange=round([ph_edges(0,i)-3.,ph_edges(1,i)-4.]) ; for j=i, Nph-1 do begin ; jrange=round([ph_edges(0,j)-3.,ph_edges(1,j)-4.]) ; IF (ph_edges(1,j) LT round(max(p.edges))) AND (ph_edges(1,i) LT round(max(p.edges))) THEN BEGIN ; Anew(i,j)=total(total(alb(irange(0):irange(1),jrange(0):jrange(1)),1) $ ; ) / anisotropy / (ph_edges(1,j)-ph_edges(0,j)) ; ;/float(irange(1)-irange(0)+1. ; ENDIF ; endfor ;endfor ;**************************************************** ;window,7 ;f1=1./p.edges(0,*)^2 ;f2=1./ph_edges(0,*)^2 ;a1=(alb##(f1)) ;a2=(anew1##(f2)) ;a3=anew##(f2) ;plot,p.edges(0,*),a1/f1,/xlog,PSYM=3 ;oplot,ph_edges(0,*),a2/f2,PSYM=4,color=2 ;oplot,ph_edges(0,*),a3/f2,line=1,color=7 ;**************************************************** Anew=transpose(Anew) one=fltarr(Nph,Nph) for k=0, Nph-1 do one(k,k)=1. ;inverted albedo correction drm_albedo=one+Anew drm_old=drm drm=drm_old # transpose(drm_albedo) ;test albedo correction in count space ;e =(ph_edges(0,*)+ph_edges(1,*))/2. ;de=transpose(ph_edges(1,*)-ph_edges(0,*)) ;ftest=e^(-2) ;old=drm#(transpose(ftest)*de) ;window,8 ;plot_oo,edges(0,*),old ;new=drm#(transpose(ftest)*de) ;oplot,edges(0,*),new,line=1 ;new2=drm_old#(transpose(drm_albedo##(ftest*de))) ;oplot,edges(0,*),new2,line=2 if verbose then begin message,/info, 'DRM without albedo correction min, max = '+trim(min(drm_old))+' '+trim(max(drm_old)) message,/info, 'DRM with albedo correction min, max = '+trim(min(drm))+' '+trim(max(drm)) message,/info, 'Photospheric albedo included into DRM' endif return, drm end