program specomp2
      implicit none
      integer ns,np,is,k,i,j,iter,istart,iend
      real*8,allocatable::s(:),m(:,:),u(:),se(:),c(:),
     1coe(:)
      real*8 wmin,wmax,dw,alpha,sum,gmin,gfree,temp,RT,glim,elim,
     1maxstep
      character*4 key
      character*80 s80
      logical lex,lwr

      write(6,5999)
5999  format(/,
     1'                  Spectral decomposion',/,/,
     1'              ref=sum(pi^2 si), penalty alpha*(pi^2-,pi^2)',/,
     1'              non-linear deviation minimization',/,/,
     1'        Input: SPECOMP.OPT options',/,
     1'               SP.LST      spectral list, first is reference',/,
     1/,
     1'       Output: test.prn    normalized reference spectrum',/,
     1'               testf.prn   fitted spectrum',/)

      alpha=0.001d0
      np=1000
      wmin=150.0d0
      wmax=1700.0d0
      iter=10
      temp=273.0d0
      glim=1.0d-3
      elim=1.0d-4
      maxstep=0.1d0
      lwr=.false.
      inquire(file='SPECOMP2.OPT',exist=lex)
      if(lex)then
       open(9,file='SPECOMP2.OPT')
991    read(9,900,end=95,err=95)key
900    format(a4)
       if(key.eq.'ALPH')read(9,*)alpha
       if(key.eq.'WMIN')read(9,*)wmin
       if(key.eq.'WMAX')read(9,*)wmax
       if(key.eq.'ITER')read(9,*)iter
       if(key.eq.'TEMP')read(9,*)temp
       if(key.eq.'GLIM')read(9,*)glim
       if(key.eq.'ELIM')read(9,*)elim
       if(key(1:3).eq.'LWR')read(9,*)lwr
       if(key(1:2).eq.'NP')read(9,*)np
       goto 991
95     close(9)
       write(6,*)' SPECOMP2.OPT read'
      endif
      inquire(file='SP.LST',exist=lex)
      if(lex)then
       ns=0
       open(9,file='SP.LST')
992    read(9,*,end=98,err=98)
       ns=ns+1
       goto 992
98     close(9)
       ns=ns-1
       write(6,*)'SP.LST found, ',ns,' spectra'
      else
       write(6,*)'SP.LST not found'
       stop
      endif
      dw=(wmax-wmin)/dble(np-1)
      RT=temp*8.314d0/4.184d0/1000.0d0
      write(6,5998)alpha,wmin,wmax,dw,np,temp,elim,glim,maxstep,lwr,
     1iter
5998  format(' penalty function prefactor          :',g12.4,/,
     1       ' frequency interval for decomposition:',/,
     1       '                                 wmin:',f12.2,/,
     1       '                                 wmax:',f12.2,/,
     1       '                                   dw:',f12.2,/,
     1       ' number of points in compared spectra:',i12,/,
     1       ' (input spectra will be interpolated)',/,
     1       '                      temperature (K):',f12.2,/,
     1       '                dev. change threshold:',g12.4,/,
     1       '               grad. change threshold:',g12.4,/,
     1       '                          max. change:',g12.4,/,
     1       '                                  lwr:',l12,/,
     1       '         maximum number of iterations:',i12,/)

      allocate(s(ns*np),m(ns,ns),u(np),c(ns),se(np),coe(ns))

      open(39,file='SP.LST')
c     read reference:
      read(39,9980)s80
9980  format(a80)
      call trim(s80,istart,iend)
      call reads(wmin,wmax,dw,se,np,s80(istart:iend))
c     normalize  sum |se_i|=1
      call qnorm(np,se)
      open(59,file='test.prn')
      do 99 i=1,np
99    write(59,990)wmin+(i-1)*dw,se(i)
990   format(2G12.4)
      close(59)
      write(6,*)' reference read,'
      write(6,*)' normalized spectrum written to test.prn'
      write(6,*)
      do 1 is=1,ns
      read(39,9980)s80
      call trim(s80,istart,iend)
      call reads(wmin,wmax,dw,u,np,s80(istart:iend))
      do 1 k=1,np
1     s(k+(is-1)*np)=u(k)
      close(39)
      write(6,*)' sub-spectra read'
c
c     normalize all 
      sum=0.0d0
      do 11 i=1,ns*np
11    sum=sum+s(i)**2
      sum=sum/dble(ns)
      do 12 i=1,ns*np
12    s(i)=s(i)/dsqrt(sum)
      write(6,*)' sub-spectra normalized'

c     correlation matrix:
      do 3 i=1,ns
      do 3 j=i,ns
      sum=0.0d0
      do 40 k=1,np
40    sum=sum+s(k+(i-1)*np)*s(k+(j-1)*np)
      m(j,i)=sum
3     m(i,j)=sum
c     scalar product:
      do 6 i=1,ns
      c(i)=0.0d0
      do 6 k=1,np
6     c(i)=c(i)+se(k)*s(k+(i-1)*np)
      write(6,*)'Correlation matrix done'

      call minimini(ns,coe,s,se,c,m,alpha,np,iter,glim,elim,maxstep)

c     write composed spectrum:
      open(9,file='testf.prn')
      do 39 i=1,np
      sum=0.0d0
      do 391 j=1,ns
391   sum=sum+coe(j)*s(i+(j-1)*np)
39    write(9,901)wmin+dble(i-1)*dw,sum
901   format(f12.2,1x,g20.4)
      close(9)
      write(6,*)' fitted spectrum written to testf.prn'

c     normalize coefficients:
      call norm(ns,coe)

c     find lowest dG
      gmin=0.0d0
      do 10 is=1,ns
      if(coe(is).gt.1.0d-200)then
       gfree=-log(coe(is))
      else
       gfree=-log(1.0d-200)
      endif
      if(is.eq.1.)gmin=gfree
10    if(gfree.lt.gmin)gmin=gfree
c     write to file, also relative -RT ln p
      open(9,file='test.txt')
      do 7 is=1,ns
      if(coe(is).gt.1.0d-200)then
       gfree=-log(coe(is))
      else
       gfree=-log(1.0d-200)
      endif
7     write(9,6001)is,coe(is),RT*(gfree-gmin)
6001  format(i5,2f12.6)
      close(9)
      write(6,*)' coefficients and free energies in test.txt'

      stop

      end

      subroutine reads(wmin,wmax,dw,u,np,fn)
      implicit none
      character*(*) fn
      integer*4 np,i,j,nmin,nmax
      real*8 u(*),wmin,wmax,dw,x,xn,y,yn,w,wn,s,sn,ws
      do 2 i=1,np
2     u(i)=0.0d0
      i=0
      xn=0.0d0
      yn=0.0d0
      open(9,file=fn)
1     read(9,*,end=99,err=99)x,y
      i=i+1
      if(i.gt.1)then
       if((x.ge.wmin.and.x.le.wmax).or.(xn.ge.wmin.and.xn.le.wmax))then
        if(xn.gt.x)then
         w=xn
         s=yn
         wn=x
         sn=y
        else
         wn=xn
         sn=yn
         w=x
         s=y
        endif
        nmin=int((wn-wmin)/dw)+1
        nmax=int((w-wmin)/dw)+2
        do 3 j=nmin,nmax
        ws=wmin+dble(j-1)*dw
3       if(ws.ge.wn.and.ws.le.w)u(j)=sn+(s-sn)*(ws-wn)/(w-wn)
       endif
      endif
      xn=x
      yn=y
      goto 1
99    close(9)
      return
      end

      subroutine norm(n,s)
      implicit none
      integer*4 n,i
      real*8 s(*),sum
      sum=0.0d0
      do 1 i=1,n
1     sum=sum+dabs(s(i))
      do 2 i=1,n
2     s(i)=s(i)/sum
      return
      end

      subroutine qnorm(n,s)
      implicit none
      integer*4 n,i
      real*8 s(*),sum
      sum=0.0d0
      do 1 i=1,n
1     sum=sum+s(i)**2
      sum=1.0d0/dsqrt(sum)
      do 2 i=1,n
2     s(i)=s(i)*sum
      return
      end

      function dd(coe,s,ns,se,c,m,g,alpha,np,ic,gn)
      implicit none
      integer*4 ns,i,k,j,np,ic
      real*8 dd,coe(*),s(*),se(*),c(*),m(ns,ns),g(*),sum,sump,
     1av,alpha,ck,sd,gn

c     coeff average:
      av=0.0d0
      do 2 j=1,ns
2     av=av+coe(j)**2
      av=av/dble(ns)

c     sum:
      sum=0.0d0
      do 7 i=1,ns
7     sum=sum+(coe(i)**2-av)**2
      sum=sum*alpha
      do 6 i=1,np
      sd=0.0d0
      do 51 j=1,ns
51    sd=sd+coe(j)**2*s(i+(j-1)*np)
6     sum=sum+(se(i)-sd)**2
      dd=sum

      if(ic.eq.0)return

c     gradient g, its norm gn:
      gn=0.0d0
      do 1 k=1,ns
      sump=0.0d0
      do 3 j=1,ns
3     sump=sump+coe(j)**2*m(j,k)
      av=av/dble(ns)
      ck=coe(k)
      g(k)=4.0d0*(sump-c(k)+alpha*(ck**2-av))*ck
1     gn=gn+g(k)**2
      gn=dsqrt(gn)
      return
      end

      subroutine asa(g0,g,ns)
      implicit none
      integer*4 i,ns
      real*8 g0(*),g(*)
      do 101 i=1,ns
101   g0(i)=g(i)
      return
      end

      subroutine minimini(ns,coe,s,se,c,m,alpha,np,nstepmax,
     1glim,elim,maxstep)
      implicit none
      integer*4 nstep,lstep,ns,i,np,nstepmax
      real*8 e0,gnorm0,glim,etot,e00,dh,maxstep,eold,h,hh0,elim,e1,
     1coe(*),gn,s(*),se(*),c(*),m(ns,ns),alpha,dd,secder,
     2gradh,enew,hh1,hh2,e2,gnorm1
      real*8,allocatable::g0(:),h0(:),x0(:),g(:)
      allocate (g0(ns),h0(ns),x0(ns),g(ns))
      do 201 i=1,ns
201   g(i)=0.0d0
      gn=0.0d0
      gnorm1=0.0d0
      gnorm0=0.0d0
      do 36 i=1,ns
36    coe(i)=dsqrt(1.0d0/dble(ns))
      nstep=0

      etot=dd(coe,s,ns,se,c,m,g,alpha,np,1,gnorm0)
      write(6,6009)nstep,etot,gnorm0
6009  format(' Iteration',i5,' sig ',g12.4,' grad ',g12.4)
      e0=etot
      call asa(g0,g,ns)
      call asa(h0,g,ns)
      call asa(x0,coe,ns)
      if(gnorm0.le.glim)goto 300

200   e00=etot
      dh=maxstep/gnorm0
      secder=-1.0d0
      nstep=nstep+1
      write(6,555)nstep
555   format(' Iteration ',i4)
      eold=etot
      e0=eold
      h=0
      hh0=0
      lstep=0

100   lstep=lstep+1
      h=h+dh
      do 101 i=1,ns
101   coe(i)=x0(i)-h*h0(i)
      etot=dd(coe,s,ns,se,c,m,g,alpha,np,0,gn)
      enew=etot
c     end of linear search:
      if(dabs(enew-eold).lt.0.5d0*elim)goto 400
      if(dabs((enew-eold)/dh).lt.0.5d0*glim*gnorm0)goto 400
      if(lstep.eq.1)then
        e1=etot
        hh1=h
      else
       if(lstep.eq.2)then
        e2=etot
        hh2=h
       else
        e0=e1
        hh0=hh1
        e1=e2
        hh1=hh2
        e2=etot
        hh2=h
       endif
       if(hh1.eq.hh0.or.hh2.eq.hh0.or.hh2.eq.hh1)then 
        secder=-1.d0 
       else
        secder=-2.0d0*((hh2-hh1)*(e1-e0)-(e2-e1)*(hh1-hh0))
     1  /((hh2-hh0)*(hh1-hh0)*(hh2-hh1))
        gradh=(e0-e1-secder/2.0d0*(hh0**2-hh1**2))/(hh0-hh1)
       endif
      endif

      if(secder.gt.0.0d0)then
       dh=-gradh/secder-h
       if(dabs(dh).gt.maxstep/gnorm0)then 
        if(dh.gt.0.0d0) then
         dh=maxstep/gnorm0 
        else 
         dh=-maxstep/gnorm0
        endif
       endif
       eold=enew
       goto 100
      endif
      if(enew.gt.eold) then 
       eold=enew
       dh=-dh/2
       goto 100
      endif
      eold=enew
      goto 100

400   write(6,4000)etot,lstep
4000  format(10x,'e = ',g12.4,' after ',i6,' linear steps')
      etot=dd(coe,s,ns,se,c,m,g,alpha,np,1,gnorm1)
      call asa(x0,coe,ns)
      do 102 i=1,ns
102   h0(i)=g(i)+(gnorm1/gnorm0)**2*h0(i)
      call asa(g0,g,ns)
      write(6,6005)gnorm1,glim
      write(6,6006)dabs(e00-etot),elim
6005  format(10x,'Gradient        ',g12.4,'; limit ',g14.4)
6006  format(10x,'Change of energy',g12.4,'; limit ',g14.4)
      if(dabs(e00-etot).lt.elim.or.gnorm1.lt.glim)goto 300
      e0=etot
      gnorm0=gnorm1
      if(nstep.gt.nstepmax)then 
       write(6,*)'Too many steps'
       goto 300
      endif
      goto 200

300   do 373 i=1,ns
373   coe(i)=coe(i)**2
      return
      end

      subroutine trim(fn,is,ie)
      implicit none
      character*80 fn
      integer*4 is,ie
      do 1 is=1,len(fn)
1     if(fn(is:is).ne.' ')goto 2
2     do 3 ie=len(fn),1,-1
3     if(fn(ie:ie).ne.' ')goto 4
4     return
      end