function bcs_spec,Chan_struct,Te6,Td6,Wave,			$
		Velocity=Vel,cnorm=cnorm,			$	; Inputs
		ioncal=ioncal,atocal=atocal,Ashift=Ashift,	$	; Inputs
		elem_num=elem_num,abun=abun,utdoppw=utdoppw,	$	; Inputs
		noline=noline,nocont=nocont,dens=dens,		$	; Inputs
		ion_mult=ion_mult,newcal=newcal,		$	; Inputs
		waveout=wave0,bcs_cont=bcs_cont,Ainfo=Ainfo,	$	; Outputs
		ion_data=ion_data,ato_data=ato_data,		$	; Outputs
		rockw=urockw,edg_waveout=wave1				; Outputs
;+
; NAME:
;  bcs_spec
; PURPOSE:
;  Calculate a synthetic Yohkoh BCS spectrum in photons cm^2 s^-1 A^-1
;  for an EM=1.e50 cm^-3
; CALLING SEQUENCE:
;  Spectrum = bcs_spec(Chan,Te6)		; Use default wavelengths
;  Spectrum = bcs_spec(Chan,Te6,Td6,Wave,waveo=waveo,Vel=Vel,Ainfo=Ainfo)
;
; INPUTS:
;  Chan	     = BSC structure or
;              channel number (1 to 4)	1=Fe XXVI, 2=Fe XXV, 3=Ca XIX, 4= S XV
;              (Must be a scalar)
;  Te6	     = Electron temperature in MK (Can be a vector)
; OPTIONAL INPUTS:
;  Td6       = Doppler  temperature in MK (Can be a vector -- must be the same
;		length as Te6).  If not present, defaults to Td6=Te6
;  Wave	     = A vector describing the start wavelength of a vector of bins.  
;              bcs_spec will return n_elements(wave)-1 bins.
;              Wave must be monotonically increasing or decreasing.
; OPTIONAL INPUT KEYWORDS:
;   Velocity = Scalar velocity of the plasma in km/s (Negative is blue-shifted)
;   Ashift   = Shift the entire spectrum by this amount
;   ioncal   = Specify the ionization balance calculation (0=default)
;   atocal   = Specify the atomic calculation file        (0=default -- see get_atomic)
;   cnorm    = Continuum normalization. Default is cnorm = 1.0
;   elem_num = To specify table of abundances to use (see get_elemabun routine)
;   abun     = Value of abundance [N(Z)/N(H)] -- to override default list returned by get_elemabun
;   noline   = If set, do not calculate the line contribution
;   nocont   = If set, do not calculate the continuum contribution
;   uTDOPPW  = User specified Gaussian Width (FWHM) in Ang.  
;	       This over-rides Td6 and detector width
;   dens     = log10(Dens cm-3).  Only available for S XV (chan=4).
;   newcal   = To force a new calculation by bcs_line
;   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 KEYWORDS:
;   waveout  = The wavelengths of the centers of the output bins = (wave1(1:*)+wave1)/2. 
;   bcs_cont = The continuum contribution only.
;   Ainfo    = Structure giving information about the calculation
;   ion_data = Data structure containing the ionization balance calculation
;   ato_data = Data structure containing the atomic data for line calculation
;   edg_waveout = The wavelength edges of the output bins, if Wave is passed in
;                 this should be the same array
;	   
; MODIFICATION HISTORY:
;    8-oct-93, J. R. Lemen (LPARL), Written
;   27-oct-93, JRL, Fixed Ainfo for the case of Te6 as a vector
;    9-mar-94, JRL, Added the dens=dens keyword
;    9-nov-94, JRL, Added the ion_mult keyword
;   16-jan-95, JRL, Added the rockw keyword
;   19-dec-96, Jmm, Added edg_waveout keyword, to pass out bin edge array, no
;                   alteration of the code needed.
;-
on_error,2 			; Return to calling routine

if n_params() lt 2 then begin
  doc_library,'bcs_spec'
  return,-1
endif

; --- Physical Constants

C_light = 3.e5			; (km/s)    Speed light km/s

; --- Check that the n_elements of Te6 and Td6 are equal

if n_elements(Td6) eq 0 then Td6 = Te6

if n_elements(Te6) ne n_elements(Td6) then begin
  message,'n_elements(Te6) must equal n_elements(Td6) ==>',/cont
  print,'n_elements(Te6) = ',n_elements(Te6)
  print,'n_elements(Td6) = ',n_elements(Td6)
  return,-1
endif 

; --- Check that the channel number is valid
Chan = gt_bsc_chan(chan_struct)			; Make sure chan_struct is a scalar
if n_elements(chan) gt 1 then message,'Chan_struct must be a scalar'

if (Chan lt 1) or (Chan gt 4) then 	$
  message,'Value of Chan is invalid  Must be between 1 and 4 ==>'+string(chan)

; --- Check that Shift and Vel are scalars

if (n_elements(Ashift) gt 1) or (n_elements(Vel) gt 1) then 	$
  message,'AShift and Velocity must be scalars ==>'+$
		string(n_elements(AShift),n_elements(Vel))

; --- Check that the wavelength vector is defined

if n_elements(wave) ne 0 then wave1 = wave else wave1 = gt_bsc_wave(chan_struct)
wave1 = wave1(where(wave1 gt 0))	; gt_bsc_wave returns some 0 wavelengths

if n_elements(wave1) eq 1 then begin
  message,'Wavelength vector must be at least two points',/continue
  message,'Wavelength vector describes start and stop wavelengths of each bin',/noname
  return,-1
endif

dwave = wave1(1:*) - wave1		; Calculate bin widths
if min(sgn(dwave)) ne max(sgn(dwave)) then 	$
  message,'Wavelength vector must monotonically increase or decrease'

wave0 = wave1 + dwave/2			; Wavelength of center of each bin
; ===================
;  wave1 describes the start and stop each bin
;  wave0 describes the wavelength of the center of each bin
;  dwave is the width (Ang) of each bin
; ===================

; --- Set up the parameters for the line broadening 

  TDOPPW = bcs_broad(chan_struct,Td6=Td6,/gauss,ROCKW=ROCKW)	    ; TDOPPW includes instrumental
  if n_elements(uROCKW) ne 0. then ROCKW = uROCKW 		    ; Pass in the rocking curve directly
  if n_elements(uTDOPPW) ne 0. then TDOPPW = uTDOPPW else uTDOPPW=0. ; Pass in the width directly

  if n_elements(cnorm) eq 0 then cnorm = 1.		; continuum normalization

; --- Set up the control parameters ---

if not keyword_set(nocont) then nocont = 0b		; 1   = no continuum calculated
if not keyword_set(noline) then noline = 0b		; 1   = no lines     calculated
if n_elements(Vel) eq 0    then vel = 0.		; vel = bulk velocity (neg is red shift)
if n_elements(Ashift) eq 0 then Ashift = 0.		; Ashift = bulk wavelength shift

; --- Compute the continuum
if not keyword_set(nocont) then begin
  cont = reform(transpose(conflx(Te6,wave0,2)))		; ph s-1 A-1
  bcs_cont = mk_bcs_spec(0., 0.,wave1,wave0,cont,vel=vel,Ashift=Ashift) * cnorm
  flux = bcs_cont
endif else flux = fltarr(n_elements(wave0),n_elements(Te6))

; --- Compute the Line Flux

if not keyword_set(noline) then begin	
   bcs_line,Chan,Te6,							$	; Inputs
		wave_line,line_intensity,dens=dens,			$	; Outputs
		ioncal=ioncal,atocal=atocal,abun=abun,elem_num=elem_num,$	; Inputs
		ion_mult=ion_mult, newcal=newcal,			$	; Inputs
		ato_data=ato_data,ion_data=ion_data,Ainfo=Ainfo			; Outputs
  flux = flux + mk_bcs_spec(TDOPPW,ROCKW,wave1,wave_line,line_intensity,$
		vel=vel,Ashift=Ashift)
endif


; --- Set up the Ainfo Information data structure ---

if n_elements(Ainfo) eq 0 then begin	; if /noline is set
  bsc_struct,fit=Ainfo			; Use the BSC .fit structure for atomic information
  if n_elements(Te6) gt 1 then Ainfo = replicate(Ainfo,n_elements(Te6))
endif

if n_elements(Te6) eq 1 then begin	; This is because [Te6] will fail if
  TTe6 = Te6(0) 			;- n_elements(Te6) = 1
  TTd6 = Td6(0)
endif else begin
  TTe6 = Te6	
  TTd6 = Td6
endelse

Ainfo.chan 	= chan			; BCS channel number
Ainfo.Te6 	= TTe6			; Electron Temperature (MK)
Ainfo.Td6 	= TTd6			; Doppler  Temperature (MK)
Ainfo.EM50 	= 1.			; Emission measure
Ainfo.Vel 	= Vel 			; Bulk Velocity (km/s)
Ainfo.Ashift 	= Ashift		; Bulk wavelength shift
Ainfo.cnorm 	= cnorm			; Continuum Normalization 
Ainfo.UTDOPPW	= UTDOPPW		; User supplied Gaussian broadening
Ainfo.nocont 	= nocont		; 1 if No continuum calculation
Ainfo.noline 	= noline		; 1 if No line calculation

return,reform(flux)
end