program specompa
c     only simple spectra comparison
      implicit none
      integer*4 ns0,np0
      parameter (ns0=2000,np0=4000)
      character*80 sn(ns0), s80
      integer*4 is,i,j,js,jc
      real*8 as(ns0),ks(ns0),tem,bs(ns0),norm1,normj,w

      integer n1,n2,ITER,nall(ns0),ns,jsc
      real*8 s1x(np0),s1y(np0),s2x(np0),s2y(np0),RTOL,RMS,a,k,b,
     1sp(ns0),sallx(ns0,np0),sally(ns0,np0),sfit(ns0,np0),rs
      COMMON/MATVAR/s1x,s1y,s2x,s2y,n1,n2,RTOL,ITER,RMS,rs,a,k,b,
     1sallx,sally,nall,ns,sp,sfit,jsc

      real*8 wmin,wmax,aa,kk,bb,al,kin,bin,nn1
      integer*4 npar,iwr,icol,iinv,ino
      common/opt/wmin,wmax,aa,kk,npar,bb,iwr,icol,iinv,al,kin,bin,
     1ino
c      
      write(6,6000)
6000  format(/,' Specomp compares spectra listed in SP.LST',/,
     1        ' to the first spectrum in the list.',/)
      call iopar
c
      open(9,file='SP.LST')
      ns=0
1     read(9,9000,end=99,err=99)s80
9000  format(a80)
      ns=ns+1
      if(ns.gt.ns0)call report('too many spectra')
      sn(ns)=s80
      goto 1
99    close(9)
      write(6,*)'Raw input spectra used'

      write(6,6003)ns
6003  format(i4,' spectra',/,/,
     1'spec        RMS')
c
      is=1
      call ls(s1x,s1y,n1,np0,sn(is))
c     
c     norm of the first spectrum
      nn1=1.0d0
      do 100 i=1,n1-1
      w=s1x(i)
      if(w.ge.wmin.and.w.le.wmax)then
       nn1=nn1+dabs(s1y(i)*(s1x(i+1)-s1x(i)))
      endif
 100  continue
      write(6,7009)nn1
 7009 format(' Norm of the first spectrum: ',g10.4)

      open(91,file='SP.OUT')
      do 2 js=1,ns
      call ls(s2x,s2y,n2,np0,sn(js))
c
c     deviation from the first spectrum:
      RMS=0.0d0
      do 111 i=1,n1-1
      w=s1x(i)
      if(w.ge.wmin.and.w.le.wmax)then
       RMS=RMS+dabs(s2y(i)-s1y(i)*(s1x(i+1)-s1x(i)))
      endif
 111  continue
      RMS=RMS/nn1
c
      write(6,6002)js,RMS
2     write(91,6002)js,RMS
6002  format(i4,f8.4)

      close(91)
      stop
      end

      subroutine ls(xs,ys,np,np0,name)
c     loads in the spectrum      
      implicit none
      integer*4 np0,np
      real*8 xs(*),ys(*),x,y
      character*(*) name
      open(9,file=name)
      np=0
 1    read(9,*,end=99,err=99)x,y
      np=np+1
      if(np.gt.np0)call report('too many points')
      xs(np)=x
      ys(np)=y
      goto 1
 99   close(9)
      return
      end

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

      subroutine norm(n,sx,sy,an,wmin,wmax,ino)
      integer*4 n,i,ino
      real*8 an,sx(*),sy(*),wmin,wmax
      an=0.0d0
      if(ino.eq.0)then
       do 1 i=1,n-1
1      if(sx(i).ge.wmin.and.sx(i).le.wmax)
     1 an=an+sy(i)**2*abs(sx(i+1)-sx(i))
       an=dsqrt(an)
      else
       do 11 i=1,n-1
11     if(sx(i).ge.wmin.and.sx(i).le.wmax)
     1 an=an+abs(sy(i))*abs(sx(i+1)-sx(i))
      endif
      if(an.gt.1.0d-24)then
       do 2 i=1,n
 2     sy(i)=sy(i)/an
      endif
      return
      end

      subroutine iopar
      character*10 s10
      real*8 wmin,wmax,aa,kk,bb,al,kin,bin
      integer*4 npar,iwr,icol,iinv,ino
      common/opt/wmin,wmax,aa,kk,npar,bb,iwr,icol,iinv,al,kin,bin,
     1ino
      logical lex
      npar=3
      wmin=-10000.0d0
      wmax=+10000.0d0
      kk=1.0d3
      aa=1.0d1
      bb=4.0d-3
      iwr=0
      icol=0
      al=0.001d0
      iinv=0
      kin=1.0d0
      bin=0.0d0
      ino=0
      inquire(file='SP.OPT',exist=lex)
      if(lex)then
       open(9,file='SP.OPT')
1      read(9,90,end=99,err=99)s10
90     format(a10)
       if(s10(1:4).eq.'wmin')read(9,*)wmin
       if(s10(1:4).eq.'alph')read(9,*)al
       if(s10(1:4).eq.'wmax')read(9,*)wmax
       if(s10(1:2).eq.'kk')read(9,*)kk
       if(s10(1:2).eq.'aa')read(9,*)aa
       if(s10(1:2).eq.'bb')read(9,*)bb
       if(s10(1:3).eq.'bin')read(9,*)bin
       if(s10(1:3).eq.'kin')read(9,*)kin
       if(s10(1:4).eq.'npar')read(9,*)npar
       if(s10(1:4).eq.'iinv')read(9,*)iinv
       if(s10(1:3).eq.'ino')read(9,*)ino
c
c      print control:
       if(s10(1:3).eq.'iwr')read(9,*)iwr
c
c      collective spectra fitting:
       if(s10(1:4).eq.'icol')read(9,*)icol

       goto 1
99     close(9)
      endif
      write(6,6000)wmin,wmax,aa,kk,bb,al,kin,bin,npar
6000  format('wmin: ',f8.1,' wmax: ',f8.1,/,
     1'parameter restrictions:',/,
     2' aa (intenzity)         : ',f8.1,/,
     3' kk (x-scale multiplier): ',f8.1,/,
     4' bb (x-scale shift)     : ',f8.4,/,/,
     6' al (1/N, decomposition): ',f8.4,/,/,
     7' ki (syst. x-multiplier): ',f8.1,/,
     8' bi (syst. x-shift)     : ',f8.1,/,
     5' npar                   : ',i8,/)
      return
      end