pro bcs_line,Chan_struct,Te6,wave,intensity,newcal=newcal,	$	
		wave_range=wave_range,				$	; Inputs
		ioncal=ioncal,atocal=atocal,ion_mult=ion_mult,	$	; Inputs
		elem_num=elem_num,abun=abun,dens=dens,		$	; Inputs
		Ainfo=Ainfo,					$	; Output
		nel=nel,ion_data=ion_data,ato_data=ato_data		; Outputs
;+
; NAME:
;   bcs_line
; PURPOSE:
;   Compute spectral line intensities in the BCS wavelength range.
;   Units are ph s-1 for EM=1.e50 cm^3 at the source (i.e., the Sun).
;
; CALLING SEQUENCE:
;   bcs_line,Chan,Te6,wave,intensity
;   bcs_line,Chan,Te6,wave,intensity,Ainfo=Ainfo
;   bcs_line,Chan,Te6,wave,intensity,Ainfo=Ainfo,dens=dens
;   bcs_line,Chan,Te6,wave,intensity,/newcal
;
; INPUTS:
;   Chan_struct = BCS Channel Number or BSC index
;   Te6	        = Vector of temperatures (MK)
;
; OUTPUTS:
;   wave      = Wavlengths of calculated lines
;   intensity = Returned line intensities (ph s-1 for EM=1.e50 cm^3 at the source).
;
; OPTIONAL INPUT PARAMETERS:
;   wave_range = Range of wavelengths to include (Not implemented yet)  ;***
;   ioncal    = The ionization balance calculation (def = 0)
;   atocal    = Specify the atomic calculation to use (see get_atomic)
;   elem_num  = The abundance table to use (check get_elemabun)
;   abun      = If defined, allows the abundance value (relative to H) to be 
;	        directly specified
;   dens      = log10(Density cm-3).  Only available for S XV (chan=4).
;   newcal    = To force a new calculation
;   ion_mult  = Multipliers for the ionization fractions.  This must be a
;	        floating vector of length 10.  If ion_mult is undefined, or 
;		is vector that is less than 10 in length, then ion_mult will 
;		be set to:  ion_mult = replicate(1.,10)  (i.e., all 1's).
; OPTIONAL OUTPUT PARAMETERS:
;   Ainfo     = Structure describing the atomic parameters data.
;   ion_data  = Ionization balance structure
;   ato_data  = Atomic data parameter structure
;   nel	      = Number of electrons in line
; COMMON BLOCKS:
;   bcs_line_comm = Holds the results from the previous call.  Will speed things
;  		    up slightly if the chan and Te6 don't change.
; RESTRICTIONS:
;   dens must be a scalar
; MODIFICATION HISTORY:
;   8-oct-93, J. R. Lemen (LPARL), Written
;  27-oct-93, JRL, Improved Ainfo code.  Made common block work for Te6 as a vector.
;   9-mar-94, JRL, No longer pre-define atocal before calling get_atomic
;		   Added the dens=dens and newcal=newcal options
;   8-nov-94, JRL, Added ion_mult keyword
;  17-dec-94, JRL, Fixed a bug introduced with the last modification
;-
on_error,2			; Return to caller

; --- Check the input parameters

if n_elements(Te6) eq 0 then message,'Te6 is undefined'
if n_elements(chan_struct) eq 0 then message,'Chan is undefined'

chan = gt_bsc_chan(chan_struct)
if n_elements(chan) gt 1 then message,'Chan must be a scalar'


; --- Set up a common block to save the results for next time

common bcs_line_comm,chan_prev,Te6_prev,ato_prev,ion_prev,	$
		wave_prev,intensity_prev,nel_prev,ion_data_prev,$
		qdens_prev,ion_mult_prev
if n_elements(Te6_prev) eq 0 then begin
   Te6_prev = -99
   chan_prev= -99
   ato_prev = -99
   ion_prev = -99
   qdens_prev = -99
   ion_mult_prev = replicate(1.,10)
endif


Atomic_Z = ([0,26,26,20,16])(chan)		; Atomic number for each channel
if n_elements(ioncal) eq 0 then ioncal = 0	; Default ionization balance

; --- Read the atomic data file
ato_data = get_atomic(Chan,atocal=atocal)
atocal1  = ato_data.cal			; Don't read if atocal not defined but chan already read
if n_elements(dens) eq 0 then dens = 0.	
  if n_elements(ion_mult) lt 10 then ion_multx = replicate(1.,10) else ion_multx = ion_mult

if (min(ion_multx - ion_mult_prev) eq 0) and	$
   (max(ion_multx - ion_mult_prev) eq 0) then old_ion = 1 else old_ion = 0

if (chan eq chan_prev) and 					$
   (n_elements(te6) eq n_elements(te6_prev)) and 		$
   (min(Te6-Te6_prev) eq 0) and	(max(Te6-Te6_prev) eq 0) and 	$
   (atocal1 eq ato_prev)  and (ioncal eq ion_prev) and		$
   (dens eq qdens_prev) and old_ion and (not keyword_set(newcal)) then begin
   wave      = wave_prev
   intensity = intensity_prev
   nel       = nel_prev
   ion_data  = ion_data_prev
   qdens     = qdens_prev
   ion_multx = ion_mult_prev
;  message,'Using previous calculation',/cont		;***
endif else begin       
;  message,'Making new calculation',/cont		;***

; --- Compute the line excitation rates (cm^3 s-1)

  collexc,ato_data,Te6,wave1,Nel1,exrate1,dens=dens,qdens=qdens
  drexc,  ato_data,Te6,wave2,Nel2,exrate2

  wave   = [wave1,wave2]
  Nel    = [Nel1, Nel2]
  exrate = [exrate1,exrate2]

; --- Set up the ionization balance
  ion_frac = transpose(get_ionbal(Atomic_Z,Te6,ioncal=ioncal,	$
			stages=stage1,ion_data=ion_data))

  Nel_b = atomic_Z - stage1 + 1
  lookup = intarr(max(Nel_b)+1) & lookup(Nel_b) = indgen(n_elements(Nel_b))
  ion_frac = ion_frac * rebin(ion_multx(nel_b),	$
		n_elements(nel_b),n_elements(Te6),/sample)

; --- Multipy excitation rates by ionization balance ---------

  intensity = exrate			; Set up the output array
  for i=0,n_elements(wave)-1 do intensity(i,*) = intensity(i,*) * ion_frac(lookup(nel(i)),*)

  intensity_prev = intensity		; Copy into the common variables to speed things up
  wave_prev      = wave			;- next time if the channel and Te6 don't change
  nel_prev       = nel			;- (This only works if Te6 has a length of 1)
  chan_prev      = chan
  Te6_prev       = Te6
  ato_prev	 = atocal1
  ion_prev	 = ioncal
  ion_data_prev  = ion_data
  qdens_prev	 = qdens
  ion_mult_prev	 = ion_multx
endelse					; (chan eq chan_prev) and ...

; --- Set up the abundances ------------
elem_abun = get_elemabun(elem_num,info=elem_info)

; Compute N(H) / Ne  (usually about .83)
N_H_N_e = 1./total(elem_abun.abun * findgen(n_elements(elem_abun.abun)))

; Did the caller specify his own abundance?
if n_elements(abun) ne 0 then begin
  L_abun = abun
  elem_info = 'User supplied abund.'
endif else begin
  L_abun = elem_abun.abun(Atomic_Z)
  elem_info = elem_abun.info
endelse

; --- Multiply by abundance and put into the correct units -----

intensity = intensity * L_abun * N_H_N_e	; Multipy by [N(Z)/N(H)] * [N(H)/Ne]

intensity = intensity * 1.e25 * 1.e25		; EM = 1.e50 cm-3

; --- Create the Ainfo information structure ---

bsc_struct,fit_bsc=Ainfo			; Use the BSC .fit structure

Ainfo.Z  		= Atomic_Z		; Atomic Number
Ainfo.st$elem 		= byte(ion_data.element); Element (string)
Ainfo.Chan 		= Chan			; BCS Channel Number
Ainfo.EM50 		= 1.			; Emission Measure
Ainfo.elem_num 		= elem_num		; Number of Abundance table used
Ainfo.abun 		= L_abun		; Abundance used [N(Z)/N(H)]
Ainfo.ion_mult		= ion_multx		; Multipliers used for ion fractions
Ainfo.density		= qdens			; Log10(Density cm-3)

nn = n_elements(Ainfo.st$abun_info) < strlen(elem_info) 
Ainfo.st$abun_info	= (byte(elem_info))(0:nn-1)	; Text describing abundance

Ainfo.ioncal		= ioncal		; Number of ionization balance calculation
Ainfo.atocal		= atocal1		; Numer of atomic calculation

nn = n_elements(Ainfo.st$ion_info) < strlen(ion_data.head)
Ainfo.st$ion_info	= (byte(ion_data.head))(0:nn-1)	; Text describing ionization balance calculation

break_file,ion_data.file,disk,dir,ion_file,ext & ion_file = ion_file + ext
nn = n_elements(Ainfo.st$ion_file) < strlen(ion_file)
Ainfo.st$ion_file	= (byte(ion_file))(0:nn-1)	; File containing ionization balance calculation

nn = n_elements(Ainfo.st$ato_info) < strlen(ato_data.head)
Ainfo.st$ato_info	= (byte(ato_data.head))(0:nn-1)	; Text describing atomic data calculation

break_file,ato_data.file,disk,dir,ato_file,ext & ato_file = ato_file + ext
nn = n_elements(Ainfo.st$ato_file) < strlen(ato_file)
Ainfo.st$ato_file	= (byte(ato_file))(0:nn-1)	; File containing atomic data parameters
     
if n_elements(Te6) gt 1 then begin
   Ainfo = replicate(Ainfo,n_elements(Te6))
   Ainfo.Te6		= Te6			; Electron Temperature
endif else Ainfo.Te6	= Te6(0)

end