pro cross,arr,rel_pos
;+
;     Routine to do simple numerical (ie non-transform) cross correlation 
;     of a set of spectra relative to one of the spectra within the set.
;
;  Input:
;          arr     - 2-d array containing the set of spectra (eg DATA from
;                    test_rd)  arr(bins,spectra)
;
;  Output:
;         rel_pos      - position of each spectrum relative to template (bins)
;
;
;  Use:
;          IDL> cross, array, posn
;
;  Prompts:
;             Which is template? >   give spectrum number which is to be used
;                                    as the template ie against which others
;                                    matched
;
;             What range of spectra to use? >  give the start and end numbers
;                                              for the spectra to be used
;
;             Which bins? >          give the range of bins to be used in
;                                    the correlation
;
;             Number of correlation steps? >   give number of shifts to be
;                                              used in match. The spectra are
;                                              shifted by +/- half this value
;                                        
;
;   CDP     March 92
;
;-
;
;  size of ouput array
;
s = size(arr)
n25 = s(2) - 1
;
;  prompt for parameters
;
read,'Which is the template? > ',tem
read,'What range of spectra to use? > ',n1,n2
read,'Which bins? > ',b1,b2
read,'Number of correlation steps? > ',nc
ans = ' '
read,'Plot as we go? > ',ans
ans = strupcase(ans)

;
;  template
;
template = arr(*,tem)

;
;  cc function
;
cc  = fltarr(nc)
;
;  only use peak +/- 2
;
xc = findgen(5)
;
;  output array
;
rel_pos = fltarr(n25+1)
;
;  range of shifts
;
nc2 = fix(nc/2)
;
;  shift and accumulate
;
for j=n1,n2 do begin
   for i=-nc2,nc2 do begin
      cc(i+nc2) = total(arr(b1:b2,j) * shift(template(b1:b2),i))/1000.
   endfor
   p = where(cc eq max(cc))
   peak = p(0)
   if peak ge 2 and peak le (n_elements(cc)-3) then begin
      cc5 = cc(peak-2:peak+2)
      coff = poly_fit(xc,cc5,2,yy)
      rel_pos(j) = -coff(1)/(2.0*coff(2)) + peak - nc2 - 2
      if ans eq 'Y' then begin
         plot,xc,cc5,psym=5,/ynoz
         oplot,xc,yy
         pp = -coff(1)/(2.0*coff(2))
         oplot,[pp,pp],[0,max(yy)],thick=2
         wait,0.25
      endif
   endif else begin
      rel_pos(j) = 0.0
   endelse
print,j,rel_pos(j)
endfor

end