PRO 	odr_kp,dev=dev,$
	x=x,$		;array of [CO2] data to be passed in
	y=y,$		;array of d13C data to be passed in	
	file=file,$	;optional file of x,y (co2,d13C) data
	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 of odr line equation with error
	text2		;text string of lsq line equation
;
;Procedure to use Orthogonal Distance Regression (from Num. Recipes)
;to fit lines to 13C v 1/CO2, given 13C and CO2 data -- Keeling Plot
;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 ??
;modified by JBM 11/3/00
;		 01/02/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_test.txt'
	ENDIF
	
	CCG_FREAD,file=file, nc=2, xydata
	xdata=1./xydata(0,*)
	ydata=xydata(1,*)
ENDIF ELSE BEGIN
	xdata=1./x
	ydata=y
ENDELSE

;fitting parameters
	;polynomial degree of fit
	n=1
	;number of fixed parameters
	p=0
	;outlier criterion in sigma factors
	s=400

	;Estimate initial natural variability uncertainty for x and y data
	sigma_x=0.00003	;for 1/CO2 in ppm
	sigma_y=0.05	;for d13C in per mil

;Fit initial straight line to C13 vs 1/CO2 


jj=INDGEN(N_ELEMENTS(xdata))       	
time=INDGEN(N_ELEMENTS(xdata))
fitcf=FINDGEN(2)	;linear fits (Num. Rec.) of C13 vs 1/CO2
errcf=FINDGEN(2)

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='~/idl/odr/xy.dat',nc=5,N_ELEMENTS(jj),$
			n,sigma_x,sigma_y,p
  		CCG_FWRITE,file='~/idl/odr/xy.dat',nc=3,/APPEND,$
  			time,xdata,ydata
  		SPAWN,'~/idl/odr/fitexy'
  		CCG_FREAD,file='~/idl/odr/results.NR',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, '   y-intercept=',fitcf(0),'  +-',errcf(0)
  		PRINT, '   slope=',fitcf(1),'  +-',errcf(1)
  		PRINT, '   chisquared=',chisquared,' chi2prob=',chi2prob
;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
  
iter=1
;re-calculate sigmas
sigma_x=STDEV(xdata-ydata/slope+fitcf(0)/slope) 
sigma_y=STDEV(ydata-xdata*slope-fitcf(0))
;stop
change=1.     ;Initial value - will be relative change in slope between iter.
diff=1.

	;WHILE iter LT 20 AND (ABS(change) GT 0.0005 OR ABS(diff) GT 0.05) DO BEGIN 


        	;'s'Sigma filter for outliers    
   		;Retained data
   		ja=WHERE(ABS(ydata-fitcf(0)-slope*xdata) LT s*sigma_y)

      		;Rejected data
   		jr=WHERE(ABS(ydata-fitcf(0)-slope*xdata) GE s*sigma_y)

    		IF N_ELEMENTS(ja) GT 2 THEN BEGIN 
      			CCG_FWRITE,file='~/idl/odr/xy.dat',nc=5,N_ELEMENTS(jj),$
      				n,sigma_x,sigma_y,p
      			CCG_FWRITE,file='~/idl/odr/xy.dat',nc=3,/append,$
              			time,xdata,ydata
                
      			
      			SPAWN,'~/idl/odr/fitexy'
      		
      			CCG_FREAD,file='~/idl/odr/results.NR',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,'   y-intercept=',fitcf(0),'    +-',errcf(0)
      			PRINT, '   slope=',fitcf(1),'    +-',errcf(1)
      			PRINT, '   chisquared=',chisquared,' chi2prob=',chi2prob
      			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
      			diff = chi2prob - 0.5
			sigma_x=sigma_x/(1+0.5*diff)
			sigma_y=sigma_y/(1+0.5*diff)
      			;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(ja)-ydata(ja)/slope+fitcf(0)/slope)
      			sigma_y=STDEV(ydata(ja)-xdata(ja)*slope-fitcf(0))
    		ENDELSE
 
	;ENDWHILE
  
	IF iter GE 19 THEN BEGIN
    		slope=-99.99
   		error=-9.99
	ENDIF

ENDIF         ;wrt 'IF N_ELEMENTS(jj) gt 2'

;calculate least squares fit coefficients
result2=POLY_FIT(xdata,ydata,1,yfit=yfit,sigma=sigma)

;interr=CCG_MEAN(xdata)*errcf(1)
interr=errcf[0]
interr2=sigma[0]
			
;
;skip plot if program is called from another program
;

IF NOT KEYWORD_SET(x) OR NOT KEYWORD_SET(y) THEN BEGIN

	;Create plot
	IF N_ELEMENTS(jj) GT 2 AND N_ELEMENTS(ja) GT 2 THEN BEGIN
	CCG_OPENDEV, dev=dev, pen=pen,portrait=0

	xtext='1/[CO!D2!N]'
  	ytext='!4d!3!U13!NC'

	IF NOT KEYWORD_SET(x) OR NOT KEYWORD_SET(y) THEN BEGIN
		titletext=file
	ENDIF ELSE BEGIN
		titletext = 'keeling plot'		
	ENDELSE

	PLOT, xdata, ydata, PSYM=4,/nodata, $
  		title=titletext, $
  		xtitle=xtext, $
  		ytitle=ytext, $
  		charsize=1.3
	OPLOT, xdata(ja), ydata(ja), PSYM=4

	;Overplot rejected data IF any
	IF jr(0) NE -1 THEN OPLOT,xdata(jr),ydata(jr), PSYM=2

		;Overplot the fitted line
		np=50   ;# of points in line
		xline=MIN(xdata(ja)) + FINDGEN(np)*(MAX(xdata(ja))-MIN(xdata(ja)))/(np-1)
		yline=fitcf(0)+fitcf(1)*xline 
		OPLOT,xline,yline,THICK=2

		OPLOT,[MIN(xdata(jj)),MAX(xdata(jj))],[0,0]          ;plot axes through (0,0)
		OPLOT,[0,0],[MIN(ydata(jj)),MAX(ydata(jj))]          ;plot axes through (0,0)

		;compare with standard least squares fit
		OPLOT, xdata,yfit
			
		;Add some text
		text=STRCOMPRESS('number of points = '+ $
            	     STRING(FORMAT='(I4)',N_ELEMENTS(xdata(ja)))) 
		XYOUTS,0.65,0.23,/NORMAL,text
		
		text=STRCOMPRESS('ODR y-intercept = '+STRING(FORMAT='(F7.2)',fitcf(0)) $
      			+ '   +-' +STRING(FORMAT='(F5.2)',interr))
		XYOUTS,0.65,0.19,/NORMAL,text

		text=STRCOMPRESS('LSq y-intercept = '+STRING(FORMAT='(F7.2)',result2(0)) $
      			+ '   +-' +STRING(FORMAT='(F5.2)',interr2))
		XYOUTS,0.65,0.15,/NORMAL,text

		CCG_CLOSEDEV, dev=dev
	ENDIF

ENDIF		;wrt 'IF NOT KEYWORD_SET ... outside call

;construct return arrays
coeff=FLTARR(2,4)
coeff(0,0)=fitcf(0)		;ODR int
coeff(1,0)=interr		;ODR int err
coeff(0,1)=fitcf(1)		;slope
coeff(1,1)=errcf(1)		;ODR slope err
coeff(0,2)=chisquared		;ODR chi2
coeff(1,2)=chi2prob		;ODR chi2prob
coeff(0,3)=result2(0)		;LSq intercept
coeff(1,3)=result2(1)		;LSq slope

yfit1=fitcf(0) + fitcf(1)*xdata

yfit2=result2(0) + result2(1)*xdata

text1=STRCOMPRESS('ODR: y='+STRING(FORMAT='(F9.1)',fitcf(1))+$
	'x+'+STRING(FORMAT='(F7.2)',fitcf(0))+'+/-'+$
	STRING(FORMAT='(F7.2)',interr),/RE)
text2=STRCOMPRESS('LSQ: y='+STRING(FORMAT='(F9.1)',coeff(1,3))+$
	'x+'+STRING(FORMAT='(F7.2)',coeff(0,3)),/RE)

;stop
END