program specomps
      implicit none
      integer*4 ns,is,nsf,npmax,iargc,ic,i,j,jj,n
      real*8 err,sf,p,yp
      integer*4,allocatable::np(:)
      real*8,allocatable::x(:,:),y(:,:),wmin(:),wmax(:)
      character*80 list,s80
      character*80,allocatable::sn(:)
c      
      write(6,6000)
6000  format(/,' Specomps - simple comparison of spectra',/,
     1         '   to the first one in the list',/,/,
     1         '   Usage: specomps <list> <sf> (n> <wmin1> <wmax1> '
     1         '<wmin2> <wmax2> ... ',/,/,
     1         '      <list> list of spectra',/,
     1         '      <sf>   scale factor for spectra >1',/,
     1         '      <n>    number of ranges',/,
     1         '      <wmin1> wmin for the first range',/,
     1         '      <wmax1> wmax for the first range',/,
     1         '      ...',/,/)

      if(iargc().lt.2)stop
      call getarg(1,list)
      call getarg(2,s80)
      read(s80,*)sf
      call getarg(3,s80)
      read(s80,*)n
      allocate(wmin(n),wmax(n))
      do 6 i=1,n
      call getarg(2+2*i,s80)
      read(s80,*)wmin(i)
      call getarg(3+2*i,s80)
 6    read(s80,*)wmax(i)
c
      ns=1
      do 2 ic=1,2
      if(ic.eq.2)allocate(sn(ns),np(ns))
      open(9,file=list)
      ns=0
1     read(9,9000,end=99,err=99)s80
9000  format(a80)
      if(s80(1:1).eq.'#')goto 1
      write(6,9000)s80
      ns=ns+1
      if(ic.eq.2)sn(ns)=s80
      goto 1
99    close(9)
2     continue
c
      npmax=0
      do 3 is=1,ns
      np(is)=nsf(sn(is))
  3   if(np(is).gt.npmax)npmax=np(is)

      write(6,6003)ns,npmax
6003  format(i4,' spectra, maxinal number of points',i6,/,/,
     1' Spectral norms:')
      allocate(x(ns,npmax),y(ns,npmax))

      do 7 is=1,ns
      call ls(x,y,npmax,np,ns,is,sn(is))
 7    call norm(wmin,wmax,x,y,n,npmax,np,ns,is)

      do 8 is=2,ns
      do 8 i=1,np(is)
 8    x(is,i)=x(is,i)*sf
      write(6,603)2,ns,sf
 603  format(/,' X of spectra ',i4,' -',i4,' scaled by ',f10.4,/,/,
     1' Spectra errors:')

      do 4 is=1,ns
      err=0.0d0
      do 5 i=1,np(is)
      p=x(is,i)
      do 11 j=1,n
 11   if(p.ge.wmin(j).and.p.le.wmax(j))goto 12
      goto 5
c     find corresponding point in S1:
 12   do 121 jj=2,np(1)
      if(x(1,jj).ge.p.and.x(1,jj-1).le.p)goto 13
 121  if(x(1,jj).le.p.and.x(1,jj-1).ge.p)goto 13
      goto 5
13    yp=y(1,jj)
      if(dabs(x(1,jj-1)-p).lt.dabs(x(1,jj)-p))yp=y(1,jj-1)
      err=err+dabs((y(is,i)-yp)*(x(is,i)-x(is,i-1)))
  5   continue
4     write(6,609)is,err
609   format(i4,f12.6)

      end

c
      function nsf(nm)
      implicit none
      integer*4 np,nsf
      real*8 x
      character*(*) nm
      open(9,file=nm)
      np=0
 1    read(9,*,end=99,err=99)x,x
      np=np+1
      goto 1
 99   close(9)
      nsf=np
      return
      end

      subroutine norm(wmin,wmax,x,y,n,npmax,np,ns,is)
      implicit none
      integer*4 npmax,np(*),i,is,ns,n,j
      real*8 x(ns,npmax),y(ns,npmax),wmin(*),wmax(*),a
      a=0.0d0
      do 1 i=2,np(is)
      do 11 j=1,n
 11   if(x(is,i).ge.wmin(j).and.x(is,i).le.wmax(j))goto 12
      goto 1
 12   a=a+dabs(y(is,i)*(x(is,i)-x(is,i-1)))
 1    continue
      write(6,600)is,a
 600  format(i6,':',e12.4)
      do 2 i=1,np(is)
 2    y(is,i)=y(is,i)/a
      return
      end

      subroutine ls(x,y,npmax,np,ns,is,nm)
c     loads in the spectrum      
      implicit none
      integer*4 np(ns),i,is,ns,npmax
      real*8 x(ns,npmax),y(ns,npmax)
      character*(*) nm
      open(9,file=nm)
      do 1 i=1,np(is)
 1    read(9,*)x(is,i),y(is,i)
      close(9)
      return
      end