program ssimil
      implicit none
      integer*4 np,nsi,nt,npn,ic,ix,np1,npall,nsn
      real*8 norm,n1,xs0,xs1,xs2,an,di,s,dev,ds,smin,smax,sm,
     1devf,wmax,wmin
      real*8,allocatable::sx(:),s1(:),sxf(:),st(:),s2(:),sf(:,:)
      logical lstable,linv
      character*80 s1s,s2s
      write(6,6000)
6000  format(/,'  ssimil: similarity index of two spectra',/,/,
     1'      Input: s1.prn .. reference spectrum',/,
     1'             s2.prn .. trial (scaled) spectrum',/,
     1'             SSIMIL.OPT .. options',/,
     1'     Output: SI.PRN .. similarity index vs scale factor',/)

      call readopt(wmin,wmax,smin,smax,nsi,xs0,xs1,xs2,lstable,linv,ic)

c     load detailed scaling table:
      if(lstable)then
       nt=nsn('STABLE.TXT')
       allocate(sf(2,nt))
       call lstab(nt,sf)
      else
       nt=1
       allocate(sf(2,nt))
      endif

      if(iargc().ge.2)then
       call getarg(1,s1s)
       call getarg(2,s2s)
      else
       s1s='S1.prn'
       s2s='S2.prn'
      endif

c     first (reference) spectrum:
      nt=npall(s1s)
      np=np1(s1s,wmin,wmax)
      allocate(sx(np),s1(np),s2(np))
      call ls(sx,s1,s1s,wmin,wmax)
      n1=norm(np,s1,ic)
      if(linv)s1=-s1

c     read second:
      npn=npall(s2s)
      allocate(sxf(npn),st(npn))
      call lsf(npn,sxf,st,s2s)
      call lsr(sx,s2,np,xs0,xs1,xs2,lstable,nt,sf,
     1sxf,st,1.0d0,npn)
      an=norm(np,s2,ic)

      write(6,598)nt,np,n1,npn,np,an,an/n1
 598  format(' s1    :',i5,' points' /,
     1       '        ',i5,' between wmin and wmax',/,
     1       '        ','s1 norm:',e12.4,/,
     1       ' s2    :',i5,' points' /,
     1       '        ',i5,' between wmin and wmax (interpol.)',/,
     1       '        ','s2 norm:',e12.4,' (',e12.4,' rel.)',/)

      call nspe(np,s1,n1)
      call nspe(np,s2,an)
      dev=devf(np,s1,s2,ic)
      write(6,76)an,dev
76    format(' Second norm deviation:',e12.4,f12.6)

c     loop over scale factor
      ds=(smax-smin)/nsi
      s=smin-ds
      sm=smin
      di=1.0d90
      if(ic.eq.1)di=1.0d0/di
      open(9,file='SI.PRN')
      do 1 ix=1,nsi
      s=s+ds
      call lsr(sx,s2,np,xs0,xs1,xs2,lstable,nt,sf,
     1sxf,st,s,npn)
      an=norm(np,s2,ic)
      call nspe(np,s2,an)
      dev=devf(np,s1,s2,ic)
      if(ic.eq.1)then
       if(ix.eq.1)then
        sm=s
        di=dev
       endif
       if(dev.gt.di)then
        sm=s
        di=dev
       endif
      else
       if(dabs(dev).lt.di)then
        sm=s
        di=dabs(dev)
       endif
      endif
1     write(9,789)s,dev
789   format(2f12.4)
      write(6,79)sm,di
79    format(/,' Best sf and error:',2f12.4,/,
     1' Full dependence at SI.PRN')
      close(9)
      end

      subroutine lsr(xs,ys,np,xs0,xs1,xs2,lstable,nt,sf,
     1sx,sy,ss,npn)
c     loads in the full spectrum, adjust number of points to np
c     with the scaled xs
      implicit none
      integer*4 np,i,npn,j,nt
      real*8 xs(*),ys(*),xi,xj,xjp,xs0,xs1,xs2,sf(2,nt),xjf,
     1sx(*),sy(*),ss
      logical lstable
      do 1 i=1,np
      xi=xs(i)
      ys(i)=0.0d0
      do 2 j=1,npn-1
c     scaled values:
      if(lstable)then
       xj =xjf(sx(j  ),nt,sf)
       xjp=xjf(sx(j+1),nt,sf)
      else
       xj =xs0+xs1*sx(j  )+xs2*sx(j  )*2
       xjp=xs0+xs1*sx(j+1)+xs2*sx(j+1)*2
      endif
c     final scaling
      xj=xj*ss
      xjp=xjp*ss
      if(xj.le.xi.and.xjp.ge.xi)then
       ys(i)=sy(j)+(sy(j+1)-sy(j))*(xi-xj)/(xjp-xj)
       goto 1
      endif
 2    continue
 1    continue
      return
      end

      subroutine lsn(ns,sn)
      implicit none
      integer*4 j,ns
      character*80 sn(ns)
      open(9,file='SP.LST')
      do 1 j=1,ns
1     read(9,9000,end=99,err=99)sn(j)
9000  format(a80)
99    close(9)
      write(6,*)'names loaded'
      return
      end


      function devf(np,s1,s2,ic)
      implicit none
      integer*4 np,ix,ic
      real*8 s1(*),s2(*),dev,devf
      dev=0.0d0
      if(ic.eq.0)then
       do 77 ix=1,np
77     dev=dev+dabs(s1(ix)-s2(ix))
      else
       if(ic.eq.1)then
        do 78 ix=1,np
78      dev=dev+s1(ix)*s2(ix)
       else
        if(ic.eq.2)then
         do 79 ix=1,np
79       dev=dev+(s1(ix)-s2(ix))**2
         dev=dsqrt(dev)
        else
         write(6,*)ic
         call report('Unknown error definition')
        endif
       endif
      endif
      devf=dev
      return
      end



      subroutine report(s)
      character*(*) s
      write(6,*)s
      stop
      end


      function nsn(fn)
      implicit none
      integer*4 nsn,ns
      character*(*) fn
      open(9,file=fn,status='old')
      ns=0
101   read(9,*,end=99,err=99)
      ns=ns+1
      goto 101
99    close(9)
      nsn=ns
      return
      end


      subroutine readopt(wmin,wmax,smin,smax,nsi,xs0,xs1,xs2,
     1lstable,linv,ic)
      implicit none
      integer*4 nsi,ic
      real*8 smin,smax,wmin,wmax,xs0,xs1,xs2
      logical lstable,lex,linv
      character*10 s10,s10c
      nsi=100
      smin=0.95d0
      smax=1.05d0
      wmin=200.0d0
      wmax=1900.0d0
      xs0=0.0d0
      xs1=1.0d0
      xs2=0.0d0
      ic=1
      lstable=.false.
      linv=.false.

      inquire(file='SSIMIL.OPT',exist=lex)
      if(lex)then
       open(9,file='SSIMIL.OPT')
1      read(9,90,end=99,err=99)s10
90     format(a10)
       write(6,90)s10
       if(s10(1:4).eq.'wmax')read(9,*)wmax
       if(s10(1:4).eq.'wmin')read(9,*)wmin
       if(s10(1:4).eq.'smax')read(9,*)smax
       if(s10(1:4).eq.'smin')read(9,*)smin
       if(s10(1:4).eq.'stab')read(9,*)lstable
       if(s10(1:4).eq.'lsta')read(9,*)lstable
       if(s10(1:4).eq.'linv')read(9,*)linv
       if(s10(1:3).eq.'xs0')read(9,*)xs0
       if(s10(1:3).eq.'xs1')read(9,*)xs1
       if(s10(1:3).eq.'xs2')read(9,*)xs2
       if(s10(1:3).eq.'nsi')read(9,*)nsi
       if(s10(1:4).eq.'norm')read(9,*)ic
       backspace 9
       read(9,90)s10c
       if(s10.eq.s10c)call report(' Unknown option '//s10)
       goto 1
99     close(9)
      endif
      write(6,6000)wmin,wmax,smin,smax,nsi,lstable,linv,
     1xs0,xs1,xs2,ic
6000  format(' wmin:',f12.1,'   wmax:',f12.1,/,
     1       ' smin:',f12.3,'   smax:',f12.3,',',i5,' points',/,
     1       ' stable:',l12,' invert:',l12,/,
     1       ' xs0 xs1 xs2:',3f12.4,/,
     1       ' norm: ',i2,/)
      if(ic.eq.0)write(6,600)
600   format(/,'                      /              ',/,
     1         '   Normalization: S / | |S| dv       ',/,
     1         '                      /              ',/)
      if(ic.eq.1.or.ic.eq.2)write(6,601)
601   format(/,'                     /  /  2   \1/2 ',/,
     1         '  Normalization: S / |  | S dv |    ',/,
     1         '                     \  /      /    ',/)
      if(ic.eq.0)write(6,603)
603   format(/,'               /                       ',/,
     1         '  Similarity: | |S1 - S2| dv    0 ... 2',/,
     1         '              /                        ',/,/)
      if(ic.eq.1)write(6,604)
604   format(/,'              /                    ',/,
     1         '  Similarity: | S1 x S2 dv  1 .. -1',/,
     1         '              /                    ',/,/)
      if(ic.eq.2)write(6,605)
605   format(/,'              / /                 \1/2       ',/,
     1         '  Similarity: | | (S1 - S2)^2 dv  |    0 .. 2',/,
     1         '              \ /                 /          ',/,/)
      return
      end

      function np1(s,wmin,wmax)
      implicit none
      integer*4 np1,n
      real*8 wmin,wmax,x
      character*(*) s
      open(9,file=s,status='old')
      n=0
101   read(9,*,end=99,err=99)x
      if(x.ge.wmin.and.x.le.wmax)n=n+1
      goto 101
99    close(9)
      np1=n
      return
      end


      subroutine lstab(ns,sf)
      implicit none
      integer*4 j,ns
      real*8 sf(2,ns)
      open(9,file='STABLE.TXT')
      do 1 j=1,ns
1     read(9,*,end=99,err=99)sf(1,j),sf(2,j)
99    close(9)
      write(6,*)ns,' assignments loaded'
      return
      end

      subroutine ls(xs,ys,name,wmin,wmax)
c     loads in the spectrum  within wmin wmax    
      implicit none
      integer*4 i
      real*8 xs(*),ys(*),wmin,wmax,x,y
      character*(*) name
      i=0
      open(9,file=name)
 11   read(9,*,end=99,err=99)x,y
      if(wmin.le.x.and.x.le.wmax)then
       i=i+1       
       xs(i)=x
       ys(i)=y
      endif
      goto 11
  99  close(9)
      return
      end

      function norm(np,s,ic)
      implicit none
      integer*4 i,np,ic
      real*8 norm,s(*)
      norm=0.0d0
      if(ic.eq.0)then
       do 1 i=1,np
   1   norm=norm+dabs(s(i))
      else
       if(ic.eq.1.or.ic.eq.2)then
        do 2 i=1,np
   2    norm=norm+s(i)**2
        norm=dsqrt(norm)
       else
        write(6,*)ic
        call report('unknown norm option')
       endif
      endif
      return
      end

      subroutine neg(np,s)
      implicit none
      integer*4 i,np
      real*8 s(*)
      do 1 i=1,np
   1  s(i)=-s(i)
      return
      end

      subroutine nspe(np,s,norm)
      implicit none
      integer*4 np,i
      real*8 norm,s(*)
      do 1 i=1,np
   1  s(i)=s(i)/norm
      return
      end

      function npall(s)
      implicit none
      integer*4 npall,n
      character*(*) s
      n=0
      open(9,file=s)
101   read(9,*,end=99,err=99)
      n=n+1
      goto 101
99    close(9)
      npall=n-1
      return
      end

      subroutine lsf(n,xs,ys,name)
c     loads in the full spectrum      
      implicit none
      integer*4 i,n
      real*8 xs(*),ys(*)
      character*(*) name
      open(9,file=name)
      do 11 i=1,n
 11   read(9,*)xs(i),ys(i)
      close(9)
      return
      end

      function xjf(xj,nt,sf)
      implicit none
      integer*4 nt,i
      real*8 xjf,xj,sf(2,nt),xe1,xe2,xc1,xc2
      do 1 i=1,nt-1
      xe1=sf(1,i )
      xe2=sf(1,i+1)
      xc1=sf(2,i )
      xc2=sf(2,i+1)
c     xe1,xe2 ... experimental values, xc1,xc2 ... corresponding calc.
      if((xj.ge.xc1.and.xj.le.xc2).or.
     1   (xj.le.xc1.and.xj.ge.xc2))then
       xjf=xe1+(xj-xc1)*(xe2-xe1)/(xc2-xc1)
       return
      endif
 1    continue
      if(dabs(xj-sf(2,nt)).lt.dabs(xj-sf(2,1)))then
       xjf=sf(1,nt)+xj-sf(2,nt)
      else
       xjf=sf(1,1)+xj-sf(2,1)
      endif
      return
      end