PRO ODR_OLD,dev=dev,$ x=x,$ ;optional array of x data to be passed in y=y,$ ;optional array of y data to be passed in file=file,$ ;optional file of x,y data noplot=noplot,$ ;cancels plot (when calling from outside) coeff,$ ;fitting coefficients and associated errors passed out yfit1,$ ;vector of y-values to be used in plotting odr line yfit2,$ ;vector of y-values to be used in plotting lsq line text1,$ ;text string for equation of odr line text2 ;text string for equation of lsq line ; ;Procedure to use Orthogonal Distance Regression (from Num. Recipes) ;to fit lines to any x and y. ;The procedure estimates x and y uncertainty in the data due to ;natural variability by relating errors to chi-sq probabilities. ;These errors allow for a realistic estimation of slope and intercept ;uncertainties. ; ;P Tans date=? ;modified by JBM 11/6/00,4/9/01 ; ; ;Read in text file containing x and y data to be fit ;and create separate arrays, only if x,y data (or ;file) not passed from outside IF NOT KEYWORD_SET(x) OR NOT KEYWORD_SET(y) THEN BEGIN IF NOT KEYWORD_SET(file) THEN BEGIN file='~john/idl/odr/kp_test2.txt' ;file='~john/idl/odr/data.txt' ENDIF CCG_FREAD,file=file,/nomessages,nc=2,xydata xdata=xydata(0,*) ydata=xydata(1,*) ENDIF ELSE BEGIN xdata=x ydata=y ENDELSE ;revert to least squares fit if only 2 data points ;and process output arrays IF N_ELEMENTS(xdata) EQ 1 THEN BEGIN coeff=FLTARR(2,4) yfit1=-999.9 yfit2=-999.9 text1='-999' text2='-999' ENDIF ELSE IF N_ELEMENTS(xdata) LT 3 THEN BEGIN result=POLY_FIT(xdata,ydata,1,yfit) ;construct return arrays coeff=FLTARR(2,4) coeff(0,0)=result(0) ;lsq int coeff(1,0)=-999.9 ;lsq int err coeff(0,1)=result(1) ;lsq slope coeff(1,1)=-999.9 ;lsq slope err coeff(0,2)=-999.9 ;lsq chi2 coeff(1,2)=-999.9 ;lsq chi2prob coeff(0,3)=result(0) ;LSq intercept coeff(1,3)=result(1) ;LSq slope yfit1=yfit yfit2=yfit ;full line equations text1=STRCOMPRESS('LSQ: y='+STRING(FORMAT='(F9.1)',result(1))+$ '('+STRING(FORMAT='(F7.2)',coeff(1,1))+')'+'x+'+ $ STRING(FORMAT='(F7.2)',result(0)),/RE) text2=STRCOMPRESS('LSQ: y='+STRING(FORMAT='(F9.1)',coeff(1,3))+$ 'x+'+STRING(FORMAT='(F7.2)',coeff(0,3)),/RE) ;just slopes and error text1=STRCOMPRESS(STRING(FORMAT='(F9.1)',result(1))+$ '('+STRING(FORMAT='(F7.2)',coeff(1,1))+')',/RE) text2=STRCOMPRESS(STRING(FORMAT='(F9.1)',coeff(1,3)),/RE) ENDIF ELSE BEGIN ;fitting parameters ;polynomial degree of fit n=1 ;number of fixed parameters p=0 ;outlier criterion in sigma factors s=3. ;Estimate initial natural variability uncertainty for x and y data sigma_x=0.1 ;for CO2 in ppm sigma_y=0.01 ;for d13C ;sigma_y=SQRT((0.01*360.)^2+(0.1*8.)^2) ;for d13C*CO2 ;Fit initial straight line to CO2 vs C13 (or any x,y) jj=INDGEN(N_ELEMENTS(xdata)) time=INDGEN(N_ELEMENTS(xdata)) fitcf=FINDGEN(2) ;linear fits (Num. Rec.) of C13 vs CO2 etc. errcf=FINDGEN(2) ;ydata=ydata*1000.D ;xd=xdata ;yd=ydata ;xdata=yd ;ydata=xd ;sx=sigma_x ;sy=sigma_y ;sigma_x=sy ;sigma_y=sx IF N_ELEMENTS(jj) GT 2 THEN BEGIN ;write file read by 'fitexy' code ;line 1:N_ELEMENTS, order of fit, sigma_x,sigma_y, # of fixed params ;rest of lines: x and y data CCG_FWRITE,file='~john/idl/odr/xy.dat', /nomessages,$ nc=5,N_ELEMENTS(jj),n,sigma_x,sigma_y,p CCG_FWRITE,file='~john/idl/odr/xy.dat', /nomessages,$ nc=3,/APPEND,time,xdata,ydata SPAWN,'~john/idl/odr/fitexy' CCG_FREAD,file='~john/idl/odr/results.NR', /nomessages,$ nc=2,result1 fitcf(0)=result1(0,0) fitcf(1)=result1(0,1) errcf(0)=result1(1,0) errcf(1)=result1(1,1) chisquared=result1(0,2) chi2prob=result1(1,2) int=fitcf(0) slope=fitcf(1) PRINT, ' slope=',fitcf(1),' +-',errcf(1) PRINT, ' chisquared=',chisquared,' chi2prob=',chi2prob PRINT, ' intercept=',fitcf(0),' +-',errcf(0) ;stop ;Iterative loop: Change sigma_x and sigma_y based on the previous ;fit (residuals in x- and in y-direction), and re-fit. If this ;iteration converges we end up with a straight line whereby the ;squared residual deviations perpendicular to the line have been ;minimized. ;Note: in that case the ratio of sigma_y to sigma_x equals the slope. ;Reject outliers if present, and re-fit curve. ;Select accepted data ;slope must converge to within 0.01 per cent iter=1 ;re-calculate sigmas sigma_x=STDEV(xdata-ydata/slope+fitcf(0)/slope) sigma_y=STDEV(ydata-xdata*slope-fitcf(0)) change=1. ;Initial value - will be relative change in slope(or int) between iterations WHILE iter LT 10 AND ABS(change) GT 0.0001 DO BEGIN ;'s'Sigma filter for outliers ;Retained data ja=WHERE(ABS(ydata-int-slope*xdata) LT s*sigma_y) ;Rejected data jr=WHERE(ABS(ydata-int-slope*xdata) GE s*sigma_y) IF N_ELEMENTS(ja) GT 2 THEN BEGIN CCG_FWRITE,file='~john/idl/odr/xy.dat',/nomessages,$ nc=5,N_ELEMENTS(jj),n,sigma_x,sigma_y,p CCG_FWRITE,file='~john/idl/odr/xy.dat',/nomessages,$ nc=3,/append,time,xdata,ydata SPAWN,'~john/idl/odr/fitexy' CCG_FREAD,file='~john/idl/odr/results.NR', /nomessages,$ nc=2,result1 fitcf(0)=result1(0,0) fitcf(1)=result1(0,1) errcf(0)=result1(1,0) errcf(1)=result1(1,1) chisquared=result1(0,2) chi2prob=result1(1,2) PRINT, iter,' slope=',fitcf(1),' +-',errcf(1) PRINT, ' chisquared=',chisquared,' chi2prob=',chi2prob PRINT, ' intercept=',fitcf(0),' +-',errcf(0) newslope=fitcf(1) ;output value newint=fitcf(0) ;output value error=errcf(1) ;output value change=(newslope-slope)/slope ;change=(newint-int)/int iter=iter+1 slope=newslope int=newint sigma_x=STDEV(xdata(ja)-ydata(ja)/slope+fitcf(0)/slope) sigma_y=STDEV(ydata(ja)-xdata(ja)*slope-fitcf(0)) ENDIF ELSE BEGIN ;wrt 'IF N_ELEMENTS(ja) gt 2' newslope=-99.99 error=-9.99 change=(newslope-slope)/slope iter=iter+1 slope=newslope sigma_x=STDEV(xdata(jr)-ydata(jr)/slope+fitcf(0)/slope) sigma_y=STDEV(ydata(jr)-xdata(jr)*slope-fitcf(0)) ENDELSE ENDWHILE IF iter GE 19 THEN BEGIN slope=-99.99 error=-9.99 ENDIF ;stop ;this iterative loop keeps the ratio of sigma-x and sigma-y constant from ;the previous loop ;diff is the difference between calculated chi squared probability and 0.5 ;(50% chi2prob). The optimization of the error bar size is set to stop when ;diff is less than 0.05 (0.45