program specomp3
      implicit none
      integer ns,np,is,k,i,j,iter,istart,iend,kmax,ng,ix,iy,npsi,
     1delfi,delpsi,psimin,psimax,fimin,fimax,nbf,n,ifi,ipsi,l,nfi,
     1nsp
      real*8,allocatable::s(:),u(:),se(:),fi(:),si(:),
     1coe(:),psi(:),dG(:),ckl(:)
      real*8 wmin,wmax,dw,sum,gmin,temp,RT,glim,elim,w,wsum,dist,gam,
     1maxstep,gkl,mexp,xs0,xs1,xs2,gamma2,dx,dy,x,y,xper,yper,xd,yd,
     1xdd
      character*4 key
      character*80 s80
      logical lex,lwr,lnorm,lodd

      write(6,5999)
5999  format(/,
     1'                  Spectral decomposion',/,/,
     1'              ref=sum(pi si)',/,
     1'              pi=exp(-dG/RT)',/,
     1'              dG=su(c_kl g_kl)',/,
     1'              g_kl ... goniometric functions',/,
     1'              non-linear deviation minimization',/,/,
     1'        Input: SPECOMP3.OPT options',/,
     1'               ANGLE.LST options',/,
     1'               SP.LST    spectral list, first is reference',/,
     1'               pes.txt PES estimate (optional) ',/,
     1/,
     1'       Output: test.prn    normalized reference spectrum',/,
     1'               testf.prn   fitted spectrum',/)

      n=0
c     read angles
      open(9,file='ANGLE.LST')
1     read(9,*,end=99,err=99)
      n=n+1
      goto 1
99    rewind 9
      allocate(fi(n),psi(n))
      do 2 i=1,n
2     read(9,*)fi(i),fi(i),psi(i)
      close(9)
      write(6,*)n,' angle values'

c     parameters for smooth plotting
      delfi=5
      delpsi=5
      psimin=-185
      fimin=-185
      fimax=185
      psimax=185

c     other parameters
c     maximum frequency:
      kmax=3
c     number of points in working spectra (all are resampled to it):
      np=1000
c     workable range:
      wmin=200.0d0
      wmax=1800.0d0
c     maximum number of iterations:
      iter=20
c     temperature in K:
      temp=273.0d0
c     gradient and energy change convergence limits:
      glim=1.0d-3
      elim=1.0d-4
      maxstep=0.1d0
c     extensive output:
      lwr=.false.
c     number of initial guesses:
      ng=1
c     x-scaling x' = xs0 + xs1 * x + xs2 * x^2
      xs0=0.0d0
      xs1=1.0d0
      xs2=0.0d0
c     normalize the spectra:
      lnorm=.true.
c     uneven distribution in the fi psi plane:
      lodd=.false.
c     idealized numbers of points in the uneven case:
      nfi=20
      npsi=20
c     distance radius for spectral weights:
      gam=10.0d0
c     periodic conditions, periodicity:
      xper=360.0d0
      yper=360.0d0
      inquire(file='SPECOMP3.OPT',exist=lex)
      if(lex)then
       open(9,file='SPECOMP3.OPT')
993    read(9,900,end=95,err=95)key
900    format(a4)
       if(key.eq.'IFI ')read(9,*)delfi
       if(key.eq.'IPSI')read(9,*)delpsi
       if(key.eq.'PSII')read(9,*)psimin
       if(key.eq.'PSIA')read(9,*)psimax
       if(key.eq.'FII ')read(9,*)fimin
       if(key.eq.'FIA ')read(9,*)fimax
       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.eq.'KMAX')read(9,*)kmax
       if(key.eq.'NGUE')read(9,*)ng
       if(key.eq.'NORM')read(9,*)lnorm
       if(key.eq.'XPER')read(9,*)xper
       if(key.eq.'YPER')read(9,*)yper
       if(key(1:3).eq.'LWR')read(9,*)lwr
       if(key(1:3).eq.'XS0')read(9,*)xs0
       if(key(1:3).eq.'XS1')read(9,*)xs1
       if(key(1:3).eq.'XS2')read(9,*)xs2
       if(key(1:2).eq.'NP')read(9,*)np
       if(key(1:3).eq.'ODD')read(9,*)lodd
       if(key(1:3).eq.'NFI')read(9,*)nfi
       if(key(1:3).eq.'NPS')read(9,*)npsi
       if(key(1:3).eq.'GAM')read(9,*)gam
       goto 993
95     close(9)
       write(6,*)' SPECOMP3.OPT read'
      endif
      inquire(file='SP.LST',exist=lex)
      if(lex)then
       ns=0
       open(9,file='SP.LST')
992    read(9,*,end=918,err=918)
       ns=ns+1
       goto 992
918    close(9)
       ns=ns-1
       write(6,*)'SP.LST found, ',ns,' spectra'
      else
       write(6,*)'SP.LST not found'
       stop
      endif
      if(ns.ne.n)call report('ns <> n')
      allocate(dG(ns))

      nbf=(kmax+kmax+1)**2

      write(6,60)kmax,nbf,delfi,delpsi,fimin,fimax,
     1psimin,psimax
60    format(' kmax ',i4,/,
     1' number of basis functions ',i4,/,
     3' smooth pes:',/,
     4' delfi  delpsi ',2i10,/,
     4' fimin  fimax  ',2i10,/,
     4' psimin psimax ',2i10,/)


      dw=(wmax-wmin)/dble(np-1)
      RT=temp*8.314d0/4.184d0/1000.0d0
      gamma2=gam**2
      write(6,5998)wmin,wmax,dw,np,temp,elim,glim,maxstep,lwr,
     1iter,xs0,xs1,xs2,lodd,nfi,npsi,gam
5998  format(' ',/,
     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,/,
     1       '                      x-scaling 0 1 2:',3f6.3,/,
     1       '             uneven plotting Nfi Npsi:',l6,2i6,/,
     1       '                                gamma:',f12.2,/)

      allocate(s(ns*np),u(np),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
      if(lnorm)call qnorm(np,se)
      open(59,file='test.prn')
      do 991 i=1,np
991   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 902 is=1,ns
      read(39,9980)s80
      call trim(s80,istart,iend)
      call readsx(xs0,xs1,xs2,wmin,wmax,dw,u,np,s80(istart:iend))
      do 902 k=1,np
902   s(k+(is-1)*np)=u(k)
      close(39)
      write(6,*)' sub-spectra read'

      if(lodd)then
c      for uneven distribution create idealized set
       nsp=nfi*npsi
       allocate(si(np*nsp))
       write(6,*)nsp,' idealized spectra allocated ..'
       dx=(fimax-fimin)/dble(nfi-1)
       dy=(psimax-psimin)/dble(nfi-1)
       x=fimin-dx
       i=0
       do 13 ix=1,nfi
       x=x+dx
       y=psimin-dy
       do 13 iy=1,npsi
       y=y+dy
       i=i+1
c      idealized spectrum at this point  si = sum (w_is s_is) / sum (w_is)
       do 14 k=1,np
14     si(k+(i-1)*np)=0.0d0
       wsum=0.0d0
       do 15 is=1,ns
       xd=xdd(x, fi(is),xper)
       yd=xdd(y,psi(is),yper)
       dist=xd**2+yd**2
c      effective weight of spectrum is at xy:
       w=1.0d0/(1.0d0+dist/gamma2)
       wsum=wsum+w
       do 15 k=1,np
15     si(k+(i-1)*np)=si(k+(i-1)*np)+w*s(k+(is-1)*np)
       do 13 k=1,np
13     si(k+(i-1)*np)=si(k+(i-1)*np)/wsum
       deallocate(s)
       ns=nsp
       n=ns
       
       allocate(s(ns*np))
       do 17 k=1,np*ns
17     s(k)=si(k)

       deallocate(si)
       deallocate(fi)
       allocate(fi(n))

       write(6,*)' ... and calculated'
c      I do not know why the following two lines do not work:
       deallocate(psi)
       allocate(psi(n))

       x=fimin-dx
       i=0
       do 18 ix=1,nfi
       x=x+dx
       y=psimin-dy
       do 18 iy=1,npsi
       y=y+dy
       i=i+1
       fi(i)=x
18     psi(i)=y
       write(6,*)' ... and used as reference'

       deallocate(coe,dG)
       allocate(coe(ns),dG(ns))
      endif
c
c     normalize all 
      if(lnorm)then
       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'
      endif


c     Gibbs energy = sum(ckl*gkl), coefficients:
      allocate(ckl(nbf))

      call minimini(ns,s,se,np,iter,glim,elim,maxstep,
     1nbf,kmax,fi,psi,RT,ckl,ng)

      write(6,*)'Minimization finished'

      do 904 j=1,ns
      dG(j)=0.0d0

      do 905 k=0,2*kmax
      do 905 l=0,2*kmax
905   dG(j)=dG(j)+ckl(l+1+(2*kmax+1)*k)*gkl(k,l,fi(j),psi(j))

      if(j.eq.1)then
       gmin=dG(j)
      else
       if(dG(j).lt.gmin)gmin=dG(j)
      endif
904   coe(j)=mexp(-dG(j)/RT)

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     write PES and coefs to file
      open(9,file='test.txt')
      do 7 is=1,ns
7     write(9,6001)is,coe(is),dG(is)-gmin
6001  format(i5,2f12.6)
      close(9)
      write(6,*)' coefficients and free energies in test.txt'

c     fine PES:
      open(9,file='pes.txt')
      do 8 ifi=fimin,fimax,delfi
      do 8 ipsi=psimin,psimax,delpsi
      sum=0.0d0
      do 906 k=0,2*kmax
      do 906 l=0,2*kmax
906   sum=sum+ckl(l+1+(2*kmax+1)*k)*gkl(k,l,dble(ifi),
     1dble(ipsi))
8     write(9,909)ifi,ipsi,sum-gmin
909   format(2i10,f12.5)
      close(9)
      write(6,*)' fine smoothed function in pes.txt'
      stop

      end

      subroutine readsx(xs0,xs1,xs2,wmin,wmax,dw,u,np,fn)
c     read spectrum with immediate scaling
      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,xs0,xs1,xs2
      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
      x=xs0+x*(xs1+xs2*x)
      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 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(corm,s,ns,se,g,np,ic,gn,ckl,kmax,RT,
     1fi,psi)
      implicit none
      integer*4 ns,i,k,j,np,ic,kmax,l
      real*8 dd,s(*),se(*),g(*),sum,sump,mexp,
     1sd,gn,fi(*),fj,psi(*),ckl(*),gkl,RT,f1,corm(*)
      real*8,allocatable:: c(:),m(:,:),dG(:),es(:)
      allocate(c(ns),m(ns,ns),dG(ns),es(ns))
c
c     Gibbs energies from goniometric functions:
      do 1 i=1,ns
      dG(i)=0.0d0
      do 1 k=0,2*kmax
      do 1 l=0,2*kmax
      j=l+1+(2*kmax+1)*k
1     dG(i)=dG(i)+ckl(j)*gkl(k,l,fi(i),psi(i))
       
c     precalculate probabilities to run faster:
      do 2 i=1,ns
2     es(i)=mexp(-dG(i)/RT)

c     correlation matrix:
      do 3 i=1,ns
      f1=es(i)
      do 3 j=i,ns
      fj=es(j)
      m(j,i)=f1*fj*corm(i+(j-1)*ns)
3     m(i,j)=f1*fj*corm(i+(j-1)*ns)
c     scalar product:
      do 6 i=1,ns
      c(i)=0.0d0
      f1=es(i)
      do 6 k=1,np
6     c(i)=c(i)+se(k)*s(k+(i-1)*np)*f1

      sum=0.0d0
      do 10 i=1,np
      sd=0.0d0
      do 51 j=1,ns
51    sd=sd+es(j)*s(i+(j-1)*np)
10    sum=sum+(se(i)-sd)**2
      dd=sum

      if(ic.eq.0)goto 99

c     gradient g, its norm gn:
      gn=0.0d0
      do 7 k=0,2*kmax
      do 7 l=0,2*kmax
      sump=0.0d0
      do 9 i=1,ns
      do 8 j=1,ns
8     sump=sump+m(i,j)*gkl(k,l,fi(j),psi(j))
9     sump=sump-c(i)*gkl(k,l,fi(i),psi(i))
      g(l+1+(2*kmax+1)*k)=-2.0d0*sump/RT
7     gn=gn+g(l+1+(2*kmax+1)*k)**2
      gn=dsqrt(gn)
99    deallocate(c,m,dG)
      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,s,se,np,nstepmax,
     1glim,elim,maxstep,nbf,kmax,fi,psi,RT,ckl,ng)
      implicit none
      integer*4 nstep,lstep,ns,i,np,nstepmax,nbf,kmax,ig,ng,k,l
      real*8 e0,gnorm0,glim,etot,e00,dh,maxstep,eold,h,hh0,elim,e1,
     1gn,s(*),se(*),dd,secder,RT,ckl(*),an,eog,fir,psir,gi,delfir,pi,
     2gradh,enew,hh1,hh2,e2,gnorm1,fi(*),psi(*),delpsir,gkl,fio,psio
      real*8,allocatable::g0(:),h0(:),x0(:),g(:),cklold(:),corm(:)
      integer*4 ix,iy,iz
      real*8 frandom
      common /randc/ ix, iy, iz
      logical lex
      pi=4.0d0*datan(1.0d0)
      ix=1
      iy=1000
      iz=99
      eog=0.0d0
      fio=0.0d0
      psio=0.0d0
      if(nbf.gt.ns)call report('nbf>ns is not realistic')
c     gradient variables:
      allocate (g0(nbf),h0(nbf),x0(nbf),g(nbf),cklold(nbf),corm(ns*ns))
c     try ng initial guesses:
      do 1000 ig=1,ng
      write(6,*)' guess ',ig
      do 201 i=1,nbf
201   g(i)=0.0d0
      gn=0.0d0
      gnorm1=0.0d0
      gnorm0=0.0d0
c     initial guess:
      if(ig.eq.1)then
       inquire(file='pes.txt',exist=lex)
       if(lex)then
        write(6,*)'first guess from pes.txt'
        do 364 i=1,nbf
364     ckl(i)=0.0d0
        i=0
        open(55,file='pes.txt')
551     read(55,*,err=552,end=552)fir,psir,gi
        do 365 k=0,2*kmax
        do 365 l=0,2*kmax
365     ckl(l+1+(2*kmax+1)*k)=ckl(l+1+(2*kmax+1)*k)+gi*gkl(k,l,fir,psir)
        i=i+1
        if(i.gt.1.and.fir .ne.fio )delfir =fir-fio
        if(i.gt.1.and.psir.ne.psio)delpsir=psir-psio
        fio=fir
        psio=psir
        goto 551
552     close(55)
        write(6,*)i,' points'
        if(delfir*delpsir.eq.0.0d0)then
         delfir =20.0d0
         delpsir=20.0d0
        endif
        do 366 i=1,nbf
366     ckl(i)=ckl(i)*dabs(delfir*delpsir*pi*pi/180.0d0/180.0d0)
       else
        do 36 i=1,nbf
36      ckl(i)=1.0d0/dsqrt(dble(nbf))
       endif
      else
       do 361 i=1,nbf
361    ckl(i)=frandom(0)
      endif

c     precalculate correlation matrix:
      do 1 k=1,ns
      do 1 l=1,ns
      an=0.0d0
      do 2 i=1,np
2     an=an+s(i+(k-1)*np)*s(i+(l-1)*np)
1     corm(k+(l-1)*ns)=an
      nstep=0

      etot=dd(corm,s,ns,se,g,np,1,gnorm0,ckl,kmax,RT,fi,psi)
      write(6,6009)nstep,etot,gnorm0
6009  format(' Iteration',i5,' sig ',g12.4,' grad ',g12.4)
      e0=etot
      call asa(g0,g,nbf)
      call asa(h0,g,nbf)
      call asa(x0,ckl,nbf)
      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,nbf
101   ckl(i)=x0(i)-h*h0(i)
      etot=dd(corm,s,ns,se,g,np,0,gn,ckl,kmax,RT,
     1fi,psi)
      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(corm,s,ns,se,g,np,1,gnorm1,ckl,kmax,RT,
     1fi,psi)
      call asa(x0,ckl,nbf)
      do 102 i=1,nbf
102   h0(i)=g(i)+(gnorm1/gnorm0)**2*h0(i)
      call asa(g0,g,nbf)
      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   continue

      if(ig.gt.1)write(6,600)ig,etot,eog
600   format('guess ',i2,' actual energy : ',g10.2,' old: ',g10.2)

      if(ig.eq.1)then
       do 103 i=1,nbf
103    cklold(i)=ckl(i)
       eog=etot
      else
       if(etot.gt.eog)then
        write(6,*)' a previous guess better'
        etot=eog
        do 104 i=1,nbf
104     ckl(i)=cklold(i)
       else
        write(6,*)' guess kept'
        do 105 i=1,nbf
105     cklold(i)=ckl(i)
        eog=etot
       endif
      endif

1000  continue

      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

      function gkl(k,l,fi,psi)
      implicit none
      integer*4 js,jw,k,l,iw,is
      real*8 gkl,gk,gl,fi,psi,wj,wi,pi,ak,al
      pi=4.0d0*datan(1.0d0)
      jw=(k+1)/2
      wj=dble(jw)
      if(mod(k,2).eq.0)then
       js=0
      else
       js=1
      endif
      if(js.eq.1)then
       gk=sin(wj*fi*pi/180.0d0)
      else
       gk=cos(wj*fi*pi/180.0d0)
      endif
      if(k.eq.0)then
       ak=1.0d0/dsqrt(2.d0*pi)
      else
       ak=1.0d0/dsqrt(pi)
      endif

      iw=(l+1)/2
      wi=dble(iw)
      if(mod(l,2).eq.0)then
       is=0
      else
       is=1
      endif
      if(is.eq.1)then
       gl=sin(wi*psi*pi/180.0d0)
      else
       gl=cos(wi*psi*pi/180.0d0)
      endif
      if(l.eq.0)then
       al=1.0d0/dsqrt(2.d0*pi)
      else
       al=1.0d0/dsqrt(pi)
      endif

      gkl=ak*al*gk*gl

      return
      end

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

      function mexp(x)
      implicit none
      real*8 x,mexp,f
      if(x.gt.20.0d0)then
       f=4.8d8
      else
       if(x.lt.-20.0d0)then
        f=2.0d-9
       else
        f=dexp(x)
       endif
      endif
      mexp=f
      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

      function xdd(x1,x2,per)
      implicit none
      real*8 d,p2,xdd,x1,x2,per
      p2=per/2.0d0
      d=x2-x1
1     if(d.gt.p2)then
       d=d-per
       goto 1
      endif
2     if(d.lt.-p2)then
       d=d+per
       goto 2
      endif
      xdd=d
      return
      end