FUNCTION unbin_spec, file=file, resp=resp, silent=silent, texp=texp, $
                     qual=qual, backscal=backscal, w_u=w_u, w_l=w_l

; ----------------------------------------------------------
;+
; NAME:
;       UNBIN_SPEC
;
; PURPOSE:
;      Load and RGS spectrum from a FITS file and 'unbin' 
;
; AUTHOR:
;      Simon Vaughan (U.Leicester)
;
; CALLING SEQUENCE:
;      spec = UNBIN_SPEC(src, resp)
;
; INPUTS:
;      file     - names of FITS file of spectra (full path)
;      resp     - name of FITS file of response matrix (full path)
;
; OPTIONAL INPUTS:
;      silent   - (logical) switch off on-screen output?
;
; OUTPUTS:
;      spec     - wavelengths of each photon in the spectrum (Angstrom)
;
; OPTIONAL OUTPUTS:
;      qual     - quality flag (0=good, 1=bad) for each channel
;      backscal - BACKSCAL value for each channel from the spectral file
;      w_l      - wavelength (Angstrom) for lower boundary of each channel
;      w_u      - wavelength (Angstrom) for upper boundary of each
;                 channel
;      texp     - Expsure time of the spectrum
;
; DETAILS:
;      Called by RGS_PLOT.
;      Performs an 'unbinning' of finely binned spectral data in a
;      FITS file ('src') given a response matrix ('resp'). Currently
;      the columns of the FITS file that are used are those suitable
;      for processing RGS data, but could be adjusted if needed. The
;      'unbinning' works by redistributing each count in the spectrum
;      within the wavelength range of the channel it was collected
;      in. The redistribution is random and uniform between w[i] and
;      w[i]+dw[i] where w[i] and dw[i] are the lower boundary and
;      width of the wavelength bin of the ith count in Angstroms. This
;      of course assumes the counts are distributed unformly within
;      each spectral channel, an assumption that should be true if the
;      input data are given with small bin widths. In the case of the
;      RGS we recommend using 3400 channels. 
;
;      Here's exactly what happens:
;       - load the data from the FITS file as channel, counts, quality
;       - load the corresponding response file FITS
;       - use channel-energy table to convert spectral channel to energy
;       - convert data wavelength: wavelength, bin width, counts, quality
;       - create 'unbinned' spectrum by assigning each count to a
;          wavelength uniformly distributed in [w, w+dw]
;       - keep track of 'quality' flag
;       
;
; EXAMPLE USAGE:
;          src = 'spec.FIT'
;          resp = 'resp.FIT'
;          spec = UNBIN_SPEC(file=src, resp=resp)
;          plot, HISTOGRAM(spec, BINSIZE=0.05), PSYM=10
;
; PROCEDURES CALLED:
;          SPEC_READFITS [READFITS, TBGET, SXPAR]
;
;
; HISTORY:
;      26/11/2009 - v1.0 - first working version
;      15/12/2009 - v1.1 - added TEXP output
;
; NOTES:
;-
; ----------------------------------------------------------


; options for compilation (recommended by RSI)

  COMPILE_OPT idl2

; watch out for errors

  on_error,0

; ----------------------------------------------------------
; Check arguments

; check number of spectra to process

  if (N_ELEMENTS(file) ne 1) then begin
      print,'** UNBIN_SPEC requies one file for FILE input.'
      return,0
  endif

; check number of response files to process

  if (N_ELEMENTS(resp) ne 1) then begin
      print,'** UNBIN_SPEC requies one file for RESP input.'
      return, 0
  endif

; ----------------------------------------------------------
; Main routine

; load the spectrum from FITS file
; Columns are 1 = channel
;             2 = counts/bin
;             3 = quality (0=good, 1=bad)
;             4 = AREASCAL (0.0 - 1.0)
;             5 = BACKSCAL (arb. units)
;
; In this case we want the counts, quality and BACKSCAL only.

  data = SPEC_READFITS(file, cols=[2,3,5], texp=texp)

; extract the vectors of counts/bin, quality and BACKSCAL

  counts = REFORM(data[*,0])
  qual = REFORM(data[*,1])
  backscal = REFORM(data[*,2])

; Load the channel-energy table from the response matrix

  data = SPEC_READFITS(resp, extension=2, cols=[2,3], /DOUBLE)

; lower/upper energies of channel boundaries

  e_l = REFORM(data[*,0])
  e_u = REFORM(data[*,1])

; convert channel boundary energies to wavelengths (Angstrom)

  apkev = 12.3984428D

  w_l = apkev/e_u
  w_u = apkev/e_l

; channel widths

  dw = w_u - w_l

; total number of counts and channels in the spectrum

  n = TOTAL(counts)
  n_chan = N_ELEMENTS(counts)
  n_chan2 = N_ELEMENTS(e_l)

; cumulative number of counts in channels <i

  cumsum = TOTAL(counts, /CUMULATIVE)

; display data parameters

  if (not KEYWORD_SET(silent)) then begin
      print,'-- Number of channels in data:     ', n_chan
      print,'-- Number of channels in response: ', n_chan2
      print,'-- Number of counts:               ', n
      print,'-- '
  endif

; prepare array for output

  spec = MAKE_ARRAY(n, /DOUBLE)

; form the 'unbinned' spectrum by assigning each count a wavelength
; uniformly distributed in [w_l[i], w_u[i]] where i is the channel
; that the count was collected in.

  for i=0, n_chan-1 do begin

      if (counts[i] eq 0) then continue

      rand = RANDOMU(seed, counts[i])

      j = cumsum[i] - counts[i] 
      k = j + counts[i] - 1
;      print,i,j,k,counts[i], cumsum[i], n
      spec[j:k] = rand * dw[i] + w_l[i]

  endfor

; ----------------------------------------------------------
; return to calling routine

  return, spec

END

; ----------------------------------------------------------
; -                                                        -
; -                                                        -
; -                                                        -
; -                                                        -
; -                                                        -
; -                                                        -
; -                                                        -
; -                                                        -
; -                                                        -
; -                                                        -
; -                                                        -
; -                                                        -
; -                                                        -
; -                                                        -
; -                                                        -
; -                                                        -
; -                                                        -
; -                                                        -
; -                                                        -
; -                                                        -
; -                                                        -
; -                                                        -
; -                                                        -
; -                                                        -
; ----------------------------------------------------------

FUNCTION unbin, x_l, x_u, y, u, du

; ----------------------------------------------------------
;+
; NAME:
;      UNBIN
;
; PURPOSE:
;      Put a binned array onto a finer, even grid
;
; AUTHOR:
;      Simon Vaughan (U.Leicester)
;
; CALLING SEQUENCE:
;      spec = UNBIN(x_l, x_u, y, u, du)
;
; INPUTS:
;      x_l      - (vector) lower boundary of input channels
;      x_u      - (vector) upper boundary of input channels
;      y        - (vector) input data
;      u        - (vector) lower boundary of output channels
;      du       - (float) width of output channels
;
; OPTIONAL INPUTS:
;      none
;
; OUTPUTS:
;      z        - the output array 
;
; OPTIONAL OUTPUTS:
;      none
;
; DETAILS:
;      Map an input array onto an finely gridded output array. The
;      input array is a vector of values y[i] in channels that cover
;      the range x_l[i] - x_u[i]. This is to be mapped onto an array
;      z whose channels cover the ranges u[j] - u[j]+du. Where a
;      single channel of the output z spans two channels of the input
;      y, the output is a simple weighted mean of the two channels of
;      y that contribute. The input array may be unevenly gridded,
;      i.e. dx = x_l[i] - x_u[i] may vary from channel to channel.
;
;      We do assume both the input and output channels (x_l, x_u and
;      u) are in acending order.
;
;       
;
; EXAMPLE USAGE:
;      spec = UNBIN(x_l, x_u, y, u, du)
;
; PROCEDURES CALLED:
;      none
;
;
; HISTORY:
;      22/12/2009 - v1.0 - first working version
;
; NOTES:
;-
; ----------------------------------------------------------

; options for compilation (recommended by RSI)

  COMPILE_OPT idl2

; watch out for errors

  on_error,0

; ----------------------------------------------------------
; Check arguments

; check number of channels to process

  n_y = N_ELEMENTS(y)
  n_x1 = N_ELEMENTS(x_l)
  n_x2 = N_ELEMENTS(x_u)

  if (n_y EQ 0) then begin
      print,'** UNBIN requies more data for input Y.'
      return,-1
  endif

; check the input arrays are all the same size

  if (n_x1 NE n_y or n_x2 NE n_y) then begin
      print,'** Arrays X_L, X_U and Y of differing size in UNBIN.'
      return,-1
  endif

; set the default output array if not given as input

  dx = min(x_u - x_l)
  yrange = x_u[n_y-1] - x_l[0]

  if KEYWORD_SET(du) then begin
      if (du GT dx) then begin
          print,'** Output array is less fine than input array in UNBIN.'
          return,-1
      endif
  endif

  n_z = N_ELEMENTS(u)
  if (n_z EQ 0) then begin

      if KEYWORD_SET(du) then begin
          n_z = yrange / du
      endif else begin
          n_z = 5 * n_y
          du = yrange / n_z
      endelse

      u = indgen(n_z) * du + x_l[0]

  endif
  
  if (not KEYWORD_SET(du)) then du = u[1] - u[0]

; ----------------------------------------------------------
; Main routine

; create output array

  z = MAKE_ARRAY(n_z)

; upper boundary of output channels

  u_u = u + du

; loop over each input channel, i

  for i = 1, n_y-1 do begin

          mask = WHERE(u_u GT x_l[i] and u_u LT x_u[i], count)
       
 ; if no output channels are included move along...

          if (count EQ 0) then CONTINUE

; find the highest and lower output channels included

          ch_1 = mask[0]
          frac_1 = (u_u[ch_1] - x_l[i]) / du

          z[ch_1] = frac_1 * y[i] + (1-frac_1) * y[i-1]

          if (count GT 1) then begin
              ch_N = mask[1:count-1]
              z[ch_N] = y[i]
          endif

  endfor

; ----------------------------------------------------------
; return to calling routine

  return, z

END

; ----------------------------------------------------------
; -                                                        -
; -                                                        -
; -                                                        -
; -                                                        -
; -                                                        -
; -                                                        -
; -                                                        -
; -                                                        -
; -                                                        -
; -                                                        -
; -                                                        -
; -                                                        -
; -                                                        -
; -                                                        -
; -                                                        -
; -                                                        -
; -                                                        -
; -                                                        -
; -                                                        -
; -                                                        -
; -                                                        -
; -                                                        -
; -                                                        -
; -                                                        -
; ----------------------------------------------------------

PRO rgs_plot, obs=obs, src=src, bkg=bkg, resp=resp, cent=cent, vlim=vlim, $
              z=z, ps=ps, dv=dv 

; ----------------------------------------------------------
;+
; NAME:
;       RGS_PLOT
;
; Purposex:
;      Plot RGS data in velocity space
;
; AUTHOR:
;      Simon Vaughan (U.Leicester)
;
; CALLING SEQUENCE:
;      rgs_plot, filename="file.qdp"
;
; INPUTS:
;      src      - list of names of FITS files of spectra (full path)
;      resp     - list of names of FITS files of response matrices (full path)
;
; OPTIONAL INPUTS:
;      bkg      - list of names of FITS files of background spectra (full path)
;      cent     - rest-frame wavelengths of lines to plot (Angstrom)
;      vlim     - velocity limits to plot (in km/s), e.g. [-5000, 2000]
;      z        - redshift of source (default = 0.0)
;      ps       - TRUE/FALSE: produce a PostScript file output (velocity.ps)
;      dv       - (scalar) velocity bin width (km/s)
;
; OUTPUTS:
;      velocity.qdp - QDP file showing final results
;
; OPTIONAL OUTPUTS:
;      none
;
; DETAILS:
;      Produce a 'composite' line profile in velocity space from raw
;      RGS spectral files (FITS). The input is a set of files for the
;      source spectra, the (optional) background spectra, and response
;      matricies. Each source spectral file should have a
;      corresponding response (and background if used), and should be
;      the non-background subtracted RGS spectrum in high resolution
;      (e.g. 3400 channels). If there are multiple spectra, e.g. one
;      for each of RGS1 and RGS2, or multiple observations, these may
;      be combined. For example:
; 
;        src = ['SRC_RGS1.FIT','SRC_RGS2.FIT']
;        bkg = ['BKG_RGS1.FIT','BKG_RGS2.FIT']
;        resp = ['MAT1.FIT','MAT2.FIT']
;        rgs_plot,src=src, bkg=bkg, resp=resp, cent=[16.006], z=0.002336, vlim=[-6000,2000]
;
;      We cannot do the more simple trick of summing the counts in
;      each channel over all spectra because the channel-energy
;      conversion for each spectrum can be (and usually is) different
;      for each RGS spectrum.  A similar problem applies to the
;      BACKSCAL and QUALITY information, which is channel
;      dependent. The quality flag shows whether the channel is good
;      or not (e.g. bad due to bad column or missing chip) and this
;      multiplied by the overall exposure gives an effective exposure
;      per channel. The BACKSCAL values are also per channel. As the
;      channel boundaries (in wavelengths) are different from spectrum
;      to spectrum, these cannot simply be averaged by averaging the
;      data in each channel. What we do is re-sample the exposure and
;      BACKSCAL vectors onto a much finer grid of wavelengths
;      (specified if necesaary using the XXX keywords) that is the
;      same for all spectral files.
;
;      Each spectrum (and background) is 'unbinned' using the
;      UNBIN_SPEC function. Briefly, this randomised the observed counts within
;      the wavelength range of each spectral channel, to produce a
;      list of unbinned 'observed' wavelengths for each count. The
;      counts for all the source spectra are then combined (likewise
;      the background if used). The channel-wavelength conversion
;      comes from the response files. 
;
;      The photon wavelengths are then de-redshifted. For each
;      rest wavelength specified with the CENT input the photon
;      wavelengths are converted to velocity shifts (conversion
;      assuming v << c), and those counts in the velocity range
;      specified by the VLIM input are collected. This is repeated for
;      each line, and for the source and background spectra
;      individually. The velocities of each photon from each line are
;      then collected together (likewise for background if used) and
;      are binned in even velocity units (specified by the DV
;      input). For example, one may wish to examine the velocity range
;      -3000,+3000 km/s around the lines at 16.006. This is done by
;      setting cent=[16.006, 18.967] and vlim=[-3000,3000] and then
;      specifiying e.g. dv=300 (km/s) for the final binning of the
;      velocity spectrum. If a background is specified the background
;      is subtracted from the source spectrum. The resulting spectrum
;      is in units of counts per velocity bin, with velocity in units
;      of km/s. Errors come from propagating the usual Poisson
;      variances (e.g. error = sqrt(N1+N2) for N1 source counts and N2
;      background counts). The result is a composite of all the counts
;      for each spectrum and for each line (shifted into velocity
;      shifts relative to that line). 
;
;      The source and background region sizes (from the BACKSCAL
;      vector) and the exposure time per channel are then
;      combined. The area is calculated for each channel by an
;      exposure time weighted average of the areas of the contributing
;      spectral files. The ratio of source and background regions
;      provides the wavelength dependent background scaling factor
;      required for the background subtraction.
;
; EXAMPLE USAGE:
;
;
; PROCEDURES CALLED:
;          SPEC_READFITS, UNBIN_SPEC, UNBIN, PS_OPEN, PS_CLOSE,
;          WRITE_TABLE, PLOT_ERR, READFITS, TBGET, SXPAR
;
;
; HISTORY:
;      12/11/2009 - v0.1 - first working version
;      26/11/2009 - v0.2 - modified to work on raw (total, background) 
;                          fits files, multiple spectra, and combine
;                          the counts from each spectrum and each line
;                          then bin to produce a single "composite"
;                          velocity profile in counts/bin against
;                          velocity. 
;      22/12/2009 - v0.3 - extended to cope with differing
;                          wavelength-channel conversions in each
;                          spectrum (using the UNBIN function).
;      27/02/2012 - v0.4 - Fixed minor bug in de-redshift conversion
;      14/03/2012 - v0.5 - extended wavelength range to 5.5 - 38 A
;                          replaced v(z) with relativistic formula
;
; NOTES:
;
;-
; ----------------------------------------------------------

; options for compilation (recommended by RSI)

  COMPILE_OPT idl2

; watch out for errors

  on_error,0


; ----------------------------------------------------------
; Use a look-up table to map the numbers in the OBS input to to the
; files of the NGC 4051 observation kept on SVs data disc. 

  if (not KEYWORD_SET(src)) then begin

      if (N_ELEMENTS(obs) eq 0) then obs=INDGEN(15)+1
      
      l_src = READ_TABLE('/data/49/sav2/xmm/ngc4051/srclist.txt', /TEXT)
      l_bkg = READ_TABLE('/data/49/sav2/xmm/ngc4051/bkglist.txt', /TEXT)
      l_resp = READ_TABLE('/data/49/sav2/xmm/ngc4051/matlist.txt', /TEXT)
      
      src = l_src[(obs-1)*2]
      bkg = l_bkg[(obs-1)*2]
      resp = l_resp[(obs-1)*2]
      
      src = [src,l_src[(obs-1)*2+1]]
      bkg = [bkg,l_bkg[(obs-1)*2+1]]
      resp = [resp,l_resp[(obs-1)*2+1]]

  endif

; ----------------------------------------------------------
; Check arguments

; set the centroid if not specified

  if (N_ELEMENTS(cent) eq 0) then cent = 21.602

; number of lines to combine

  n_line = N_ELEMENTS(cent)

; set the velocity limit if not specified

  if (N_ELEMENTS(vlim) eq 0) then begin
      print,'** VLIM parameter not given in RGS_PLOT. Using default setting.'
      vlim = [-5000.0, 5000.0]
  endif
  if (N_ELEMENTS(vlim) eq 1) then begin
      print,'** Only one value of VLIM parameter given in RGS_PLOT.'
      print,'** Will use [-VLIM, +VLIM] as suitable range.'
      vlim = [-vlim, vlim]
  endif

; set the redshift if not specified

  if (N_ELEMENTS(z) eq 0) then z = 0.0

; set the velocity bin width (km/s) if not specified

  if (N_ELEMENTS(dv) ne 1) then dv = 200.0

; number of spectra to process

  n_file = N_ELEMENTS(src)
  if (n_file eq 0) then begin
      filename=''
      read,'-- Enter file name (ENTER to list current directory): ',filename
      if (filename eq '') then begin
          list = findfile()
          print, list
          read,'-- Enter file name: ',src
      endif
      n_file = N_ELEMENTS(src)
  endif

; make sure number of background spectra matches number of spectra 

  if KEYWORD_SET(bkg) then begin
      n_bkg = N_ELEMENTS(bkg)
      if (n_bkg ne n_file) then begin
          print,'** Number of background files \(BKG\) is different from'
          print,'** number of spectral files \(SRC\) in RGS_PLOT.'
          return
      endif
  endif

; make sure number of response matrices matches number of spectra 

  n_resp = N_ELEMENTS(resp)
  if (n_resp ne n_file) then begin
      print,'** Number of response files \(RESP\) is different from'
      print,'** number of spectral files \(SRC\) in RGS_PLOT.'
      return
  endif

; ----------------------------------------------------------
; speed of light (km/s)

  c = 299792.458D

; ----------------------------------------------------------
; Loop over each spectrum file (FITS). For each one 'unbin' the
; spectrum using the UNBIN_SPEC function (see above). The result is a
; list of wavelengths for each count in the original (binned)
; spectrum. 
;
; Also calculte the total exposure time, summed over all spectra,
; along with the exposure time per channel (exposure multiplied by
; quality flag), exposure time-weighted area. These are mapped onto a
; fine grid of wavelengths.

; contruct output wavelength grid, u
; XXX specify options XXX

  u_min = 5.5
  u_max = 38.0
  u_range = (u_max - u_min)
  du = 0.002
  n_u = round(u_range/du)
  u = indgen(n_u)*du + u_min

; loop over the source files

  totime = MAKE_ARRAY(n_u)
  src_area = MAKE_ARRAY(n_u)

  src_t = 0.0
  for i = 0, n_file-1 do begin
      print,'-- Processing ',src[i]
      spec_i = UNBIN_SPEC(file=src[i], resp=resp[i], /SILENT, w_l=w_l, w_u=w_u, $
                          texp=texp, qual=qual, backscal=backscal)

     if (i eq 0) then begin
          src_spec = spec_i
      endif else begin
          src_spec = [src_spec,spec_i]
      endelse

      qual = round(1 - qual)
      exptime = qual * texp
      src_t = src_t + texp

      u_time = UNBIN(w_l, w_u, exptime, u, du)
      u_area = UNBIN(w_l, w_u, backscal, u, du)

      src_area = src_area + u_area * u_time
      totime = totime + u_time

; check things are accumulating correctly
;
;     plot, u, totime, psym=10
;     oplot, u, u_time*1.1, psym=10, linestyle=3
;     oplot, w_l, exptime, psym=10, linestyle=4

  endfor

; Expsure time and (time weighted) area per channel for source

  src_exp = totime
  mask = WHERE(src_exp gt 0.0, COMPLEMENT=cmask, count)
  if (count GT 0) then src_area[mask] = src_area[mask] / src_exp[mask]
  if (count LT n_u) then src_area[cmask] = 0.0

; ----------------------------------------------------------
; Do the same for the background (if required)

  totime = MAKE_ARRAY(n_u)
  bkg_area = MAKE_ARRAY(n_u)

  bkg_t = 0.0
  if KEYWORD_SET(bkg) then begin
      for i = 0, n_file-1 do begin
          print,'-- Processing ',bkg[i]
          spec_i = UNBIN_SPEC(file=bkg[i], resp=resp[i], /SILENT, w_l=w_l, w_u=w_u, $
                              texp=texp, qual=qual, backscal=backscal)

          if (i eq 0) then begin
              bkg_spec = spec_i
          endif else begin
              bkg_spec = [bkg_spec,spec_i]
          endelse
          
          qual = round(1 - qual)
          exptime = qual * texp
          bkg_t = bkg_t + texp
          
          u_time = UNBIN(w_l, w_u, exptime, u, du)
          u_area = UNBIN(w_l, w_u, backscal, u, du)
          
          bkg_area = bkg_area + u_area * u_time
          totime = totime + u_time
          
      endfor

; Expsure time and (time weighted) area per channel for background

      bkg_exp = totime
      mask = WHERE(bkg_exp gt 0.0, COMPLEMENT=cmask, count)
      if (count GT 0) then bkg_area[mask] = bkg_area[mask] / bkg_exp[mask]
      if (count LT n_u) then bkg_area[cmask] = 0.0

; The source/background region scaling factor, per channel

      scale = MAKE_ARRAY(N_ELEMENTS(src_area))
      mask = WHERE(bkg_area gt 0.0, COMPLEMENT=cmask)
      scale[mask] = src_area[mask] / bkg_area[mask]
      scale[cmask] = 0.0

  endif

; ----------------------------------------------------------
; correct wavelengths for redshift

  src_spec = src_spec / (z + 1.0D)

  if KEYWORD_SET(bkg) then  bkg_spec = bkg_spec / (z + 1.0D)
  
; open the PS device if requested

  if keyword_set(ps) then ps_open, "velocity.ps", /COLOUR

; ----------------------------------------------------------
; make a few quick diagnostic plots

  dw = 0.1
  hist = HISTOGRAM(src_spec, BINSIZE=dw, locations=binw)
  plot, binw, hist, PSYM=10, xtitle="Wavelength (!6!sA!r!u!9 %!6!n)", $
    ytitle="counts/bin"

  if KEYWORD_SET(bkg) then begin
      hist = HISTOGRAM(bkg_spec, BINSIZE=dw, locations=binw)
      oplot, binw, hist, PSYM=10
  endif

  plot, u, src_area, PSYM=10, xtitle="Wavelength (!6!sA!r!u!9 %!6!n)", ytitle="Area"

  if KEYWORD_SET(bkg) then begin
      oplot, u, bkg_area, PSYM=10
  endif

  if KEYWORD_SET(bkg) then begin
      plot, u, scale, PSYM=10, xtitle="Wavelength (!6!sA!r!u!9 %!6!n)", ytitle="src/bkg area"
  endif

  plot, u, src_exp, PSYM=10, xtitle="Wavelength (!6!sA!r!u!9 %!6!n)", ytitle="Exposure"


; ----------------------------------------------------------
; loop over each line to consider
; For each one collect the counts within vlim[0] <= v <= vlim[1]
; for both source and background. 

  for i = 0, n_line-1 do begin

; calcuate Doppler shift z = (w[obs] - w[lab]) / w[lab] 
; 14/03/2012 - replaced with relativistic formula

     delta = src_spec/cent[i]
     delta2 = delta * delta
     src_v = (delta2 - 1.0D) / (delta2 + 1.0D) * c

     if KEYWORD_SET(bkg) then begin
;        bkg_v = (bkg_spec/cent[i] - 1.0D) * c
         delta = bkg_spec/cent[i]
         delta2 = delta * delta
         bkg_v = (delta2 - 1.0D) / (delta2 + 1.0D) * c
     endif

; extract src and bkg data data from velocities in range vlim[0] -
; vlim[1] only

      mask = WHERE(src_v LE vlim[1] and src_v GE vlim[0], count)
      src_i = src_v[mask]

      if KEYWORD_SET(bkg) then begin
          mask = WHERE(bkg_v LE vlim[1] and bkg_v GE vlim[0], count)
          bkg_i = bkg_v[mask]
      endif

; include the current line in the composite array

      if (i eq 0) THEN BEGIN
          src_n = src_i
          if KEYWORD_SET(bkg) then  bkg_n = bkg_i
      endif else begin
          src_n = [src_n,src_i]
          if KEYWORD_SET(bkg) then bkg_n = [bkg_n,bkg_i]
      endelse

  endfor

; sort counts from all lines into velocity order

  indx = SORT(src_n)
  src_n = src_n[indx]

  if KEYWORD_SET(bkg) then begin
      indx = SORT(bkg_n)
      bkg_n = bkg_n[indx]
  endif

; ----------------------------------------------------------
; Bin the counts by wavelengths to get a binned spectrum
; in ct/bin units

  hist = HISTOGRAM(src_n, BINSIZE=dv, locations=binv, min=vlim[0], max=vlim[1])
  src_biny = hist

  if KEYWORD_SET(bkg) then begin
      hist = HISTOGRAM(bkg_n, BINSIZE=dv, min=vlim[0], max=vlim[1])
      bkg_biny = hist
  endif

  n_bin = N_ELEMENTS(binv)
   
; ----------------------------------------------------------
; loop over each line to consider. Find channels that 
; correspond to vlim[0] <= v <= vlim[1] for each line.
; Calculate the mean exposure time and src/bkg area within
; each velocity bin. (The mean is the average over all channels, from
; all lines, contributing to a given velocity bin. Where only a
; factional of a channel falls with a given velocity bin, the
; contribution to the mean is weighted accordingly.)

  src_exp_vbin = MAKE_ARRAY(n_bin)
  src_area_vbin = MAKE_ARRAY(n_bin)
  bkg_exp_vbin = MAKE_ARRAY(n_bin)
  bkg_area_vbin = MAKE_ARRAY(n_bin)

; Fixed bug in these two lines - 27/02/2012

  u_l = u / (1.0D + z)
  u_u = (u+du) / (1.0D + z)

  for i = 0, n_line-1 do begin

; convert channel boundary wavelengths to velocity in units (km/s)
; NB: because the RGS channels are ordered in decreasing energy,
; i.e. increasing wavelength, the channel-velocity conversion
; produces velocities in acending order.

; 14/03/2012 - replaced with relativistic formula

      delta = u_l/cent[i]
      delta2 = delta * delta
      v_l = (delta2 - 1.0D) / (delta2 + 1.0D) * c
      delta = u_u/cent[i]
      delta2 = delta * delta
      v_u = (delta2 - 1.0D) / (delta2 + 1.0D) * c

;      v_l = (u_l/cent[i] - 1.0D) * c
;      v_u = (u_u/cent[i] - 1.0D) * c

; inner loop: over each velocity bin within each line

      for j = 0, n_bin-1 do begin

; extract channels from velocities in range of Jth velocity bin

          lo_v = binv[j]
          hi_v = binv[j] + dv
          mask = WHERE(v_u GT lo_v and v_l LT hi_v, count)

; how much of the first and last of these channels are actually
; included in the velocity range of Jth velocity bin?

          if (count EQ 0) then CONTINUE

          ch_1 = mask[0]
          frac_1 = (v_u[ch_1] - lo_v) / (v_u[ch_1] - v_l[ch_1])

          ch_N = mask[count-1]
          frac_N = (hi_v - v_l[ch_N]) / (v_u[ch_N] - v_l[ch_N])

          n_channels = count - 2 + frac_1 + frac_N
 
          src_exp_vbin_i = TOTAL(src_exp[mask]) - (1-frac_1)*src_exp[ch_1] - $
            (1-frac_N)*src_exp[ch_N] 

          src_area_vbin_i = TOTAL(src_area[mask]) - (1-frac_1)*src_area[ch_1] - $
            (1-frac_N)*src_area[ch_N] 
 
; check the channels are being assigned correctly
;          print,i,j,count,lo_v,ch_1,hi_v,ch_N,frac_1,frac_N,n_channels,src_exp_vbin[j]

; XXX need to average properly over lines (i...) XXX

          src_exp_vbin[j]  = src_exp_vbin[j] + (src_exp_vbin_i / n_channels)
          src_area_vbin[j] = src_area_vbin[j] + src_area_vbin_i

          if KEYWORD_SET(bkg) then begin
              bkg_exp_vbin_i = TOTAL(bkg_exp[mask]) - (1-frac_1)*bkg_exp[ch_1] - $
                (1-frac_N)*bkg_exp[ch_N] 

              bkg_area_vbin_i = TOTAL(bkg_area[mask]) - (1-frac_1)*bkg_area[ch_1] - $
                (1-frac_N)*bkg_area[ch_N] 

              bkg_exp_vbin[j]  = bkg_exp_vbin[j] + (bkg_exp_vbin_i / n_channels)
              bkg_area_vbin[j] = bkg_area_vbin[j] + bkg_area_vbin_i
          endif

      endfor

  endfor

; ----------------------------------------------------------
; convert to count rate and then subtract the (area scale corrected)
; background spectrum 
; Calculate the error by propogating the Poisson [sqrt(N)] error on the
; scaled source and background counts


  net_biny = MAKE_ARRAY(n_bin)
  bkg_rate = MAKE_ARRAY(n_bin)
  err      = MAKE_ARRAY(n_bin)
  backscal = MAKE_ARRAY(n_bin)

  if KEYWORD_SET(bkg) then begin
      mask = WHERE(src_exp_vbin GT 0 and bkg_area_vbin GT 0 and bkg_exp_vbin GT 0)

      backscal[mask] = src_area_vbin[mask] / bkg_area_vbin[mask]
      bkg_rate[mask] = bkg_biny[mask]/bkg_exp_vbin[mask] * backscal[mask]
      net_biny[mask] = src_biny[mask]/src_exp_vbin[mask] - bkg_rate[mask]

      err[mask] = sqrt(src_biny[mask]/src_exp_vbin[mask]^2 + $
                       bkg_biny[mask]/bkg_exp_vbin[mask]^2*backscal[mask]^2)

  endif else begin
      mask = WHERE(src_exp_vbin GT 0)
      net_biny[mask] = src_biny[mask] / src_exp_vbin[mask]
      err[mask] = sqrt(src_biny[mask]) / src_exp_vbin[mask]
  endelse

; ----------------------------------------------------------
; plot the resulting histogram(s)
;

  title = "Wavelengths"
  for i = 0, n_line-1 do begin
      title = title + STRING(cent[i], FORMAT='(f7.3)')
      title = title +"!6!sA!r!u!9 %!6!n "
  endfor

  yrange = [0.0, 1.2*max(src_biny)]

  plot, binv, src_biny, PSYM=10, xrange=vlim, /xstyle, yrange=yrange, $ 
    /ystyle, xtitle="Velocity (km s!U-1!N)", ytitle="counts/bin", title=title
    
  yrange = [0.0, 1.2*max(net_biny)]

  plot, binv, net_biny, PSYM=10, xrange=vlim, /xstyle, yrange=yrange, $ 
    /ystyle, xtitle="Velocity (km s!U-1!N)", ytitle="counts/s/bin", title=title
    
  plot_err, binv, net_biny, err

; plot the exposure time per channel from zero to max
;
  y = src_exp_vbin/max(src_exp_vbin)*yrange[1]*0.975
  oplot, binv, y, PSYM=10, THICK=10, COLOR=100
;
;  if KEYWORD_SET(bkg) then begin
;      y =  backs/max(backs)*yrange[1]*0.95
;      oplot, v_l, y, PSYM=10, LINESTYLE=3
;  endif

; plot the area scaling per channel from zero to max

  y = src_area_vbin/max(src_area_vbin)*yrange[1]*0.95
  oplot, binv, y, PSYM=10, THICK=10, COLOR=160

  if KEYWORD_SET(bkg) then begin
      oplot, binv, bkg_rate, PSYM=10, linestyle=3
  endif

; close the PS device if requested

  if keyword_set(ps) then ps_close

; ----------------------------------------------------------
; Make a 'pretty' plot

; xxxxxxxxxxxxxxxxxxxxxx

; open the PS device if requested

  if keyword_set(ps) then ps_open, "velocity-pretty.ps", charsize=1.4, position=[0.1,0.1,0.9,0.9]

  title = "Wavelengths"
  for i = 0, n_line-1 do begin
      title = title + STRING(cent[i], FORMAT='(f7.3)')
      title = title +"!6!sA!r!u!9 %!6!n "
  endfor

  yrange = [0.0, 1.2*max(net_biny)]

  plot, binv, net_biny, PSYM=10, xrange=vlim, /xstyle, yrange=yrange, $ 
    /ystyle, xtitle="Velocity (km s!U-1!N)", ytitle="counts s!U-1!N bin!U-1!N", title=title
    
  plot_err, binv, net_biny, err

  if KEYWORD_SET(bkg) then begin
      oplot, binv, bkg_rate, PSYM=10, linestyle=3
  endif

; close the PS device if requested

  if keyword_set(ps) then ps_close

; ----------------------------------------------------------
; open a file for the output
; and copy the QDP header information

  openw,lun,"velocity.qdp",/get_lun
    printf, lun, "READ SERR 2"
    printf, lun, "!    velocity            flux          error"
    printf, lun, "!       (km/s)       (counts)       (counts)"
  free_lun,lun

; store the output in a file
; after clipping to required velocity range
; note, the velocity is the CENTRE of the bin

  v = binv + 0.5*dv
  y = net_biny
  dy = err

  data = TRANSPOSE([[v],[y],[dy]])

  WRITE_TABLE, data, "velocity.dat"

  SPAWN, 'cat velocity.dat >> velocity.qdp'

; ----------------------------------------------------------
; open a file for the output
;
; ----------------------------------------------------------

END