;****************************************************************************
pro combine,infile,inrec,outfile,outrec,average=average,noscale=noscale, $
            sumup=sumup,examine=examine,helpme=helpme,stp=stp,edit=edit, $
            nosave=nosave             ;,etype=etype,sigma=sigma
if keyword_set(helpme) then begin
   print,' '
   print,'* COMBINE - combine 3 or more spectra with median filtering'
   print,'*           averages 2 spectra'
   print,'*    calling sequence: COMBINE,infile,inrec,outfile,outrec'
   print,'*       INFILE:  input .ICD file'
   print,'*       INREC:   array of record numbers to combine'
   print,'*       OUTFILE: output .ICD file, default=INFILE'
   print,'*       OUTREC:  output record, default=next free record'
   print,'*'
   print,'*    KEYWORDS:'
   print,'*       AVERAGE: use averaging, not median filtering'
;   print,'*       ETYPE:   set to force ETYPE when .ICD file is created'
   print,'*       EXAMINE: examine data at all stages.'
   print,'*       NOSCALE: do not scale flux vectors'
;   print,'*       SIGMA:   sigma limit for median filtering
   print,'*       SUMUP:   sum all flux, no scaling - weights by exp. time'
   print,' '
   return
   endif
;
if keyword_set(edit) then examine=1
if n_elements(infile) eq 0 then begin
   infile=''
   read,' Enter name of .ICD file: ',infile
   endif
if n_elements(inrec) eq 0 then begin
   inrec=-1 & ii=0
   print,' Enter list of record numbers, -1 to end'
   while ii ge 0 do begin
      read,'>>',ii
      inrec=[inrec,ii]
      endwhile
   k=where(inrec ge 0,nk)
   if nk lt 2 then begin
      print,' Fewer than 2 records specified - returning'
      return
      endif
   inrec=inrec(k)
   endif
if n_elements(outfile) eq 0 then outfile=infile
if n_elements(outrec) eq 0 then outrec=-1
nrec=n_elements(inrec)
if nrec lt 2 then begin
   print,' Fewer than 2 records specified - returning'
   return
   endif 
;
gdat,infile,h,w,f,e,inrec(0)
if keyword_set(examine) then plot,w,f,/ynoz
ft=f*h(5)                   ;flux
np=h(7)
kp=indgen(np/2)+np/4            ;middle half of vector for scaling
tf=total(f(kp))
h0=h
h0(9)=1
epstype=h(33)
darr=fltarr(np,nrec) & earr=darr
darr(0,0)=f & earr(0,0)=e
scalef=fltarr(nrec)+1.
;
for i=1,nrec-1 do begin
   gdat,infile,h,w,f,e,inrec(i)
   ft=ft+f*h(5)                   ;total flux
   scalef(i)=tf/total(f(kp))
   darr(0,i)=f & earr(0,i)=e
   h0(5)=h0(5)+h(5)                   ;increment time
   h0(9)=h0(9)+1
   if keyword_set(examine) then begin
      print,' Record',i,':  scale factor= ',scalef(i),' Accum time= ',h0(5)
      if not keyword_set(noscale) then sf=scalef(i) else sf=1.
      oplot,w,f*sf
      endif
   endfor   
if keyword_set(edit) then stop,' Edit data now'
;
if (epstype eq 0) and (not keyword_set(sumup)) then average=1
if (nrec eq 2) and (not keyword_set(sumup)) then average=1
;
if keyword_set(noscale) then scale=0 else scale=1         ;scaling turned off
;
if scale eq 1 then for i=1,nrec-1 do begin                ;scale flux vector
   f=darr(*,i)*scalef(i)
   darr(0,i)=f
   endfor
fbar=sum(darr,1)/nrec     ;average flux
;
if keyword_set(examine) then begin
   if not keyword_set(sumup) then sumup=0
   if not keyword_set(average) then average=0
   if not keyword_set(noscale) then noscale=0
   print,'Average:',average,' Sumup:',sumup,' Noscale:',noscale,' Etype:',epstype
   endif
;
case 1 of
   keyword_set(sumup):   f=ft/h0(5)             ;total flux
   keyword_set(average): f=fbar                 ;average flux
   else: begin                                  ;median filter
      f=medarr(darr)                           ;median array
;      if not keyword_set(sigma) then sigma=5.      
; csig3,sigma,nrec,darr,earr,f,expo
;k=where(expo eq 0,nzero)
;      carr=narr/(carr > .001)            ;sigma
;      narr=narr/float((expo > 1))      ;counts per second
;      if nzero gt 0 then begin         ;no exposure - average adjacent points
;      print,'nzero,element =',nzero,k
;         for i=0,nzero-1 do begin
;            narr(k(i))=(narr(k(i)-1)+narr(k(i)+1))/2.
;            carr(k(i))=(carr(k(i)-1)+carr(k(i)+1))/2.
      end
   endcase
;
case 1 of             ;coadd errors
   epstype eq 30: begin    
      e=sqrt(sum(earr*earr,1))         ;s/n vector
      end
   else: begin
      e=fltarr(np)
      for i=0,np-1 do e(i)=f(i)/(stddev(darr(i,*))>1.e-25)   ;standard deviation
      h0(33)=30
      end
   endcase
;
if keyword_set(examine) then begin
   gc,11
   oplot,w,f,color=55
   erbar,2,f,f/e,w,color=15
   wshow
;   stp=1
   endif
;
if not keyword_set(nosave) then kdat,outfile,h0,w,f,e,outrec    ;,epstype=etype
;
if keyword_set(stp) then stop,'COMBINE>>>'
return
end