Climate Audit

michel

There is also the free version, GDL:-

http://en.wikipedia.org/wiki/GNU_Data_Language

• Hugo M
  Posted Sep 26, 2011 at 5:08 AM | Permalink

  michel,

  regrettably, IDL implements more functions than GDL currently does. “JULDAY” is probably only one of them.

  GDL> main
  % Compiled module: MAIN.
  % Compiled module: OBS.
  % Compiled module: FILE_LINES.
  % Compiled module: STRSPLIT.
  % Compiled module: UNIQ.
  % Compiled module: MEAN.
  % OBS: Function not found: JULDAY
  % Execution halted at: OBS 45 obs.pro
  % MAIN 15 main.pro
  % $MAIN$
  GDL>

  • Joel
    Posted Sep 26, 2011 at 12:14 PM | Permalink

    You can get the source code for JULDAY at:

    http://idlastro.gsfc.nasa.gov/ftp/ittvislib/julday.pro

    • P. Solar
      Posted Sep 27, 2011 at 2:19 AM | Permalink

      Thanks for filling that gap. However the next whole is TS_SMOOTH. I guess it won’t be the last. GDL is clearly less well equipt in terms of its library of fucntions at this stage and probably would require more coding effort than porting to R.

      Unfortunately neither approach allows identical replication of L&C but replication by similar means in another language should confirm the same result and if there’s a significant difference (not due to porting errors) that would be significant as well.

    • Hugo M
      Posted Sep 27, 2011 at 4:20 PM | Permalink

      P. Solar,

      I got the L&C-11 scripts running with gdl version 0.9.1 and closely reproduced at least “Fig.jpg”. Here is a list of urls for all needed IDL library routines:

      
h ttp://idlastro.gsfc.nasa.gov/ftp/ittvislib/julday.pro
h ttp://idlastro.gsfc.nasa.gov/ftp/ittvislib/label_date.pro
h ttp://idlastro.gsfc.nasa.gov/ftp/pro/plot/legend.pro
h ttp://idlastro.gsfc.nasa.gov/ftp/ittvislib/regress.pro
h ttp://idlastro.gsfc.nasa.gov/ftp/ittvislib/timegen.pro
h ttp://idlastro.gsfc.nasa.gov/ftp/ittvislib/ts_coef.pro
h ttp://idlastro.gsfc.nasa.gov/ftp/ittvislib/ts_fcast.pro
h ttp://idlastro.gsfc.nasa.gov/ftp/ittvislib/ts_smooth.pro
h ttp://astro.uni-tuebingen.de/software/idl/textoidl/textoidl.pro
h ttp://coco.atmos.washington.edu/ion_script/gamap2/graphics/screen2jpg.pro
h ttp://idl-coyote.googlecode.com/svn-history/r293/trunk/coyote/retired/tvread.pro


      Put the above files into the L&C-2011 scripts directory; cd into this directory and put the patch below into a file named, say, “LC-11.diff”. Then apply the lc11.diff using the shell command patch -p3 < LC-11.diff

      One last thing: My gdl version behaves very strangely (spurious syntax error messages) when gdl was called with main.pro on the command line. This errors did not show up when “main” was called from the gdl prompt.

      
diff -ur a/lindzen_choi_2011/scripts/amip.pro b/lindzen_choi_2011/scripts/amip.pro
--- a/lindzen_choi_2011/scripts/amip.pro 2011-09-25 18:40:42.000000000 +0200
+++ b/lindzen_choi_2011/scripts/amip.pro 2011-09-27 20:06:10.000000000 +0200
@@ -3,7 +3,7 @@
ObsName='ERBE&CERES'
clim_period=59
tot_period=288
-path=['./amip_olr/','./amip_osr/']
+path=['../amip_olr/','../amip_osr/']
file=['ncar_ccsm3_0_20n20s.dat',$
'mpi_echam51_20n20s.dat',$
'iap_fgoals1_0_g1_20n20s.dat',$
diff -ur a/lindzen_choi_2011/scripts/legend.pro b/lindzen_choi_2011/scripts/legend.pro
--- a/lindzen_choi_2011/scripts/legend.pro 2011-08-21 16:37:24.000000000 +0200
+++ b/lindzen_choi_2011/scripts/legend.pro 2011-09-27 21:09:00.000000000 +0200
@@ -442,7 +442,7 @@
if vectorfont[i] ne '' then begin
; if (num eq 1) and vertical then xp = x + xt/2 ; IF 1, CENTERED.
xyouts,xp,yp,vectorfont[i],width=width,color=colors[i] $
- ,charsize=charsize,align=xalign,charthick = charthick,/norm,font=font
+ ,charsize=charsize,align=xalign,charthick = charthick,/norm;,font=font
xt = xt > width
xp = xp + width/2.
endif else begin
@@ -456,11 +456,11 @@
TEXT_ONLY:
if vertical and (vectorfont[i] eq '') and symline and (linestyle[i] eq -99) then x=x0 + xspacing
xyouts,x,y,delimiter,width=width,/norm,color=textcolors[i], $
- charsize=charsize,align=xalign,charthick = charthick,font=font
+ charsize=charsize,align=xalign,charthick = charthick;,font=font
x = x + width*xsign
if width ne 0 then x = x + 0.5*xspacing
xyouts,x,y,items[i],width=width,/norm,color=textcolors[i],charsize=charsize, $
- align=xalign,charthick=charthick,font=font
+ align=xalign,charthick=charthick;,font=font
x = x + width*xsign
if not vertical and (i lt (n-1)) then x = x+2*xspacing; ADD INTER-ITEM SPACE
xfinal = (x + xspacing*margin)
diff -ur a/lindzen_choi_2011/scripts/main.pro b/lindzen_choi_2011/scripts/main.pro
--- a/lindzen_choi_2011/scripts/main.pro 2011-09-25 18:40:40.000000000 +0200
+++ b/lindzen_choi_2011/scripts/main.pro 2011-09-27 22:50:22.000000000 +0200
@@ -1,4 +1,5 @@
pro main
+myctblack=0

      ; -----------------------------------------------------------------------------
      ; This is the main program calling all procedures, calculating feedbacks, and
      @@ -24,7 +25,7 @@
      ; endfor

      ; ---------------
      -; Display 1 - Timeseries
      +print,"Display 1 - Timeseries"
      ; ---------------
      window,1,xs=600,ys=500
      loadct,5
      @@ -32,19 +33,19 @@
      !p.charsize=1.3
      !P.background=255
      dummy=label_date(date_format=('%y'))
      - plot,jul,sst,color=!myct.black,$
      + plot,jul,sst,color=myctblack,$
      xtitle='Year',ytitle='T (K)',$
      xtickunits=['time'],xtickformat='label_date',xstyle=1,xminor=5
      - plot,jul,replicate(0,258),color=!myct.black,yr=[-4,8],$
      + plot,jul,replicate(0,258),color=myctblack,yr=[-4,8],$
      xtitle='Year',ytitle=textoidl('R (Wm^{-2})'),$
      xtickunits=['time'],xtickformat='label_date',xstyle=1,xminor=5
      oplot,jul,olr,color=112 ;red
      oplot,jul,osr,color=47 ;blue
      - legend,label=['LW','SW'],line=[0,0],lcolor=[112,47],$
      - charsize=1.3,boxcolor=255,spacing=1.3
      + legend,['LW','SW'],line=[0,0],colors=[112,47],$
      + charsize=1.3,outline_color=255,spacing=1.3

      ; ---------------
      -; Feedback analysis for ERBE and CERES;
      +print, "Feedback analysis for ERBE and CERES"
      ; Call lc11.pro procedure, which include both LC11 method and simple regression.
      ; ---------------
      ; Call procedure for calculating climate feedback
      @@ -55,8 +56,12 @@
      LagDim,SlopDim,ConstDim,CorrDim,SlopDim2,ConstDim2,CorrDim2

      ; ---------------
      -; Display 2 - Timeseries
      - goto,jump
      +print, "Display 2 - Timeseries"
      +; goto,jump
      +; goto not available in GDL outside of procedures, hence:
      +if 0 ne 0 then begin
      +
      +
      ; ---------------
      ; SST timeseries graph with warming and cooling intervals
      window,2,xs=600,ys=300
      @@ -78,7 +83,7 @@
      endfor
      ;
      ; ---------------
      -; Display 3 - Feedback analysis for ERBE and CERES
      +print,"Display 3 - Feedback analysis for ERBE and CERES"
      ; ---------------
      window,3,xs=600,ys=600
      !p.multi=[0,2,2,0,0]
      @@ -90,12 +95,12 @@
      xtitle='T (K)',ytitle='R (W/m2)' ;textoidl('R (Wm^{-2})')
      oplot,findgen(10)-5,SlopDim2(5)*(findgen(10)-5)+ConstDim2(5)
      if ConstDim2(5) ge 0 then begin
      - legend,label=['y='+string(SlopDim2(5),format='(f4.1)')+'x+'$
      - +strtrim(string(ConstDim2(5),format='(f4.1)'),1)],boxcolor=-5,charsize=1.3
      + legend,['y='+string(SlopDim2(5),format='(f4.1)')+'x+'$
      + +strtrim(string(ConstDim2(5),format='(f4.1)'),1)],outline_color=-5,charsize=1.3
      endif
      if ConstDim2(5) lt 0 then begin
      - legend,label=['y='+string(SlopDim2(5),format='(f4.1)')+'x'$
      - +strtrim(string(ConstDim2(5),format='(f4.1)'),1)],boxcolor=-5,charsize=1.3
      + legend,['y='+string(SlopDim2(5),format='(f4.1)')+'x'$
      + +strtrim(string(ConstDim2(5),format='(f4.1)'),1)],outline_color=-5,charsize=1.3
      endif
      f_sim= SlopDim2(5) ; Feedback from simple regression
      ; Scatter plots for Lindzen-Choi Method
      @@ -105,14 +110,14 @@
      xtitle='T (K)',ytitle='R (W/m2)';textoidl('R (Wm^{-2})')
      oplot,findgen(10)-5,SlopDim(pos(0))*(findgen(10)-5)+ConstDim(pos(0))
      if ConstDim(pos(0)) ge 0 then begin
      - legend,label=['y='+string(SlopDim(pos(0)),format='(f4.1)')+'x+'$
      + legend,['y='+string(SlopDim(pos(0)),format='(f4.1)')+'x+'$
      +strtrim(string(ConstDim(pos(0)),format='(f4.1)'),1)],$
      - boxcolor=-5,charsize=1.3
      + outline_color=-5,charsize=1.3
      endif
      if ConstDim(pos(0)) lt 0 then begin
      - legend,label=['y='+string(SlopDim(pos(0)),format='(f4.1)')+'x'$
      + legend,['y='+string(SlopDim(pos(0)),format='(f4.1)')+'x'$
      +strtrim(string(ConstDim(pos(0)),format='(f4.1)'),1)],$
      - boxcolor=-5,charsize=1.3
      + outline_color=-5,charsize=1.3
      endif
      f_lc=SlopDim(pos(0)) ; Feedback from LC11
      ; Lead-lag slopes
      @@ -120,26 +125,28 @@
      xtitle='Lagged month',xminor=2, $
      ytitle='R/T (W/m2/K)';textoidl('R/T (Wm^{-2}K^{-1})')
      oplot,LagDim,SlopDim2,linestyle=1
      - legend,label=['LC','Simple'],line=[0,1],boxcolor=-5,charsize=1.3,spacing=2
      + legend,['LC','Simple'],line=[0,1],outline_color=-5,charsize=1.3,spacing=2
      ; Lead-lag correlations
      plot,LagDim,CorrDim,xr=[-5,5],xs=1,yr=[-1,1],$
      title='Lag Correlation',xtitle='Lagged month',xminor=2
      oplot,LagDim,CorrDim2,linestyle=1
      - legend,label=['LC','Simple'],line=[0,1],boxcolor=-5,charsize=1.3
      + legend,['LC','Simple'],line=[0,1],outline_color=-5,charsize=1.3

      - jump: print,'Skip Display 2 and 3'
      +;jump:
      + endif
      + print,'Skip Display 2 and 3'

      ; ---------------
      -; Open AMIP data
      +print, "Open AMIP data"
      ; Call amip.pro procedure, which include 11 AMIP model outputs
      ; ---------------
      -; Call proceedure for reading AMIP model outputs
      +print,"Call proceedure for reading AMIP model outputs"
      amip,olr_amip,osr_amip
      -; Generate julian day corresponding to AMIP outputs
      +print,"Generate julian day corresponding to AMIP outputs"
      julmo=timegen(start=julday(1,1,1985),final=julday(12,1,2008),units='M')

      ; ---------------
      -; Display 4 - Timeseries for AMIP data
      +print, "Display 4 - Timeseries for AMIP data -1"
      ; ---------------
      title=['CCSM3','ECHAM5/MIP-OM','FGOALS-g1.0',$
      'GFDL-CM2.1','GISS-ER','INM-CM3.0','IPSL-CM4','MRI-CGCM2.3.2',$
      @@ -148,76 +155,86 @@
      window,4,xs=800,ys=800
      loadct,5
      !p.multi=[0,2,7,0,0]
      +
      + print, "Display 4 - Timeseries for AMIP data -2"
      +
      dummy=label_date(date_format=('%y'))
      - plot,julmo,replicate(0,288),color=!myct.black,yr=[-8,8],$
      + plot,julmo,replicate(0,288),color=myctblack,yr=[-8,8],$
      ytitle='R (W/m2)',$ ;textoidl('R (Wm^{-2})'),$
      xticklen=0.1,xtickunits=['time'],xtickformat='label_date',$
      xstyle=1,xminor=5,charsize=2,xmargin=[6,2],$
      ystyle=1,yminor=2,yticks=2,ymargin=[1.9,0.3]
      - legend,label=['ERBE&CERES'],$
      - boxcolor=-5,charsize=1.3,spacing=1.4,halign=0.03,valign=0.91
      + legend,['ERBE&CERES'],$
      + outline_color=-5,charsize=1.3,spacing=1.4 ;,halign=0.03,valign=0.91
      +
      +print, "Display 4 - Timeseries for AMIP data -3"
      oplot,jul,olr,color=112 ;red
      oplot,jul,osr,color=47 ;blue
      - legend,label=['LW','SW'],line=[0,0],lcolor=[112,47],$
      - charsize=1,boxcolor=255,spacing=1
      + legend,['LW','SW'],line=[0,0],colors=[112,47],$
      + charsize=1,outline_color=255,spacing=1
      + print, "Display 4 - Timeseries for AMIP data -4"
      for m=0,10 do begin
      - plot,julmo,replicate(0,288),color=!myct.black,yr=[-8,8],$
      + print, "Display 4 - Timeseries for AMIP data -5"
      + plot,julmo,replicate(0,288),color=myctblack,yr=[-8,8],$
      ytitle='R (W/m2)',$ ;textoidl('R (Wm^{-2})'),$
      xticklen=0.1,xtickunits=['time'],xtickformat='label_date',$
      xstyle=1,xminor=5,charsize=2,xmargin=[6,2],$
      ystyle=1,yminor=2,yticks=2,ymargin=[1.9,0.3]
      - legend,label=[title(m)],$
      - boxcolor=-5,charsize=1.3,spacing=1.4,halign=0.03,valign=0.91
      + legend,[title(m)],$
      + outline_color=-5,charsize=1.3,spacing=1.4 ;,halign=0.03,valign=0.91
      oplot,julmo,ts_smooth(olr_amip(*,m),3),color=112 ;red
      oplot,julmo,ts_smooth(osr_amip(*,m),3),color=47 ;blue
      endfor

      ; ---------------
      -; Display 5 - Feedback analysis for AMIP 11 models
      +print, "Display 5 - Feedback analysis for AMIP 11 models"
      ; ---------------
      window,5,xs=800,ys=600
      !p.multi=[0,4,3,0,0]
      !x.margin=[7,2]
      !y.margin=[3,2]

      -; ERBE&CERES results
      +print, "ERBE&CERES results"
      +
      pos=where(CorrDim eq max(CorrDim(5:10)))
      plot,LagDim,SlopDim,xr=[-5,5],xs=1,yr=[-4,6],title ='ERBE&CERES',$
      xtitle='Lagged month',xminor=2,charsize=2, $
      ytitle='R/T (W/m2/K)';textoidl('R/T (Wm^{-2}K^{-1})')
      oplot,LagDim,SlopDim2,linestyle=1
      - legend,label=['LC','Simple'],line=[0,1],boxcolor=-5,charsize=1.3,spacing=1.4
      - legend,label=['F_LC ='+string(SlopDim(pos(0)),format='(f4.1)'),$
      + legend,['LC','Simple'],line=[0,1],outline_color=-5,charsize=1.3,spacing=1.4
      + legend,['F_LC ='+string(SlopDim(pos(0)),format='(f4.1)'),$
      'F_sim='+string(SlopDim2(5),format='(f4.1)')],$
      - boxcolor=-5,charsize=1.3,spacing=1.4,halign=0.05,valign=0.95
      -; 11 AMIP model results
      + outline_color=-5,charsize=1.3,spacing=1.4
      + ;,halign=0.05,valign=0.95
      +
      +print, '11 AMIP model results'
      for m=0,10 do begin ; 11 models

      olrmo=olr_amip(*,m) ; AMIP monthly LW flux from 1985
      osrmo=osr_amip(*,m) ; AMIP monthly SW flux from 1985

      -; Smoothing and scaling
      +print, "Smoothing and scaling"
      sstmo=ts_smooth(sst*0.78,3)
      olrmo=ts_smooth(olrmo,3)
      osrmo=ts_smooth(osrmo,3)

      -; Call LC11.pro procedure
      +print, "Call LC11.pro procedure"
      ; put 'olrmo' for LW, put 'osrmo' for SW, put 'olrmo+osrmo' for Net
      fluxmo=olrmo
      lc11,julmo,sstmo,fluxmo,int1,int2,dT,dTmax,dR,dRmax,$
      LagDim,SlopDim,ConstDim,CorrDim,SlopDim2,ConstDim2,CorrDim2
      ;
      - ; Lead-lag slopes
      +print, "Lead-lag slopes"
      pos=where(CorrDim eq max(CorrDim(5:10)))
      plot,LagDim,SlopDim,xr=[-5,5],xs=1,yr=[-4,6],title =title(m),$
      xtitle='Lagged month',xminor=2,charsize=2, $
      ytitle='R/T (W/m2/K)';textoidl('R/T (Wm^{-2}K^{-1})')
      oplot,LagDim,SlopDim2,linestyle=1
      - ;legend,label=['LC','Simple'],line=[0,1],boxcolor=-5,charsize=1.3,spacing=1.4
      - legend,label=['F_LC ='+string(SlopDim(pos(0)),format='(f4.1)'),$
      + ;legend,['LC','Simple'],line=[0,1],outline_color=-5,charsize=1.3,spacing=1.4
      + legend,['F_LC ='+string(SlopDim(pos(0)),format='(f4.1)'),$
      'F_sim='+string(SlopDim2(5),format='(f4.1)')],$
      - boxcolor=-5,charsize=1.3,spacing=1.4,halign=0.05,valign=0.95
      + outline_color=-5,charsize=1.3,spacing=1.4;,halign=0.05,valign=0.95

      endfor
      - screen2jpg,'Fig.jpg'
      + screen2jpg,'Fig-gdl.jpg'
      end
      diff -ur a/lindzen_choi_2011/scripts/obs.pro b/lindzen_choi_2011/scripts/obs.pro
      --- a/lindzen_choi_2011/scripts/obs.pro 2011-09-25 18:40:42.000000000 +0200
      +++ b/lindzen_choi_2011/scripts/obs.pro 2011-09-27 20:05:35.000000000 +0200
      @@ -5,7 +5,7 @@

      pro obs,jul,sst,sstmo,olr,osr

      -path1='./obs/'
+path1='../obs/'
file1='oisst_36_20n20s_1985_1999.txt'
file2='amipsst_mon_20n20s_1985_2008.txt'
file3='erbe_s10n_36_20n20s_1985_1999.txt'


    • Hugo M
      Posted Sep 27, 2011 at 10:55 AM | Permalink

      The last missing function is regress(). The L&C-2011 script seems to use an IDL-7.x version of regress(), which isn’t easily availabe on the net, but is contained as IDL source code within IDL installations. All others missing functions are available:

      
label_date.pro
legend.pro
textoidl.pro
ts_fcast.pro
julday.pro
ts_coef.pro
ts_smooth.pro


      Since these files are mostly from older IDL versions, some minor editing of L&C scripts is needed (such as adapting named parameter to the older definitions or replacing color=”!myct.black” by color=0). If anyone understands the regress() call wrt the parameters “Correlation” and “Const” so far that the call could be adapted to use the older version available on the net, then a GDL version of the L&C scripts could become fully functional. It may make it easier to develop and test the R-version.

      Actual 7.x – version of regress()

      Result = REGRESS( X, Y, [, CHISQ=variable] [, CONST=variable] [, CORRELATION=variable] [, /DOUBLE] [, FTEST=variable] [, MCORRELATION=variable] [, MEASURE_ERRORS=vector] [, SIGMA=variable] [, STATUS=variable] [, YFIT=variable] )

      http://star.pst.qub.ac.uk/idl/REGRESS.html

      Obsoleted 5.4 Version of regress()

      REGRESS(X, Y, Weights, Yfit, Const, Sigma, Ftest, R, Rmul, Chisq)

      http://star.pst.qub.ac.uk/idl/Obs_in_54.html

    • jgc
      Posted Sep 28, 2011 at 1:53 AM | Permalink

      Here is the IDL 8 version, which should be compatible with IDL 7x

      ; $Id: //depot/idl/releases/IDL_80/idldir/lib/regress.pro#1 $
      ;
      ; Copyright (c) 1982-2010, ITT Visual Information Solutions. All
      ;       rights reserved. Unauthorized reproduction is prohibited.
      
      ;+
      ; NAME:
      ;	REGRESS
      ;
      ; PURPOSE:
      ;	Perform a multiple linear regression fit.
      ;
      ;	REGRESS fits the function:
      ;		Y[i] = Const + A[0]*X[0,i] + A[1]*X[1,i] + ... +
      ;                      A[Nterms-1]*X[Nterms-1,i]
      ;
      ; CATEGORY:
      ;       G2 - Correlation and regression analysis.
      ;
      ; CALLING SEQUENCE:
      ;
      ;	Result = REGRESS(X, Y)
      ;
      ; INPUTS:
      ;
      ;       X:	The array of independent variable data.  X must
      ;		be dimensioned as an array of Nterms by Npoints, where
      ;		there are Nterms coefficients (independent variables) to be
      ;		found and Npoints of samples.
      ;
      ;       Y:	The vector of dependent variable points.  Y must have Npoints
      ;		elements.
      ;
      ; OUTPUTS:
      ;
      ;	REGRESS returns a column vector of coefficients that has Nterms
      ;	elements.
      ;
      ; KEYWORDS:
      ;
      ;   CHISQ:   Sum of squared errors divided by MEASURE_ERRORS if specified.
      ;
      ;   CONST:   Constant term. (A0)
      ;
      ;	CORRELATION: Vector of linear correlation coefficients.
      ;
      ;   COVAR:   Covariance matrix of the coefficients.
      ;
      ;	DOUBLE:  if set, force computations to be in double precision.
      ;
      ;	FTEST:	The value of F for goodness-of-fit test.
      ;
      ;	MCORRELATION:   The multiple linear correlation coefficient.
      ;
      ;   MEASURE_ERRORS: Set this keyword to a vector containing standard
      ;       measurement errors for each point Y[i].  This vector must be the same
      ;       length as X and Y.
      ;
      ;     Note - For Gaussian errors (e.g. instrumental uncertainties),
      ;        MEASURE_ERRORS should be set to the standard
      ; 	     deviations of each point in Y. For Poisson or statistical weighting
      ; 	     MEASURE_ERRORS should be set to sqrt(Y).
      ;
      ;   SIGMA:   The 1-sigma error estimates of the returned parameters.
      ;
      ;     Note: if MEASURE_ERRORS is omitted, then you are assuming that the
      ;           regression is the correct model. In this case,
      ;           SIGMA is multiplied by SQRT(CHISQ/(N-M)), where N is the
      ;           number of points in X. See section 15.2 of Numerical Recipes
      ;           in C (Second Edition) for details.
      ;
      ;   STATUS = Set this keyword to a named variable to receive the status
      ;          of the operation. Possible status values are:
      ;          0 for successful completion, 1 for a singular array (which
      ;          indicates that the inversion is invalid), and 2 which is a
      ;          warning that a small pivot element was used and that significant
      ;          accuracy was probably lost.
      ;
      ;    Note: if STATUS is not specified then any error messages will be output
      ;          to the screen.
      ;
      ;   YFIT:   Vector of calculated Y's.
      ;
      ; PROCEDURE:
      ;	Adapted from the program REGRES, Page 172,
      ;	Bevington, Data Reduction and Error Analysis for the
      ;	Physical Sciences, 1969.
      ;
      ; MODIFICATION HISTORY:
      ;	Written, DMS, RSI, September, 1982.
      ;	Added RELATIVE_WEIGHT keyword    W. Landsman   August 1991
      ;       Fixed bug in invert  Bobby Candey 1991 April 22
      ;       Added STATUS argument.  GGS, RSI, August 1996
      ;   CT, RSI, March 2000: Fixed status flag. Cleaned up help.
      ;   CT, RSI, July 2000: Change arguments to keywords.
      ;         Add MEASURE_ERRORS [equivalent to 1/sqrt(Weights)],
      ;         removes need for RELATIVE_WEIGHT.
      ;-
      ;
      FUNCTION REGRESS,X,Y, $
      	weights_old,yfit_old,const_old,sigma_old, $       ; obsolete arguments
      	ftest_old,r_old,rmul_old,chisq_old,status_old, $  ; obsolete arguments
      	RELATIVE_WEIGHT=relative_weight, $                ; obsolete keyword
      	CHISQ=chisq, $
      	CONST=const, $
      	DOUBLE=double, $
      	FTEST=ftest, $
      	CORRELATION=correlation, $
      	MCORRELATION=mcorrelation, $
      	MEASURE_ERRORS=measure_errors, $
      	SIGMA=sigma, $
      	STATUS=status, $
      	YFIT=yfit
      
      COMPILE_OPT idl2
      
      ON_ERROR,2              ;Return to caller if an error occurs
      sy = SIZE(Y)            ;Get dimensions of x and y.
      sx = SIZE(X)
      
      ndimX = sx[0]
      nterm = (ndimX EQ 1) ? 1 : sx[1]   ; # of terms (coefficients)
      nptsX = sx[ndimX]       ;# of observations (samples)
      npts = sy[1]            ;# of observations (samples)
      
      IF (nptsX NE npts) THEN MESSAGE, $
      	'X and Y have incompatible dimensions.'
      IF (ndimX EQ 1) THEN x = REFORM(x, 1, npts, /OVER)   ; change X to a 2D array
      
      double = (N_ELEMENTS(double) GT 0) ? KEYWORD_SET(double) : $
      	((SIZE(x,/TNAME) EQ 'DOUBLE') OR (SIZE(y,/TNAME) EQ 'DOUBLE'))
      one = double ? 1d : 1.0
      two = double ? 2d : 2.0
      
      
      IF (N_PARAMS() GT 2) THEN BEGIN      ; old-style REGRESS (with arguments)
      	IF (N_ELEMENTS(weights_old) NE npts) THEN MESSAGE, $
      		'Weight has been replaced by MEASURE_ERRORS. See Online Help.'
      	no_weights = KEYWORD_SET(relative_weight)
      	weights = no_weights ? one : weights_old  ; note there is no ^2
      ENDIF ELSE BEGIN      ; new style (with keywords)
      	IF KEYWORD_SET(relative_weight) THEN $
      		MESSAGE,/INFO,'Keyword RELATIVE_WEIGHT is obsolete.'
      	no_weights = N_ELEMENTS(measure_errors) NE npts
          ; note different meaning from Weights argument
      	weights = no_weights ? one : 1/measure_errors^two
      ENDELSE
      
      sw = no_weights ? npts : TOTAL(weights, DOUBLE=double) ;sum of weights
      ymean = TOTAL(y*weights, DOUBLE=double)/sw   ;y mean
      wgt = no_weights ? one : REBIN(TRANSPOSE(weights),nterm,npts)
      xmean = TOTAL(wgt*x,2, DOUBLE=double)/sw
      wmean = sw/npts
      wgt = TEMPORARY(wgt)/wmean
      ww = weights/wmean
      
      nfree = npts-1          ;degs of freedom
      sigmay = SQRT(TOTAL(ww * (y-ymean)^2)/nfree) ;weights*(y(i)-ymean)
      xx = x - REBIN(xmean,nterm,npts)     ;x(j,i) - xmean(i)
      wx = TEMPORARY(wgt) * xx             ;weights(i)*(x(j,i)-xmean(i))
      sigmax = SQRT(TOTAL(xx*wx,2)/nfree)  ;weights(i)*(x(j,i)-xm)*(x(k,i)-xm)
      array = (wx # TRANSPOSE(xx))/(nfree * sigmax #sigmax)
      
      array = INVERT(array, status)
      IF (status EQ 1L && ~ARG_PRESENT(status)) THEN BEGIN
      	IF (ndimX EQ 1) THEN x = REFORM(x, /OVER)    ; change X back to a vector
      	MESSAGE, "Inversion failed due to singular array."
      END
      IF (status EQ 2L && ~ARG_PRESENT(status)) THEN BEGIN
      	MESSAGE, /INFO, $
      	    "Inversion used a small pivot element. Results may be inaccurate."
      endif
      
      correlation = (TEMPORARY(wx) # (y - ymean)) / (sigmax * sigmay * nfree)
      a = (correlation # array)*(sigmay/sigmax)         ;get coefficients
      yfit = a # x                            ;compute fit
      const = ymean - TOTAL(a*xmean)             ;constant term
      yfit = yfit + const                        ;add it in
      freen = npts-nterm-1 > 1                ;degs of freedom, at least 1.
      chisq = TOTAL(ww*(y - yfit)^2)*wmean ;weighted chi squared
      
      
      ; correction factor for no weights (see Num.Rec. sec 15-2)
      varnce = no_weights ? chisq/freen : one/wmean
      
      ; error terms
      sigma = SQRT(array[LINDGEN(nterm)*(nterm+1)]*varnce/(nfree*sigmax^2))
      
      mult_covariance = TOTAL(a*correlation*sigmax/sigmay)
      ftest = (mult_covariance GE 1.) ? 1e6 : $
      	mult_covariance/nterm / ((1.-mult_covariance)/freen)
      mcorrelation = SQRT(mult_covariance)
      IF (ndimX EQ 1) THEN x = REFORM(x, /OVER)    ; change X back to a vector
      
      
      ; fill in obsolete arguments, if necessary
      IF (N_PARAMS() GT 3) THEN BEGIN
      	status_old = status
      	chisq_old = chisq/freen    ; reduced chi-square
      	rmul_old = mcorrelation
      	r_old = correlation
      	ftest_old = ftest
      	sigma_old = sigma
      	const_old = const
      	yfit_old = yfit
      ENDIF
      
      RETURN, a
      END