function mk_bcs_spec,TDOPPW,ROCKW,wave1,wave0,intensity,velocity=vel,Ashift=Ashift
;+
; NAME:
;   mk_bcs_spec
; PURPOSE:
;   Called by bcs_spec to compute a BCS spectrum (ph cm-2 s-1 Ang-1 at 1AU)
; CALLING SEQUENCE:
;   spectrum = mk_bcs_spec(TDOPPW,ROCKW,wave1,wave0,intensity)
;   spectrum = mk_bcs_spec(TDOPPW,ROCKW,wave1,wave0,intensity,vel=vel)
;
; INPUTS:
;   TDOPPW	= FWHM of Gaussian line broadening term     (Ang)
;   ROCKW	= FWHM of Lorentzian line broadening term   (Ang)
;   Wave1	= A vector describing the start wavelength of a vector of bins. (Ang)
;         	  bcs_spec will return n_elements(wave)-1 bins.
; 		  Wave must be monotonically increasing or decreasing.
;   wave0	= Vector of wavelengths corresponding to intensity. (Ang)
;   intensity	= Vector of line or continuum intensities (photons s-1)
;		  If TDOPPW=ROCKW=0, then assume intensity is ph s-1 Ang-1
;
; OPTIONAL INPUT KEYWORDS:
;   velocity	= Line-of-sight velocity in km/s (negative is blue-shifted)
;   Ashift	= Bulk wavelength shift
;
; OUTPUTS:
;   Note:  This routine normalizes to per Ang and puts the source at 1 AU
;   Output spectrum is in photons/cm^2/sec/Ang at 1 AU
;
; MODIFICATION HISTORY:
;  17-sep-93, J. R. Lemen (LPARL), Written
;  29-oct-93, JRL, Fixed a bug with the Vel= option
;-
AU = 1.495979e13	; cm
c_light = 3.e5		; km/s


wave_cent = (wave1(1:*)+wave1)/2.	; Get central wavelength of each bin

num = n_elements(Intensity(0,*))	; Number of spectra (Temperatures)
outflux = fltarr(n_elements(wave_cent),num)

if n_elements(vel) eq 0 then vel = 0.		; Set up the velocity shift
if n_elements(Ashift) eq 0 then Ashift = 0.	; Set up the wavelength shift
swave0 = wave0 * (1. + vel / c_light) + Ashift	; Shifted wavelengths

; If TDOPPW=0 and ROCKW=0 then simply linearly interpolate (no broadening):
; (Also, assume that intensity is already ph cm-2 A-1)

if (total(TDOPPW) eq 0) and (total(ROCKW) eq 0) then begin
    for j=0,num-1 do 		$	; Loop over temperature
	outflux(0,j) = dspline(swave0,intensity(*,j),wave_cent,interp=0)
endif else begin

; --  Apply Line Broadening  --

; Set up TDOPPW as a function of wave0 and num
    if n_elements(TDOPPW) eq 1 then FF = replicate(TDOPPW(0),num) else FF = TDOPPW

     for j=0,num-1 do 		$	; Loop over temperature:    
	outflux(0,j) = wvoigt(wave_cent,swave0,FF(j),ROCKW) # intensity(*,j)
endelse

; Place the source at 1 AU
outflux = outflux / (4. * !pi * AU^2) 	; Photons/cm^2/s/Ang at 1 AU

return,outflux
end