!------------------------------------------------------------------------- ! This is the latest version of the BRAGG script ! A "better" flat top is achieved by carefully choosing the fitting range, ! which is now done automatically. !------------------------------------------------------------------------- defaults clear datafile = `c:\extrema\Scripts\bragg.dat' read\-messages datafile xdata ydata set tension 0.001 npts = 100 generate x1 xdata[1],,xdata[#] npts y1 = smooth(xdata,ydata,x1) statistics\-messages y1 ylev\max iymax\imax xupr = x1[iymax] ! window 5 set plotsymbol -1 set %xnumberheight 5 set %ynumberheight 5 set %xloweraxis 20 set ylabel 'Normalized Dose' set ylabelon 1 set %ylabelheight 5 set xlabel 'Depth (mm)' set xlabelon 1 set %xlabelheight 5 graph xdata ydata set plotsymbol 0 graph\overlay x1 y1 set textcolor blue set textalign 1 set textinteractive 0 set %xtextlocation 24 set %ytextlocation 80 set %textheight 8 text 'raw Bragg peak' ! dx = 1 nstep = 5 inquire 'shift amount ( '//rchar(dx)//' ) & #functions ( '//rchar(nstep)//' ) >> ' dx nstep destroy m m[1:npts,1] = y1 ! destroy a1 scalar\fit 'a1' a1 = 1 params = '[a1;' do i = [2:nstep] m[1:npts,i] = step(y1,-abs(dx)/(x1[2]-x1[1])*(i-1)) destroy\expand 'a'//rchar(i) scalar\fit 'a'//rchar(i) 'a'//rchar(i)=1 params = params//'a'//rchar(i)//';' enddo paramLen = 1 + 3*min(nstep,9) + 4*(nstep>9)*(nstep-9) params[paramLen:paramLen] = ']' xlow = xupr-(nstep-1)*abs(dx) ! thresh = 0.01 ! lower limit of fit parameters itmax = 10 ! max. # of fit iterations yf[1:npts] = ylev fit_again: fit\weight\itmax\-messages 10*(x1>xlow)*(x1<xupr) itmax yf=m<>params pmin=1.e30 ipmin=0 do i = [1:nstep] temp = 'a'//rchar(i) value = evaluate(temp) if ( value<pmin & value~=0 ) then ipmin = i pmin = value endif enddo if ( pmin>=thresh ) then goto done_fit destroy\expand 'a'//rchar(ipmin) scalar 'a'//rchar(ipmin) txtemp='a'//rchar(ipmin) txtemp=0 goto fit_again done_fit: window 7 set %xnumberheight 5 set %ynumberheight 5 set %xloweraxis 15 set ylabel 'Normalized Dose' set ylabelon 1 set %ylabelheight 5 set xlabel 'Depth (mm)' set xlabelon 1 set %xlabelheight 5 graph x1 m<>params set textcolor blue set textalign 1 set textinteractive 0 set %xtextlocation 18 set %ytextlocation 20 set %textheight 8 text 'spread out Bragg peak' ! window 4 set %xnumberheight 5 set %ynumberheight 5 graph x1 m[*,1]*a1 vector output nstep output[1] = a1 do i = [2:nstep] temp = 'a'//rchar(i) graph\overlay x1 m[*,i]*temp output[i] = evaluate(temp) enddo set %textheight 8 get %xloweraxis xlaxis get %yupperaxis yuaxis get %textheight txthit set %xtextlocation xlaxis+3 set %ytextlocation yuaxis-txthit set textalign 1 set textinteractive 0 text 'shift = '//rchar(dx) set %ytextlocation yuaxis-2.5*txthit text '#functions = '//rchar(nstep) defaults window 0