pro xgaussfit_pd_event, event common xgaussfit_block, widget, xg par=xg.par xtitle=xg.xtitle ytitle=xg.ytitle case event.value of 'Done' : begin widget_control, event.top, /destroy return end 'Save': begin file=pickfile(file='xgaussfit.res', path=!sumer_config.work_dir,/write) save, xg, filename=file return end 'Display' : begin widget_control, widget.text, $ set_value='Displaying the profiles with specfied parameters...' par1=par nline=n_elements(par)/3-1 plot, xg.x, xg.y, pos=xg.pos_dev, /dev, $ xr=xg.xr, xst=1, yr=xg.yr, yst=1, psym=10, xtitle=xtitle, ytitle=ytitle x1=findgen(101)*(xg.xr(1)-xg.xr(0))/100.+xg.xr(0) y1=mgaussfit(x1, x1, x1*0+1, par1, /value) oplot, x1, y1, thick=2 par1(0:2)=0. bg=y1-mgaussfit(x1, x1, x1*0+1, par1, /value) for l=1, nline do begin par1=[par(0:2)*0, par(l*3), par(l*3+1), par(l*3+2)] y1=mgaussfit(x1, x1, x1+1, par1, /value)+bg oplot, x1, y1, linest=1, thick=2 endfor return end 'Hardcopy' : begin xps_setup, ps_setup ps_setup.printer = 'lpr -h' cd, current=current defsysv, '!sumer_config', exist=exist if exist eq 1 then tmp=execute('dir = !sumer_config.work_dir') $ else dir=current ps_setup.filename= 'idl.ps' ps_setup.hard=0 ps_setup.portrait=1 xps_setup, ps_setup, status=accept if not accept then return set_plot, 'ps' device, filename=concat_dir(dir,ps_setup.filename), $ portrait=ps_setup.portrait, $ encapsulated = ps_setup.encapsulated, color=ps_setup.color device, /times, /bold, /isolatin1 case ps_setup.portrait of 1: begin xsize=18 factor=xsize /float(xg.d_win.x_size) ysize=factor*xg.d_win.y_size xoff=(21.5-xsize-1)/2. yoff=(28.-ysize-1)/2. end 0 : begin xsize=18 factor=xsize /float(xg.d_win.x_size) ysize=factor*xg.d_win.y_size xoff=(21.5-ysize-1)/2. yoff=28.-(28.-xsize-1)/2. end endcase device, set_character_size=xg.d_win.x_ch_size*factor*1000*[1, 12./8.] device, xoff=xoff, yoff=yoff, xs=xsize, ys=ysize par1=par nline=n_elements(par)/3-1 nline=n_elements(par)/3-1 x1=findgen(101)*(xg.xr(1)-xg.xr(0))/100.+xg.xr(0) y1=mgaussfit(x1, x1, x1*0+1, par1, /value) par1(0:2)=par1(0:2)*0 bg = y1-mgaussfit(x1, x1, x1*0+1, par1, /value) plot, xg.x, xg.y, pos=xg.pos, /norm, $ xr=xg.xr, xst=1, yr=xg.yr, yst=1, psym=10, thick=4, $ xtitle=xg.xtitle, ytitle=xg.ytitle, font=0, xthick=4, ythick=4 oplot, x1, y1, thick=4 for l=1, nline do begin par1=[par(0:2)*0, par(l*3), par(l*3+1), par(l*3+2)] y1=mgaussfit(x1, x1, x1*0+1, par1, /value)+bg oplot, x1, y1, linest=1, thick=4 endfor if not ps_setup.encapsulated then begin xyouts, xsize*0.5*1000, -0.5*1000, /dev, font=0, $ align=0.5, charsize=0.8,$ 'Fitting Results: '+ $ 'Background level = '+strtrim(string(par(0), format='(g12.4)'),2) + $ '; Background slope = '+strtrim(string(par(1),format='(g12.4)'),2) for l=1, nline do begin xyouts, xsize*0.5*1000, -(0.5+0.5*l)*1000, /dev, font=0, align=0.5, $ charsize=0.8, $ 'Line '+strtrim(string(l), 2)+': '+ $ 'Peak Intensity = '+strtrim(string(par(l*3), format='(g14.7)'),2) +$ '; Center Position = '+strtrim(string(par(l*3+1), format='(g14.7)'),2)+$ '; Standard width = '+strtrim(string(par(l*3+2), format='(g14.7)'),2) endfor endif device, /close set_plot, 'x' if ps_setup.hard then spawn, $ ps_setup.printer+ ' '+ concat_dir(dir, ps_setup.filename) end 'Fit' : begin widget_control, widget.text, $ set_value='Fitting and Displaying the profiles with fitted parameters...' par1=par nline=n_elements(par)/3-1 free = [0] if xg.slope_free then free=[free, 1] for l=1, nline do free = [free, l*3, l*3+1, l*3+2] y1=mgaussfit(xg.x, xg.y, xg.y*0+1, par1, free_par=free, /value) w = 1./(xg.noise(0)+xg.noise(1)*(y1>0.)+xg.noise(2)*y1^2) npar=n_elements(par1) weight = fltarr(npar) weight(0) =1./(xg.noise(0)+xg.noise(1)* abs(par1(0))+xg.noise(2)*par1(0)^2)   weight(1)=weight(0) for l=1, nline do begin weight(l*3-1) = 1./0.5^2 weight(l*3) = 1./(par1(3*l+2)*0.5)^2 weight(l*3+1) = 1./ (par1(3*l+2)*0.25)^2 endfor y1=mgaussfit(xg.x, xg.y, w, par1, free_par=free, $ expected=par1(free), weight=weight(free)) par=par1 xg.par=par for k=0, n_elements(free)-1 do $ widget_control, widget.field(free(k)), $ set_value=strtrim(string(xg.par(free(k)), format='(g14.7)'),2) widget_control, widget.draw, get_value=id wset, id plot, xg.x, xg.y, pos=xg.pos_dev, /dev, $ xr=xg.xr, xst=1, yr=xg.yr, yst=1, psym=10, $ xtitle=xtitle, ytitle=ytitle x1=findgen(101)*(xg.xr(1)-xg.xr(0))/100.+xg.xr(0) y1=mgaussfit(x1, x1, x1*0+1, par1, /value) oplot, x1, y1, thick=2 par1(0:2)=0. bg=y1-mgaussfit(x1, x1, x1*0+1, par1, /value) for l=1, nline do begin par1=[par(0:2)*0, par(l*3), par(l*3+1), par(l*3+2)] y1=mgaussfit(x1, x1, x1+1, par1, /value)+bg oplot, x1, y1, linest=1, thick=2 endfor return end endcase end pro xgaussfit_event, event common xgaussfit_block, widget, xg ;par=xg.par ;print, event.id xtitle=xg.xtitle ytitle=xg.ytitle if event.id eq widget.pdmenu then xgaussfit_pd_event, event s=where(widget.field eq event.id, count) if count eq 1 then begin widget_control, widget.field(s), get_value=value xg.par(s) = float(value) return endif s=where(widget.noise eq event.id, count) if count eq 1 then begin s=s(0) widget_control, widget.noise(s), get_value=value xg.noise(s) = value return endif s=where(widget.pick eq event.id, count) if count eq 1 then begin s=s(0) widget_control, widget.draw, get_value=id wset, id plot, xg.x, xg.y, pos=xg.pos_dev, /dev, $ xr=xg.xr, xst=1, yr=xg.yr, yst=1, xtitle=xtitle, ytitle=ytitle case s of 0 : begin widget_control, widget.text, $ set_value='Select Two points on the background level!' xg.npoint = 2 end 1 : begin widget_control, widget.text, $ set_value='Select 2 points on the background level ' xg.npoint = 2 end 2 : begin widget_control, event.id, get_value=ss xg.slope_free=ss end else : begin case s mod 3 of 0 : begin widget_control, widget.text, $ set_value= $ 'Select 1 point on the base level and 1 point on the peak.' xg.npoint = 2 end 1 : begin widget_control, widget.text, $ set_value='Select 1 point on the center of the line!' xg.npoint=1 end 2 : begin widget_control, widget.text, $ set_value='Select 2 points on the the half maximum!' xg.npoint=2 end endcase end endcase xg.button=s return endif if event.id eq widget.draw then begin ; print, 'draw' if xg.npoint gt 0 and event.press gt 0 then begin xmark = event.x ymark = event.y x = (xmark-xg.pos_dev(0))/float(xg.pos_dev(2)-xg.pos_dev(0)) $ *(xg.xr(1)-xg.xr(0))+xg.xr(0) y = (ymark-xg.pos_dev(1))/float(xg.pos_dev(3)-xg.pos_dev(1)) $ *(xg.yr(1)-xg.yr(0))+xg.yr(0) ; print, 'x, y=', x, y widget_control, widget.draw, get_value=id wset, id mark_on_image, xmark, ymark, color=!d.table_size-1, /mark_only xg.point(*, xg.npoint-1)=[x,y] xg.npoint = xg.npoint-1 endif if xg.npoint eq 0 then begin if xg.button eq 0 then begin xg.par(0) = 0.5*(xg.point(1,0)+xg.point(1,1)) ; widget_control, widget.field(0), $ set_value=strtrim(string(xg.par(0,0), format='(g14.7)'),2) widget_control, widget.pick(3), sensitive=1 endif if xg.button eq 1 then begin xg.par(1) = $ (xg.point(1,0)-xg.point(1,1))/(xg.point(0,0)-xg.point(0,1))*xg.par(5) widget_control, widget.field(1), $ set_value=strtrim(string(xg.par(1,0), format='(g14.7)'),2) endif if xg.button mod 3 eq 0 and xg.button gt 2 then begin xg.par(xg.button) = xg.point(1,0)-xg.point(1,1) xg.par(xg.button+1) = xg.point(0,1) widget_control, widget.field(xg.button), $ set_value=strtrim(string(xg.par(xg.button), format='(g14.7)'),2) widget_control, widget.field(xg.button+1), $ set_value=strtrim(string(xg.par(xg.button+1), format='(g14.7)'),2) for k=1, 2 do $ widget_control, $ widget.pick((xg.button+k)<( n_elements(widget.pick)-1)), sensitive=1 if xg.button eq n_elements(widget.pick)-1 then widget_control, $ widget.pdmenu, sensitive=1 endif if xg.button mod 3 eq 1 and xg.button gt 2 then begin xg.par(xg.button) = xg.point(0,0) widget_control, widget.field(xg.button), $ set_value=strtrim(string(xg.par(xg.button), format='(g14.7)'),2) endif if xg.button mod 3 eq 2 and xg.button gt 2 then begin xg.par(xg.button) = 0.4246*abs(xg.point(0,0)-xg.point(0,1)) widget_control, widget.field(xg.button), $ set_value=strtrim(string(xg.par(xg.button), format='(g14.7)'),2) if xg.button eq 5 then widget_control, widget.pick(1), sensitive=1 if xg.button lt n_elements(widget.pick)-1 then $ widget_control, widget.pick(xg.button +1), sensitive=1 endif endif return endif end ;+ ; NAME : ; XGAUSSFIT ; PURPOSE : ; Interactive Gaussian Fitting ; CALLING SEQUENCE : ; ; XGAUSSFIT, X, Y, Par, Ncomp=Ncomp, Noise=Noise ; ; INPUTS: ; ; X independent variables ; Y dependent variables ; ; OUTPUT : ; ; PAR fitted gaussian parameters ; ; Y = par(0) + par(1)*X + par(2)*X^2 ; + par(3)*exp(-0.5*((X-par(4))/par(5))^2) ; + par(6)*exp(-0.5*((X-par(7))/par(8))^2) + ... ; under the positive constraints, i.e. par(3)>0, par(6)>, ... ; ; KEYWORDS ; ; XTITLE x-axis title ; YTITLE y-axis title ; NCOMP number of gaussian components (default=1) ; NOISE a 3-element array defining noise characteristics ; (default=[1,0 0]) ; ; RMS NOISE = sqrt(NOISE(0)+NOISE(1)*SIGNAL ; +NOISE*SIGNAL^2) ; HISTORY: ; 1998 December 2 J. Chae updated hardcopy part. ;- pro xgaussfit, x, y, par1, ncomp=ncomp, noise=noise_arr, $ xtitle=xtitle, ytitle=ytitle, slope_free=slope_free common xgaussfit_block, widget,xg common mark_on_image_block, mark if n_elements(slope_free) eq 0 then slope_free=0 if n_elements(ncomp) eq 0 then ncomp=1 if n_elements(mark) ne 0 then mark_old=mark device, set_character_size=[8, 12] base=widget_base(title='Tool For Interactive Gaussian Fitting', col=1) junk={CW_PDMENU_S, flags:0, name:''} menu = [$ { CW_PDMENU_S, 0, 'Done'}, $ { CW_PDMENU_S, 0, 'Display'}, $ { CW_PDMENU_S, 0, 'Hardcopy'}, $ { CW_PDMENU_S, 0, 'Save'}, $ { CW_PDMENU_S, 2, 'Fit'} $ ] wd_xsize=600 wd_ysize=450 base1 = widget_base(base, col=2) pdmenu = cw_pdmenu(base1, menu, /return_full_name) base2 = widget_base(base1, col=3) noise=lonarr(3) if n_elements(noise_arr) eq 0 then noise_arr=[1., 0,0] ;print, noise_arr noise(0) = cw_field(base2, title=' NOISE^2 = ', value=noise_arr(0), $ xsize=10, /all_events , float=1) noise(1) = cw_field(base2, title=' + SIGNAL x ', value=noise_arr(1), $ xsize=10, /all_events, float=1) noise(2) =cw_field(base2,title='+ SIGNAL^2 x ', value=noise_arr(2), $ xsize=10, /all_events, float=1) text = widget_text(base, xsize=60, ysize=1) draw = widget_draw(base, xsize=wd_xsize, ysize=wd_ysize, button_events=1) par_base = widget_base(base, col=6 , map=1, frame=1) field=lonarr(3,ncomp+1) par=fltarr(3*ncomp+3) if ncomp eq 0 then par=[par, 1.0e-10, 1., 1.0] pick = lonarr(3, ncomp+1) initial=bytarr(3, ncomp+1) pick(0,0) = widget_button(par_base, value='Background level') field(0,0) = widget_text(par_base, value='0.', xsize=14, /edit, /all_events) pick(1,0) =widget_button(par_base, value='Slope') field(1,0) = widget_text(par_base, value='0.',xsize=14, /edit, /all_events) pick(2,0) =cw_bgroup(par_base, ['fixed slope', 'free slope'], $ exclusive=1, row=1, set_value=slope_free) ;field(2,0) = widget_text(par_base, value='0.', xsize=14, /edit, /all_events) if ncomp ge 1 then for nline=1, ncomp do begin line = 'Line '+string(nline, format='(I1)') par_base = widget_base(base, col=6 , map=1, frame=1) pick(0,nline) = widget_button(par_base, value=line+' Peak') field(0,nline) =widget_text(par_base, value='0.', xsize=14, /edit, /all_events) pick(1,nline) = widget_button(par_base, value=' Center') field(1, nline)= widget_text(par_base, value='0.', xsize=14, /edit, /all_events) pick(2,nline) = widget_button(par_base, value=' Standard Width') field(2, nline) = widget_text(par_base, value='0.', xsize=14, /edit, /all_events) endfor widget_control, base, /realize widget = {pdmenu:pdmenu,draw:draw, text:text, par_base:par_base, $ field:field, pick:pick, noise:noise} widget_control, draw, get_value=draw_id wset, draw_id pos=[ 0.15, 0.15, 0.95, 0.95] pos_dev = [pos(0)*wd_xsize, pos(1)*wd_ysize, $ pos(2)*wd_xsize, pos(3)*wd_ysize] xr=[min(x, max=m), m] ymax=max(y, min=m) yr=[m-0.5*(ymax-m), ymax+0.5*(ymax-m)] if n_elements(xtitle) eq 0 then xtitle='' if n_elements(ytitle) eq 0 then ytitle='' xg = {x:x, y:y, pos:pos, pos_dev:pos_dev, $ xr:xr, yr:yr, par:par, npoint:0, button:-1,$ initial :initial, point:fltarr(2,3), noise:noise_arr, $ d_win:!d, xtitle:xtitle, ytitle:ytitle, slope_free:0 , $ title:''} plot, x, y, pos=pos, /norm, $ xr=xr, xst=1, yr=yr, yst=1, xtitle=xtitle, ytitle=ytitle ;oplot, x, y widget_control, text, set_value='Welcome to Tool for Gaussian Fitting' widget_control, widget.pdmenu, sensitive=1 for k=1, n_elements(widget.pick)-1 do $ widget_control, widget.pick(k), sensitive=0 widget_control, widget.pick(1,0), sensitive=1 widget_control, widget.pick(2,0), sensitive=1 device, set_character_size=[8, 12] xmanager, 'xgaussfit', base, /modal par1=xg.par noise_arr=xg.noise if n_elements(mark_old) ne 0 then mark=mark_old slope_free=xg.slope_free end