program specompG
      implicit none
      integer*4 ns,nsn,N,np1,np,nt,js,ix,nk,nn,maxit,ndim,nd,ipp1,ipp2
      real*8 wmin,wmax,na,norm,xs0,xs1,xs2,an,dev,n1,g,
     1xmin,xmax,ymin,ymax,temp
      real*8,allocatable::sf(:,:),sx(:),sy(:),sally(:,:),aa(:),syt(:),
     1c(:),sfit(:),xy(:,:),salln(:,:),xyn(:,:),pxy(:),ran(:,:)
      integer*4,allocatable::ii(:)
      character*80,allocatable::sn(:)
      logical lstable,lcon,lwr,leq,lran
      write(6,6000)
6000  format(/,'  SpecompG compares spectra listed in SP.LST to the',/,
     1        ' first spectrum in the list, makes a free-energy (G)',/,
     1        ' map smoothed by plane waves,',/,
     1        '            S1=sum (i>1) pi Si, pi=exp(-Gi/kT',/,
     1        '            Gi = cik gk, gk ... plane waves',/,/,
     1        '      Input: SP.LST .. list of spectra',/,
     1        '             XY.LST .. xy coordinates',/,
     1        '             SP.OPT .. options',/,
     1        '     Output: COEF.TXT .. coefficients (pi)',/,
     1        '             EMAP.TXT .. map of energies and prob.')

      call readopt(wmin,wmax,lstable,xs0,xs1,xs2,lcon,g,nk,lwr,
     1nn,xmin,xmax,ymin,ymax,maxit,temp,ndim,leq,nd,ipp1,ipp2)
      allocate(pxy(ndim),ran(ndim,2))
      call readopt2(pxy,ndim,ran,lran)

c     getn number of spectra:
      ns=nsn('SP.LST')
      write(6,*)ns,' spectra in SP.LST'
      N=ns-1

c     load coordinates
      allocate(xy(N,ndim))
      call lcn(N,xy,ndim)

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

c     get number of points within wmin wmax in first
      write(6,*)wmin,wmax
      np=np1(sn(1),wmin,wmax)
      write(6,*)np,' points in the first spectrum'
      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)
      else
       nt=1
       allocate(sf(2,nt))
      endif

c     read decomposed spectrum
      allocate(sx(np),sy(np))
      call ls(sx,sy,sn(1),wmin,wmax)
      n1=norm(np,sy)
      write(6,599)n1
 599  format(' Decomposed spectrum norm:',e12.4)
c     normalize
      call nspe(np,sy,n1)

c     read rest, adjust number of points and scale:
      write(6,761)
761   format(' Loading subspectra')
      allocate(sally(N,np),syt(np),ii(np),aa(np))
      ii=0
      aa=0.0d0
      na=0.0d0
      do 7 js=2,ns
      call lsr(ns,sx,syt,np,js,sn,xs0,xs1,xs2,lstable,nt,sf,lcon,g,
     1ii,aa)
      if(js.eq.2.and.lwr)write(6,76)
76    format(' Spectrum    norm  coords   deviation from the first:')
      an=norm(np,syt)
      na=na+an
      dev=0.0d0
      do 77 ix=1,np
77    dev=dev+dabs(sy(ix)-syt(ix)/an)
      if(lwr)then
       write(6,78)js,an
78     format(i6,e12.4,$)
       do 80 ix=1,ndim
80     write(6,481)xy(js-1,ix)
481    format(f12.6,$)
       write(6,482)dev
482    format(f12.6)
      endif
c     save it to the buffer:
7     call sav(js-1,np,syt,sally,N)
      na=na/dble(N)
      write(6,598)N,na
 598  format(i6,' subspectra loaded, average norm:',E12.4)

c     scale subspectra by everage norm
      do 79 js=1,N
      do 79 ix=1,np
79    sally(js,ix)=sally(js,ix)/na
c
c     simulate spectra equivalently spaced in a coordinate
      if(leq)then
       allocate(salln(nd**ndim,np),xyn(nd**ndim,ndim))
       deallocate(sn)
       allocate(sn(1+nd**ndim))
       call seteq(N,ndim,np,nd,nk,xy,lwr,sally,salln,xyn,lran,ran,
     1 sn,pxy)
c      re-define N:
       N=nd**ndim
       allocate(c(N),sfit(np))
       call cg(N,np,c,sy,sfit,salln,n1,nk,sx,xyn,pxy,
     1 nn,xmin,xmax,ymin,ymax,lwr,maxit,temp,ndim,ipp1,ipp2)
      else
       allocate(c(N),sfit(np))
       call cg(N,np,c,sy,sfit,sally,n1,nk,sx,xy,pxy,
     1 nn,xmin,xmax,ymin,ymax,lwr,maxit,temp,ndim,ipp1,ipp2)
      endif
      call sort(c,N,sn)
      call wrc(N,c,sn)
      end

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

      function ddp(x,y,p)
      real*8 ddp,x,y,p,a
      a=dabs(x-y)
      if(dabs(x-y-p).lt.a)a=dabs(x-y-p)
      if(dabs(x-y+p).lt.a)a=dabs(x-y+p)
      ddp=a
      return
      end

      subroutine seteq(N,ndim,np,nd,nk,xy,lwr,sally,salln,xyn,lran,ran,
     1sn,pxy)
      implicit none
      integer*4 N,ndim,np,nd,nk,ip,ic,i,j,k,nst,ii
      real*8 xy(N,ndim),q,dn,salln(nd**ndim,np),pxy(ndim),xt,
     1sally(N,np),xyn(nd**ndim,ndim),cc,ran(ndim,2),dx,ddp
      real*8,allocatable::delta(:),xa(:),xi(:),x1(:),x2(:),xc(:)
      integer*4,allocatable::kk(:)
      logical lwr,lran
      character*80 sn(*)
      allocate(delta(ndim),xa(ndim),xi(ndim),x1(ndim),x2(ndim),
     1xc(ndim))
      write(6,61)nd**ndim
61    format(i6,' equivalently-spaced subspectra requested')
      do 20 i=1,nd**ndim
      do 20 j=1,80
20    sn(i)(j:j)=' '
      if(nd.lt.nk)call report('nd < nk!')
      do 1 i=1,ndim
      xi(i)=xy(1,i)
      xa(i)=xy(1,i)
      do 2 j=1,N
      if(xy(j,i).lt.xi(i))xi(i)=xy(j,i)
2     if(xy(j,i).gt.xa(i))xa(i)=xy(j,i)
      delta(i)=(xa(i)-xi(i))/dble(nd-1)
1     write(6,60)i,xi(i),xa(i),delta(i)
60    format(' Dimension',i2,', min:',e12.4,', max:',e12.4,
     1', step:',e12.4)
      if(lran)then
       write(6,*)'Replaced by pre-defined ranges:'
       do 19 i=1,ndim
       xi(i)=ran(i,1)
       xa(i)=ran(i,2)
       delta(i)=(xa(i)-xi(i))/dble(nd-1)
19     write(6,60)i,xi(i),xa(i),delta(i)
      endif
      salln=0.0d0
      allocate(kk(ndim))
c     loop over the new spectra:
      do 5 ii=1,nd**ndim
      call setkk(ndim,nd,ii,kk)
      write(sn(1+ii)(1:25),611)ii
 611  format('Arbitrary spectrum ',i6)
c     xc,x1,x2: center and borders of ndim dimensional cube:
c     xyn: new coordinates of spectrum ii
      do 6 j=1,ndim
      xc(j)=xi(j)+delta(j)*dble(kk(j)-1)
      x1(j)=xc(j)-delta(j)/2.0d0
      x2(j)=xc(j)+delta(j)/2.0d0
      write(6,612)x1(j),x2(j)
 612  format(f7.1,'..',f7.1,$)
6     xyn(ii,j)=xc(j)
      write(6,*)
      write(6,609)ii
 609  format(i6,$)
      do 7 j=1,ndim
7     write(6,610)kk(j)
 610  format(i4,$)
c     nst: number of supspectra in a brick around xc:
      nst=0
      q=0.0d0
      do 8 j=1,N
      do 9 k=1,ndim
      xt=xy(j,k)
9     if((xt       .lt.x1(k).or.xt       .gt.x2(k)).and.
     1   (xt+pxy(k).lt.x1(k).or.xt+pxy(k).gt.x2(k)).and.
     1   (xt-pxy(k).lt.x1(k).or.xt-pxy(k).gt.x2(k)))goto 8
      nst=nst+1
      dn=0.0d0
      do 10 k=1,ndim
      dx=ddp(xy(j,k),xc(k),pxy(k))
10    dn=dn+(dx/delta(k))**2
      dn=dsqrt(dn)
      q=q+1.0d0/(dn+0.25d0)
8     continue
      if(nst.gt.0)then
       if(lwr)write(6,62)nst
62     format(' from',i5,' originals')
       do 11 j=1,N
       do 12 k=1,ndim
       xt=xy(j,k)
12     if((xt       .lt.x1(k).or.xt       .gt.x2(k)).and.
     1    (xt+pxy(k).lt.x1(k).or.xt+pxy(k).gt.x2(k)).and.
     1    (xt-pxy(k).lt.x1(k).or.xt-pxy(k).gt.x2(k)))goto 11
       dn=0.0d0
       do 13 k=1,ndim
       dx=ddp(xy(j,k),xc(k),pxy(k))
13     dn=dn+(dx/delta(k))**2
       dn=dsqrt(dn)
       do 14 ip=1,np
14     salln(ii,ip)=salln(ii,ip)+1.0d0/(dn+0.25d0)*sally(j,ip)/q
c      salln .. average spectrum in the rectangle around xc
11     continue
      else
c      if nst=0, find closest spectrum from all:
       ic=1
       cc=0.0d0
       do 15 k=1,ndim
       dx=ddp(xy(1,k),xc(k),pxy(k))
15     cc=cc+dx**2
       do 17 j=1,N
       dn=0.0d0
       do 16 k=1,ndim
       dx=ddp(xy(j,k),xc(k),pxy(k))
16     dn=dn+dx**2
       if(dn.lt.cc)then
        ic=j
        cc=dn
       endif
17     continue
       if(lwr)write(6,63)ic
63     format(' from number',i5)
       do 18 ip=1,np
18     salln(ii,ip)=sally(ic,ip)
      endif
5     continue
      return
      end

      subroutine  setkk(ndim,nk,k,kk)
      implicit none
      integer*4 ndim,nk,k,kk(ndim),ki,up,ii
      kk=0
      do 101 ki=ndim,1,-1
      up=k+nk**(ki-1)-1
      do 100 ii=ki,ndim-1
100   up=up-nk**ii*(kk(ii+1)-1)
101   kk(ki)=up/nk**(ki-1)
c     kk(1) .. fastest index, kk(ndim) .. slowest
      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 lcn(nc,xy,ndim)
      implicit none
      integer*4 j,nc,ndim,ix
      real*8 xy(nc,ndim)
      open(9,file='XY.LST')
      do 1 j=1,nc
1     read(9,*)(xy(j,ix),ix=1,ndim)
      close(9)
      write(6,*)'coordinates loaded'
      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 readopt(wmin,wmax,lstable,xs0,xs1,xs2,lcon,g,nk,
     1lwr,nn,xmin,xmax,ymin,ymax,maxit,temp,ndim,leq,nd,ipp1,ipp2)
      implicit none
      integer*4 nk,nn,maxit,ndim,nd,ipp1,ipp2,i
      real*8 wmin,wmax,xs0,xs1,xs2,g,xmin,xmax,ymin,ymax,temp,num
      logical lex,lstable,lcon,lwr,leq
      character*10 s10,s10c
      ipp1=1
      ipp2=2
      ndim=2
      nn=20
      lstable=.false.
      leq=.false.
      wmin=0.0d0
      wmax=4000.0d0
      lcon=.false.
      xs0=0.0d0
      xs1=1.0d0
      xs2=0.0d0
      g=20.0d0
      nk=5
      nd=25
      lwr=.false.
      xmin=0.0d0
      ymin=0.0d0
      xmax=360.0d0
      ymax=360.0d0
      maxit=100
      temp=293.15d0
c
c      Option (default)      
c      ----------------
c      maxi (maxit=100)      maximum number of iterations
c      wmin (wmin=0 cm-1)    start of active spectral interval 
c      wmax (wmax=4000 cm-1) end of active spectral interval 
c      xmin (xmin=0)         xmin in degrees. for plotting
c      xmax (xmax=360)       xmax 
c      ymin (ymin=0)         ymin
c      ymax (ymax=360)       ymax
c      ipp1 (ipp1=1)         coordinate for 2D plot
c      ipp2 (ipp1=2)         coordinate for 2D plot
c      nn   (nn=20)          number of values in one dimension, for plotting
c      stab (lstable=false)  scaling in STABLE.TXT
c      equi (leq=false)      simulate equivalent spectra distance in coords
c      ndd  (nd=25)          number of spectra for equivalent distance
c      conv (lcon=false)     convolute spectra with Gaussians (see "ggg")
c      ggg  (g=20 cm-1)      bandwidth for convolution (conv)
c      temp (temp=293.15)    temperature in K
c      xs0  (xs0=0)          x-scaling x'=xs0 + xs1*x + xs2*x^2
c      xs1  (xs1=1) 
c      xs2  (xs2=0)
c      lwr  (lwr=false)      extended output
c      nk   (nk=5)           number of plane waves in one dimension
c      ndim (ndim=2)         dimension of the coordinate space

      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)
       write(6,90)s10
       if(s10(1:4).eq.'maxi')read(9,*)maxit
       if(s10(1:4).eq.'wmax')read(9,*)wmax
       if(s10(1:4).eq.'wmin')read(9,*)wmin
       if(s10(1:4).eq.'xmax')read(9,*)xmax
       if(s10(1:4).eq.'xmin')read(9,*)xmin
       if(s10(1:4).eq.'ymax')read(9,*)ymax
       if(s10(1:4).eq.'ymin')read(9,*)ymin
       if(s10(1:4).eq.'stab')read(9,*)lstable
       if(s10(1:4).eq.'lsta')read(9,*)lstable
       if(s10(1:4).eq.'conv')read(9,*)lcon
       if(s10(1:4).eq.'equi')read(9,*)leq
       if(s10(1:4).eq.'temp')read(9,*)temp
       if(s10(1:4).eq.'ndim')read(9,*)ndim
       if(s10(1:4).eq.'ipp1')read(9,*)ipp1
       if(s10(1:4).eq.'ipp2')read(9,*)ipp2
       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.'ggg')read(9,*)g
       if(s10(1:3).eq.'lwr')read(9,*)lwr
       if(s10(1:3).eq.'ndd')read(9,*)nd
       if(s10(1:2).eq.'nk')read(9,*)nk
       if(s10(1:2).eq.'nn')read(9,*)nn

c      to be read by readopt2:
       if(s10(1:4).eq.'peri')read(9,*)
       if(s10(1:4).eq.'lran')read(9,*)
       if(s10(1:4).eq.'rang')then
        do 66 i=1,ndim
66      read(9,*)
       endif

       backspace 9
       read(9,90)s10c
       if(s10.eq.s10c)call report(' Unknown option '//s10)
       goto 1
99     close(9)
      endif
      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


      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)
      implicit none
      integer*4 i,np
      real*8 norm,s(*)
      norm=0.0d0
      do 1 i=1,np
   1  norm=norm+dabs(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

      subroutine lsr(ns,xs,ys,np,js,sn,xs0,xs1,xs2,lstable,nt,sf,
     1lcon,g,ii,aa)
c     loads in the full spectrum, adjust number of points to np
c     with the scaled xs
      implicit none
      integer*4 np,i,js,npn,j,nt,npall,ns,ii(np)
      real*8 xs(*),ys(*),xi,xj,xjp,xs0,xs1,xs2,sf(2,nt),xjf,g,
     1ef,pi,aa(np)
      real*8 ,allocatable::sx(:),sy(:),sp(:)
      character*80 sn(ns)
      logical lstable,lcon
      pi=4.0d0*datan(1.0d0)
      npn=npall(sn(js))
      allocate(sx(npn),sy(npn),sp(npn))
c     loads in the full spectrum:
      call lsf(npn,sx,sy,sn(js))
c     convolute
      if(lcon)then
       do 3 i=1,npn
       sp(i)=0.0d0
       do 31 j=1,npn-1
       ef=((sx(i)-sx(j))/g)**2
31     if(ef.lt.20.0d0)sp(i)=sp(i)+exp(-ef)*sy(j)*dabs(sx(j)-sx(j+1))
3      sp(i)=sp(i)/dsqrt(pi)/g
       sy=sp
      endif
c     make np points only:
      if(js.eq.2)then
       write(6,*)' Resampling x-scale for 2nd spectrum'
       do 1 i=1,np
       xi=xs(i)
       ys(i)=0.0d0
       ii(i)=0
       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
        ii(i)=j
        aa(i)=(xi-xj)/(xjp-xj)
        ys(i)=sy(j)+(sy(j+1)-sy(j))*aa(i)
        goto 1
       endif
 2     continue
 1     continue
      else
       do 4 i=1,np
       j=ii(i)
       if(j.ne.0)then
        ys(i)=sy(j)+(sy(j+1)-sy(j))*aa(i)
       else
        ys(i)=0.0d0
       endif
 4     continue
      endif
      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


      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

      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

      subroutine inic(N,c)
      implicit none
      integer*4 N,i
      real*8 c(N)
      do 1 i=1,N
1     c(i)=1.0d0/dble(N)
      return
      end

      subroutine mkfit(N,np,sfit,p,sally)
      implicit none
      integer*4 ii,np,i,N
      real*8 sfit(*),p(*),sally(N,np)
      do 3 ii=1,np
      sfit(ii)=0.0d0
      do 3 i=1,N
3     sfit(ii)=sfit(ii)+p(i)*sally(i,ii)
      return
      end

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

      subroutine setg(m,N,np,g,sy,sfit,sally,p,gf,gnorm)
      implicit none
      integer*4 N,i,j,m,np
      real*8 g(m),gnorm,s,sy(np),sfit(np),p(N),gf(m,N),sally(N,np)
      real*8,allocatable::intg(:)
      allocate(intg(N))
      intg=0.0d0
      do 3 j=1,N
      do 3 i=1,np
3     intg(j)=intg(j)+(sfit(i)-sy(i))*sally(j,i)
      gnorm=0.0d0
      do 1 i=1,m
      s=0.0d0
      do 2 j=1,N
   2  s=s+p(j)*gf(i,j)*intg(j)
      s=-2.0d0*s
      gnorm=gnorm+s**2
  1   g(i)=s
      gnorm=dsqrt(gnorm)
      return
      end

      function ggf(k1,x,px)
      implicit none
      integer*4 k1,kk1
      real*8 x,px,ggf,pi
      pi=4.0d0*datan(1.0d0)
      kk1=k1/2
      if(mod(k1,2).eq.0)then
       ggf=sin(dble(kk1)*x/px*2.0d0*pi)
      else
       ggf=cos(dble(kk1)*x/px*2.0d0*pi)
      endif
      return
      end

      subroutine fb(gf,m,N,xy,pxy,nk,ndim)
      implicit none
      integer*4 N,m,i,nk,k,ndim,ki
      real*8 gf(m,N),pxy(ndim),xy(N,ndim),ggf,pr
      integer*4 ,allocatable::kk(:)
      allocate(kk(ndim))
      do 2 k=1,m
      call setkk(ndim,nk,k,kk)
      do 2 i=1,N
      pr=1.0d0
      do 102 ki=1,ndim
102   pr=pr*ggf(kk(ki),xy(i,ki),pxy(ki))
2     gf(k,i)=pr       
      return
      end

      function gfx(nk,k,x,y,pxy,ndim,ipp1,ipp2)
      implicit none
      integer*4 nk,k,ndim,ipp1,ipp2
      real*8 gfx,x,y,pxy(ndim),ggf
      integer*4 ,allocatable::kk(:)
      allocate(kk(ndim))
      call setkk(ndim,nk,k,kk)
      if(ndim.eq.1)then
       gfx=ggf(kk(1),x,pxy(1))
      else
       gfx=ggf(kk(ipp1),x,pxy(ipp1))*ggf(kk(ipp2),y,pxy(ipp2))
      endif
      return
      end

      subroutine mkp(N,m,c,gf,p)
      implicit none
      integer*4 N,m,i,k
      real*8 c(m),gf(m,N),p(N),dG
      do 1 i=1,N
      dG=0.0d0
      do 2 k=1,m
2     dG=dG+c(k)*gf(k,i)
1     p(i)=exp(-dG)
      return
      end

      subroutine cg(N,np,p,sy,sfit,sally,n1,nk,
     1sx,xy,pxy,nn,xmin,xmax,ymin,ymax,lwr,maxit,temp,ndim,ipp1,ipp2)
      implicit none
      integer*4 N,nstep,lstep,i,maxit,np,nk,m,nn,ndim,ipp1,ipp2
      real*8 p(N),maxstep,sy(np),sfit(np),sx(np),ef,temp,
     1e0,gnorm0,e00,dh,secder,eold,h,hh0,sally(N,np),hh1,e2,hh2,
     1e1,gradh,n1,etot,enew,gnorm1,xy(N,ndim),pxy(ndim),xmin,xmax,
     1ymin,ymax
      real*8,allocatable::p0(:),g0(:),h0(:),g(:),c0(:),gf(:,:),ck(:)
      logical lwr
      maxstep=0.1d0
      nstep=0
c     number of parameters to optimize:
      m=nk**ndim
      allocate(p0(N),g0(m),h0(m),g(m),c0(m),gf(m,N),ck(m))
c     fill basis
      call fb(gf,m,N,xy,pxy,nk,ndim)
c     initial values of probabilities:
      call inic(N,p0)
c     initial values of pw coefficients:
      c0=0.0d0
      c0(1)=-log(p0(1))
c     fitted spectrum Sfit=sum(i) p0(i) Si
      call mkfit(N,np,sfit,p0,sally)
c     write initial average, scaled:
      call nspe(np,sfit,1.0d0/n1)
      call wrs('aves.prn',np,sx,sfit)
      call nspe(np,sfit,n1)

c     "energy" to be minimized"
      etot=ef(sy,sfit,np)
c     gradient:
      call setg(m,N,np,g0,sy,sfit,sally,p0,gf,gnorm0)
      e0=etot
      h0=g0

222   e00=etot
      dh=maxstep/gnorm0
      secder=-1.0d0
      nstep=nstep+1
      if(lwr)write(6,600)nstep,etot
 600  format(' Linear search at point ',i6,', E = ',e12.4)
      eold=etot
      e0=etot
      h=0.0d0
      hh0=0.0d0
      lstep=0

111   lstep=lstep+1
      h=h+dh
      do 1 i=1,m
1     ck(i)=max(0.0d0,c0(i)-h*h0(i))
      call mkp(N,m,ck,gf,p)
      call mkfit(N,np,sfit,p,sally)
      etot=ef(sy,sfit,np)
      enew=etot
c     possibly end of linear search?:
      if(lstep.gt.100)then
       write(6,*)'Too many linear steps'
       goto 444
      endif
      if(dabs(enew-eold)<1.0d-5*eold)goto 444
      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.0d0
       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 111
      endif
      if(enew.gt.eold)then
       eold=enew
       dh=-dh/2
       goto 111
      endif
      eold=enew
      goto 111

444   if(lwr)write(6,601)lstep,etot
 601  format(i6,' linear steps ',e12.4)
      call mkp(N,m,ck,gf,p)
      call mkfit(N,np,sfit,p,sally)
      etot=ef(sy,sfit,np)
      call setg(m,N,np,g,sy,sfit,sally,p,gf,gnorm1)
      c0=ck
      do 2 i=1,m
      h0(i)=g(i)+(gnorm1/gnorm0)**2*h0(i)
   2  g0(i)=g(i)
      write(6,602)gnorm1,etot,etot-e00
 602  format(' Grad:',e12.4,' E:',e12.4,' dE:',e12.4)
      call nspe(np,sfit,1.0d0/n1)
      call wrs('fit.prn',np,sx,sfit)
      if(dabs(e00-etot).lt.1.0d-5*e00)then
       write(6,*)'Optimization converged'
       goto 555
      endif
      if(nstep.gt.maxit)then
       write(6,*)'too many steps'
       goto 555
      endif
      e0=etot
      gnorm0=gnorm1
      goto 222
555   write(6,605)etot
 605  format(' Error:',e12.4)
      call wmap(m,nn,nk,ck,pxy,xmin,xmax,ymin,ymax,temp,ndim,ipp1,ipp2)
      return
      end

      subroutine wmap(m,nn,nk,ck,pxy,xmin,xmax,ymin,ymax,temp,ndim,
     1ipp1,ipp2)
      implicit none
      integer*4 nn,i,j,k,m,nk,ndim,ipp1,ipp2
      real*8 ck(m),x,y,dx,dy,xmin,ymin,xmax,ymax,dG,p,gfx,Gmin,
     1temp,kT,pxy(ndim)
      kT=temp/503.228d0
      dx=(xmax-xmin)/dble(nn)
      open(9,file='EMAP.TXT')
      open(91,file='EP.PRN')
      if(ndim.eq.1)then
       Gmin=0.0d0
       do 421 k=1,m
421    Gmin=Gmin+ck(k)*gfx(nk,k,xmin,0.0d0,pxy,ndim,ipp1,ipp2)

       x=xmin-dx
       do 4 i=1,nn+1
       x=x+dx
       dG=0.0d0
       do 42 k=1,m
42     dG=dG+ck(k)*gfx(nk,k,x,0.0d0,pxy,ndim,ipp1,ipp2)
4      Gmin=min(dG,Gmin)

       x=xmin-dx
       do 5 i=1,nn+1
       x=x+dx
       dG=0.0d0
       do 6 k=1,m
6      dG=dG+ck(k)*gfx(nk,k,x,0.0d0,pxy,ndim,ipp1,ipp2)
c      Free Energy in kcal/mol:
       dG=(dG-Gmin)*kT
       p=exp(-dG)
       write(91,941)x,dG
941    format(2e12.4)
5      write(9,94)i,x,dG,p
94     format(i6,3e12.4)
      else
       dy=(ymax-ymin)/dble(nn)

       Gmin=0.0d0
       do 321 k=1,m
321    Gmin=Gmin+ck(k)*gfx(nk,k,xmin,ymin,pxy,ndim,ipp1,ipp2)

       x=xmin-dx
       do 3 i=1,nn+1
       x=x+dx
       y=ymin-dy
       do 3 j=1,nn+1 
       y=y+dy
       dG=0.0d0
       do 32 k=1,m
32     dG=dG+ck(k)*gfx(nk,k,x,y,pxy,ndim,ipp1,ipp2)
3      Gmin=min(dG,Gmin)

       x=xmin-dx
       do 1 i=1,nn+1
       x=x+dx
       y=ymin-dy
       do 1 j=1,nn+1 
       y=y+dy
       dG=0.0d0
       do 2 k=1,m
2      dG=dG+ck(k)*gfx(nk,k,x,y,pxy,ndim,ipp1,ipp2)
c      Free Energy in kcal/mol:
       dG=(dG-Gmin)*kT
       p=exp(-dG)
1      write(9,91)i,j,x,y,dG,p
91     format(2i6,4e12.4)
      endif
      close(9)

      return
      end

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

      subroutine sort(co,N,sn)
      implicit none
      real*8 co(*),ci,cmax,s
      integer*4 n,i,imax,j,it,is
      integer*4 ,allocatable::ind(:)
      character*80 sn(*)
      real*8,allocatable::c(:)
      allocate(c(N))
      s=0.0d0
      do 5 i=1,N
      s=s+co(i)
5     c(i)=co(i)
      allocate(ind(N))
      do 10 i=1,N
  10  ind(i)=i
      do 1 i=1,N-1
      ci=c(i)
      imax=i
      cmax=dabs(ci)
      it=ind(i)
      do 2 j=i+1,N
      if(dabs(c(j)).gt.cmax)then
       imax=j
       cmax=dabs(c(j))
      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 cofficients absulute/normalized: '
      do 3 i=1,min(10,N)
      is=ind(i)+1
3     write(6,679)i,ind(i),c(i),c(i)/s,sn(is)(1:50)
679   format(2i4,E12.4,' /',f12.6,1x,a50) 
      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)
      return
      end

      subroutine readopt2(pxy,ndim,ran,lran)
      implicit none
      integer*4 ndim,i
      real*8 pxy(ndim),ran(ndim,2)
      character*10 s10
      logical lex,lran
c     pxy (all 360)         periodicity of each coordinate, in degrees
      pxy=360.0d0
      lran=.true.
      do 3 i=1,ndim
      ran(i,1)=  0.0d0
3     ran(i,2)=360.0d0
      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.'peri')read(9,*)pxy
       if(s10(1:4).eq.'lran')read(9,*)lran
       if(s10(1:4).eq.'rang')then
        do 2 i=1,ndim
2       read(9,*)ran(i,1),ran(i,2) 
       endif
       goto 1
99     close(9)
      endif
      return
      end