PRO lambdacolfit
;NAME:lambdacolfit
;
;PURPOSE:calibrate the wavelength for the SED mode -- fit a line to the 
;row number and the wavelength - for each of the 24 columns
;
;CALLING SEQUENCE: lambdacolfit
;
;INPUTS:name of file with wavelength in microns and the col #
;"LambdaCal.txt"
;
;OPTIONAL KEYWORD:
;
;OUTPUTS: LambdaCols.txt - 
;
;SIDE EFFECTS:
;
;RESTRICTIONS:??
;
;EXAMPLES:
;
;MODIFICATION HISTORY:started on 3/3/03
;
;------------------------------------

name='LambdaCal.txt'
Close,10
OPENR,10, name

Close,20
OPENW,20, 'LambdaCols.txt'
PRINTF, 20, 'Col #         Lambda(0)  +/-          Dispersion      +/-'
;spawn,'wc '+infile,wc_out
spawn,'wc '+name,wc_out
nlines = (strsplit(strtrim(wc_out[0]),/extract))[0]

;colnum=dindgen(24)
colnum=fltarr(24)
dump=strarr(nlines)
lambda=strarr(nlines)
maxcol=strarr(nlines)
pound=strarr(nlines)
startcol=strarr(nlines)
endcol=strarr(nlines)

READF,10,dump
pound=strmid(dump,0,1)
lambda=(strmid(dump,6,7))
maxcol=(strmid(dump,19,7))
startcol=(strmid(dump,32,7))
endcol=(strmid(dump,45,7))
print, pound


For j=0,23 do begin
;FOr j=0,2 do begin
k=0

FOR i=0,nlines-1 do begin
 IF (pound(i) ne '#') THEN BEGIN
  IF(j ge FIX(startcol(i)) AND j le FIX(endcol(i))) THEN BEGIN
   k=k+1
 ENDIF
 ENdif
ENDFOR

;ncol=where( fix(startcol) le j AND fix(endcol) ge j,k)
;print,'k=',k,'ncol=',ncol
xs=fltarr(k)
ys=fltarr(k)
xys=fltarr(k)
xsquare=fltarr(k)
;sigma2=fltarr(k)
;intsigma2=fltarr(k)
;slopesigma2=fltarr(k)
k=0

FOR i=0,nlines-1 do begin
 IF (pound(i) ne '#') THEN BEGIN 
   IF(j ge FIX(startcol(i)) AND j le FIX(endcol(i))) THEN BEGIN
   ys(k)=FLOAT(lambda(i))
; ys=FLOAT(lambda(ncol))
;print,'ys=',ys
   xs(k)=FLOAT(maxcol(i))
; xs=Float(maxcol(ncol))
;print,'xs=',xs
   k=k+1
 ENDIF
 ENDIF
ENDFOR


for h= 0,k-1 DO Begin
   xys(h)=xs(h)*ys(h)
   xsquare(h)=xs(h)*xs(h)
;  xys(h)=lambda(h)*maxcol(h)
;  xsquare(h)=lambda(h)*lambda(h)
endfor

delta=(k*TOTAL(xsquare)-(TOTAL(xs))^2)
;Slope
 m=(k*TOTAL(xys)-TOTAL(xs)*TOTAL(ys))/delta
;y intercept
 b=(TOTAL(xsquare)*TOTAL(ys) - TOTAL(xs)*TOTAL(xys))/delta

;UNCERTAINTIES
;sigma^2
sigma2=(TOTAL((ys-(m*xs+b))^2))/(k-2)
;y-intercept uncertainty squared
 intsigma2=(sigma2/delta)*TOTAL(xs^2)
;slope uncertainty squared
slopesigma2=k*sigma2/delta

PRINT, 'For col #',j
PRINT, 'Slope= ',m,'  +/-',SQRT(slopesigma2)
PRINT, 'Y-Int= ',b,'  +/-',SQRT(intsigma2)
PRINTF,20, j, b, SQRT(intsigma2), m, SQRT(slopesigma2)
;PLOT, xs, (m*xs+b),ytit='wavelength in microns', xtit='row #'


;stop

endfor ; each col - j
;close,20

;*****************************************************
k=0
FOR i=0,nlines-1 do begin
  IF (pound(i) ne '#') THEN BEGIN
   k=k+1
  ENdif
ENDFOR
xs=fltarr(k)
ys=fltarr(k)
xys=fltarr(k)
xsquare=fltarr(k)
k=0

FOR i=0,nlines-1 do begin
  IF (pound(i) ne '#') THEN BEGIN
   ys(k)=FLOAT(lambda(i))
   xs(k)=FLOAT(maxcol(i))
   k=k+1
  ENDIF
ENDFOR

for h= 0,k-1 DO Begin
   xys(h)=xs(h)*ys(h)
   xsquare(h)=xs(h)*xs(h)
;  xys(h)=lambda(h)*maxcol(h)
;  xsquare(h)=lambda(h)*lambda(h)
endfor

delta=(k*TOTAL(xsquare)-(TOTAL(xs))^2)
;Slope
 m=(k*TOTAL(xys)-TOTAL(xs)*TOTAL(ys))/delta
;y intercept
 b=(TOTAL(xsquare)*TOTAL(ys) - TOTAL(xs)*TOTAL(xys))/delta


;UNCERTAINTIES
;sigma^2
sigma2=(TOTAL((ys-(m*xs+b))^2))/(k-2)
;y-intercept uncertainty squared
 intsigma2=(sigma2/delta)*TOTAL(xs^2)
;slope uncertainty squared
slopesigma2=k*sigma2/delta

Print, ' '
PRINT, 'For all cols combined:'
PRINT, 'Slope= ',m,'  +/-', SQRT(slopesigma2)
PRINT, 'Y-Int= ',b,'  +/-',SQRT(intsigma2)
PRINTF,20, '  '
PRINTF,20, 'For all cols combined:'
PRINTF,20, 'Lam(0)= ',b,'  +/-',SQRT(intsigma2)
PRINTF,20, 'Dispersion=',m,'  +/-',SQRT(slopesigma2)
;PLOT, xs, (m*xs+b),ytit='wavelength in microns', xtit='row #'
close,20


END