;+
; NAME:
;	RD_BCSCAL
; PURPOSE:
;       to obtain sensitivity factor, central wavelength, and
;       physical position of right-hand boundaries of BCS bins.
; CALLING SEQUENCE:
;	RD_BCSCAL,CHAN
; INPUTS:
;	CHAN         - BCS channel
; KEYWORD OUTPUTS:
;       W0           - start wavelength of detector (A)
;       EDGES        - physical position of right-hand edges of bin
;       SENS         - sensitivity factor for bin
;       DISP         - dispersion (A/bin)
;       EFAREA       - detector effective area(cm^2)
;       BINARC       - bins per arcmin
;       ROCKW        - rocking width (A)
;       GAUSSW       - gaussian detector width (A)
;       BSTART       - first valid bin number
;       BEND         - last valid bin number
;       VERSION      - version number of CAL file
;	HEAD	     - Header string of the CAL file
; RESTRICTIONS
;       Will reject old version calibration files that do not contain
;       the version number and effective area value.
; PROCEDURE:
;       Reads latest version of Sin(theta) vs Sensitivity file for
;       specified channel.
;       Determines first and last valid bin numbers (bstart,bend) from
;       first and last non-zero sin(theta) values.
;       Computes physical positions of RHS edges of each bin
;       and start wavelengths of each bin using dispersion relation.
; RESTRICTIONS:
;       The bins in the BCS calibration files correspond to ungrouped bins.
;       The physical position, wavelength, and sensitivity of grouped data
;       will depend on the grouping plan.
; HISTORY:
;       Written Oct. 1992 by D. Zarro (Applied Research Corp).
;	22-oct-93, J. R. Lemen (LPARL), Added Head keyword. Read file using rd_tfile to
;					avoid comment lines which begin with ";"
;	 3-nov-93, JRL, Updated to handle Format Ver 3 files (double column)
;-
pro rd_bcscal,chan,w0=w0,edges=edges,sens=sens,$
    disp=disp,efarea=efarea,binarc=binarc,rockw=rockw,bstart=bstart,$
    bend=bend,sinth=sinth,version=version,gaussw=gaussw,head=head

on_error,1

if n_params() eq 0 then return

twod=[4.,4.,4.,6.532] & maxb=256
rock=[14.,15.,50.,68.]             ;crystal rocking widths (arcsec)
rock=[0.2085e-3, 0.2063e-3, 0.3889e-3, 1.371e-3]      ;-- FWHM  in A 
gaussd=[2.71E-4, 5.62E-4, 2.94E-4, 10.0E-4]           ;-- FWHM (gaussian detector widths)                                                    


;
;-- Find the BCS calibration file for this channel
;

get_bcscal,chan,calfil,version=version
;
;-- Read DISPERSION, EFFECTIVE AREA and 256 SIN(THETA) and SENSITIVITY
;   values for this channel  -- Use rd_tfile to avoid comment lines
;

buff = rd_tfile(calfil,nocomm=';')	; Will treat ";" as the comment character
head = buff(0)				; Header 
if strpos(buff(1),'VERSION') eq -1 then 		$
   message,'CALFIL is out of date; please check setup environment' $
   else begin
      File_version = 0. & reads,buff(1),file_version
   endelse
disp=0.   & reads,buff(2),disp		; Dispersion
efarea=0. & reads,buff(3),efarea	; Effective area (cm-2)

if file_version le 2 then begin
  if n_elements(buff) ne 512 + 4 then message,'Few calibration bins than expected'
  sinth = float(buff(4:4+maxb-1))
  sens  = float(buff(4+maxb:4+2*maxb-1))	; Sensitivity
endif else begin
; Read file_version 3 files (double column)
  if n_elements(buff) ne 256 + 4 then message,'Few calibration bins than expected'
  sinth = fltarr(maxb)
  sens  = fltarr(maxb)
  xx = 0. & yy = 0.
  for i=0,maxb-1 do begin
    reads,buff(4+i),xx,yy
    sinth(i) = xx & sens(i) = yy
  endfor
endelse
  
;
;-- Find non-zero SINTH and SENSIT values
;

find=where( (sinth gt 0.) ,count)
bstart=0 & bend=maxb-1
if count gt 0 then begin
 bstart=find(0)
 bend=find(count-1)
endif

;
;-- Compute bin central wavelengths 
;

twodt=twod(chan-1)

;
;-- Assume unit sensitivity where not defined
;

find=where(sens le 0,count)
if count gt 0 then sens(find)=1.

;-- Compute left hand edges of first valid bin

fsinth=sinth(bstart)-(sinth(bstart+1)-sinth(bstart))/2.

;
;-- Compute RHS edges boundaries of each bin
;

edges=findgen(maxb)+1
for i=bstart,bend-1 do $
 edges(i)=edges(bstart-1)+(twodt/disp)*( (sinth(i+1)+sinth(i))/2. - fsinth)

;-- Last valid bin

edges(bend)=edges(bstart-1)+(twodt/disp)*(sinth(bend)-fsinth+$
           (sinth(bend)-sinth(bend-1))/2.)

for i=bend+1,maxb-1 do edges(i)=edges(i-1)+1.

; 
;-- compute wavelength offset 
;

w0=twodt*sinth(bstart)-disp*(edges(bstart)+edges(bstart-1))/2.

;-- angle to bin conversion

msinth=(sinth(bstart)+sinth(bend))/2.
binarc=twodt*(1./60./!radeg/disp)*sqrt(1.- msinth^2)

rockw=rock(chan-1)
gaussw=gaussd(chan-1)
return & end