function poly4fnc,x,a ; y=a[0]+a[1]*x+a[2]*x^2+a[3]*x^3 yder0=dblarr(n_elements(x))+1.0d0 yder1=x yder2=x^2 yder3=x^3 return,[ [y], [yder0], [yder1], [yder2] ,[yder3] ] end ; ; ; pro polyfit4deg,x,y,params,meas_err=meas_err,winn=winn ; if(!version.os_family eq 'unix') then windev='X' else windev='win' origdev=!d.name tvlct,rcol_orig,gcol_orig,bcol_orig,/get set_plot,windev ; rcol=[0.,1.,0.,1.]*255 gcol=[0.,1.,1.,0.]*255 bcol=[0.,1.,0.,0.]*255 tvlct,rcol,gcol,bcol ; device,dec=0 if(n_elements(meas_err) eq 0) then meas_err=sqrt(y) if(n_elements(winn) gt 0) then begin window,winn endif else begin window,/free winn=!window endelse a=[1.d0,1.d0,1.d0,1.d0] ftng_agn_ln: plot,x,y,/nodata oplot,x,y,psym=1,color=3 tmpstr=' ' print,'' print,'' print,'function: a0+a1*x+a2*x^2+a3*x^3' print,'a0=',a[0],'(leave blank for no change)' read,tmpstr if(strlen(tmpstr) gt 0) then a[0]=double(tmpstr) print,'' print,'a1=',a[1],'(leave blank for no change)' read,tmpstr if(strlen(tmpstr) gt 0) then a[1]=double(tmpstr) print,' ' print,'a2=',a[2],'(leave blank for no change)' read,tmpstr if(strlen(tmpstr) gt 0) then a[2]=double(tmpstr) print,' ' print,'a3=',a[3],'(leave blank for no change)' read,tmpstr if(strlen(tmpstr) gt 0) then a[3]=double(tmpstr) print,' ' print,' ' print,'Choose whether a parameter is taken in the fitting' print,' Otherwise it is kept constant' print,' ' pmtft=[0,0,0,0] tmpstr=' ' print,'Take a[0] into fitting? (Y/n)' read,tmpstr tmpstr=strupcase(tmpstr) if(tmpstr ne 'N') then pmtft[0]=1 print,' ' tmpstr=' ' print,'Take a[1] into fitting? (Y/n)' read,tmpstr tmpstr=strupcase(tmpstr) if(tmpstr ne 'N') then pmtft[1]=1 print,' ' tmpstr=' ' print,'Take a[2] into fitting? (Y/n)' read,tmpstr tmpstr=strupcase(tmpstr) if(tmpstr ne 'N') then pmtft[2]=1 print,' ' tmpstr=' ' print,'Take a[3] into fitting? (Y/n)' read,tmpstr tmpstr=strupcase(tmpstr) if(tmpstr ne 'N') then pmtft[3]=1 print,' ' fitxvec=lmfit(x,y,a,chisq=ch_sq,/double,fita=pmtft,measure_errors=meas_err,$ function_name='poly4fnc',sigma=sgm_a) print,'function: a0+a1*x+a2*x^2+a3*x^3' for i=0,3 do begin aistr='a['+strcompress(string(i),/remove_all)+']='+strcompress(string(a[i]),/remove_all) if(pmtft[i] eq 1) then print,aistr,' +/-',abs(sgm_a[i]/a[i]*100.),' %' $ else print,aistr,' +/- 0.0 %' endfor xxstep=(max(x)-min(x))/100. xxx=findgen(101)*xxstep+min(x) fitvec=a[0]+a[1]*xxx+a[2]*xxx^2+a[3]*xxx^3 window,winn plot,x,y,/nodata oplot,x,y,psym=1,color=3 oplot,xxx,fitvec,color=2 print,'' print,'--------------------------------' print,'Continue with fitting? (Y/n)' print,' ' sw_cft=' ' read,sw_cft sw_cft=strupcase(sw_cft) if(sw_cft ne 'N') then goto,ftng_agn_ln ; params=a tvlct,rcol_orig,gcol_orig,bcol_orig set_plot,origdev end