c     Program for the 'Propagation of short range order' in Chapter 10
c     Radial disribution function up to the seventh neighbor distance will
be found.
c     Set of 21 nonlinear equations is solved by iteration method with some
c     approximate initial values, included in the program, of 21 unknown
parameters.
      implicit real*8(a-h,o-z)
      dimension g(21),r(6),x(21)
      dimension mset(4),b(21,6),q(21,3),wx(21,21),wcut(21)
      logical cnv
      external fcts
      common/parm/iref
      common/rh/rho
      common/core/h
      open(unit=17,file="tempout.dat",status="unknown")
c      open(unit=2,file="input.data",status="unknown")
      db=1.d-5
      db=1.d-2
      gmin=1.d-10
      id=21
      mset(1)=id
      mset(2)=id
      mset(3)=200
      mset(4)=1
      rho=1.d0/6.d0
      ddh=dsqrt(10.d0)
      ddh=dsqrt(ddh)
      data x/
     -0.50991390300635D-01,0.22901297737811D-01,0.40449521427772D-01,
     -0.23407688100452D-01,0.37947172218659D-01,0.21957563197631D-01,
     -0.51461491708616D-02,0.58426965743075D-02,0.14982968752074D-01,
     -0.63060799284091D-02,0.56268683306135D-02,0.25421783020549D-02,
     -0.24307262742896D-02,0.47994411214763D-02,0.44772424310040D-02,
     -0.23235067934634D-02,0.46116724018749D-02,0.25474065730193D-02,
     -0.42993642937286D-02,0.22203879597745D-02,0.87238590078792D-03/
c       close(unit=2,status="keep")
c     This program is capable of finding accurate values of the correlation
functions,
c     even though the initial values are approximate.
      x(1)=0.05d0
      x(2)=0.02d0
      x(3)=0.04d0
      x(4)=0.02d0
      x(5)=0.03d0
      x(6)=0.022d0
      x(7)=0.051d0
      x(8)=0.d0
      x(9)=0.d0
      x(10)=0.d0
      x(11)=0.d0
      x(12)=0.d0
      x(13)=0.d0
      x(14)=0.d0
      x(15)=0.d0
      x(16)=0.d0
      x(17)=0.d0
      x(18)=0.d0
      x(19)=0.0d0
      x(20)=0.d0
c      x(21)=0.d0
      do 10 i=1,id
          b(i,2)=x(i)
          b(i,3)=b(i,2)*db+1.d-14
          b(i,4)=0.d0
         b(i,5)=1.d0
         b(i,6)=1.d0
         q(i,2)=-gmin
         q(i,3)=gmin
         wcut(i)=1.d0
         wx(i,i)=1.d-06
      do 10 j=1,id
         if(i.ne.j) wx(i,j)=0.d0
   10 continue
c
      call minmax(mset,cnv,b,fcts,q,wx,id,wcut,ll)
      do 40 i=1,id
           x(i)=b(i,1)
           g(i)=q(i,1)
           b(i,2)=b(i,1)
           b(i,3)=b(i,1)*db+1.d-14
   40 continue
c      ss=entr(x)
      write(17,100) iref,cnv
      write(17,110) ll
      write(17,200)  x
      write(17,400)g
      do 60 i=1,6
      r(i)=x(i)*36.d0
 60   continue
c      call prt(x)
      write(17,300) r
      if(ll.eq.mset(3)) stop
 110  format(1h ,1x,'ite =',i6)
 200  format(7(1h,4x,'x =  ',1pd20.14,3x,d20.14,3x,d20.14/))
 100  format(1h ,'iref=',i10 ,2x,' converg = ',l2)
 300  format(1x,'rad:: ',6f10.5)
 400  format(3(1x,'g :',1p7d10.1/))
 600  format(5x,d20.10)
      close(unit=17,status="keep")
      stop
      end
c
      subroutine fcts(mset,g,x)
      implicit real*8(a-h,o-z)
      dimension g(21),c2(6),x(21),g5(4),d(23,20),c3(23)
      dimension mset(4)
      common /parm/iref
      common /rh/rho
      common/core/h
c
      data c2/1.d0,4.d0,2.d0,2*4.d0,8.d0/
      data c3/6.d0,2*12.d0,5*24.d0,48.d0,12.d0,2*24.d0,6.d0,12.d0,
     -4.d0,24.d0,24.d0,48.d0,36.d0,2*24.d0,48.d0,8.d0/
c
      iref=0
      do 5 i=1,21
      g(i)=0.d0
   5  continue
c
      do 20 i=1,23
      do 20 j=1,20
      d(i,j)=0.d0
  20  continue
c
c      d(1,1)=dtr1(x(1))
c      d(2,2)=dtr1(x(2))
      d(3,1)=dtr2(x(1),x(4),1)
      d(3,4)=dtr2(x(1),x(4),2)
      d(4,2)=dtr2(x(2),x(3),1)
      d(4,3)=dtr2(x(2),x(3),2)
      d(5,2)=dtr2(x(2),x(4),1)
      d(5,4)=dtr2(x(2),x(4),2)
      d(6,2)=dtr2(x(2),x(5),1)
      d(6,5)=dtr2(x(2),x(5),2)
      d(7,3)=dtr2(x(3),x(6),1)
      d(7,6)=dtr2(x(3),x(6),2)
      d(8,4)=dtr2(x(4),x(6),1)
      d(8,6)=dtr2(x(4),x(6),2)
      d(9,5)=dtr2(x(5),x(6),1)
      d(9,6)=dtr2(x(5),x(6),2)
c      d(10,1)=dtr3(x(1),x(2),x(2),x(7),1)
c      d(10,2)=dtr3(x(1),x(2),x(2),x(7),2)*2.d0
c      d(10,7)=dtr3(x(1),x(2),x(2),x(7),4)
      call cubic(x(1),x(2),x(6),x(8))
      d(11,1)=dtr3(x(1),x(2),x(6),x(8),1)
      d(11,2)=dtr3(x(1),x(2),x(6),x(8),2)
      d(11,6)=dtr3(x(1),x(2),x(6),x(8),3)
      d(11,8)=dtr3(x(1),x(2),x(6),x(8),4)
      call cubic(x(1),x(3),x(5),x(9))
      d(12,1)=dtr3(x(1),x(3),x(5),x(9),1)
      d(12,3)=dtr3(x(1),x(3),x(5),x(9),2)
      d(12,5)=dtr3(x(1),x(3),x(5),x(9),3)
      d(12,9)=dtr3(x(1),x(3),x(5),x(9),4)
c
      call cubic(x(1),x(4),x(4),x(10))
      d(13,1)=dtr3(x(1),x(4),x(4),x(10),1)
      d(13,4)=dtr3(x(1),x(4),x(4),x(10),2)*2.d0
      d(13,10)=dtr3(x(1),x(4),x(4),x(10),4)
c
      call cubic(x(1),x(6),x(6),x(11))
      d(14,1)=dtr3(x(1),x(6),x(6),x(11),1)
      d(14,6)=dtr3(x(1),x(6),x(6),x(11),2)*2.d0
      d(14,11)=dtr3(x(1),x(6),x(6),x(11),4)
c
      call cubic(x(2),x(2),x(2),x(12))
      d(15,2)=dtr3(x(2),x(2),x(2),x(12),1)*3.d0
      d(15,12)=dtr3(x(2),x(2),x(2),x(12),4)
c
      call cubic(x(2),x(2),x(6),x(13))
      d(16,2)=dtr3(x(2),x(2),x(6),x(13),1)*2.d0
      d(16,6)=dtr3(x(2),x(2),x(6),x(13),3)
      d(16,13)=dtr3(x(2),x(2),x(6),x(13),4)
c
      call cubic(x(2),x(3),x(4),x(14))
      d(17,2)=dtr3(x(2),x(3),x(4),x(14),1)
      d(17,3)=dtr3(x(2),x(3),x(4),x(14),2)
      d(17,4)=dtr3(x(2),x(3),x(4),x(14),3)
      d(17,14)=dtr3(x(2),x(3),x(4),x(14),4)
c
      call cubic(x(2),x(4),x(5),x(15))
      d(18,2)=dtr3(x(2),x(4),x(5),x(15),1)
      d(18,4)=dtr3(x(2),x(4),x(5),x(15),2)
      d(18,5)=dtr3(x(2),x(4),x(5),x(15),3)
      d(18,15)=dtr3(x(2),x(4),x(5),x(15),4)
c
      call cubic(x(2),x(6),x(6),x(16))
      d(19,2)=dtr3(x(2),x(6),x(6),x(16),1)
      d(19,6)=dtr3(x(2),x(6),x(6),x(16),2)*2.d0
      d(19,16)=dtr3(x(2),x(6),x(6),x(16),4)
c
      call cubic(x(3),x(4),x(6),x(17))
      d(20,3)=dtr3(x(3),x(4),x(6),x(17),1)
      d(20,4)=dtr3(x(3),x(4),x(6),x(17),2)
      d(20,6)=dtr3(x(3),x(4),x(6),x(17),3)
      d(20,17)=dtr3(x(3),x(4),x(6),x(17),4)
c      write(17,130) iref
      call cubic(x(4),x(4),x(6),x(18))
      d(21,4)=dtr3(x(4),x(4),x(6),x(18),1)*2.d0
      d(21,6)=dtr3(x(4),x(4),x(6),x(18),3)
      d(21,18)=dtr3(x(4),x(4),x(6),x(18),4)
c
      call cubic(x(4),x(5),x(6),x(19))
      d(22,4)=dtr3(x(4),x(5),x(6),x(19),1)
      d(22,5)=dtr3(x(4),x(5),x(6),x(19),2)
      d(22,6)=dtr3(x(4),x(5),x(6),x(19),3)
      d(22,19)=dtr3(x(4),x(5),x(6),x(19),4)
c
      call cubic(x(6),x(6),x(6),x(20))
      d(23,6)=dtr3(x(6),x(6),x(6),x(20),1)*3.d0
      d(23,20)=dtr3(x(6),x(6),x(6),x(20),4)
c
      call dfiv(x(1),x(2),x(7),x(21),g5)
c
      do 50 j=1,6
      do 50 i=3,23
      g(j)=g(j)+c3(i)*d(i,j)/3.d0
   50 continue
      do 70 i=3,6
      g(i)=g(i)+c2(i)*dpa(x(i))
70    continue
c
      do 60 i=7,20
      g(i)=g(i)+d(i+3,i)*c3(i+3)/3.d0
   60 continue
c
c    Free energy :   g4(1234)+4*g4(0123)+2*g3(001)+4*g3(002)+
c                          +4*g3(122)+g2(1)+4*g2(2)
c                          +[G5(01234)-g4(1234)-4*g4(0123)-2*g3(001)-4*g3(002)
c                              -4*g3(122)-2*g2(1)-4*g2(2)]
c    g4(1234),g4(0123),g3(001),g3(002),g3(122) do not contribute
c

      g(1)=g(1)+g5(1)-dpa(x(1))
      g(2)=g(2)+g5(2)
      g(7)=g5(3)
      g(21)=g5(4)
c      call dfour1(x(1),x(4),x(6),x(11),g41)
c      g(1)=g(1)+g41(1)-d(14,1)-dtr1(x(1))-dpa(x(1))
c      g(4)=g(4)+g41(2)-2.d0*d(8,4)-dpa(x(4))
c      g(6)=g(6)+g41(3)-2.d0*d(8,6)-d(14,6)-2.d0*dpa(x(6))
c      g(11)=g(11)+g41(4)-d(14,11)
      return
      end
c
      subroutine dfiv(x1,x2,x7,x21,g)
      implicit real*8(a-h,o-z)
      dimension d(6),r(6),w(6),g(4)
      common /rh/rho
      data w/1.d0,4.d0,2.d0,4.d0,4.d0,1.d0/
      r(1)=x21
      r(2)=x7-x21
      r(3)=x1-2.d0*x7+x21
      r(4)=x2-2.d0*x7+x21
      r(5)=rho-x1-2.d0*x2+3.d0*x7-x21
      r(6)=1.d0-5.d0*rho+2.d0*x1+4.d0*x2-4.d0*x7+x21
      do 10 i=1,6
           d(i)=w(i)*dlg(r(i))
   10 continue
      g(1)=d(3)-d(5)+2.d0*d(6)
      g(2)=d(4)-2.d0*d(5)+4.d0*d(6)
      g(3)=d(2)-2.d0*(d(3)+d(4))+3.d0*d(5)-4.d0*d(6)
      g(4)=d(1)-d(2)+d(3)+d(4)-d(5)+d(6)
      return
      end
c
      subroutine cubic(a1,a2,a3,tri)
      implicit real*8(a-h,o-z)
      common/rh/rho
      a=rho-a1-a2
      b=rho-a2-a3
      c=rho-a3-a1
      a4=1.d0-3.d0*rho+a1+a2+a3
      b2=a*b+b*c+c*a-(a1*a2+a1*a3+a1*a4+a2*a3+a2*a4+a3*a4)
      b1=a*b*c+a1*a2*a3+a1*a2*a4+a1*a3*a4+a2*a3*a4
      b0=-a1*a2*a3*a4
      q=b1/3.d0-b2**2/9.d0
      r=(b1*b2-3.d0*b0)/6.d0-b2**3/27.d0
      s1=(r+dsqrt(q*q*q+r*r))**(1.d0/3.d0)
      s2=-(-r+dsqrt(q*q*q+r*r))**(1.d0/3.d0)
      tri=s1+s2-b2/3.d0
      return
      end
c
c      real*8 Function dtr1(a)
c       implicit real*8(a-h,o-z)
c       common/rh/rho
c      dtr1=dlg((1.d0-3.d0*rho+a)/(1.d0-2.d0*rho+a))
c      return
c      end
c
      real*8 Function dtr2(a1,a2,n)
       implicit real*8(a-h,o-z)
       common/rh/rho
      if(n.eq.1) then
      dtr2=dlg((rho-a1)*(1.d0-3.d0*rho+a1+a2)/((rho-a1-a2)
     - *(1.d0-2.d0*rho+a1)))
      return
      else
      dtr2=dlg((rho-a2)*(1.d0-3.d0*rho+a1+a2)/((rho-a1-a2)
     - *(1.d0-2.d0*rho+a2)))
      return
      endif
      end
c
      real*8 Function dtr3(a1,a2,a3,t,n)
      implicit real*8(a-h,o-z)
      dimension r(8)
      common/rh/rho
       r(1)=t
       r(2)=a2-t
       r(3)=a3-t
       r(4)=a1-t
       r(5)=rho-a2-a3+t
       r(6)=rho-a1-a2+t
       r(7)=rho-a1-a3+t
       r(8)=1.d0-3.d0*rho+a1+a2+a3-t
       if(n.eq.1) dtr3=dlg(r(4)*r(8)/(r(6)*r(7)))-dpa(a1)
       if(n.eq.2) dtr3=dlg(r(2)*r(8)/(r(5)*r(6)))-dpa(a2)
       if(n.eq.3) dtr3=dlg(r(3)*r(8)/(r(5)*r(7)))-dpa(a3)
c       if(n.eq.4) dtr3=d(1)-d(2)-d(3)-d(4)+d(5)+d(6)+d(7)-d(8)
       if(n.eq.4) dtr3=dlg(r(1)*r(5)*r(6)*r(7)/(r(2)*r(3)*r(4)*r(8)))
       return
       end
c
      real*8 function dpa(z)
      implicit real*8(a-h,o-z)
      common /rh/rho
      dpa=dlg(z*(1.d0-2.d0*rho+z)/((rho-z)*(rho-z)))
      return
      end
c
      real*8 function dlg(x)
      implicit real*8 (a-h,o-z)
      common/parm/iref
c
      if(x.lt.1.d-14) go to 2
      dlg=dlog(x)
      return
c
    2 dlg=-14.d0
      iref=iref+1
      return
      end
c     Effects of pairs up to the 13th neighbor and all triangles containing
c     up to the 13th neighbor distance can be included with this program.
c     x(1) is the second neighbor and x(12) is the 13th neighbor because
c     the nearest neighbor pair correlation function is set equal to zero.
      real*8 function entr(x)
      implicit real*8(a-h,o-z)
      dimension x(97),c3(105),c2(13),g2(13),g3(105)
      common/rh/rho
      data c2/2.d0,1.d0,4.d0,2.d0,2*4.d0,8.d0,3.d0,6.d0,3*4.d0,12.d0/
      data c3/6.d0,2*12.d0,6*24.d0,12.d0,
     -48.d0,24.d0,48.d0,4*24.d0,36.d0,36.d0,24.d0,24.d0,
     -3.d0,12.d0,2*24.d0,6.d0,2*12.d0,24.d0,12.d0,2*24.d0,12.d0,
     -24.d0,18.d0,4.d0,24.d0,3*12.d0,24.d0,2*24.d0,48.d0,
     -2*24.d0,3*48.d0,36.d0,48.d0,48.d0,5*24.d0,
     -48.d0,36.d0,12.d0,24.d0,24.d0,3*48.d0,24.d0,24.d0,
     -24.d0,2*12.d0,48.d0,24.d0,48.d0,2*24.d0,48.d0,2*24.d0,
     -36.d0,24.d0,36.d0,48.d0,96.d0,24.d0,2*48.d0,8.d0,24.d0,36.d0,
     -48.d0,24.d0,72.d0,24.d0,72.d0,72.d0,48.d0,2*12.d0,24.d0,12.d0,
     -60.d0,24.d0,24.d0,48.d0,24.d0/
c
      g1=rho*dlog(rho)+(1.d0-rho)*dlog(1.d0-rho)
      g2(1)=2.d0*rho*dlg(rho)+(1.d0-2.d0*rho)*dlg(1.d0-2.d0*rho)-2.d0*g1
      do 10 i=2,13
      g2(i)=pa(x(i-1))-2.d0*g1
 10    continue
c
       g3(1)=tr1(x(1))
       g3(2)=tr1(x(2))
c
       g3(3)=tr2(x(1),x(4))
       g3(4)=tr2(x(2),x(3))
       g3(5)=tr2(x(2),x(4))
       g3(6)=tr2(x(2),x(5))
       g3(7)=tr2(x(3),x(6))
       g3(8)=tr2(x(4),x(6))
       g3(9)=tr2(x(4),x(7))
       g3(10)=tr2(x(4),x(9))
       g3(11)=tr2(x(5),x(6))
       g3(12)=tr2(x(5),x(10))
       g3(13)=tr2(x(6),x(8))
       g3(14)=tr2(x(6),x(11))
       g3(15)=tr2(x(6),x(12))
       g3(16)=tr2(x(7),x(8))
       g3(17)=tr2(x(7),x(12))
       g3(18)=tr2(x(8),x(9))
       g3(19)=tr2(x(9),x(12))
       g3(20)=tr2(x(10),x(11))
       g3(21)=tr2(x(10),x(12))
       g3(22)=tr3(x(1),x(1),x(7),x(13))
       g3(23)=tr3(x(1),x(2),x(2),x(14))
       g3(24)=tr3(x(1),x(2),x(6),x(15))
       g3(25)=tr3(x(1),x(3),x(5),x(16))
       g3(26)=tr3(x(1),x(4),x(4),x(17))
       g3(27)=tr3(x(1),x(4),x(8),x(18))
       g3(28)=tr3(x(1),x(4),x(12),x(19))
       g3(29)=tr3(x(1),x(5),x(11),x(20))
       g3(30)=tr3(x(1),x(6),x(6),x(21))
       g3(31)=tr3(x(1),x(6),x(10),x(22))
       g3(32)=tr3(x(1),x(7),x(9),x(23))
       g3(33)=tr3(x(1),x(8),x(8),x(24))
       g3(34)=tr3(x(1),x(8),x(12),x(25))
       g3(35)=tr3(x(1),x(12),x(12),x(26))
       g3(36)=tr3(x(2),x(2),x(2),x(27))
       g3(37)=tr3(x(2),x(2),x(6),x(28))
       g3(38)=tr3(x(2),x(2),x(7),x(29))
       g3(39)=tr3(x(2),x(2),x(9),x(30))
       g3(40)=tr3(x(2),x(2),x(10),x(31))
       g3(41)=tr3(x(2),x(3),x(4),x(32))
       g3(42)=tr3(x(2),x(3),x(8),x(33))
       g3(43)=tr3(x(2),x(3),x(12),x(34))
       g3(44)=tr3(x(2),x(4),x(5),x(35))
       g3(45)=tr3(x(2),x(4),x(8),x(36))
       g3(46)=tr3(x(2),x(4),x(11),x(37))
       g3(47)=tr3(x(2),x(4),x(12),x(38))
       g3(48)=tr3(x(2),x(5),x(8),x(39))
       g3(49)=tr3(x(2),x(5),x(12),x(40))
       g3(50)=tr3(x(2),x(6),x(6),x(41))
       g3(51)=tr3(x(2),x(6),x(7),x(42))
       g3(52)=tr3(x(2),x(6),x(9),x(43))
       g3(53)=tr3(x(2),x(6),x(10),x(44))
       g3(54)=tr3(x(2),x(8),x(8),x(45))
       g3(55)=tr3(x(2),x(8),x(11),x(46))
       g3(56)=tr3(x(2),x(8),x(12),x(47))
       g3(57)=tr3(x(2),x(9),x(10),x(48))
       g3(58)=tr3(x(2),x(11),x(12),x(49))
       g3(59)=tr3(x(2),x(12),x(12),x(50))
       g3(60)=tr3(x(3),x(3),x(7),x(51))
       g3(61)=tr3(x(3),x(4),x(6),x(52))
       g3(62)=tr3(x(3),x(4),x(10),x(53))
       g3(63)=tr3(x(3),x(5),x(9),x(54))
       g3(64)=tr3(x(3),x(6),x(8),x(55))
       g3(65)=tr3(x(3),x(6),x(12),x(56))
       g3(66)=tr3(x(3),x(7),x(11),x(57))
       g3(67)=tr3(x(3),x(10),x(12),x(58))
       g3(68)=tr3(x(4),x(4),x(6),x(59))
       g3(69)=tr3(x(4),x(4),x(9),x(60))
       g3(70)=tr3(x(4),x(4),x(10),x(61))
       g3(71)=tr3(x(4),x(5),x(6),x(62))
       g3(72)=tr3(x(4),x(6),x(8),x(63))
       g3(73)=tr3(x(4),x(6),x(11),x(64))
       g3(74)=tr3(x(4),x(6),x(12),x(65))
       g3(75)=tr3(x(4),x(7),x(8),x(66))
       g3(76)=tr3(x(4),x(7),x(12),x(67))
       g3(77)=tr3(x(4),x(8),x(9),x(68))
       g3(78)=tr3(x(4),x(8),x(10),x(69))
       g3(79)=tr3(x(4),x(9),x(12),x(70))
       g3(80)=tr3(x(4),x(10),x(12),x(71))
       g3(81)=tr3(x(5),x(5),x(7),x(72))
       g3(82)=tr3(x(5),x(6),x(8),x(73))
       g3(83)=tr3(x(5),x(6),x(12),x(74))
       g3(84)=tr3(x(5),x(8),x(10),x(75))
       g3(85)=tr3(x(5),x(9),x(11),x(76))
       g3(86)=tr3(x(5),x(10),x(12),x(77))
       g3(87)=tr3(x(6),x(6),x(6),x(78))
       g3(88)=tr3(x(6),x(6),x(7),x(79))
       g3(89)=tr3(x(6),x(6),x(9),x(80))
       g3(90)=tr3(x(6),x(7),x(10),x(81))
       g3(91)=tr3(x(6),x(8),x(11),x(82))
       g3(92)=tr3(x(6),x(8),x(12),x(83))
       g3(93)=tr3(x(6),x(9),x(10),x(84))
       g3(94)=tr3(x(6),x(11),x(12),x(85))
       g3(95)=tr3(x(6),x(12),x(12),x(86))
       g3(96)=tr3(x(7),x(8),x(12),x(87))
       g3(97)=tr3(x(7),x(9),x(9),x(88))
       g3(98)=tr3(x(7),x(10),x(10),x(89))
       g3(99)=tr3(x(7),x(11),x(11),x(90))
       g3(100)=tr3(x(8),x(8),x(10),x(91))
       g3(101)=tr3(x(8),x(9),x(12),x(92))
       g3(102)=tr3(x(8),x(10),x(11),x(93))
       g3(103)=tr3(x(8),x(10),x(12),x(94))
       g3(104)=tr3(x(9),x(12),x(12),x(95))
       g3(105)=tr3(x(10),x(11),x(12),x(96))
       return
       end
c
      subroutine fiv(x1,x2,x7,x21,ent5)
      implicit real*8(a-h,o-z)
      dimension r(7),w(7)
      common /rh/rho
      data w/1.d0,4.d0,2.d0,4.d0,4.d0,1.d0,1.d0/
      r(1)=x21
      r(2)=x7-x21
      r(3)=x1-2.d0*x7+x21
      r(4)=x2-2.d0*x7+x21
      r(5)=rho-x1-2.d0*x2+3.d0*x7-x21
      r(6)=rho
      r(7)=1.d0-5.d0*rho+2.d0*x1+4.d0*x2-4.d0*x7+x21
      ent5=0.d0
      do 10 i=1,7
      ent5=ent5+w(i)*r(i)*dlg(r(i))
10    continue
      return
      end
c
      real*8 function tr1(a)
      implicit real*8(a-h,o-z)
      dimension r(5)
      common /rh/rho
      r(1)=a
      r(2)=rho
      r(3)=rho-a
      r(4)=rho-a
      r(5)=0.5d0+a
      g=0.d0
      do 10 i=1,5
      g=g+r(i)*dlg(r(i))
10    continue
      g1=rho*dlg(rho)+(1.d0-rho)*dlg(1.d0-rho)
      g21=2.d0*rho*dlg(rho)+(1.d0-2.d0*rho)*dlg(1.d0-2.d0*rho)
      g22=a*dlg(a)+2.d0*(rho-a)*dlg(rho-a)
     -      +(1.d0-2.d0*rho+a)*dlg(1.d0-2.d0*rho+a)
      tr1=g-2.d0*g21-g22+3.d0*g1
      return
      end
c
      real*8 function tr2(a1,a2)
      implicit real*8(a-h,o-z)
      dimension r(6)
      common /rh/rho
      r(1)=a2
      r(2)=a1
      r(3)=rho-a2
      r(4)=rho-a1-a2
      r(5)=rho-a1
      r(6)=0.5d0+a1+a2
      g=0.d0
      do 10 i=1,6
      g=g+r(i)*dlg(r(i))
10    continue
      g1=rho*dlg(rho)+(1.d0-rho)*dlg(1.d0-rho)
      g21=2.d0*rho*dlg(rho)+(1.d0-2.d0*rho)*dlg(1.d0-2.d0*rho)
      g22=a1*dlg(a1)+2.d0*(rho-a1)*dlg(rho-a1)
     -      +(1.d0-2.d0*rho+a1)*dlg(1.d0-2.d0*rho+a1)
      g23=a2*dlg(a2)+2.d0*(rho-a2)*dlg(rho-a2)
     -      +(1.d0-2.d0*rho+a2)*dlg(1.d0-2.d0*rho+a2)
      tr2=g-2.d0*g21-g22-g23+3.d0*g1
      return
      end
c
      real*8 function tr3(a1,a2,a3,t)
      implicit real*8(a-h,o-z)
      dimension r(8)
      common /rh/rho
       r(1)=t
       r(2)=a2-t
       r(3)=a3-t
       r(4)=a1-t
      r(5)=rho-a2-a3+t
      r(6)=rho-a1-a2+t
      r(7)=rho-a1-a3+t
      r(8)=1.d0-3.d0*rho+a1+a2+a3-t
      g=0.d0
      do 10 i=1,8
      g=g+r(i)*dlg(r(i))
10    continue
      g1=rho*dlg(rho)+(1.d0-rho)*dlg(1.d0-rho)
      g21=a1*dlg(a1)+2.d0*(rho-a1)*dlg(rho-a1)
     -      +(1.d0-2.d0*rho+a1)*dlg(1.d0-2.d0*rho+a1)
      g22=a2*dlg(a2)+2.d0*(rho-a2)*dlg(rho-a2)
     -      +(1.d0-2.d0*rho+a2)*dlg(1.d0-2.d0*rho+a2)
      g23=a3*dlg(a3)+2.d0*(rho-a3)*dlg(rho-a3)
     -      +(1.d0-2.d0*rho+a3)*dlg(1.d0-2.d0*rho+a3)
      tr3=g-g21-g22-g23+3.d0*g1
      return
      end
c
      real*8 function pa(z)
      implicit real*8(a-h,o-z)
      common/rh/rho
      pa=z*dlg(z)+2.d0*(rho-z)*dlg(rho-z)
     -  +(1.d0-2.d0*rho+z)*dlg(1.d0-2.d0*rho+z)
      return
      end
c
c     subroutine minmax(mset,convrg,b,fcts,q,wx,nd,wxcut,l)
c     program is adapted to on june 1,1983
      subroutine minmax(mset,convrg,b,fcts,q,wx,nd,wcut,l,xlamda)
      implicit real*8(a-h,o-z)
      dimension b(nd,6),q(nd,3),wx(nd,nd),w(97,97),wcut(nd)
      logical feedbk, first, convrg,nostop
      dimension qfdbk(97),mfb(97),nofb(97),dxpm(10000)
      dimension c(97)
      dimension beta(97)
      dimension res(97),wres(97),qnom(97),mset(4),pm(97,97)
      dimension xpm(97,97),xr(97),db(97),xrt(97),xpmt(97,97)
      equivalence (dxpm(1),xpm(1))
      data nc4/10000000/
      data clamp/700.d0/
      twom28=1.d0/2.d0**48
      avrage=0.d0
      reduce=0.01d0
      xlamda=1.d-9
      kb=mset(1)
      kq=mset(2)
      kstop=mset(3)
      assign 505 to jgo2
      rlold=0.d0
      moptm = 0
      do 10 m=1,kq
      if( wx(m,m) ) 5, 6, 6
    5 moptm=m
      go to 10
    6 c(m)=wx(m,m)
      rlold=rlold+c(m)
   10 continue
      if(moptm) 11, 12, 11
   11 fkq=dfloat((kq-1)*128)
      avrage=rlold/fkq
      c(moptm)=avrage*0.001d0
      wx(moptm,moptm)=c(moptm)
   12 continue
      l=0
      kount=0
      first = .true.
      feedbk=.false.
   20 do 21 m=1,kb
   21 b(m,1)=b(m,2)
      nostop=.true.
      go to 86
    4 kount=kount+1
  701 rlnew=0.d0
      convrg=.true.
      do 23 m=1,kq
      do 23 k=1,kq
   23 w(m,k)=wx(m,k)
      do 27 m=1,kq
      res(m)=q(m,1)-0.5d0*(q(m,2)+q(m,3))
      if(q(m,1).ge.q(m,2).and.q(m,1).le.q(m,3)) go to 26
      if ( m - moptm ) 24, 27, 24
   24 convrg=.false.
      go to 27
   26 w(m,m)=w(m,m)*wcut(m)
   27 continue
      if(.not.convrg) go to 30
      nostop=.false.
      reduce=1.d0
      go to 86
30    do 29 m=1,kq
      wres(m)=0.d0
      do 28 k=1,kq
28    wres(m)=wres(m)+w(m,k)*res(k)
29    rlnew=rlnew+wres(m)*res(m)
      if(first) go to 35
      go to jgo2, (505, 510)
505   if(.not.convrg) go to 515
      if(moptm) 507, 525, 507
507   reduce=1.d0
      assign 510 to jgo2
      c(moptm)=avrage
      kstop=nc4
      do 509 m=1,kq
      nofb(m)=0
509   wx(m,m)=c(m)
      first=.true.
      go to 701
  510 if(rlold-rlnew) 57, 57, 240
  515 if(rlnew-rlold) 40, 57, 57
   35 first=.false.
      lsize=1
      go to 250
   40 if(feedbk) go to 235
  209 jfb = 1
      go to 411
235   if(lsize) 410, 666, 410
666   kfb=kfb+1
      if(kfb-6) 410, 400, 400
400   feedbk=.false.
      do 405 m=1,jfb
      j=mfb(m)
      k=nofb(m)
405   q(j,k)=qfdbk(m)
      first=.true.
      go to 701
410   jfb=jfb+1
411   do 450 m=1,kq
      if(wcut(m)) 450, 415, 450
415   if(q(m,1).lt.q(m,2).or.q(m,1).gt.q(m,3)) go to 450
      temp=q(m,2)-q(m,3)
      if(temp) 420, 450, 420
420   if(dabs((q(m,1)-q(m,2))/temp) - 0.002d0) 455, 425, 425
425   if(dabs((q(m,1)-q(m,3))/temp) - 0.002d0) 460, 450, 450
450
continue
      jfb=jfb-1
      go to 255
455   nofb(jfb)=2
      go to 461
460   nofb(jfb)=3
461   j=nofb(jfb)
      qfdbk(jfb)=q(m,j)
      mfb(jfb)=m
      k=5-j
      q(m,j)=q(m,k)
      feedbk=.true.
      l=l+1
      kfb=0
      first=.true.
      go to 701
240   do 933 m=1,kq
      if(wcut(m)) 933, 241, 933
241   temp=(q(m,2)-q(m,3))
      if(temp) 242, 901, 242
242   if(q(m,1)-q(m,2)) 933, 245, 243
243   if(q(m,1)-q(m,3)) 244, 247, 933
244   if(dabs((q(m,1)-q(m,2))/temp)-0.002d0) 245, 246, 246
245   nofb(m)=1
931   qfdbk(m)=q(m,3)
      q(m,3)=q(m,2)
      first=.true.
      go to 933
246   if(dabs((q(m,1)-q(m,3))/temp)-0.002d0) 247, 933, 933
247   nofb(m)=1
932   qfdbk(m)=q(m,2)
      q(m,2)=q(m,3)
      first=.true.
      go to 933
901   if(m-moptm) 910, 933, 910
902   if(qfdbk(m)-q(m,2)) 903, 933, 905
903   if(q(m,1)-q(m,2)) 904, 933, 933
904   q(m,2)=qfdbk(m)
      first=.true.
      go to 933
905   if(q(m,1)-q(m,2)) 933, 933, 906
906   q(m,3)=qfdbk(m)
      first=.true.
      go to 933
910   if(nofb(m)) 911, 902, 911
911   if(qfdbk(m)-q(m,2)) 912, 933, 914
912   if(q(m,1)-q(m,2)) 913, 930, 930
913   nofb(m)=nofb(m)+1
      if(nofb(m).lt.5) go to 933
      q(m,2)=qfdbk(m)
      nofb(m)=0
      go to 931
914   if(q(m,1)-q(m,2)) 930, 930, 915
915   nofb(m)=nofb(m)+1
      if(nofb(m).lt.5) go to 933
      q(m,3)=qfdbk(m)
      nofb(m)=0
      go to 932
930   nofb(m)=1
933   continue
      if(first) go to 701
250   continue
255   rlold=rlnew
      if(lsize) 610, 600, 620
600   lsize=-1
      if(xlamda-1.d-15) 620, 602, 602
602   xlamda=xlamda*5.d-1
603   do 604 m=1,kb
604   b(m,1)=beta(m)
      go to 70
610   lsize=1
      go to 602
620   lsize=0
      if(.not.feedbk) l=1+l
      if(l-kstop) 260, 110, 110
260   if(xlamda.ge.1.d-15) xlamda=xlamda*reduce
      do 41 m=1,kb
41    beta(m)=b(m,1)
      do 50 m=1,kb
      b(m,1)=b(m,1)+b(m,3)
      call fcts(mset,qnom,b)
44    do 45 k=1,kq
45    pm(k,m)=(q(k,1)-qnom(k))/b(m,3)
50    b(m,1)=beta(m)
60    do 67 m=1,kq
67    xrt(m)=wres(m)
      do 75 j=1,kq
      do 75 k=1,kb
      xpmt(k,j)=0.d0
      do 75 m=1,kq
75    xpmt(k,j)=xpmt(k,j)+pm(m,k)*w(m,j)
      go to 70
57    if(lsize) 640, 630, 640
630   do 58 m=1,kb
58    b(m,1)=beta(m)
59    xlamda=xlamda*8.d0
      if(xlamda.gt.1.d35) go to 300
      go to 70
640   xlamda=xlamda*2.d0
      first=.true.
      go to 603
70    n=0
      do 80 j=1,kb
      do 78 k=1,kb
      dtemp=0.d0
      do 77 m=1,kq
77    dtemp=dtemp+xpmt(k,m)*pm(m,j)
      n=n+1
      if(j.eq.k) dtemp=dtemp+xlamda*b(k,6)
78    dxpm(n)=dtemp
      dtemp=0.d0
      do 79 m= 1,kq
79    dtemp=dtemp+pm(m,j)*xrt(m)
80    xr(j)=dtemp
      call matrix(xpm,kb,kb,xr,1,determ,id)
      if(id.eq.1) go to 1019
      go to 57
1019  do 2023 i=1,kb
2023  db(i)=xr(i)
      nostop=.false.
      do 82 m=1,kb
      if(dabs(db(m)).gt.dabs(clamp*b(m,3))) go to 57
      if(dabs(db(m))-dabs(b(m,1)*twom28)) 82, 81, 81
81    nostop=.true.
      b(m,1)=b(m,1)+db(m)
      if(b(m,1).le.b(m,5)) go to 83
      b(m,1)=b(m,5)
      go to 82
83    if(b(m,1).lt.b(m,4)) b(m,1)=b(m,4)
82    continue
      if(.not.nostop) go to 302
86    call fcts(mset,q,b)
87    if(mset(4).eq.0) go to 84
      write(6,2062) l,(b(m,1),m=1,kb)
2062  format(1h0,i5,(1h ,1p5d10.2))
      write(6,2063) (q(m,1),m=1,kq)
2063  format(6x,(1h ,1p5d10.2))
84    if(xlamda.le.1.d35) go to 109
85    nostop=.false.
      write(6,3005)
3005  format(28h lambda has become too large)
109   if(nostop) go to 4
      if(reduce-1.d0) 550, 525, 550
525   convrg=.true.
      go to 299
110   write(6,3000)
3000  format(44h0the number of iterations have been exceeded)
      go to 300
550   convrg=.false.
299   return
300   nostop=.false.
302   reduce=1.d0
      call fcts(mset,q,b)
      do 301 m=1,kq
      if(q(m,1).ge.q(m,2).and.q(m,1).le.q(m,3)) go to 301
      reduce=0.01d0
      go to 87
301   continue
      go to 87
      end
*
      subroutine matrix(a,n1,i1,b,m1,determ,id)
c
      implicit real*8(a-h,o-z)
      dimension a(i1,i1),b(i1,m1),index(50,3)
      equivalence (irow,jrow),(icolum,jcolum),(amax,t,swap)
      data eps4/1.d-7/
c
c     initialization
      bmax=0.d0
      id=1
      m=m1
      n=n1
10    determ=1.d0
15    do 20 j=1,n
20    index(j,3)=0
30    do 550 i=1,n
c
c     search for pivot element
c
40    amax=0.d0
45    do 105 j=1,n
      if(index(j,3).eq.1) go to 105
60    do 100 k=1,n
      if ( index(k,3)-1) 80, 100, 100
80    if(amax.ge.dabs(a(j,k))) go to 100
85    irow=j
90    icolum=k
      amax=dabs(a(j,k))
100   continue
105   continue
      index(icolum,3)=index(icolum,3)+1
260   index(i,1)=irow
270   index(i,2)=icolum
      if(bmax.eq.amax) id=2
c
c     interchange rows to put pivot element on diagonal
c
130   if(irow.eq.icolum) go to 310
140   determ=-determ
150   do 200 l=1,n
160   swap=a(irow,l)
170   a(irow,l)=a(icolum,l)
200   a(icolum,l)=swap
      if(m) 310, 310, 210
210   do 250 l=1,m
220   swap=b(irow,l)
230   b(irow,l)=b(icolum,l)
250   b(icolum,l)=swap
c
c     divide pivot row by pivot element
c
  310 pivot =a(icolum,icolum)
c
c     eliminate overflow and underflow due to either pivot or determ
c
      if(pivot) 311, 312, 313
311   dabpvt=-pivot
      go to 314
312   pivot=1.d-35
      go to 315
313   dabpvt=pivot
314   if(dabpvt.ge.1.d+35)pivot=sign(1.d+35,pivot)
      if(dabpvt.le.1.d-35)pivot=sign(1.d-35,pivot)
315   if(determ) 316, 317, 318
316   dabdet=-determ
      go to 319
317   determ=1.d-35
      go to 320
318   dabdet=determ
319   if(dabdet.ge.1.d+35) determ=sign(1.d+35,determ)
      if(dabdet.le.1.d-35) determ=sign(1.d-35,determ)
c
c     end of elimination of overflow
c
320   determ=determ*pivot
330   a(icolum,icolum)=1.d0
340   do 352 l=1,n
c
c     elimination of overflow
c
c     if(pivot .ne. 0.d0) go to 350
      if(dabpvt .gt. 1.d-35) go to 350
c
c     end of elimination of overflow
c
      a(icolum,l)=0.d0
      go to 352
350   a(icolum,l)=a(icolum,l)/pivot
352   continue
355   if(m) 380, 380, 360
360   do 379 l=1,m
c
c     elimination of overflow
c
c     if(pivot .ne. 0.d0) goto 370
      if(dabpvt .gt. 1.d-35) go to 370
c
c     end of elimination of overflow
c
      b(icolum,l)=0.d0
      go to 379
370   continue
c 370 b(icolum,l)=b(icolum,l)/pivot
      b1=b(icolum,l)
      if(b1) 372, 373, 374
372   dab1=-b1
      go to 375
373   b(icolum,l)=1.d-35
      go to 376
374   dab1=b1
375   if(dab1.ge.1.d+35) b(icolum,l)=sign(1.d+35,b1)
      if(dab1.le.1.d-35) b(icolum,l)=sign(1.d-35,b1)
376   continue
      b(icolum,l)=b(icolum,l)/pivot
379   continue
c
c     reduce non-pivot rows
c
      bmax=0.d0
380   do 550 l1=1,n
390   if(l1.eq.icolum) go to 550
400   t=a(l1,icolum)
420   a(l1,icolum)=0.d0
430   do 450 l=1,n
      if(t) 431, 432, 433
431   dabt=-t
      go to 434
432   t=1.d-35
      go to 435
433   dabt=t
434   if(dabt.ge.1.d+35) t=sign(1.d+35,t)
      if(dabt.le.1.d-35) t=sign(1.d-35,t)
435   continue
      sub=a(icolum,l)*t
      a(l1,l)=a(l1,l)-sub
      if(index(l1,3).eq.1) go to 450
      if(dabs(a(l1,l))-eps4*dabs(sub)) 449, 449, 450
449   bmax=dmax1(bmax, dabs(a(l1,l)))
450   continue
455   if(m) 550, 550, 460
460   do 500 l=1,m
500   b(l1,l)=b(l1,l)-b(icolum,l)*t
550   continue
c
c     interchange columns
c
600   do 710 i=1,n
610   l=n+1-i
620   if(index(l,1) .eq. index(l,2)) go to 710
630   jrow=index(l,1)
640   jcolum=index(l,2)
650   do 705 k=1,n
660   swap=a(k,jrow)
670   a(k,jrow)=a(k,jcolum)
700   a(k,jcolum)=swap
705   continue
710   continue
      do 730 k=1,n
      if(index(k,3) .ne. 1) go to 715
730   continue
      return
715   id=2
      return
      end


