program specompr
c     minimize (S-piSi) by random combinations
      implicit none
      integer*4 ns,nsn,np1,np,N,js,maxit,nt,i1,i2,i3,ii,i
      real*8 xs1,xs0,xs2,wmin,wmax,rn,
     1norm,s0,eb,er,err,tol,o
      character*10  s10
      real*8,allocatable::sx(:),sy(:),sally(:,:),c(:),cb(:),
     1sfit(:),syt(:),sf(:,:),c1(:),c2(:),c3(:)
      character*80,allocatable::sn(:)
      logical lstable,unscaled
      integer*4 ix,iy,iz
      common /randc/ ix, iy, iz
      ix=1
      iy=1000
      iz=99
      tol=1.0d-2
c      
      write(6,6000)
6000  format(/,' Specompr randomly combines spectra in SP.LST',
     1   /,   '   to reproduce the first spectrum in the list.',/)
      call iopar(xs0,xs1,xs2,wmin,wmax,maxit,lstable,unscaled)
c
c     getn number of spectra:
      ns=nsn('SP.LST')
      write(6,*)ns,' spectra in SP.LST'
      N=ns-1

      allocate(sn(ns))

c     load spectral names:
      call lsn(ns,sn)

c     get number of points in first
      np=np1(sn(1),wmin,wmax)
      write(6,*)np,' points in the first spectrum within wmin...wmax'
      if(np.lt.ns)write(6,*)' Warning: np < ns!'

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

      allocate(sx(np),sy(np),sally(N,np),c(ns),cb(ns),
     1sfit(np),syt(np),c1(ns),c2(ns),c3(ns))

c     read decomposed spectrum, calculate its norm:
      call ls(sx,sy,sn(1),wmin,wmax)
      s0=norm(sy,np)
      write(6,601)s0
601   format(' Norm0:',g12.4)

c     read rest, adjust number of points and scale:
      write(6,*)'Loading subspectra'
      do 7 js=2,ns
      call lsr(sx,syt,np,sn(js),xs0,xs1,xs2,lstable,nt,sf)
c     save it to the buffer:
7     call sav(js-1,np,syt,sally,N)
      write(6,*)ns-1,' subspectra loaded'

      call setc(c,N)
      call mkfit(N,np,sfit,c,sally,s0)
      eb=err(np,sfit,sy)
      write(6,619)eb
 619  format(' Initial error:',e12.4)

      ii=0
      do 3 i1=1,maxit
c     random combination of coefficients, first level
      call setc(c1,N)
      do 3 i2=1,maxit
c     random combination of coefficients, second level
      call modi(N,c1,c2,0.5d0)
      do 3 i3=1,maxit
c     random combination of coefficients, third level
      call modi(N,c2,c3,0.25d0)
c     fitted spectrum
      call mkfit(N,np,sfit,c3,sally,s0)
      er=err(np,sfit,sy)
      if(er.lt.eb)then
       eb=er
       cb=c3
       write(6,609)i1,i2,i3,er,(c3(i),i=1,10) 
 609   format(' Match at ',3i5,e12.4,10f6.3)
       ii=ii+1
       write(s10,80)ii
 80    format(i10)
       do 5 i=1,len(s10)
5      if(s10(i:i).ne.' ')goto 55
 55    call wrs(s10(i:10)//'fit.prn',np,sx,sfit)
      endif
3     continue
     
c     fine iterations
      o=0.5d0
99    call modi(N,cb,c,o)
      er=err(np,sfit,sy)
      if(er.lt.eb)then
       eb=er
       cb=c
       write(6,629)o,er
 629   format(' Fine natch at ',2e12.4)
       goto 99
      else
       o=o*0.5d0
       if(o.gt.tol)goto 99
      endif

      rn=0.0d0
      do 1 ii=1,N
   1  rn=rn+c(ii)
      do 2 ii=1,N
   2  c(ii)=c(ii)/rn
      call sort(c,N,sn)
      open(61,file='COEF.TXT')
      call wrc(61,N,c,sn)
      close(61)

c     fitted spectrum
      call mkfit(N,np,sfit,c,sally,s0)
      call wrs('fit.prn',np,sx,sfit)

      if(unscaled)then
       write(6,*)'Loading unscaled subspectra ...'
       do 8 js=2,ns
       call lsr(sx,syt,np,sn(js),0.0d0,1.0d0,0.0d0,.false.,0,sf)
8      call sav(js-1,np,syt,sally,N)
       write(6,*)ns-1,' .. loaded'
       call mkfit(N,np,sfit,c,sally,s0)
       call wrs('fitu.prn',np,sx,sfit)
      endif

      end

      function err(np,sfit,sy)
      implicit none
      integer*4 np,i
      real*8 sfit(*),s,err,sy(np)
      s=0.0d0
      do 1 i=1,np
  1   s=s+dabs(sy(i)-sfit(i))
      err=s
      return
      end

      subroutine wrs(s,np,x,y)
      implicit none
      character*(*) s
      integer*4 i,np
      real*8 x(*),y(*)
      open(67,file=s)
      do 1 i=1,np
1     write(67,900)x(i),y(i)
900   format(f9.2,g13.5)
      close(67)
      write(6,*)s
      return
      end

      function norm(y,np)
c     norm(s) = int |s| dx
      implicit none
      integer*4 i,np
      real*8 y(*),an,norm
      an=0.0d0
      do 1 i=1,np
1     an = an + dabs(y(i))
      norm=an
      return
      end

      subroutine sav(is,np,sy,sally,N)
      integer*4 is,i,np,N
      real*8 sy(*),sally(N,np)
      do 1 i=1,np
 1    sally(is,i)=sy(i)
      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

      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.xe1.and.xj.le.xe2).or.
     1   (xj.le.xe1.and.xj.ge.xe2))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

      subroutine lsr(xs,ys,np,name,xs0,xs1,xs2,lstable,nt,sf)
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,npall,nt
      real*8 xs(*),ys(*),xi,xj,xjp,yj,yjp,xs0,xs1,xs2,sf(2,nt),xjf
      real*8 ,allocatable::sx(:),sy(:)
      character*(*) name
      logical lstable
      npn=npall(name)
      allocate(sx(npn),sy(npn))
c     loads in the full spectrum:
      call lsf(npn,sx,sy,name)
      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
      if(xj.le.xi.and.xjp.ge.xi)then
       yj =sy(j  )
       yjp=sy(j+1)
       ys(i)=yj+(yjp-yj)/(xjp-xj)*(xi-xj)
       goto 1
      endif
 2    continue
 1    continue
      return
      end

      function np1(s,wmin,wmax)
      implicit none
      integer*4 np1,n
      real*8 wmin,wmax,x
      character*(*) s
      open(9,file=s)
      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

      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
      return
      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 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

      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 iopar(xs0,xs1,xs2,wmin,wmax,maxit,
     1lstable,unscaled)
      implicit none
      integer*4 maxit
      real*8 xs0,xs1,xs2,wmin,wmax
      logical lex,lstable,unscaled
      character*10 s10
      wmin=0.0d0
      wmax=4000.0d0
      xs0=0.0d0
      xs1=1.0d0
      xs2=0.0d0
      maxit=100
      lstable=.false.
      unscaled=.true.
      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.'wmax')read(9,*)wmax
       if(s10(1:4).eq.'wmin')read(9,*)wmin
       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:4).eq.'lsta')read(9,*)lstable
       if(s10(1:4).eq.'maxi')read(9,*)maxit
       if(s10(1:4).eq.'unsc')read(9,*)unscaled
       goto 1
99     close(9)
      endif
      write(6,6000) wmin,wmax,maxit,xs0,xs1,xs2,
     1lstable,unscaled
6000  format(
     1' wmin:                    : ',f8.1,/,
     1' wmax:                    : ',f8.1,/,
     1' maxit:                   : ',i8,/,
     1' xs0                      : ',f8.3,/,
     1' xs1                      : ',f8.3,/,
     1' xs2                      : ',f8.2,/,
     1' scaling table            : ',l8,/,
     1' unscaling at the end     : ',l8,/)

       return
      end

      subroutine wrc(io,N,c,sn)
      implicit none
      integer*4 io,N,i
      real*8 c(*)
      character*80 sn(*)
      do 3 i=1,N
3     write(io,678)i,c(i),sn(i+1)(1:40)
678   format(i4,f19.4,1x,a40) 
      return
      end

      subroutine sort(co,N,sn)
      implicit none
      real*8 co(*),ci,cmax
      integer*4 n,i,imax,j,it,is
      integer*4 ,allocatable::ind(:)
      character*80 sn(*)
      real*8,allocatable::c(:)
      allocate(c(N))
      do 5 i=1,N
5     c(i)=co(i)
      allocate(ind(N))
      do 10 i=1,N
  10  ind(i)=i
      write(6,*)'sorting'
      do 1 i=1,N-1
      ci=c(i)
      imax=i
      cmax=ci**2
      it=ind(i)
      do 2 j=i+1,N
      if(c(j)**2.gt.cmax)then
       imax=j
       cmax=c(j)**2
      endif
2     continue
      if(i.ne.imax)then
       c(i)=c(imax)
       c(imax)=ci
       ind(i)=ind(imax)
       ind(imax)=it
      endif
1     continue
      write(6,*)'Largest  '
      do 3 i=1,min(10,N)
      is=ind(i)+1
3     write(6,679)i,ind(i),c(i),sn(is)(1:50)
679   format(2i4,f19.4,1x,a50) 
      return
      end      

      subroutine mkfit(N,np,sfit,c,sally,s0)
      implicit none
      integer*4 ii,np,i,N
      real*8 sfit(*),c(*),sally(N,np),s0,s,sn
      sn=0.0d0
      do 1 ii=1,np
      s=0.0d0
      do 11 i=1,N
11    s=s+c(i)*sally(i,ii)
      sn=sn+dabs(s)
1     sfit(ii)=s
      s=s0/sn
      do 2 ii=1,np
 2    sfit(ii)=sfit(ii)*s
      return
      end

      function frandom(i)
      implicit none
      real*8 frandom
      integer*4 i
c
c     Algorithm AS 183 Appl. Statist. (1982) vol.31, no.2
c
c     Returns a pseudo-random number rectangularly distributed
c     between 0 and 1.   The cycle length is 6.95E+12 (See page 123
c     of Applied Statistics (1984) vol.33), not as claimed in the
c     original article.
c
c     IX, IY and IZ should be set to integer values between 1 and
c     30000 before the first entry.
c
      integer*4 ix, iy, iz
      common /randc/ ix, iy, iz
c
      ix = mod(171 * ix, 30269)
      iy = mod(172 * iy, 30307)
      iz = mod(170 * iz, 30323)
      frandom = mod(dble(ix) / 30269.0d0 + dble(iy) / 30307.0d0 +
     1                        dble(iz) / 30323.0d0, 1.0d0)
      return
      end

      subroutine setc(c,N)
      implicit none
      integer*4 N,i
      real*8 c(*),frandom
      do 1 i=1,N
1     c(i)=frandom(i)
      return
      end

      subroutine modi(N,c1,c2,o)
      implicit none
      integer*4 N,i
      real*8 c1(*),c2(*),o,frandom
      do 4 i=1,N
   4  c2(i)=c1(i)+o*frandom(i)
      return
      end