      parameter (maxpix=8192,nwmax=200,ncolmax=10,ncor=50,ncmax=32)
      parameter (nfmax=12)
      real*4 xbuf(maxpix*nwmax),wcoefs(2,nwmax),zcoefs(2,nwmax)    
      real*4 ybuf(maxpix*nwmax),skyspec(maxpix),skycord(maxpix)
      real*4 xcord(maxpix),ycord(maxpix),zcord(maxpix),skymask(maxpix)
      real*4 skyline(maxpix),skycont(maxpix),objcont(maxpix)
      real*4 xp(2),yp(2),objline(maxpix),work1(maxpix),work2(maxpix)
      real*4 clip(nwmax),presky(maxpix),work(maxpix,nwmax)
      real*4 wcord(-ncor:ncor),vcord(-ncor:ncor)
      real*4 veloc(nwmax),wcmax(nwmax),contflux(nwmax),contnoise(nwmax)
      real*4 nullval,skytempl(maxpix),objcorr(maxpix)
      logical anynull
      integer status,inunitref,inunit,outunit,blocksize,hdutype
      integer*4 skynum(nwmax),objnum(nwmax),pronum(nwmax),calnum(nwmax)
      integer naxes(2),fps(nwmax)
      character*80 infile,outfile,errmsg,card
      character*80 listname,argv,title,filenames(nfmax)
      character*256 line
      character*40 tabnam
      character*30 object(nwmax)
      real*4 rmag(nwmax)
      character*3 specnum(nwmax),recnum(nwmax)
c *** FITSIO_SPCSKYSUB does sky subtraction for fibre data from FLAMES
      anynull=.false.
      nullval=0.0
      print *,' '

c *** try for presence of command line arguments
      m=iargc()
c *** if correct number of arguments try and interpret them
      if(m.eq.5)then
        call getarg(1,argv)
        if(index(argv(1:1),'@').gt.0)then
          listname=argv(2:80)
          ioerr=0
          open(unit=4,name=listname,status='old',iostat=ioerr)
          k=0
          do while(ioerr.eq.0.and.k.lt.nfmax)
          k=k+1
          read(4,'(a)',iostat=ioerr) filenames(k)
          enddo
          if(ioerr.ne.0)then
            nf=k-1
          else
            nf=k
          endif
          print *,'No. of files in list =',nf
        else
          nf=1
          filenames(1)=argv
        endif
        call getarg(2,argv)
        outfile=argv
        call getarg(3,argv)
        read(argv,*) wmin                      ! eg. 8250.0
        call getarg(4,argv)  
        read(argv,*) wmax                      !     8950.0
        call getarg(5,argv)
        read(argv,*) clipmax
      else
        print *,'Arguments: infile outfile wmin wmax clipmax'
        print *,' '
        print *,' infile    standard fits output from pipeline'
        print *,' outfile   output file for sky-subtracted spectra'
        print *,' wmin wmax defines analysis wavelength range (A)'
        print *,'           clipped by range of wavelengths in data'
        print *,' clipmax   overrides default plotting range'
        print *,' '
        stop
      endif

c *** open reference file readonly
      status=0
      inunitref=1
      infile=filenames(1)
      call ftopen(inunitref,infile,0,blocksize,status)
      if(status.ne.0)then
        call ftgerr(status,errmsg)
        print *,'FITSIO Error status =',status,': ',errmsg
        stop
      endif
c *** find out how many file extensions
      call ftthdu(inunitref,nhdus,status)
      nmefs=nhdus-1

c *** and set up output file based on reference file
      morekeys=0
      outunit=2
      call clobber(outunit,outfile)            ! remove if exists already
      call ftinit(outunit,outfile,1,status)    ! create new file std blocksize
c *** record PHU header records
      call ftcopy(inunitref,outunit,morekeys,status)
      if(status.ne.0)then
        call ftgerr(status,errmsg)
        print *,'FITSIO Error status =',status,': ',errmsg
        stop
      endif

c *** and get reference data buffer
      call rdbuf_fitsio(inunitref,xbuf,wcoefs,specnum,nxout,nyout)
      npix=nxout*nyout
      if(nyout.gt.maxpix)stop 'Too many pixels'

c *** get table information and flag spectra
      call fitsio_rtable(inunitref,nwmax,ncmax,
     1                   object,fps,rmag,nrows,ier)
      nsky=0
      nobj=0
      ncal=0
      do i=1,nrows
      if(index(object(i),'_Sky').gt.0.or.
     1   index(object(i),'_s').gt.0.or.
     2   index(object(i),'sky').gt.0)then
        nsky=nsky+1
        skynum(nsky)=fps(i)
      else
        if(index(object(i),'CALSIM').gt.0)then
          ncal=ncal+1
          calnum(ncal)=fps(i)
        else
          nobj=nobj+1
          objnum(nobj)=fps(i)
        endif
      endif
      enddo
      print *,'No. of sky and object fibres =',nsky,nobj,ncal

c *** get rest of files spectra and add them in
      inunit=8
      status=0
      do k=2,nf
      infile=filenames(k)
      call ftopen(inunit,infile,0,blocksize,status)
      if(status.ne.0)then
        call ftgerr(status,errmsg)
        print *,'FITSIO Error status =',status,': ',errmsg
        stop
      endif
      call rdbuf_fitsio(inunit,ybuf,zcoefs,recnum,nxoutc,nyoutc)
      if(nxoutc.ne.nxout.or.nyoutc.ne.nyout)stop 'wrong size'
      call ftclos(inunit,status)
      if(status.ne.0)then
        call ftgerr(status,errmsg)
        print *,'FITSIO Error status =',status,': ',errmsg
        stop
      endif
      do i=1,nxout
      do j=1,nyout
      jj=(j-1)*nxout
      xbuf(i+jj)=xbuf(i+jj)+ybuf(i+jj)
      enddo
      enddo
      enddo

c *** if no sky fibres prompt for alternative
      if(nsky.eq.0)then
        print *,' '
        call getchar('Ascii sky spectrum file ',title)
        ioerr=0
        open(unit=13,name=title,status='old',iostat=ioerr)
        j=0
        do while(ioerr.eq.0.and.j.lt.nyout)
        read(13,*,iostat=ioerr) d1,d2
        if(ioerr.eq.0)then
          j=j+1
          skycord(j)=d1
          skyspec(j)=d2
        endif
        enddo
        print *,'No. of points read in =',j
        kref=1
        wcoefs(1,kref)=skycord(1)                   ! bodge the wcoefs
        wcoefs(2,kref)=(skycord(j)-skycord(1))/float(j-1)
      else
        do i=1,nsky
        kk=skynum(i)
        if(i.eq.1)then
          do j=1,nyout
          jj=(j-1)*nxout
          skyspec(j)=xbuf(kk+jj)
          work(j,i)=skyspec(j)
          skycord(j)=wcoefs(1,kk)+wcoefs(2,kk)*(j-1)
          kref=kk
          enddo
        else
          do j=1,nyout
          jj=int((skycord(j)-wcoefs(1,kk))/wcoefs(2,kk))+1
          jj=min(nyout-1,max(1,jj))
          delx=(skycord(j)-wcoefs(1,kk))/wcoefs(2,kk)-jj+1    ! linear interp
          jjj=(jj-1)*nxout
          jp1=jj*nxout
          skyval=xbuf(kk+jjj)*(1.0-delx)+xbuf(kk+jp1)*delx
          skyspec(j)=skyspec(j)+skyval               
          work(j,i)=skyval                          ! save for k-sigma clipping
          enddo
        endif
        enddo

        do j=1,nyout
        skyspec(j)=skyspec(j)/float(nsky)           ! initial guess
        presky(j)=skyspec(j)
        skymask(j)=0.0                              ! initialise skymask
        skytempl(j)=0.0                             ! and template
        enddo

c *** get an estimate of generic noise
        nlin=25
        nmed=201
        niter=5
        call skyfilt(skyspec,nyout,skycont,skyline,maxpix,sigma,
     1               skymask,nlin,nmed,niter)
        do j=1,nyout
        do i=1,nsky
        clip(i)=work(j,i)
        enddo
        call sort1(clip,nsky)
        skymed=0.5*(clip(nsky/2)+clip(nsky/2+1))    ! conservative
        if(skyspec(j)-skymed.gt.3.0*sigma.and.      ! toss out cosmics
     1    (skyspec(j)-skymed)/max(1.0,skymed).gt.0.025)then
          sum=0.0
          do i=1,nsky-1
          sum=sum+clip(i)
          enddo
          skyspec(j)=sum/float(nsky-1)
        endif
        presky(j)=skyspec(j)-presky(j)
        enddo

      endif                                         ! endif for no sky

c *** could clip either side of it too !!!!!!
      ymin=1.0e6
      ymax=-1.0e6
      do j=1,nyout
      ymax=max(ymax,skyspec(j))
      ymin=min(ymin,skyspec(j))
      enddo
      yrange=min(clipmax,ymax)-ymin
      wmin=max(wmin,skycord(1))
      wmax=min(wmax,skycord(nyout))
c *** plots first up the combined sky spectrum
      call pgintro
      call pgscr(1,0.0,0.0,0.0)
      call pgscr(0,1.0,1.0,1.0)
      call pgenv(wmin,wmax,max(-75.0,ymin-0.1*yrange),
     1           min(clipmax,ymax+0.01*yrange),0,0)
      call pglabel('Wavelength','Flux','Sky spectrum')
      call pgline(nyout,skycord,skyspec)
c *** make a copy of the sky spectrum for use in cross
      write(3,'(2f10.3)') (skycord(j),skyspec(j),j=1,nyout)
      call pgline(nyout,skycord,skycont)
c      call pgline(nyout,skycord,presky)
c *** decompose sky into lines and continuum
      nlin=25
      nmed=201
      niter=5
      call skyfilt(skyspec,nyout,skycont,skyline,maxpix,sigma,
     1             skymask,nlin,nmed,niter)

c *** and make a sky mask where significant features are
      do i=1,nyout
      if(abs(skyline(i)).gt.4.0*sigma)then
        skymask(i)=50.0
      else
        skymask(i)=0.0
      endif
      enddo
c *** and grow it to include adjacent pixels
      do i=3,nyout-2
c      if(abs(skyline(i)).gt.4.0*sigma)then
      if(skymask(i).gt.40.0)then
        if(skymask(i-2).eq.0.0)skymask(i-2)=25.0
        if(skymask(i-1).eq.0.0)skymask(i-1)=25.0
        if(skymask(i+1).eq.0.0)skymask(i+1)=25.0
        if(skymask(i+2).eq.0.0)skymask(i+2)=25.0
      endif
      if(skymask(i).gt.0.0)skytempl(i)=skyline(i)     ! template for cross-corr
      enddo
      call pgsci(2)
      call pgline(nyout,skycord,skyline)
      call pgsci(3)
      call pgline(nyout,skycord,skycont)
      call pgsci(4)
      call pgline(nyout,skycord,skymask)
      call pgsci(1)
      call pgextro

c *** flag objects for processing
      do i=1,nxout
      pronum(i)=0
      enddo
      do i=1,nobj
      iarg=objnum(i)
      pronum(iarg)=1
      enddo
      do i=1,nsky
      iarg=skynum(i)
      pronum(iarg)=-1
      enddo
      do i=1,nsky
      iarg=calnum(i)
      pronum(iarg)=-2
      enddo

      wmins=wmin
      wmaxs=wmax
c *** process and plot all the data
      call pgintrom(1,2)
      call pgscr(1,0.0,0.0,0.0)
      call pgscr(0,1.0,1.0,1.0)
      ii=0
      do i=1,nxout
c      if(pronum(i).eq.1)then
        ymin=1.0e6
        ymax=-1.0e6
        wmin=max(wmins,wcoefs(1,i))
        wmax=min(wmaxs,wcoefs(1,i)+wcoefs(2,i)*(nyout-1))
        do j=1,nyout
        jj=(j-1)*nxout
        xcord(j)=wcoefs(1,i)+wcoefs(2,i)*(j-1)         
        if(xcord(j).lt.wmin)ilow=j
        if(xcord(j).lt.wmax)ihih=j
        ycord(j)=xbuf(i+jj)
        ymax=max(ymax,ycord(j))
        ymin=min(ymin,ycord(j))
        enddo
        yrange=min(clipmax,ymax)-ymin                    ! to clip cosmics
        ymin=ymin-0.05*yrange
        ymax=min(clipmax,ymax+0.05*yrange)

c *** first find continuum
        nlin=25
        nmed=201
        niter=5
        call skyfilt(ycord,nyout,objcont,objline,maxpix,sigma,skymask,
     1               nlin,nmed,niter)
c *** check for wavelength offset
        print *,'Object:',i
        call cross(skycord,xcord,skytempl,objline,wmin,wmax,
     1                        wcoefs(2,kref),nyout,nyout,off)
c *** and fix if sub-pixel
        do j=1,nyout
          objcorr(j)=objline(j)
        enddo
        if(abs(off).lt.1.0)then
          do j=2,nyout-1
          if(off.lt.0.0)then
            objcorr(j)=(1.0+off)*objline(j)-off*objline(j-1)
          else
            objcorr(j)=(1.0-off)*objline(j)+off*objline(j+1)
          endif
          enddo
        endif
c *** then do sky subtraction
        call skysub(xcord,objcorr,zcord,skyline,nyout,kref,
     1              wcoefs(1,kref),wmin,wmax,scalemin,skymask,maxpix)
c *** Hanning filter to get rid of pixellation crap
c        call hanning(zcord,nyout)

c        write(3,'(2i8,f10.3)') i,pronum(i),scalemin
c *** some plots
        do j=1,nyout
        zcord(j)=zcord(j)+objcont(j)-scalemin*skycont(j) ! put continua back on
        jj=(j-1)*nxout
        xbuf(i+jj)=zcord(j)                              ! and save for later
        ycord(j)=ycord(j)+1000.0
        enddo
        call pgenv(wmin,wmax,ymin,ymax+1000.0,0,0)
        call pglabel('Wavelength','Flux',specnum(i))
        call pgline(nyout,xcord,zcord)
        call pgline(nyout,xcord,ycord)
c        call pgsci(2)
c        call pgline(nyout,xcord,skyline)
        call pgsci(4)
        call pgline(nyout,skycord,skymask)
        call pgsci(1)

c      endif
      enddo                                              ! end of object loop
      call pgextro

c *** add keyword to header
      card='HISTORY  sky fibre subtraction done'
      call ftprec(outunit,card,status)  
c *** and write out data parts
      print *,' '
      print *,'Writing frames ......'
      print *,' '
      call ftppre(outunit,1,1,npix,xbuf,status)
      if(status.ne.0.and.status.ne.412)then
        call ftgerr(status,errmsg)
        print *,'FITSIO Error status =',status,': ',errmsg
      endif

c *** copy any remaining extensions
      if(nhdus.gt.1)then
        do i=2,nhdus       
        print *,'Copying HDU ......',i
        call ftmahd(inunitref,i,hdutype,status)
        call ftcrhd(outunit,status)
        call ftcopy(inunitref,outunit,0,status)
        if(status.ne.0)then
          call ftgerr(status,errmsg)
          print *,'FITSIO Error status =',status,': ',errmsg
        endif
        enddo
      endif

c *** finish up
      status=0
      call ftclos(outunit,status)
      call ftclos(inunitref,status)
c *** report any errors
      if(status.ne.0.and.status.ne.412)then
        call ftgerr(status,errmsg)
        print *,'FITSIO Error status =',status,': ',errmsg
      endif

      end

c *** ------------------------------------------------------------------------

      subroutine rdbuf_fitsio(iunit,map,wcoefs,specnum,nxout,nyout)
      real*4 map(*),wcoefs(2,*),nullval
      integer status,iunit,naxes(2)
      character*80 errmsg
      character*72 comment
      character*3 dum,specnum(*)
      logical anynull
c *** FITSIO_RDBUF reads in entire fits files using fitsio subroutines
c *** map is the 2-D array of r*4 data stored as a 1-D array 
c *** nxout is the size of the axis most slowly varying in this array
c *** nyout is the size of the axis most rapidly varying.
      lt=6
      lr=5
      status=0
      nullval=0.0
      anynull=.false.
c *** read assorted keywords to find array dimensions etc..
      call ftgknj(iunit,'NAXIS',1,2,naxes,nfound,status)
      if(nfound.ne.2)stop 'Not enough NAXIS keywords'
      nxout=naxes(1)
      nyout=naxes(2)
      npix=nyout*nxout
      if(npix.gt.4096*4096)stop ' Too many pixels for buffer'
      write(lt,1500) nxout,nyout 
 1500 format(/'FITS header information:',2i8)
c *** get the wavelength info
      call ftgkye(iunit,'CRVAL2',crval2,comment,status)
      call ftgkye(iunit,'CDELT2',cdelt2,comment,status)
c *** coeffs
      do k=1,nxout
      if(k.lt.10)then
        write(dum,'(i1)') k
      endif
      if(k.gt.9.and.k.lt.100)then
        write(dum,'(i2)') k
      endif
      if(k.gt.99)then
        write(dum,'(i3)') k
      endif
      wcoefs(1,k)=crval2*10.0                          ! convert to Angstroms
      wcoefs(2,k)=cdelt2*10.0
      write(specnum(k),'(i3)') k
      enddo
c *** ok now read all the data
      print *,' '
      call ftgpve(iunit,1,1,npix,nullval,map,anynull,status)
c *** report any errors
      if(status.ne.0)then
        call ftgerr(status,errmsg)
        print *,'FITSIO Error status =',status,': ',errmsg
      endif
      return
      end

c *** ------------------------------------------------------------------------

      subroutine skyfilt(skyspec,nyout,skycont,skyline,maxpix,sigma,
     1                   skymask,nlin,nmed,niter)
      real*4 skyspec(*),skycont(*),skyline(*),skymask(*)
      real*4 work(maxpix)
c *** estimate continuum for sky and set up sky line mask by iteratively 
c *** filtering 1D spectrum
      do iter=1,niter
      irej=0
      if(iter.eq.1)then
        do i=1,nyout
        skycont(i)=skyspec(i)
        enddo
      else
c *** k-sigma clip and try again
c        nmed=101             ! since clipped out ~1/3 pts should be << but !!
        do i=1,nyout
        if(abs(skyline(i)).gt.3.0*sigma.or.skymask(i).gt.40.0)then
          skycont(i)=-9999.0
          irej=irej+1
        else
          skycont(i)=skyspec(i)
        endif
        enddo
      endif
      call filt1d(skycont,nyout,nmed,nlin)
      do i=1,nyout
      skyline(i)=skyspec(i)-skycont(i)
      work(i)=abs(skyspec(i)-skycont(i))
      enddo
      call sort1(work,nyout)
      sigma=1.48*work(nyout/2)
c      write(6,'(2i6,f10.3,i6)') iter,irej,sigma,niter
      enddo

c *** ok finish off
      do i=1,nyout
      skyline(i)=skyspec(i)-skycont(i)
      enddo
      return
      end

c *** ------------------------------------------------------------------------

      subroutine skysub(xcord,ycord,zcord,skyspec,nyout,kref,wcoefs,
     1                  wmin,wmax,scalemin,skymask,maxpix)
      real*4 wcoefs(2),xcord(*),ycord(*),zcord(*),skyspec(*),skymask(*)
      real*4 sum(3,3),skysmth(maxpix),objsmth(maxpix)
c *** do simple sky subtraction with iterative derivation of scale factor
      scalecurr=1.0
      gap=1.0
c *** and also try slightly blurred versions to match MTF differences
      do i=1,nyout
      objsmth(i)=ycord(i)
      skysmth(i)=skyspec(i)
      enddo
      call hamming(objsmth,nyout)
      call hamming(skysmth,nyout)

      do iter=1,3
      gap=gap/2.0
      do j=1,3
      do i=1,3
      sum(i,j)=0.0
      enddo
      enddo
      do j=1,nyout
      if(xcord(j).ge.wmin.and.xcord(j).le.wmax
     1                                .and.skymask(j).gt.0.0)then
        jj=int((xcord(j)-wcoefs(1))/wcoefs(2))+1
        jj=min(nyout-1,max(1,jj))
        delx=(xcord(j)-wcoefs(1))/wcoefs(2)-jj+1              ! linear interp
        skyval=skyspec(jj)*(1.0-delx)+skyspec(jj+1)*delx
        sum(1,1)=sum(1,1)+abs(ycord(j)-(scalecurr-gap)*skyval)
        sum(2,1)=sum(2,1)+abs(ycord(j)-scalecurr*skyval)
        sum(3,1)=sum(3,1)+abs(ycord(j)-(scalecurr+gap)*skyval)
        sum(1,2)=sum(1,2)+abs(objsmth(j)-(scalecurr-gap)*skyval)
        sum(2,2)=sum(2,2)+abs(objsmth(j)-scalecurr*skyval)
        sum(3,2)=sum(3,2)+abs(objsmth(j)-(scalecurr+gap)*skyval)
        skyval=skysmth(jj)*(1.0-delx)+skysmth(jj+1)*delx
        sum(1,3)=sum(1,3)+abs(ycord(j)-(scalecurr-gap)*skyval)
        sum(2,3)=sum(2,3)+abs(ycord(j)-scalecurr*skyval)
        sum(3,3)=sum(3,3)+abs(ycord(j)-(scalecurr+gap)*skyval)
      endif
      enddo
c *** interpolate by local parabolic fit to get best solution
      sumfound=1.0e10
      do j=1,3
      a=sum(2,j)
      b=0.5*(sum(3,j)-sum(1,j))
      c=0.5*(sum(3,j)+sum(1,j)-2.0*sum(2,j))
      offset=-0.5*b/c
      do i=1,3
      if(sum(i,j).lt.sumfound)then
        sumfound=sum(i,j)
        scalefound=scalecurr+(i-2)*gap
        method=j                                ! flags which smearing to use
        aa=a
        bb=b
        cc=c
      endif
      enddo
      offset=-0.5*b/c
c *** clip to stay in range of scale factors tried
      offset=max(min(1.0,offset),-1.0)
      summin=a+b*offset+c*offset**2
      scalemin=scalecurr+offset*gap
c      write(6,'(4f10.1,f10.3)') (sum(i,j),i=1,3),summin,scalemin
      enddo
      print *,'Method =',method,'    Gap = ',gap
      offset=-0.5*bb/cc
      offset=max(min(1.0,offset),-1.0)
      summin=aa+bb*offset+cc*offset**2
      scalemin=scalecurr+offset*gap
c *** make sure not going for a maximum by mistake
      if(summin.gt.sumfound)then
        spreadmin=sumfound
        scalemin=scalefound
      endif
      write(6,'(4f10.1,f10.3,i6)') (sum(i,method),i=1,3),summin,
     1                              scalemin,iter
      scalecurr=scalemin
      enddo
      print *,' '

c *** apply best scale factor
      do j=1,nyout
      jj=int((xcord(j)-wcoefs(1))/wcoefs(2))+1
      jj=min(nyout-1,max(1,jj))
      delx=(xcord(j)-wcoefs(1))/wcoefs(2)-jj+1              ! linear interp
      if(method.eq.1)then
        skyval=skyspec(jj)*(1.0-delx)+skyspec(jj+1)*delx
        zcord(j)=ycord(j)-scalemin*skyval
      endif
      if(method.eq.2)then
        skyval=skyspec(jj)*(1.0-delx)+skyspec(jj+1)*delx
        zcord(j)=objsmth(j)-scalemin*skyval
      endif
      if(method.eq.3)then
        skyval=skysmth(jj)*(1.0-delx)+skysmth(jj+1)*delx
        zcord(j)=ycord(j)-scalemin*skyval
      endif

      enddo
      return
      end

c *** ------------------------------------------------------------------------

      subroutine sort1(ia,n)
      real*4 ia(n)
      real*4 it
      int=2
   10 int=2*int
      if(int.lt.n)goto 10
      int=min0(n,(3*int)/4-1)
   20 int=int/2
      ifin=n-int
      do 70 ii=1,ifin
      i=ii
      j=i+int
      if(ia(i).le.ia(j))goto 70
      it=ia(j)
   40 ia(j)=ia(i)
      j=i
      i=i-int
      if(i.le.0)goto 60
      if(ia(i).gt.it)goto 40
   60 ia(j)=it
   70 continue
      if(int.gt.1)goto 20
      return
      end

c *** ------------------------------------------------------------------------

      subroutine clobber(iunit,name)
c *** CLOBBER checks to see if fort.iunit exists and either deletes or
c ***         renames it if present
      integer ier,stat,statb(13),status,system
      character*(*) name
      character*80 filename
      character*3 del /'rm '/
      filename=name
      i=1
      do while (filename(i:i).ne.' '.and.i.lt.80)
      i=i+1
      enddo
      open(unit=iunit,file=filename)
      ier=stat(filename,statb)
      if(statb(8).eq.0.or.ier.ne.0)then
        status=system(del//filename)
      else
        print *,'File: '//filename(1:i)//' already exists - clobbering'
        print *,' '
        status=system(del//filename)
      endif
      close(unit=iunit)
      return
      end

c *** ------------------------------------------------------------------------

      subroutine filt1d(xlb,mpt,nmed,nlin)
      real*4 wbuf(16384),xlb(mpt)
      irc=0
      do 450 i=1,mpt
      if(xlb(i).lt.-999.0)goto 450
      irc=irc+1
      wbuf(irc)=xlb(i)
  450 continue
      mrc=irc
      if(mrc.eq.0)return
      call median(wbuf,mrc,nmed)
      irc=0
      do 475 i=1,mpt
      if(xlb(i).lt.-999.0)goto 475
      irc=irc+1
      xlb(i)=wbuf(irc)
  475 continue
      call padext(xlb,mpt)      ! linear interpolation to pad out missing bits
      call linear(xlb,mpt,nlin)
      return
      end

c *** ------------------------------------------------------------------------

      subroutine linear(xbuf,npt,nfilt)
      real*4 xbuf(npt),ybuf(200512)
c *** LINEAR  performs boxcar filtering on array xbuf
      if(nfilt.gt.511)stop ' too large a filter'
      if(npt.gt.200512)stop ' too many points in data array'
      if((nfilt/2)*2.eq.nfilt)nfilt=nfilt+1
      if(npt.le.nfilt)return
c *** set first and last edges equal 
      il=nfilt/2
      ilow=max(3,nfilt/4)
      ilow=(ilow/2)*2+1
      sum=0.0
      do 100 i=1,ilow
  100 sum=sum+xbuf(i)
      xmns=sum/float(ilow)
      sum=0.0
      do 150 i=1,ilow
  150 sum=sum+xbuf(npt+1-i)
      xmnf=sum/float(ilow)
c *** reflect edges before filtering
      do 200 i=1,il
      ybuf(i)=2.0*xmns-xbuf(il+ilow+1-i)
  200 ybuf(npt+i+il)=2.0*xmnf-xbuf(npt-i-ilow+1)
      do 250 i=1,npt
  250 ybuf(i+il)=xbuf(i)
c *** do linear filtering on rest
      fnfilt=float(nfilt)
      sum=0.0
      do 300 i=1,nfilt
  300 sum=sum+ybuf(i)
      xbuf(1)=sum/fnfilt
      do 400 i=2,npt
      sum=sum-ybuf(i-1)
      sum=sum+ybuf(i+nfilt-1)
  400 xbuf(i)=sum/fnfilt
      return
      end

c *** ------------------------------------------------------------------------

      subroutine median(xbuf,npt,nfilt)
      real*4 xbuf(npt),ybuf(200512),array(512)
      integer*2 point(512)
c *** MEDIAN  performs median filtering on array xbuf
      if(nfilt.gt.511)stop ' too large a filter'
      if(npt.gt.200512)stop ' too many points in data array'
      if((nfilt/2)*2.eq.nfilt)nfilt=nfilt+1
      if(npt.le.nfilt)return
      nfo2p1=nfilt/2+1
c *** set first and last edges equal 
      il=nfilt/2
      ilow=max(3,nfilt/4)
      ilow=(ilow/2)*2+1
      do 100 i=1,ilow
  100 array(i)=xbuf(i)
      call sortm(array,point,ilow)
      xmns=array(ilow/2+1)
      do 150 i=1,ilow
  150 array(i)=xbuf(npt+1-i)
      call sortm(array,point,ilow)
      xmnf=array(ilow/2+1)
c *** reflect edges before filtering
      do 200 i=1,il
      ybuf(i)=2.0*xmns-xbuf(il+ilow+1-i)
  200 ybuf(npt+i+il)=2.0*xmnf-xbuf(npt-i-ilow+1)
      do 250 i=1,npt
  250 ybuf(i+il)=xbuf(i)
c *** do median filtering on rest
      do 260 i=1,nfilt
      array(i)=ybuf(i)
  260 point(i)=i
      call sortm(array,point,nfilt)
      xbuf(1)=array(nfo2p1)
      jl=nfilt+1
      jh=nfilt+npt-1
      do 400 j=jl,jh
      do 300 i=1,nfilt
      if(point(i).eq.1)goto 290
      point(i)=point(i)-1
      goto 300
  290 point(i)=nfilt
      array(i)=ybuf(j)
      l=i
  300 continue
      call qsort(array,point,l,nfilt)
  400 xbuf(j-jl+2)=array(nfo2p1)
      return
      end
      subroutine sortm(ia,ib,n)
      real*4 ia(n)
      integer*2 ib(n),iu
      real*4 it
      int=2
   10 int=2*int 
      if(int.lt.n)goto 10
      int=min0(n,(3*int)/4-1)
   20 int=int/2
      ifin=n-int
      do 70 ii=1,ifin
      i=ii
      j=i+int
      if(ia(i).le.ia(j))goto 70
      it=ia(j)
      iu=ib(j)
   40 ia(j)=ia(i)
      ib(j)=ib(i)
      j=i
      i=i-int
      if(i.le.0)goto 60
      if(ia(i).gt.it)goto 40
   60 ia(j)=it
      ib(j)=iu
   70 continue
      if(int.gt.1)goto 20
      return
      end
      subroutine qsort(x,point,l,nfilt)
      real*4 x(nfilt)
      integer*2 it,point(nfilt)
      test=x(l)
      do 100 i=1,nfilt
      if(i.eq.l)goto 100
      if(test.gt.x(i))goto 100
      j=i
      goto 200
  100 continue
      j=nfilt+1
  200 if(j-1.eq.l)return
      if(j-l)300,500,400
  300 temp=x(l)
      it=point(l)
      npt=l-j
      do 350 i=1,npt
      ii=l-i
      x(ii+1)=x(ii)
  350 point(ii+1)=point(ii)
      x(j)=temp
      point(j)=it
      return
  400 temp=x(l)
      it=point(l)
      j=j-1
      npt=j-l
      if(npt.eq.0)goto 475
      do 450 i=1,npt
      ii=l+i
      x(ii-1)=x(ii)
  450 point(ii-1)=point(ii)
  475 x(j)=temp
      point(j)=it
  500 return
      end

c *** ------------------------------------------------------------------------

      subroutine padext(x,n)
      real*4 x(n)
c *** PADEXT  pads out array with missing points
c ***         and linearly extrapolates the ends
c ***         points .le. 0 are treated as missing
      i=0
   10 i=i+1
      if(x(i).gt.-999.0)goto 40
      goto 10
   40 il=i+1
      ilow=il-1
      do 300 i=il,n
      if(x(i).gt.-999.0)goto 250
      ic=1
   50 if(x(i+ic).gt.-999.0)goto 100
      ic=ic+1
      if(ic+i.ge.n)goto 500
      goto 50
  100 xlow=x(i-1)
      xhih=x(i+ic)
      npt=ic
      do 200 j=1,npt
      t2=float(j)/float(ic+1)
      t1=1.0-t2
  200 x(i+j-1)=t1*xlow+t2*xhih
      goto 300
  250 ihih=i
  300 continue
c *** linear extrapolation of ends
  500 if(ilow.eq.1)goto 600
      slope=x(ilow+1)-x(ilow)
      do 550 i=1,ilow-1
  550 x(i)=x(ilow)-slope*(ilow-i)
  600 if(ihih.eq.n)return
      slope=x(ihih)-x(ihih-1)
      do 650 i=ihih+1,n
  650 x(i)=x(ihih)+slope*(i-ihih)
      return
      end

c *** ------------------------------------------------------------------------

      subroutine hamming(xbuf,npt)
      real*4 xbuf(npt),ybuf(16896)
c *** HAMMING  performs linear filtering on array xbuf
      if(npt.gt.16384)stop ' too many points in data array'
      nfilt=3
      if(npt.le.nfilt)return
c *** set first and last edges equal 
      il=nfilt/2
      ilow=max(3,nfilt/4)
      ilow=(ilow/2)*2+1
      sum=0.0
      do 100 i=1,ilow
  100 sum=sum+xbuf(i)
      xmns=sum/float(ilow)
      sum=0.0
      do 150 i=1,ilow
  150 sum=sum+xbuf(npt+1-i)
      xmnf=sum/float(ilow)
c *** reflect edges before filtering
      do 200 i=1,il
      ybuf(i)=2.0*xmns-xbuf(il+ilow+1-i)
  200 ybuf(npt+i+il)=2.0*xmnf-xbuf(npt-i-ilow+1)
      do 250 i=1,npt
  250 ybuf(i+il)=xbuf(i)
c *** do linear filtering on rest
      do i=1,npt
      xbuf(i)=(ybuf(i)+4.0*ybuf(i+1)+ybuf(i+2))/6.0  ! 1-4-1 weighting
      enddo
      return
      end

c *** ----------------------------------------------------------------------

      subroutine pgintro
      ido=0
      call posei4('Do you want to alter default line width ? [N]',ido)
      if(ido.eq.1)then
        call geti4('Line width (1?)',ilw)
      endif
      call pgbegin(7,'?',1,1)
      if(ido.eq.1)call pgslw(ilw)
      return
      end

      subroutine pgextro
      call pgend
      return
      end

      subroutine pgintrom(nx,ny)
      ido=0
      call posei4('Do you want to alter default line width ? [N]',ido)
      if(ido.eq.1)then
        call geti4('Line width (1?)',ilw)
      endif
      call pgbegin(7,'?',nx,ny)
      if(ido.eq.1)call pgslw(ilw)
      return
      end

      subroutine pgenvm(xmin,xmax,ymin,ymax,xw1,xw2,yw1,yw2)
      call pgadvance
      call pgsvp(xw1,xw2,yw1,yw2)
      call pgswin(xmin,xmax,ymin,ymax)
      call pgbox('BCNST',0.0,0,'BCNST',0.0,0)
      return
      end

      subroutine pgenvnl(xmin,xmax,ymin,ymax,xw1,xw2,yw1,yw2)
      call pgadvance
      call pgsvp(xw1,xw2,yw1,yw2)
      call pgswin(xmin,xmax,ymin,ymax)
      call pgbox('BCST',0.0,0,'BCNST',0.0,0)
      return
      end

c *** ------------------------------------------------------------------------

      subroutine fitsio_rtable(iunit,nrmax,ncmax,
     1                         object,fps,rmag,nrows,ier)
      real*4 rmag(nrmax)
      integer status,hdutype,naxes(2),fps(nrmax)
      character*30 object(nrmax)
      character*80 errmsg
      character*72 comment
      character*68 value
      character*16 ttype(ncmax),tunit(ncmax),tform(ncmax)
      logical anynull
c *** FITSIO_RTABLE - reads fits binary tables
      enullval=0.0      
      inullval=0
      status=0
      ier=0
c *** find out what sort of table extension is
      call ftmahd(iunit,4,hdutype,status)
      if(hdutype.eq.1)then
        print *,'Extension is an ASCII table'
        stop ' - not implemented yet'
      else if(hdutype.eq.2)then
        print *,'Extension is a BINARY table'
      endif

c *** read assorted keywords to find table dimensions etc.
      call ftgknj(iunit,'NAXIS',1,2,naxes,nfound,status)
      if(nfound.ne.2)stop 'Not enough NAXIS keywords'
      nrows=naxes(2)
      call ftgkey(iunit,'TFIELDS',value,comment,status)
      read(value,*) ncols
      print *,'No. of columns =',ncols
      print *,'No. of rows    =',nrows

c *** now get rest of table info
      call ftgkns(iunit,'TTYPE',1,ncols,ttype,nfound,status)
      call ftgkns(iunit,'TUNIT',1,ncols,tunit,nfound,status)
      call ftgkns(iunit,'TFORM',1,ncols,tform,nfound,status)
      if(ncols.ne.nfound)stop 'Something screwy with header'

c      write(6,'(a,a,a)') (ttype(i),tunit(i),tform(i),i=1,ncols)

c *** now read data by column
      call ftgcvj(iunit,1,1,1,nrows,inullval,fps,anynull,status)
      call ftgcvs(iunit,2,1,1,nrows,' ',object,anynull,status)
      call ftgcve(iunit,8,1,1,nrows,enullval,rmag,anynull,status)
      write(6,'(a,i8,f8.2)') (object(i),fps(i),rmag(i),i=1,nrows)

c *** finish off
      if(status.ne.0)then
        call ftgerr(status,errmsg)
        print *,'FITSIO Error status =',status,': ',errmsg
        ier=1
      endif
      return
      end


