C-------------------------------------------------------------------------------
C
C   MERGE  -  Listing, unification and merging of two column ASCII spectra
C
C-------------------------------------------------------------------------------
C
C  The program does the following:
c
c  1/ Writes frequency sorted lists of basic properties of all spectra in 
c     the current directory
c
c  2/ Unifies the spectra to a common frequency grid (defined by FSTEP)
c
c  3/ Merges spectra into two files containing all A spectra and all
c     B spectra, respectively
C
C  Notes: 
C  
C  1. MERGE works only on the current directory and spectral files are
c     to have extension .SPE
C  2. no spaces allowed in file names  
C  3. data points have to be equidistant in frequency 
C  4. missing parts are filled with zero, overlapping parts are connected at
C     the middle of the overlap region
c  5. This is step 4 in producing a broadband spectrum from scanned multipage
c     paper copies   TRACE -> SPLICE -> FZERO -> MERGE
c  
c
C   Ver 24.XI.2008                                ----- Zbigniew Kisiel -----
C                          __________________________________________________
C                         | Institute of Physics, Polish Academy of Sciences |
C                         | Al.Lotnikow 32/46, Warszawa, POLAND              |
C                         |                             kisiel@ifpan.edu.pl  |
C                         |     http://info.ifpan.edu.pl/~kisiel/prospe.htm  |
C_________________________/--------------------------------------------------
C
C
      USE DFLIB
c      
      implicit real*8 (a-h,o-z)
      implicit integer*4 (i-n)
      PARAMETER (MAXPTS=30000,maxspe=150,fstep=0.5d0,maxall=1000000)
C     
      INTEGER*4 npoint(maxspe),iwk(maxspe),npconv(maxspe)
      REAL*8 FSTART(maxspe),FEND(maxspe),FINCR(maxspe)
      real*8 fspe(maxspe,maxpts),yspe(maxspe,maxpts)
      real*8 fconv(maxpts),yconv(maxpts)
      real*8 fa(maxall),ya(maxall),fb(maxall),yb(maxall)     
      CHARACTER*100 line,filin,name(maxspe),filout
c
      COMMON /FRBLK/fstart
      COMMON /LNUMS/IWK
c
C
C...WRITE THE HEADER INFORMATION
C
      WRITE(*,3344)
3344  FORMAT(1X//'  ',76('_')/' |',T79,'|'/
     * ' |     M E R G E  -  Listing, unification and merging ', 
     * T79,'|'/
     * ' |                   of two column ASCII spectra ',T79,'|'/
     * ' |',76('_'),'|'/'  version 24.XI.2008',T64,'Zbigniew KISIEL'/)
c
C...list the directory
c
      fsys=systemqq('dir *.spe/Od>spec.lst')
C
C
C...Open and go through the directory listing 
C
      open(6,file='spec.lst',status='old')        
c
c - - - - - - -
      nspec=0
3     read(6,'(a)',end=4)line
      if(line(1:1).ne.'2')goto 3
      i=len_trim(line)-1
5     if(line(i:i).ne.' ')then
        i=i-1
        goto 5
      endif
      filin=line(i+1:len_trim(line))
      if(filin(1:2).eq.'u_')goto 3
c
c...go through the contents of the spectrum
c    
      npts=0
8     open(4,file=filin(1:len_trim(filin)),status='old')
7     read(4,'(a)',err=6,end=6)line
      if(line(1:1).eq.'!')goto 7
      read(line,*,end=6,err=6)ff,yy   
      npts=npts+1    
c
      if(npts.eq.1)then
        nspec=nspec+1
        fstart(nspec)=ff
      endif
c     write(*,*)nspec,npts
      fspe(nspec,npts)=ff
      yspe(nspec,npts)=yy
      goto 7
c
6     if(npts.gt.0)then
        npoint(nspec)=npts
        fend(nspec)=ff
        fincr(nspec)=(fend(nspec)-fstart(nspec))/real(npts-1)
        name(nspec)=filin(1:len_trim(filin))
        iwk(nspec)=nspec
        write(*,514)name(nspec)(1:15),fstart(nspec),fend(nspec),
     *    npoint(nspec),fincr(nspec)
514     format(1x,a,5x,f10.3,' -- ',f10.3,i8,' points  @',
     *    f9.5)
      endif
 
      close(4)
c
c...convert spectrum to the predefined grid, such that first frequency is 
c   at an exact multiple of grid spacing
c
      nsp=1
      n=0
      freq=int(fspe(nspec,1)/fstep)*fstep
c
659   freq=freq+fstep
      if(freq.ge.fspe(nspec,npts))goto 660
      n=n+1
662   if(fspe(nspec,nsp+1).le.freq)then
        nsp=nsp+1
        goto 662
      endif
      yconv(n)=yspe(nspec,nsp)+
     *   (yspe(nspec,nsp+1)-yspe(nspec,nsp))*
     *   (freq-fspe(nspec,nsp))/(fspe(nspec,nsp+1)-fspe(nspec,nsp))
      fconv(n)=freq
      goto 659    
c
660   nconv=n
      npconv(nspec)=nconv
      do 661 n=1,nconv
        fspe(nspec,n)=fconv(n)
        yspe(nspec,n)=yconv(n)
661   continue
c
c...write out the converted spectrum
c
      filout='u_'//filin(1:len_trim(filin))
      open (2,file=filout(1:len_trim(filout)),status='unknown')     
      do 663 n=1,nconv
        write(2,'(f15.5,f12.3)')fconv(n),yconv(n)
663   continue
      close(2)
c
c      
      goto 3
c
c - - - - - - -
c
4     write(*,'(1x/i5,'' spectra found in this directory''//)')nspec
      close(6)
      fsys=systemqq('del spec.lst')      
c
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c
c...Sort the spectra in order of increasing frequency
c
      i=1
      if(nspec.gt.1)call sorth(i,nspec)                                 <-----
c
c
c...Write the LIST file of all spectra
c
      open(4,file='LIST',STATUS='UNKNOWN')
      do 510 j=1,nspec
        i=iwk(j)
        write(4,504)name(i)(1:15),fstart(j),fend(i),npoint(i),
     * fincr(i)
        if((j/10)*10.eq.j)write(4,'(1x)')
510   continue
      write(4,'(1x/i5,'' spectra listed''/)')nspec
c
      close(4)      
504   format(1x,a,5x,f9.2,' -- ',f9.2,' MHz',i9,' pts',
     *  f12.6,' MHz/pt')
c
c
c...Write the LISTa file of all A spectra
c
      open(4,file='LISTa',STATUS='UNKNOWN')
      nn=0
      do 511 j=1,nspec
        i=iwk(j)
        ii=len_trim(name(i))-4
        if(name(i)(ii:ii).ne.'a'.and.name(i)(ii:ii).ne.'A')goto 511
        write(4,504)name(i)(1:15),fstart(j),fend(i),npoint(i),
     * fincr(i)
        nn=nn+1
        if((nn/10)*10.eq.nn)write(4,'(1x)')
511   continue
      write(4,'(1x/i5,'' spectra listed''/)')nn
c
      close(4)      
c
c
c...Write the LISTb file of all B spectra
c
      open(4,file='LISTb',STATUS='UNKNOWN')
      nn=0
      do 512 j=1,nspec
        i=iwk(j)
        ii=len_trim(name(i))-4
        if(name(i)(ii:ii).ne.'b'.and.name(i)(ii:ii).ne.'B')goto 512
        write(4,504)name(i)(1:15),fstart(j),fend(i),npoint(i),
     * fincr(i)
        nn=nn+1
        if((nn/10)*10.eq.nn)write(4,'(1x)')
512   continue
      write(4,'(1x/i5,'' spectra listed''/)')nn
c
      close(4)      
c
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c
c...Merge all A spectra into a single spectral file and
c         all B spectra into a single spectral file
c
c
c   Pointers to spectra in order of frequency are contained in IWK
c   Spectra are contained in FSPE(nspe,npt)
c                            YSPE(nspe,npt)
c   Name of spectra are in   NAME(nspe)
c   Numbers of points are in NPCONV(nspe)
c
      write(*,'('' Merging .....'')')
      npta=0
      nptb=0
      do 700 j=1,nspec
        i=iwk(j)
        ii=len_trim(name(i))-4
c
c...spectrum A
c
        if(name(i)(ii:ii).eq.'a'.or.name(i)(ii:ii).eq.'A')then
c
c__first spectrum
          if(npta.eq.0)then
            nn=npconv(i)
            do 701 n=1,nn
              fa(n)=fspe(i,n)
              ya(n)=yspe(i,n)
701         continue                        
            npta=nn
            goto 700
          endif
c
c__insert points
c
          if(fspe(i,1).gt.fa(npta))then
703         ffa=fa(npta)+fstep
            if(ffa.lt.fspe(i,1))then
              npta=npta+1
              fa(npta)=ffa
              ya(npta)=0.d0
              goto 703
            endif
            do 704 n=1,npconv(i)
              npta=npta+1
              fa(npta)=fspe(i,n)
              ya(npta)=yspe(i,n)
704         continue
            goto 700
          endif
c
c__deal with overlap
          if(fspe(i,1).lt.fa(npta))then
            nover=0.5d0*(fa(npta)-fspe(i,1))/fstep
            npta=npta-nover-5
707         if(fa(npta).lt.fspe(i,nover))then
              npta=npta+1
              goto 707
            endif
            npta=npta-1
            do 708 n=nover,npconv(i)
              npta=npta+1
              fa(npta)=fspe(i,n)
              ya(npta)=yspe(i,n)
708         continue             
            goto 700
          endif
c
        endif
c
c...spectrum B
c
        if(name(i)(ii:ii).eq.'b'.or.name(i)(ii:ii).eq.'B')then
          if(nptb.eq.0)then
            nn=npconv(i)
            do 702 n=1,nn
              fb(n)=fspe(i,n)
              yb(n)=yspe(i,n)
702         continue                        
            nptb=nn
            goto 700
          endif
c
          if(fspe(i,1).gt.fb(nptb))then
705         ffb=fb(nptb)+fstep
            if(ffb.lt.fspe(i,1))then
              nptb=nptb+1
              fb(nptb)=ffb
              yb(nptb)=0.d0
              goto 705
            endif
            do 706 n=1,npconv(i)
              nptb=nptb+1
              fb(nptb)=fspe(i,n)
              yb(nptb)=yspe(i,n)
706         continue
            goto 700
          endif
c
          if(fspe(i,1).lt.fb(nptb))then
            nover=0.5d0*(fb(nptb)-fspe(i,1))/fstep
            nptb=nptb-nover-5
709         if(fb(nptb).lt.fspe(i,nover))then
              nptb=nptb+1
              goto 709
            endif
            nptb=nptb-1
            do 710 n=nover,npconv(i)
              nptb=nptb+1
              fb(nptb)=fspe(i,n)
              yb(nptb)=yspe(i,n)
710         continue             
            goto 700
          endif
c
        endif
c
700   continue
c
c...write out the two merged spectra
c
      write(*,'(1x/'' Merged spectrum written to: '',a)')'u_A.spe'
      open(2,file='u_A.spe',status='unknown')
      do 800 n=1,npta
        write(2,'(f15.5,f12.3)')fa(n),ya(n)
800   continue       
      close(2) 
c
      write(*,'('' Merged spectrum written to: '',a)')'u_B.spe'
      open(2,file='u_B.spe',status='unknown')
      do 801 n=1,nptb
        write(2,'(f15.5,f12.3)')fb(n),yb(n)
801   continue       
      close(2)
c
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c
530   stop
      end
C
C------------------------------------------------------------------------
c
      SUBROUTINE SORTH(NSTART,N)
c
c   This routine is based on the SORT2 'heapsort' routine from Numerical
c   Recipes and sorts the quantities in vector WK from WK(NSTART) to WK(N)
C   in ascending order of magnitude and also accordingly rearranges vector
C   IPT of pointers to original positions of sorted quantities.
c
      PARAMETER (MAXspe=150)
c
      COMMON /FRBLK/WK
      COMMON /LNUMS/IPT
      INTEGER IPT(maxspe),IIPT,L,N,NSTART,I,J,IR
      REAL*8 WK(maxspe),WWK
C
      L=N/2+1
      IR=N
10    CONTINUE
        IF(L.GT.NSTART)THEN
          L=L-1
          WWK=WK(L)
          IIPT=IPT(L)
        ELSE
          WWK=WK(IR)
          IIPT=IPT(IR)
          WK(IR)=WK(1)
          IPT(IR)=IPT(1)
          IR=IR-1
          IF(IR.EQ.NSTART)THEN
            WK(1)=WWK
            IPT(1)=IIPT
            RETURN
          ENDIF
        ENDIF
        I=L
        J=L+L
20      IF(J.LE.IR)THEN
          IF(J.LT.IR)THEN
            IF(WK(J).LT.WK(J+1))J=J+1
          ENDIF
          IF(WWK.LT.WK(J))THEN
            WK(I)=WK(J)
            IPT(I)=IPT(J)
            I=J
            J=J+J
          ELSE
            J=IR+1
          ENDIF
        GO TO 20
        ENDIF
        WK(I)=WWK
        IPT(I)=IIPT
      GO TO 10
c
      RETURN
      END
C
C------------------------------------------------------------------------
C------------------------------------------------------------------------