pro imspec_event, event ;+ ; this is the main event handler, which handles all but time-selection ; events. (For technical reasons, these are treated in "time_setps.pro") ;- widget_control, event.id, get_uvalue=uv widget_control, event.top, get_uvalue=data, /no_copy case uv.name of "help" : show_help,event "image": begin ; region selection in the RGB image define_img_w,Nxd,Nyd,pos1,pos2 if (event.x gt pos1[0]) and (event.x lt pos2[0]) and $ (event.y gt pos1[1]) and (event.y lt pos2[1]) then begin if event.type eq 0 then begin ; start point (DOWN) ; FROM HERE ON NEW (PixMap BUFFER VERSION) Widget_Control, event.id, Draw_Motion_Events=1 Window, /Free, /Pixmap, XSize=Nxd, YSize=Nyd uv.pixID = !D.Window Device, Copy=[0, 0, Nxd, Nyd, 0, 0, data.imgw_id] ; END OF NEW x1 = event.x y1 = event.y uv.cursor[0] = x1 uv.cursor[1] = y1 uv.cursor[2] = x1 uv.cursor[3] = y1 widget_control,event.id,set_uvalue=uv endif ; if event.type eq 2 and uv.cursor[0] ne -1 then begin ; motion (OLD XOR VERSION) ; widget_control, event.id, get_value=index ; wset, index ; device, set_graphics = 6 ; x1 = uv.cursor[0] ; y1 = uv.cursor[1] ; x2 = uv.cursor[2] ; y2 = uv.cursor[3] ; plots, [x1,x2,x2,x1,x1],[y1,y1,y2,y2,y1],/device ; x2 = event.x ; y2 = event.y ; plots, [x1,x2,x2,x1,x1],[y1,y1,y2,y2,y1],/device ; uv.cursor[2] = x2 ; uv.cursor[3] = y2 ; widget_control,event.id,set_uvalue=uv ; device, set_graphics=3 ; endif if event.type eq 2 then begin ; motion (PixMap BUFFER VERSION) WSet, data.imgw_id Device, Copy=[0, 0, Nxd, Nyd, 0, 0, uv.pixID] uv.cursor[2] = event.x uv.cursor[3] = event.y get_region, uv.cursor[0],uv.cursor[1],event.x,event.y,data,region drawbox,region widget_control,event.id,set_uvalue=uv endif if event.type eq 1 then begin ; end point (UP) ; FROM HERE ON NEW (PixMap BUFFER VERSION) WSet, data.imgw_id Device, Copy=[0, 0, Nxd, Nyd, 0, 0, data.imgw_id] WDelete, uv.pixID Widget_Control, event.id, Draw_Motion_Events=0, Clear_Events=1 ; END OF NEW x2 = event.x y2 = event.y x1 = (uv.cursor)[0] y1 = (uv.cursor)[1] ;if (x1 ne x2) and (y1 ne y2) then begin source_lab = data.src_nr data.src_nr = data.src_nr + 1 widget_control, event.id, get_value=index wset, index get_region,x1,y1,x2,y2,data,region drawbox,region,source_lab=source_lab get_spec, region,data,spec mk_spec_w, event,data,spec,source_lab,region,base_id show_spec,base_id,data.time_index,data.info uv.cursor[0]=-1 widget_control,event.id,set_uvalue=uv ;endif endif endif end "region_style": begin case event.index of 0 : data.region_style = 'rectangle' 1 : data.region_style = 'ellipse' 2 : data.region_style = 'contour' endcase end "spec_style": begin ; linlog/loglog selection in the spectrum widget_control,event.id,get_uvalue=uv_tmp old_style = uv_tmp.style if old_style eq 'loglog' then new_style='linlog' if old_style eq 'linlog' then new_style='loglog' uv_tmp.style = new_style widget_control,event.id,set_uvalue=uv_tmp,set_value=old_style widget_control,uv.base_id,get_uvalue=uv_base uv_base.style=new_style widget_control,uv.base_id,set_uvalue=uv_base ifit_old = uv_base.this_fit for ifit=0,1 do begin if uv_base.fitted[ifit] then begin uv_base.this_fit = ifit+1 widget_control,uv.base_id,set_uvalue = uv_base mk_fit, uv.base_id,data.time_index,data.info widget_control,uv.base_id,get_uvalue = uv_base endif endfor uv_base.this_fit = ifit_old widget_control,uv.base_id,set_uvalue=uv_base show_spec,uv.base_id,data.time_index,data.info end "new fit":begin widget_control,uv.base_id,get_uvalue=uv_base uv_base.this_fit = (uv_base.this_fit+1) < 2 uv_base.fitted[uv_base.this_fit-1] = 1 widget_control,uv.base_id,set_uvalue=uv_base mk_fit, uv.base_id,data.time_index,data.info show_spec,uv.base_id,data.time_index,data.info end "delete fit": begin widget_control,uv.base_id,get_uvalue=uv_base if uv_base.this_fit gt 0 then uv_base.fitted[uv_base.this_fit-1] = 0 uv_base.this_fit = max(where(uv_base.fitted ne 0))+1 widget_control,uv.base_id,set_uvalue=uv_base show_spec,uv.base_id,data.time_index,data.info end "ShowEvolution": begin if data.Nt gt 1 then begin get_spec_w,event,spec_ids if spec_ids[0] ne -1 then begin get_evol,event,data,fit_results if fit_results[0].source_lab ne 0 then begin get_evol_w,event,evol_id if evol_id eq -1 then mk_evol_w,event ; get_evol,event,data,fit_results plot_evol,event,data,fit_results endif else print, 'MAKE A SPECTRUM FIT FIRST!' endif else print, 'SELECT A REGION FIRST!' endif else print, 'NO SPECTRUM EVOLUTION POSSIBLE - THERE IS A SINGLE TIME ONLY!' end "ShowLightCurves": begin if data.Nt gt 1 then begin get_spec_w, event,spec_ids if spec_ids[0] ne -1 then begin get_lightcurve_w,event,lightcurve_id if lightcurve_id eq -1 then mk_lightcurve_w,event get_lightcurves,event,data,lightcurves plot_lightcurves,event,data,lightcurves endif else print, 'SELECT A REGION FIRST!' endif else print, 'NO LIGHTCURVE POSSIBLE - THERE IS A SINGLE TIME ONLY! end "RGB": begin widget_control,event.top,set_uvalue=data if n_elements(data.E) ge 3 then begin data.image_style='RGB' widget_control, widget_info(event.id,/parent), get_uvalue=uv_parent elab0 = string(min(data.E-data.dE/2),max(data.E+data.dE/2),format='(i0,"..",i0," keV")') widget_control,uv_parent.energy_label_id,set_value=elab0 redraw,event,data,/no_contour endif end "E_slider" : begin widget_control, event.id, get_value=v widget_control, widget_info(event.id,/parent), get_uvalue=uv_parent data.image_style='SLICE' data.energy_index = v E1 = data.E[v] - 0.5*data.dE[v] E2 = data.E[v] + 0.5*data.dE[v] widget_control,uv_parent.energy_label_id,set_value=string(E1,E2,format='(i0,"-",i0," keV")') redraw,event,data,/no_contour widget_control, event.top, set_uvalue=data, /no_copy end "MakeASCII": mk_ASCII,event,data "MakePS": mk_PS, event,data "MakeIDLSAVE": mk_IDLSAVE,event,data "redraw": redraw,event,data,/no_contour "spec_drw": begin ; interactive spectrum fitting widget_control,uv.base_id,get_uvalue=uv_base if uv_base.this_fit ne 0 then begin ; convention: this_fit = (0,1,2) for (none,1st,2nd) fit define_spec_w,Nxd,Nyd,pos1,pos2 if (event.x gt pos1[0]) and (event.x lt pos2[0]) and $ (event.y gt pos1[1]) and (event.y lt pos2[1]) then begin if event.type eq 0 then begin ; button pressed speccoord2device,uv_base.Ef1,uv_base.Sf1[*,data.time_index],uv_base,xf1,yf1 speccoord2device,uv_base.Ef2,uv_base.Sf2[*,data.time_index],uv_base,xf2,yf2 ifit = -1 for this_fit=0,1 do begin if uv_base.fitted[this_fit] eq 1 then begin if norm([event.x-xf1[this_fit],event.y-yf1[this_fit]]) le 8. $ or norm([event.x-xf2[this_fit],event.y-yf2[this_fit]]) le 8. then ifit=this_fit endif endfor if ifit ne -1 then begin uv_base.this_fit = ifit+1 widget_control,uv.base_id,set_uvalue=uv_base if norm([event.x-xf1[ifit],event.y-yf1[ifit]]) le 8. then uv.changing = 1 if norm([event.x-xf2[ifit],event.y-yf2[ifit]]) le 8. then uv.changing = 2 uv.cursor=[event.x,event.y] widget_control,event.id,set_uvalue=uv endif endif if event.type eq 2 and uv.changing ne 0 then begin ; button shifted (drag) widget_control, event.id, get_value=index wset, index device, set_graphics = 6 ; XOR mode plots, uv.cursor[0],uv.cursor[1],/device,psym=4,symsize=2,thick=2 plots, event.x,event.y,/device,psym=4,symsize=2,thick=2 uv.cursor[0] = event.x uv.cursor[1] = event.y widget_control,event.id,set_uvalue=uv device, set_graphics=3 ; overwrite mode endif if event.type eq 1 then begin ; button released (drop) if uv.changing ne 0 then begin ifit = uv_base.this_fit-1 device2speccoord,uv_base,data.time_index,event.x,event.y,xc,yc if uv.changing eq 1 then uv_base.Ef1[ifit]=xc if uv.changing eq 2 then uv_base.Ef2[ifit]=xc if uv_base.Ef2[ifit] lt uv_base.Ef1[ifit] then begin Ef1 = uv_base.Ef1[ifit] uv_base.Ef1[ifit] = uv_base.Ef2[ifit] uv_base.Ef2[ifit] = Ef1 endif widget_control,uv.base_id,set_uvalue=uv_base mk_fit,uv.base_id,data.time_index,data.info show_spec,uv.base_id,data.time_index,data.info uv.changing = 0 widget_control,event.id,set_uvalue=uv endif endif endif endif end "squit": begin widget_control,uv.base_id,/destroy redraw,event,data,/no_contour end "evquit": begin widget_control,uv.base_id,/destroy redraw,event,data,/no_contour end "lcquit": begin widget_control,uv.base_id,/destroy redraw,event,data,/no_contour end "hquit": begin widget_control,event.top,/destroy if widget_info(uv.event,/active) then redraw,uv.event,data,/no_contour end "quit" : widget_control,event.top,/destroy endcase if xalive(event.top) and keyword_set(data) then widget_control, event.top, set_uvalue=data, /no_copy end pro time_steps,event ;+ ; handles forward/backward time stepping events ;- widget_control, event.top, get_uvalue=data, /no_copy widget_control, event.id, get_value=v widget_control, widget_info(event.id,/parent), get_uvalue=uv_parent case v of "<": if data.time_index gt 0 then data.time_index=data.time_index-1 ">": if data.time_index lt data.Nt-1 then data.time_index=data.time_index+1 endcase widget_control,uv_parent.time_label_id,set_value=string(data.time_index+1,format='(i0)')$ +'/'+string(data.Nt,format='(i0)') redraw,event,data,/no_contour widget_control, event.top, set_uvalue=data, /no_copy end pro movie,event ;+ ; handles movie events. Uses program-driven events (interrupts) in conjecture ; with a logical flip-flop constructed from two flags: "movie_flag" and ; "movie_killed", which are stored in the event.top data structure (because ; the program-driven events can not carry user values). The possible state ; of [movie_flag,movie_killed] are [0,0] (initial), [1,0] (running), ; [1,1] (running, to be killed by next intrerrupt) and [0,1] (killed). ;- widget_control, event.top, get_uvalue=data, /no_copy widget_control, event.id, get_value=v widget_control, widget_info(event.id,/parent), get_uvalue=uv_parent case v of "<<": begin if data.time_index gt 0 then begin data.time_index=data.time_index-1 if data.movie_killed ne 1 then data.movie_flag = 1 if data.movie_killed eq 1 then data.movie_killed = 0 endif else data.movie_flag = 0 end ">>": begin if data.time_index lt data.Nt-1 then begin data.time_index=data.time_index+1 if data.movie_killed ne 1 then data.movie_flag = 1 if data.movie_killed eq 1 then data.movie_killed = 0 endif else data.movie_flag = 0 end "STOP": begin data.movie_flag = 0 data.movie_killed = 1 end endcase if data.movie_flag eq 1 and data.movie_killed ne 1 then begin widget_control,uv_parent.time_label_id,set_value=string(data.time_index+1,format='(i0)')$ +'/'+string(data.Nt,format='(i0)') redraw,event,data,/no_contour widget_control,event.id,timer=0.01 endif widget_control, event.top, set_uvalue=data, /no_copy end pro show_help,parent_event ;+ ; user instructions - needs the file `imspec_help.txt' ;- file = hsi_loc_file('imspec_help.txt', path='$HSI_DOC', count=count,/no_dialog) if count gt 0 then begin text = rd_ascii(file[0], error=error) if not error then dummy = dialog_message(text,title='HELP for imspec',/info) endif else error = 1 if error then dummy=dialog_message('Error opening or reading help file (imspec_help.txt)') ;text = rd_ascii('imspec_help.txt') ;dummy = dialog_message(text,title='HELP for imspec') end pro mk_fit, base_id,time_index,info ;+ ; least-square fit to a straight line in either linlog or loglog coorinates ; of the energy spectrum S(E). ;- widget_control,base_id,get_uvalue=uv if uv.this_fit ne 0 then begin ifit = uv.this_fit-1 Ef1 = uv.Ef1[ifit] ; fit boundaries Ef2 = uv.Ef2[ifit] hilfs = min(abs(Ef1-uv.E),i1) hilfs = min(abs(Ef2-uv.E),i2) if i1 ne i2 then begin Ef = (uv.E)[i1:i2] Sf = (uv.spec)[i1:i2,time_index] valid= where(Sf gt 0) ; ignore negative fluxes Ef = Ef[valid] Sf = Sf[valid] ; measure_errors = sqrt(Sf) ; errors of ln(Sf), with Poissonian Sf if uv.style eq 'loglog' then begin ; result = linfit(alog(Ef),alog(Sf),measure_errors=measure_errors,sigma=sigma) result = linfit(alog(Ef),alog(Sf),sigma=sigma) if min(finite([result,sigma])) ne 0 then begin ; catch up IEEE NaN and Inf exceptions Sf1 = exp(result[0] + result[1]*alog(Ef1)) Sf2 = exp(result[0] + result[1]*alog(Ef2)) expon = fix(alog10(exp(result[0])) ) mantissa = exp(result[0]) / 10.^expon A_exp = string(expon,format='(i0)') A_mant= string(mantissa,format='(f3.1)') comment = 'S(E)='+A_mant+'!9.!310!a'+A_exp+'!nE!a'+string(result[1],format='(f4.1)')+'!n' endif else print,'***** fit failed (insufficient counts)! *****' endif if uv.style eq 'linlog' then begin ; result = linfit(Ef,alog(Sf),measure_errors=measure_errors,sigma=sigma) result = linfit(Ef,alog(Sf),sigma=sigma) if min(finite([result,sigma])) ne 0 then begin ; catch up IEEE NaN and Inf exceptions Sf1 = exp(result[0] + result[1]*Ef1) Sf2 = exp(result[0] + result[1]*Ef2) if result[1] le 0 then sgn='-' else sgn = '' expon = fix(alog10(exp(result[0])) ) mantissa = exp(result[0]) / 10.^expon A_exp = string(expon,format='(i0)') A_mant= string(mantissa,format='(f3.1)') comment = 'S(E)='+A_mant+'!9.!310!a'+A_exp+'!ne!a'+sgn+'E/'+string(1./abs(result[1]),format='(f4.1)')+'!n' endif else print,'***** fit failed (insufficient counts)! *****' endif if min(finite([result,sigma])) ne 0 then begin if abs(sigma[1]) le 10 then begin uv.comment[ifit] = comment uv.Ef1[ifit] = Ef1 uv.Ef2[ifit] = Ef2 uv.Sf1[ifit,time_index] = Sf1 uv.Sf2[ifit,time_index] = Sf2 endif else print, '***** WARNING: sigma_fit too large ! *****' uv.fit_results[ifit,0,time_index] = exp(result[0]) ; Normalization A: S(E)=A*E^alpha or S(E)=A*exp(beta*E) uv.fit_results[ifit,1,time_index] = result[1] ; either power law exponent or (minus) inverse temperature, endif endif else print,'***** spectrum boundaries coincide - chose again! *****' endif widget_control,base_id,set_uvalue=uv ; depending on uv.style end pro follow,F,x0,y0,x_,y_ ;+ ; follow a closed contour F=F(x0,y0). Follows the array ; boundary,if hit. ; INPUT: F = fltarr(Nx,Ny), x0,y0 = float ; OUTPUT: x_,y_ = fltarr(N_steps) ;- Nx = (size(F))[1] Ny = (size(F))[2] N1 = fltarr(Nx,Ny) N1(where(F ge F(x0,y0))) = 1 N = fltarr(Nx+2,Ny+2) N[1:Nx,1:Ny] = N1 x_ = 0 y_ = 0 F0 = F[x0,y0] if F[x0-1,y0] lt F0 then cxo = [x0,y0] + [1.5,1] if F[x0+1,y0] lt F0 then cxo = [x0,y0] + [0.5,1] if F[x0,y0-1] lt F0 then cxo = [x0,y0] + [1,1.5] if F[x0,y0+1] lt F0 then cxo = [x0,y0] + [1,0.5] if keyword_set(cxo) then begin x = [x0,y0]+[1,1] dx1 = x-cxo dx2 = [-dx1(1),dx1(0)] dx3 = -dx2 cxn = x + dx1 xl = x+dx2 xr = x+dx3 xl1 = cxn+dx1+dx2 xr1 = cxn+dx1+dx3 Nl = n(xl(0),xl(1)) Nr = n(xr(0),xr(1)) Nl1 = n(xl1(0),xl1(1)) Nr1 = n(xr1(0),xr1(1)) if (Nl eq Nl1) and (Nr eq Nr1) then x = cxn + dx1 if (Nl ne Nl1) and (Nr eq Nr) then x = cxn + dx2 if (Nr ne Nr1) and (Nl eq Nl1) then x = cxn + dx3 if (Nl ne Nl1) and (Nr ne Nr1) then x = cxn + dx2 x00 = x cxo = cxn repeat begin dx1 = x-cxo dx2 = [-dx1(1),dx1(0)] dx3 = -dx2 cxn = x + dx1 xl = x+dx2 xr = x+dx3 xl1 = cxn+dx1+dx2 xr1 = cxn+dx1+dx3 Nl = n[xl(0),xl(1)] Nr = n[xr(0),xr(1)] Nl1 = n[xl1(0),xl1(1)] Nr1 = n[xr1(0),xr1(1)] if (Nl eq Nl1) and (Nr eq Nr1) then x = cxn + dx1 if (Nl ne Nl1) and (Nr eq Nr) then x = cxn + dx2 if (Nr ne Nr1) and (Nl eq Nl1) then x = cxn + dx3 if (Nl ne Nl1) and (Nr ne Nr1) then x = cxn + dx2 cxo = cxn x_ = [x_,x[0]] y_ = [y_,x[1]] endrep until (x(0) eq x00[0]) and (x[1] eq x00[1]) x_ = fix(x_[1:n_elements(x_)-1]-1) y_ = fix(y_[1:n_elements(y_)-1]-1) x_ = x_ > 0 x_ = x_ < (Nx-1) y_ = y_ > 0 y_ = y_ < (Ny-1) endif x_ = [x_,x_[0]] y_ = [y_,y_[0]] end pro get_region, x1,y1,x2,y2,data,region ;+ ; defines the (closed) boundary of a simply connected region in which ; the average spectrum is determined. The boundary is defined in terms ; of its x/y coordinates. (x1,y1,x2,y2) and (x_r,y_r) are in device ; coordinates. A minimum size of 8 arc sec is required. ;- common img_cache, img_, intmap_, img_known case data.region_style of "rectangle" : begin x1a = x1 x2a = x2 y1a = y1 y2a = y2 if abs(x2a-x1a) le 8 then begin x1a = 0.5*(x1a+x2a) - 4. x2a = 0.5*(x1a+x2a) + 4. endif if abs(y2a-y1a) le 8 then begin y1a = 0.5*(y1a+y2a) - 4. y2a = 0.5*(y1a+y2a) + 4. endif x_r = [x1a,x2a,x2a,x1a,x1a] y_r = [y1a,y1a,y2a,y2a,y1a] end "ellipse" : begin Nphi = 50. phi = findgen(Nphi)/(Nphi-1)*2*!pi x_r = 0.5*(x1+x2) + 0.5*(abs(x2-x1)>8.)*cos(phi) y_r = 0.5*(y1+y2) + 0.5*(abs(y2-y1)>8.)*sin(phi) end "contour" : begin intensity = (total(data.imc,3))[*,*,data.time_index] device2imgindx, x1,y1,data.x,data.y,ix,iy follow,intensity,ix,iy,ix_,iy_ imgindx2device,ix_,iy_,data.x,data.y,x_r,y_r end endcase region = {x:x_r,y:y_r} end pro drawbox,region,source_lab=source_lab ;+ ; draws and annotates a selected region ;- plots, region.x,region.y,/device if keyword_set(source_lab) then $ xyouts,mean(region.x),mean(region.y),string(source_lab,format='(i0)'),/device end pro redraw, event,data,no_contour=no_contour ;+ ; redraws region boxes, spectra and spectral and fits ;- wset, data.imgw_id show_img, data,no_contour=no_contour get_spec_w, event,spec_ids if spec_ids[0] ne -1 then begin for i = 0,n_elements(spec_ids)-1 do begin this_spec = spec_ids[i] widget_control,this_spec,get_uvalue = uv wset,data.imgw_id drawbox,uv.region,source_lab=uv.source_lab ifit_old = uv.this_fit for ifit=0,1 do begin if uv.fitted[ifit] then begin uv.this_fit = ifit+1 widget_control,this_spec,set_uvalue = uv mk_fit, this_spec,data.time_index,data.info widget_control,this_spec,get_uvalue = uv endif endfor uv.this_fit = ifit_old ; restore active fit widget_control,this_spec,set_uvalue = uv show_spec,this_spec,data.time_index,data.info endfor endif get_evol_w, event,evol_id if evol_id ne -1 then begin get_evol,event,data,fit_results plot_evol,event,data,fit_results endif get_lightcurve_w, event,lightcurve_id if lightcurve_id ne -1 then begin get_lightcurves,event,data,lightcurves plot_lightcurves,event,data,lightcurves endif end pro get_spec_w, event,spec_ids ;+ ; gets the identities of all spectrum widgets [id1,id2,...] and ; returns -1 if no valid spectrum widget is found. The spectrum ; widgets must be childs of event.top. ;- spec_ids= -1 next_spec = widget_info(event.top,find_by_uname='spec') while next_spec ne 0 do begin widget_control,next_spec,get_uvalue = uv if uv.type eq 'spectrum' then spec_ids = [spec_ids,next_spec] next_spec = widget_info(next_spec,/sibling) endwhile if n_elements(spec_ids) gt 1 then spec_ids = spec_ids[1:n_elements(spec_ids)-1] end pro get_evol_w, event,evol_id ;+ ; gets the identity of the spectral evolution widget, which must ; be a child of "event.top". Returns -1 if no evolution widget exists. ; evol_id = -1 next_spec = widget_info(event.top,find_by_uname='evol') if next_spec ne 0 then evol_id = next_spec end pro get_lightcurve_w, event,lightcurve_id ;+ ; gets the identity of the lightcurve widget, which must ; be a child of "event.top". Returns -1 if no evolution widget exists. ; lightcurve_id = -1 next_spec = widget_info(event.top,find_by_uname='lightcurve') if next_spec ne 0 then lightcurve_id = next_spec end pro mk_lightcurve_w,event ;+ ; creates a new lightcurve widget ;- define_ev_w,Nxd,Nyd lcbase= widget_base(event.top,/column,uname='lightcurve',title='Light Curves') menu = widget_base(lcbase,/row) button1 = widget_button(menu,uvalue={name:'lcquit',base_id:lcbase},value='Dismiss') lcdraw = widget_draw(lcbase,xsize=Nxd,ysize=Nyd) widget_control, lcbase,/realize xmanager,"Lightcurve", lcbase, event_handler="imspec_event" widget_control,lcbase,set_uvalue={type:'lightcurve',lc_drw:lcdraw} end pro get_lightcurves, event,data,lightcurves ;+ ; PURPOSE: Extracts the lightcurves of all selected regions. ; INPUT: a valid event rooting to the main (base) widget. ; OUTPUT: lightcurves = array of structures, each with tags ; source label and intensity (total counts, all energies). ;- get_spec_w, event,spec_ids lightcurves = {source_lab:0,intensity:fltarr(data.Nt)} if spec_ids[0] ne -1 then begin for i=0,n_elements(spec_ids)-1 do begin this_spec = spec_ids[i] widget_control,this_spec,get_uvalue = uv device2imgindx, (uv.region.x),(uv.region).y,data.x,data.y,ix,iy N_x = n_elements(data.x) N_y = n_elements(data.y) N_E = n_elements(data.E) inside = polyfillv(ix,iy,N_x,N_y) ; indices of points inside the countor (region.x,region.y) intensity = total(total((reform(data.imc,N_x*N_y,N_E,data.Nt))[inside,*,*],1),1) lightcurves = [lightcurves,{source_lab:uv.source_lab,intensity:intensity}] widget_control,this_spec,set_uvalue = uv endfor endif if n_elements(lightcurves) gt 1 then lightcurves = lightcurves[1:n_elements(lightcurves)-1] end pro plot_lightcurves, event,data,lightcurves ;+ ; shows the time light curves of the selected regions ;- get_lightcurve_w,event,lightcurve_id widget_control, lightcurve_id, get_uvalue = uv widget_control, uv.lc_drw, get_value = windex wset, windex N_curves = n_elements(lightcurves) t = 0.5*total(data.time_bands,1) ; center times of time_bands !p.multi=[0,1,N_curves] for i=0,N_curves-1 do begin utplot,t,lightcurves[i].intensity,anytim(t,/ex),xtit='time [UT]',title='Region ' $ +string(lightcurves[i].source_lab,format='(i0)')+': Light Curve',ytit='img cube units',ymarg=[4,2],psym=10 endfor !p.multi=0 end pro mk_PS, event,data ;+ ; Creates a PS file of image and spectra ;- common img_cache, img_,intmap_,img_known define_img_w,Nxd,Nyd,pos1,pos2 Nxpixel = pos2[0]-pos1[0] Nypixel = pos2[1]-pos1[1] x1 = congrid(data.x,Nxpixel,/interp) y1 = congrid(data.y,Nypixel,/interp) if img_known[data.time_index] eq 0 then mk_img_cache,data intmap = intmap_[*,*,data.time_index] psfile = 'idl.ps' print,'Program creates '+psfile+'...' dname = !d.name set_plot,'ps' device, file=psfile,/color,xsize=17,ysize=22,xoff=1,yoff=1 !p.multi=[0,2,3] start_time = anytim(data.time_bands[0,data.time_index], /ex) stop_time = anytim(data.time_bands[1,data.time_index], /ex) tme_str = string(start_time[4:6],start_time[0:2],stop_time[0:2], $ format='(i0,".",i0,".",i4,", ",i2.2,":",i2.2,":",i2.2,"-",i2.2,":",i2.2,":",i2.2)') contour, intmap_[*,*,data.time_index],x1,y1,xst=1,yst=1,c_linest=1, $ xtit='heliospheric x [arcsec]',ytit='heliospheric y [arcsec]',tit='RHESSI intensity' get_spec_w, event,spec_ids if spec_ids[0] ne -1 then begin for i=0,n_elements(spec_ids)-1 do begin widget_control,spec_ids[i],get_uvalue = uv device2imgcoord,(uv.region).x,(uv.region).y,data.x,data.y,xc,yc plots,xc,yc,thick=2 xyouts,mean(xc),mean(yc),string(uv.source_lab,format='(i0)') endfor for i=0,n_elements(spec_ids)-1 do begin widget_control,spec_ids[i],get_uvalue = uv case uv.style of 'loglog' : xlog = 1 'linlog' : xlog = 0 endcase if xlog then xr = [0.5*min(uv.E),2.5*max(uv.E)] $ else xr = [0.2*min(uv.E),1.2*max(uv.E)] plot,uv.E,uv.spec[*,data.time_index],psym=10,xst=1,yst=1, $ xtitle='E [keV]',xlog=xlog,ylog=1, $ title='Region '+string(uv.source_lab,format='(i0)')+': Spectrum', $ xrange=xr,yrange=[0.5*min(uv.spec),2.5*max(uv.spec)] if max(uv.fitted) ne 0 then begin ifit= 0 Ef1 = uv.Ef1[ifit] Ef2 = uv.Ef2[ifit] Sf1 = uv.Sf1[ifit,data.time_index] Sf2 = uv.Sf2[ifit,data.time_index] plots, [Ef1,Ef2],[Sf1,Sf2],linest=2,thick=2 plots, [Ef1,Ef2],[Sf1,Sf2],psym=4,symsize=2,thick=2 if xlog then xlab = sqrt(Ef1*Ef2) else xlab = 0.5*(Ef1+Ef2) ylab = 1.5*sqrt(Sf1*Sf2) xyouts,xlab,ylab,uv.comment[ifit],size=1.0 endif endfor endif xyouts, 0.2,1.05,/normal,'DATA FROM: '+tme_str device, /close set_plot,dname !p.multi=0 end pro mk_IDLSAVE,event,data ;+ ; creates IDL SAVE file with spectrum data - see file for an explantion. ;- save_file = 'idl.dat' get_spec_w, event,spec_ids if spec_ids[0] ne -1 then begin print, 'Program creates '+save_file time_bands = data.time_bands E1 = data.E - data.dE/2. E2 = data.E + data.dE/2. energy_bands = transpose([[E1],[E2]]) i = 0 this_spec = spec_ids[i] widget_control,this_spec,get_uvalue = uv device2imgcoord,uv.region.x,uv.region.y,data.x,data.y,xc,yc spectrum_data = {source_lab:uv.source_lab,source_region:{x:xc,y:yc},spec:uv.spec} for i=1,n_elements(spec_ids)-1 do begin this_spec = spec_ids[i] widget_control,this_spec,get_uvalue = uv device2imgcoord,uv.region.x,uv.region.y,data.x,data.y,xc,yc spectrum_data = [spectrum_data,{source_lab:uv.source_lab,source_region:{x:xc,y:yc},spec:uv.spec}] endfor get_evol, event,data,fit_results README = [ 'The variables in imspec`s idl.dat are: TIME_BANDS, ENERGY_BANDS ', $ 'SPECTRUM_DATA and FIT_RESULTS. The units are keV, seconds, and ', $ 'photons/(cm^2*s*asec^2*keV). The fit models are S(E)=A*exp(-E/kT) if ', $ 'fit_results.style=linlog, and S(E)=A*E^alpha if fit_results.style=loglog. ', $ 'The meaning of fit_results.coefs is (A,kT) if fit_results.style=linlog, ', $ 'and (A,alpha) if fit_results.style=loglog. The contour surrounding the ', $ 'source region is stored in spectrum_data.source_region (heliocentric x/y ', $ '[arc sec]). ' ] save,file=save_file,time_bands,energy_bands,spectrum_data,fit_results,README endif else print, 'NO FIT DATA TO SAVE!' end pro mk_ASCII,event,data ;+ ; prints the fit results into an ASCII file - see file for an explantion. ;- ascii_file = 'idl.asc' get_evol, event,data,fit_results Nfits = n_elements(fit_results) if fit_results[0].style ne '' then begin start_time = anytim(data.time_bands[0,data.time_index], /ex) stop_time = anytim(data.time_bands[1,data.time_index], /ex) tme_str = string(start_time[4:6],start_time[0:2],stop_time[0:2], $ format='(i0,".",i0,".",i4,", ",i2.2,":",i2.2,":",i2.2,"-",i2.2,":",i2.2,":",i2.2)') print, 'Program creates '+ascii_file close,1 openw,1,ascii_file printf,1,'DATA FROM: '+tme_str printf,1,'**************************************' for i=0,Nfits-1 do begin printf,1,'FIT STYLE: '+fit_results[i].style printf,1,'SOURCE LABEL: '+string(fit_results[i].source_lab,format='(i0)') printf,1,fit_results[i].source_region,format='("SOURCE REGION: (min(x),min(y),max(x),max(y)) = (",i0,3(",",i0),")")' printf,1,'1st/2nd FIT:',fit_results[i].ifit+1 printf,1,'ENERGY BAND [keV]: ',fit_results[i].energy_band if fit_results[i].style eq 'linlog' then begin printf,1,'FIT MODEL: S(E)=A*exp(-E/kT)' printf,1,'(rel. start time [s], rel. stop time [s], A [phot/s/cm2/asec2/keV], kT [keV]):' endif if fit_results[i].style eq 'loglog' then begin printf,1,'FIT MODEL: S(E)=A*(-E/1keV)^(-alpha)' printf,1,'(start time [s], stop time [s], A [phot/s/cm2/asec2/keV], alpha):' endif printf,1,[data.time_bands-data.time_bands[0,0],fit_results[i].coefs] printf,1,'**********************************************' endfor close,1 endif else print, 'No fit data to save!' end pro read_image_cube_fits, fitsfile,imc,x,y,energy_bands,time_bands,info ;+ ; Reads a HESSI fits file with an image cube of format [N_x,N_y,N_E,N_t], ; as provided by Kim. ;- imc = hsi_image_fitsread(fitsfile=fitsfile,control=info,ebands_arr=energy_bands,times_arr=time_bands) N_x = (size(imc))[1] N_y = (size(imc))[2] x = info.xyoffset[0] + (findgen(N_x)-N_x/2)*info.pixel_size[0] y = info.xyoffset[1] + (findgen(N_y)-N_y/2)*info.pixel_size[1] ;imc = imc[*,*,*,0] ; TEST ONLY: check the case where a ;time_bands = time_bands[*,0] ; single time interval is present end pro get_rgb_filters,E,im,rfilter,gfilter,bfilter ;+ ; creates 3 filter curves, from which the (r,g,b) components ; of the true color image are computed. ; N_E = n_elements(E) S = total(total(im,1),1) ; S=fltarr(N_E) ...... spectrum, integrated over the ; whole image at the actual time. if N_E eq 3 then begin rfilter = [1,0,0] gfilter = [0,1,0] bfilter = [0,0,1] endif if N_E gt 3 then begin i = 0 ; divide S into three parts of equal counts while (total(S[0:i]) le 0.333*total(S)) do i = i+1 i_rg = i>1 while (total(S[0:i]) le 0.666*total(S)) do i = i+1 i_gb = (i>i_rg+1)<(N_E-2) rfilter = fltarr(N_E) & rfilter[0:i_rg] = 1. gfilter = fltarr(N_E) & gfilter[i_rg+1:i_gb] = 1. bfilter = fltarr(N_E) & bfilter[i_gb:N_E-1] = 1. endif rfilter = rfilter/total(rfilter) gfilter = gfilter/total(gfilter) bfilter = bfilter/total(bfilter) end pro spec2rgb, im,E,img ;+ ; PURPOSE: reduces the 3D image cube at fixed time to an RGB image of ; the format bytarr(N_x,N_y,3), where the last index indicates red, green ; and blue components. ;- N_x = (size(im))(1) N_y = (size(im))(2) N_E = (size(im))(3) get_rgb_filters,E,im,rfilter,gfilter,bfilter r = reform(reform(im,N_x*N_y,N_E)#rfilter,N_x,N_y) ; projection on RGB filters g = reform(reform(im,N_x*N_y,N_E)#gfilter,N_x,N_y) b = reform(reform(im,N_x*N_y,N_E)#bfilter,N_x,N_y) r = r/total(r) g = g/total(g) b = b/total(b) r = bytscl(r) g = bytscl(g) b = bytscl(b) img = [[[r]],[[g]],[[b]]] end pro define_img_w,Nxd,Nyd,pos1,pos2 ;+ ; defines size/positioning of RGB image draw widget (device coordinates) ;- Nxd = 360 Nyd = 360 pos1 = [0.15*Nxd,0.15*Nyd] ; image position within draw widget pos2 = [0.9*Nxd,0.9*Nyd] end pro define_spec_w,Nxd,Nyd,pos1,pos2 ;+ ; defines size/positioning of spectrum draw widget (device coordinates) ;- Nxd = 360 Nyd = 320 pos1 = [0.20*Nxd,0.18*Nyd] ; image position within draw widget pos2 = [0.93*Nxd,0.90*Nyd] end pro define_ev_w,Nxd,Nyd ;+ ; defines size/positioning of time evolution draw widget (device coordinates) ;- Nxd = 360 Nyd = 320 end pro mk_img_cache,data ;+ ; creates color image after reading of new fits file ;- common img_cache, img_, intmap_, img_known define_img_w,Nxd,Nyd,pos1,pos2 Nxpixel = pos2[0]-pos1[0] Nypixel = pos2[1]-pos1[1] if n_elements(data.E) ge 3 then begin spec2rgb,data.imc[*,*,*,data.time_index],data.E,img img_[*,*,*,data.time_index] = congrid(img,Nxpixel,Nypixel,3) intmap = congrid(total(data.imc[*,*,*,data.time_index],3),Nxpixel,Nypixel,/interp) endif else begin data.image_style = 'SLICE' intmap = congrid((total(data.imc,3))[*,*,data.time_index],Nxpixel,Nypixel,/interp) endelse intmap_[*,*,data.time_index] = intmap img_known[data.time_index] = 1 end pro show_img, data,no_contour=no_contour ;+ ; displays and annotates HESSI image on TRUE COLOR device ;- common img_cache, img_,intmap_,img_known if not keyword_set(no_contour) then no_contour=0 define_img_w,Nxd,Nyd,pos1,pos2 Nxpixel = pos2[0]-pos1[0] Nypixel = pos2[1]-pos1[1] x1 = congrid(data.x,Nxpixel,/interp) y1 = congrid(data.y,Nypixel,/interp) if data.image_style eq 'RGB' then begin if img_known[data.time_index] eq 0 then mk_img_cache,data img1 = img_[*,*,*,data.time_index] erase tv, img1,pos1[0],pos1[1],true=3 endif if data.image_style eq 'SLICE' then begin img1 = congrid(data.imc[*,*,data.energy_index,data.time_index],Nxpixel,Nypixel,/interp) erase tvscl,img1,pos1[0],pos1[1] endif intmap = intmap_[*,*,data.time_index] start_time = anytim(data.time_bands[0,data.time_index], /ex) stop_time = anytim(data.time_bands[1,data.time_index], /ex) tme_str = string(start_time[4:6],start_time[0:2],stop_time[0:2], $ format='(i0,".",i0,".",i4,", ",i2.2,":",i2.2,":",i2.2,"-",i2.2,":",i2.2,":",i2.2)') contour, intmap,x1,y1,nlevels=5,/noerase,nodata=no_contour,/device, $ xst=1,yst=1,position=[pos1,pos2],xtitle='heliocentric x [arcsec]', $ ytitle='heliocentric y [arcsec]',title=tme_str end pro device2imgcoord,xd,yd,x,y,xc,yc ;+ ; converts device to physical units ;- define_img_w,Nxd,Nyd,pos1,pos2 Nx = n_elements(x) Ny = n_elements(y) Nxpixel = pos2[0]-pos1[0] Nypixel = pos2[1]-pos1[1] xc = x[0] + (xd-pos1[0])*(x[Nx-1]-x[0])/float(Nxpixel) yc = y[0] + (yd-pos1[1])*(y[Ny-1]-y[0])/float(Nypixel) end pro device2speccoord,uv,time_index,xd,yd,xc,yc ;+ ; converts device to physical units ;- define_spec_w,Nxd,Nyd,pos1,pos2 x = uv.E y = uv.spec[*,time_index] wset,uv.drw_id if uv.style eq 'linlog' then begin x1 = !x.crange[0] x2 = !x.crange[1] endif if uv.style eq 'loglog' then begin x1 = 10^!x.crange[0] x2 = 10^!x.crange[1] endif y1 = 10^!y.crange[0] y2 = 10^!y.crange[1] Nxpixel = float(pos2[0]-pos1[0]) Nypixel = float(pos2[1]-pos1[1]) if uv.style eq 'linlog' then xc = x1+((xd-pos1[0])*(x2-x1)/Nxpixel) if uv.style eq 'loglog' then xc = x1*exp((xd-pos1[0])*(alog(x2/x1))/Nxpixel) yc = y1*exp((yd-pos1[1])*(alog(y2/y1))/Nypixel) end pro speccoord2device,xc,yc,uv,xd,yd ;+ ; converts physical to device units ;- define_spec_w,Nxd,Nyd,pos1,pos2 x = uv.E y = uv.spec wset,uv.drw_id if uv.style eq 'linlog' then begin x1 = !x.crange[0] x2 = !x.crange[1] endif if uv.style eq 'loglog' then begin x1 = 10^!x.crange[0] x2 = 10^!x.crange[1] endif y1 = 10^!y.crange[0] y2 = 10^!y.crange[1] Nxpixel = float(pos2[0]-pos1[0]) Nypixel = float(pos2[1]-pos1[1]) yd = alog(yc/y1)*Nypixel/alog(y2/y1)+pos1[1] if uv.style eq 'linlog' then xd = (xc-x1)*Nxpixel/(x2-x1)+pos1[0] if uv.style eq 'loglog' then xd = (alog(xc)-alog(x1))*Nxpixel/alog(x2/x1)+pos1[0] end pro device2imgindx, xd,yd,x,y,ix,iy ;+ ; converts spatial image cube indices to device coordinates ;- define_img_w,Nxd,Nyd,pos1,pos2 Nx = n_elements(x) Ny = n_elements(y) Nxpixel = pos2[0]-pos1[0] Nypixel = pos2[1]-pos1[1] xc = x[0] + (xd-pos1[0])*(x[Nx-1]-x[0])/float(Nxpixel) yc = y[0] + (yd-pos1[1])*(y[Ny-1]-y[0])/float(Nypixel) ix = round(Nx*(xc-min(x))/(max(x)-min(x))) iy = round(Ny*(yc-min(y))/(max(y)-min(y))) ix = ix > 1 ix = ix < (Nx-2) iy = iy > 1 iy = iy < (Ny-2) end pro imgindx2device, ix,iy,x,y,xd,yd ;+ ; converts spatial image cube indices to device coordinates ;- define_img_w,Nxd,Nyd,pos1,pos2 Nx = n_elements(x) Ny = n_elements(y) Nxpixel = pos2[0]-pos1[0] Nypixel = pos2[1]-pos1[1] xc = ix*(max(x)-min(x))/float(Nx) + min(x) yc = iy*(max(y)-min(y))/float(Ny) + min(y) xd = (xc - x[0])*float(Nxpixel)/(x[Nx-1]-x[0]) + pos1[0] yd = (yc - y[0])*float(Nypixel)/(y[Ny-1]-y[0]) + pos1[1] end pro get_spec,region,data,spec ;+ ; extracts spectrum from selected region. Input: (region.x,region.y) = closed ; countour surrounding the selected region ;- ; ***** rectangular enclosing box ***** ;x1 = min(region.x) ;y1 = min(region.y) ;x2 = max(region.x) ;y2 = max(region.y) ;device2imgcoord,[x1,x2],[y1,y2],data.x,data.y,xc,yc ;hilfs = min(abs(data.x-xc[0]),ix1) ;hilfs = min(abs(data.y-yc[0]),iy1) ;hilfs = min(abs(data.x-xc[1]),ix2) ;hilfs = min(abs(data.y-yc[1]),iy2) ;spec = total(total(data.imc[ix1:ix2,iy1:iy2,*,*],1),1)/(data.dE#replicate(1.,data.Nt)) ; ***** arbitrary region bounded by x[i],y[i], i=0,1,2,.. ***** device2imgindx, region.x,region.y,data.x,data.y,ix,iy N_x = n_elements(data.x) N_y = n_elements(data.y) N_E = n_elements(data.E) inside = polyfillv(ix,iy,N_x,N_y) ; indices of points inside the countor (region.x,region.y) ; tested by eye - ok. spec = total((reform(data.imc,N_x*N_y,N_E,data.Nt))[inside,*,*],1)/(data.dE#replicate(1.,data.Nt)) end pro mk_spec_w,event,data,spec,source_lab,region,sbase ;+ ; creates and adds a spectrum widget ;- E = data.E N_E = (size(spec))[1] N_t = (size(spec))[2] define_spec_w,Nxd,Nyd,pos1,pos2 sbase = widget_base(event.top,uname='spec',/column,/no_copy) smenu = widget_base(sbase,/row) sdraw = widget_draw(sbase,xsize=Nxd,ysize=Nyd, $ button_events=1,motion_events=1, $ uvalue={name:'spec_drw',base_id:sbase, $ changing:0,cursor:[-1,-1]}) button1 = widget_button(smenu,value='Dismiss', $ uvalue={name:'squit',base_id:sbase}) button2 = widget_button(smenu,value='linlog', $ uvalue={name:'spec_style',style:'loglog',base_id:sbase}) button3 = widget_button(smenu,value='new fit', $ uvalue={name:'new fit',base_id:sbase}) button4 = widget_button(smenu,value='delete fit', $ uvalue={name:'delete fit',base_id:sbase}) widget_control, sbase,/realize widget_control, sdraw, get_value=index wset, index iEf1 = -1 repeat iEf1 = iEf1+1 until (min(spec[iEf1,*] ne 0)) or (iEf1 ge N_E-1) iEf2 = n_elements(E) repeat iEf2 = iEf2-1 until (min(spec[iEf2,*] ne 0)) or (iEf2 le 0) if iEf2 le iEf1 then iEf2=iEf1 Ef1 = E[iEf1]*[1,1] Sf1 = [1,10]#reform(spec[iEf1,*]) Ef2 = E[iEf2]*[1,1] Sf2 = [1,10]#reform(spec[iEf2,*]) widget_control, sbase, set_uvalue={E:E,spec:spec,type:'spectrum', $ source_lab:source_lab,drw_id:index,region:region, $ style:'loglog',this_fit:0,fitted:[0,0],comment:['',''], $ Ef1:Ef1,Sf1:Sf1,Ef2:Ef2,Sf2:Sf2, $ fit_results:fltarr(2,2,data.Nt)} end pro show_spec,base_id,time_index,info ;+ ; plots spectrum in selected position box ;- widget_control,base_id,get_uvalue=uv wset, uv.drw_id define_spec_w,Nxd,Nyd,pos1,pos2 if n_elements(uv.E) gt 1 then begin case uv.style of 'loglog' : xlog = 1 'linlog' : xlog = 0 endcase if xlog then xr = [0.5*min(uv.E),2.8*max(uv.E)] $ else xr = [0.2*min(uv.E),1.4*max(uv.E)] valid = where(uv.spec gt 0) yr = [0.1*min(uv.spec[valid]),2.5*max(uv.spec[valid])] plot,uv.E,uv.spec[*,time_index],pos=[pos1,pos2],/device,psym=10,xst=1,yst=1, $ xtitle='E [keV]',ytitle='photons cm!e-2!n s!e-1!n asec!e-2!n keV!e-1!n ', $ xlog=xlog,ylog=1,title='Region '+string(uv.source_lab,format='(i0)')+': Spectrum', $ xrange=xr,yrange=yr col = [-1,200] for ifit =0,1 do begin if uv.fitted[ifit] eq 1 then begin Ef1 = uv.Ef1[ifit] Ef2 = uv.Ef2[ifit] Sf1 = uv.Sf1[ifit,time_index] Sf2 = uv.Sf2[ifit,time_index] plots, [Ef1,Ef2],[Sf1,Sf2],linest=2,thick=2,color=col[ifit] plots, [Ef1,Ef2],[Sf1,Sf2],psym=4,symsize=2,thick=2,color=col[ifit] if xlog then xlab = sqrt(Ef1*Ef2) else xlab = 0.5*(Ef1+Ef2) ylab = 1.5*sqrt(Sf1*Sf2) xyouts,xlab,ylab,uv.comment[ifit],size=1.2,color=col[ifit] endif endfor endif else xyouts,/device,0.25*Nxd,0.5*Nyd,'NO SPECTRUM - SINGLE ENERGY BIN' end pro mk_evol_w,event ;+ ; creates a new time evolution widget ;- define_ev_w,Nxd,Nyd evbase= widget_base(event.top,/column,uname='evol',title='Fit Results') menu = widget_base(evbase,/row) button1 = widget_button(menu,uvalue={name:'evquit',base_id:evbase},value='Dismiss') ;title = widget_label(menu,value='Fitted Coefficients as a function of Time') evdraw = widget_draw(evbase,xsize=Nxd,ysize=Nyd) widget_control, evbase,/realize xmanager,"Evol", evbase, event_handler="imspec_event" widget_control,evbase,set_uvalue={type:'evolution',evol_drw:evdraw} end pro get_evol, event,data,fit_results ;+ ; Extracts all available results of spectral fits. Fills in zeros at ; non-computed times. ; INPUT: a valid event rooting to the main (base) widget. ; OUTPUT: fit_results = array of structures, each with tags style ; (loglog or linlog), source label, and the fitted coefficients ; coefs = fltarr(2,Nt). ;- get_spec_w, event,spec_ids fit_results = {style:'',source_lab:0,source_region:fltarr(4),ifit:0, $ energy_band:[0.,0.],coefs:fltarr(2,data.Nt)} if spec_ids[0] ne -1 then begin for i=0,n_elements(spec_ids)-1 do begin this_spec = spec_ids[i] widget_control,this_spec,get_uvalue = uv ifit_old = uv.this_fit for ifit=0,1 do begin if uv.fitted[ifit] eq 1 then begin coefs = fltarr(2,data.Nt) for time_index=0,data.Nt-1 do begin coefs[*,time_index] = uv.fit_results[ifit,*,time_index] endfor device2imgcoord,(uv.region).x,(uv.region).y,data.x,data.y,xc,yc fit_results = [fit_results,{style:uv.style,source_lab:uv.source_lab, $ source_region:[min(xc),min(yc),max(xc),max(yc)],ifit:ifit, $ energy_band:float([uv.Ef1[ifit],uv.Ef2[ifit]]),coefs:coefs}] endif endfor uv.this_fit = ifit_old widget_control,this_spec,set_uvalue = uv endfor endif if n_elements(fit_results) gt 1 then fit_results = fit_results[1:n_elements(fit_results)-1] end pro plot_evol, event,data,fit_results ;+ ; shows the time evolution of the fitted parameters alpha, T, and A. ;- get_evol_w,event,evol_id widget_control, evol_id, get_uvalue = uv widget_control, uv.evol_drw, get_value = windex wset, windex col = [-1,200] t = 0.5*total(data.time_bands,1) ; center times of time_bands ;if n_elements(t) gt 1 then begin N_fits = n_elements(fit_results) if N_fits ge 1 and (fit_results[0].source_lab) ne 0 then begin loglog_fits = where(fit_results.style eq 'loglog') linlog_fits = where(fit_results.style eq 'linlog') n_plots = 2 if (loglog_fits[0] ne -1) and (linlog_fits[0] ne -1) then n_plots = 3 !p.multi= [0,1,n_plots] if loglog_fits[0] ne -1 then begin alpha = -reform((fit_results[loglog_fits].coefs)[1,*,*]) valid = where((alpha ne 'Inf') and (alpha ne '-Inf') and (alpha ne 0)) if valid[0] ne -1 then begin utplot,t,[min(alpha[valid])-0.5,max(alpha[valid])+0.5],anytim(t,/ex), $ xtit='',ytit='!4a!3',/nodata,charsize=1.2,ymarg=[2,1],yst=1 for i=0,n_elements(loglog_fits)-1 do begin i_valid = where(alpha[*,i] ne 0) oplot,t[i_valid],alpha[i_valid,i],color=col[fit_results[i].ifit] xyouts,t[i+1],alpha[i+1,i],string((fit_results[loglog_fits].source_lab)[i], $ format='(i0)'),charthick=2,size=1.2 endfor endif xyouts, !x.crange[0]+0.4*(!x.crange[1]-!x.crange[0]),!y.crange[0]+0.75*(!y.crange[1]-!y.crange[0]),$ 'S(E)=A(E/1keV)!a-!4a!3!n',size=1.4,charthick=1.5 endif if linlog_fits[0] ne -1 then begin beta = reform((fit_results[linlog_fits].coefs)[1,*,*]) valid = where((beta ne 'Inf') and (beta ne '-Inf') and (beta ne 0)) Temp = -1./beta utplot,t,[min(Temp[valid]),max(Temp[valid])],anytim(t,/ex),xtit='',ytit='kT [keV]', $ yst=2,/nodata,charsize=1.2,ymarg=[2,1] for i=0,n_elements(linlog_fits)-1 do begin i_valid = where(Temp[*,i] ne 'Inf') oplot,t[i_valid],Temp[i_valid,i],color=col[fit_results[i].ifit] xyouts,t[i+1],Temp[i+1,i],string((fit_results[linlog_fits].source_lab)[i],format='(i0)'),$ charthick=2,size=1.2 endfor xyouts, !x.crange[0]+0.4*(!x.crange[1]-!x.crange[0]),!y.crange[0]+0.75*(!y.crange[1]-!y.crange[0]),$ 'S(E)=Ae!a-E/kT!n',size=1.4,charthick=1.5 endif Amp = reform((fit_results.coefs)[0,*,*]) valid = where((Amp ne 'Inf') and (Amp ne 0)) utplot,t,[min(Amp[valid]),max(Amp[valid])],anytim(t,/ex),ytit='A [photons s!e-1!ncm!e-2!nasec!e-2!nkeV!e-1!n]',$ yst=2,/nodata,charsize=1.2,ymarg=[3,1] for i=0,N_fits-1 do begin i_valid = where((Amp[*,i] ne 0) and (Amp[*,i] ne 'Inf')) oplot,t[i_valid],Amp[i_valid,i],color=col[fit_results[i].ifit] xyouts,t[i+1],Amp[i+1,i],string((fit_results.source_lab)[i],format='(i0)'), $ size=1.2,charthick=2 endfor endif !p.multi=0 ;endif end ; ************** MAIN PROCEDURE "impsec.pro " ************************************** ;+ ; ; PURPOSE: ; "imspec" is a widget-based user interface which graphically displays ; the position-dependent spectrum of HESSI HXR/Gamma ray iamges. ; The program shows the image in a color code, where hard (soft) ; spectra correspond to blue (red) colors. The color lookup table is ; is defined in terms of RGB filters (see subroutine "spec2rgb"). ; Rectangular regions may be selected interactively, and the average spectrum ; within each reagion is displayed using either linlog or loglog scales, ; to which interactive curve fits for exponentoals/powerlaws can be made. ; Each region allows for up to two curve fits, such as to represent "broken ; powerlaws" or thermal spectra. The number of selectable regions is unlimited. ; If the image cube contains several time bins, the user is allowed to ; step (or flow) through a movie, and compute the temporal evolution of the ; fitted values such as powerlaw exponents or intensities. ; ; INPUT: an image cube of the format fltarr(N_x,N_y,N_E) or (N_x,N_y,N_E,N_t). ; The image cube can be obtained (1) from a HESSI style FITS file ; (e.g., 'hsi_imagecube_flare8_12x10.fits'), or from an ad hoc simulation ; (for testing purposes). ; ; OUTPUT: the primary output devise is the screen, with optional output to a ; PostScript, ASCII and IDL SAVE files. ; ; SOFTWARE/HARDWARE REQUIREMENTS: the program makes use of several SSW hessi solar soft ; routines, and should be started from the same environment as the usual hessi ; software. The color output is designed for TRUE COLOR devices. ; ; SIDE EFFECTS: ; "imspec" starts its own Xmanager, which might interfer with pre-existing Xmanagers ; and the HESSI GUI when ruinning in the same session. The "decomposed" keyword ; (transiently) set to commode with WIN devices. ; ; HISTORY: ; -written for simulated data by KA, Summer 2001. ; -Adapted for real data, March/April 2002. ; -New functionality added on behalf of Brian Dennis: broken powerlaw fits, ; movie capability, KA, April/Mai, 2002. ; -Added IDL SAVE output, and general debugging, KA, June 2002. ; -Improved drawing graphics (PixMap buffer instead of XOR mode) and more general ; integration domains (rectangles, ellipses, simply connected iso-intensity contours) ; added, KA, August 2002. ; -Color (RGB) table improved, KA, Sept. 2002. ;  -Catching up single energy band case, KA, October 2002. ;- pro imspec, fitsfile=fitsfile ;+ ; Creates main level widget and starts xmanager ;- common img_cache, img_,intmap_,img_known ; RGB/intensity image cache (speedup) device,get_decomposed=decomposed_old device,decomposed=1 loadct,0 ; ***** chose input data **** ;if not keyword_set(fitsfile) then filename = 'hsi_imagecube_flare8_12x10.fits' ;if not keyword_set(fitsfile) then filename = 'hsi_imagecube_20020421_011512.fits' if not keyword_set(fitsfile) then begin filename = dialog_pickfile (filter='*.fits',title = 'Select input file',get_path=path) if filename eq path then begin message, 'No input file selected - aborting imspec',/cont return endif fitsfile = filename endif read_image_cube_fits,fitsfile,imc,x,y,energy_bands,time_bands,info imc = imc > 0 ; just a quick-and-dirty solution if (size(imc))(0) eq 4 then Nt = (size(imc))(4) else Nt = 1 ; number of time bins in image cube E = total(energy_bands,1)/2. ; center energies dE = reform(energy_bands[1,*]-energy_bands[0,*]) ; energy band widths define_img_w,Nxd,Nyd,pos1,pos2 ; get image size and position ; ***** definition of the default widget tree contained in the main menu ***** base = widget_base(title="imspec for RHESSI",column=3);,/no_copy) ; /no_copy: pointers passed, not data menu = widget_base(base,/column) button1 = widget_button(menu,value='Help',uvalue={name:'help'}) button2 = widget_button(menu,value='Redraw',uvalue={name:'redraw'}) menu1= widget_base(menu,row=4,frame=3) tit1 = widget_label(menu1,value='Time Selection') menu11 = widget_base(menu1,column=3,event_pro='time_steps') button111 = widget_button(menu11,value='<') tme_lab = widget_label(menu11,value='1/'+string(Nt,format='(i0)'),xsize=74) button112 = widget_button(menu11,value='>') menu12 = widget_base(menu1,column=3,event_pro='movie') button121 = widget_button(menu12,value='<<') button122 = widget_button(menu12,value='STOP') button123 = widget_button(menu12,value='>>') menu2= widget_base(menu,column=1,frame=3) tit2 = widget_label(menu2,value='Image Options') menu21 = widget_base(menu2,column=1) button211 = widget_button(menu21,value='RGB',uvalue={name:'RGB'}) menu22 = widget_base(menu21,column=1) if size(energy_bands,/n_dim) gt 1 then $ slider=widget_slider(menu22,minimum=0,maximum=n_elements(E)-1, $ uvalue={name:'E_slider'},/suppress_value) elab0 = string(min(energy_bands),max(energy_bands),format='(i0,"..",i0," keV")') energy_lab = widget_label(menu22,xsize=70,value=elab0) menu3= widget_base(menu,column=1,frame=3) tit3 = widget_label(menu3,value='Integration Domain') drplist = widget_droplist(menu3,value=['rectangle','ellipse','intensity contour'],uvalue={name:'region_style'}) button8 = widget_button(menu,value='Show Spec. Evolution',uvalue={name:'ShowEvolution'}) button9 = widget_button(menu,value='Show Light Curves',uvalue={name:'ShowLightCurves'}) button10 = widget_button(menu,value='Make PS',uvalue={name:'MakePS'}) button11 = widget_button(menu,value='Make ASCII',uvalue={name:'MakeASCII'}) button12 = widget_button(menu,value='Make IDL SAVE',uvalue={name:'MakeIDLSAVE'}) button13 = widget_button(menu,value='Quit',uvalue={name:'quit'}) widget_control, base, /realize widget_control, menu, /realize widget_control, menu11,set_uvalue={time_label_id:tme_lab} widget_control, menu12,set_uvalue={bwd:button111,fwd:button112,time_label_id:tme_lab} widget_control, menu21,set_uvalue={energy_label_id:energy_lab} widget_control, menu22,set_uvalue={energy_label_id:energy_lab} draw = widget_draw(base, xsize=Nxd,ysize=Nyd, $ uvalue={name:'image',cursor:[-1,-1,-1,-1],pixID:-1}, $ button_events=1,motion_events=0) widget_control, draw, get_value=windex ; get and set draw window index wset, windex data = {imc:imc,x:x,y:y,E:E,dE:dE,imgw_id:windex,src_nr:1, $ Nt:Nt,time_index:0,time_bands:time_bands,info:info, $ movie_flag:0,movie_killed:0,image_style:'RGB', $ region_style:'rectangle',energy_index:0} widget_control, base, set_uvalue=data Nxpixel = pos2[0]-pos1[0] ; initialize image cache Nypixel = pos2[1]-pos1[1] img_known = fltarr(Nt) img_ = fltarr(Nxpixel,Nypixel,3,data.Nt) intmap_ = fltarr(Nxpixel,Nypixel,data.Nt) mk_img_cache,data ; compute initial (t0) rgb image show_img,data,/no_contour xmanager, "imspec",base,event_handler="imspec_event",/no_block device,decomposed=decomposed_old end