implicit double precision(a-h,l,m,o-z)
      parameter(np=4096)
      dimension vr(np),vk(np),pot(np),epot(np),expy(np)
      dimension yr(np),yk(np),yn(np),cr(np),ck(np),alfa(np)
     +         ,sk(np),gr(np)
      save
     
      pi=dacos(-1.d0)

C ---> External Parameters
      Write(6,*) "rho?"
      read(5,*) rho
      Write(6,*) "temp?"
      read(5,*) temp
      Write(6,*) "nclos?"
      read(5,*) nclos
      Write(6,*) "istart?"
      read(5,*) istart
      eta=pi*rho/6.d0
      beta=1.d0/temp
      amix=0.50d0
      tol=1.d-10
      itermax=10000

C ---> Setup of the Grid Mesh
      dr=0.005d0
      dk=pi/dble(np)/dr
      do i=1,np
         vr(i)=dble(i)*dr
         vk(i)=dble(i)*dk
      enddo

C ---> Build the potential
      sigma=1.d0
      pot=0.d0
      do i=1,np
         if (vr(i).le.sigma) then
         pot(i)=0.500E+35
         epot(i)=0.d0
         else
         pot(i)=0.d0
         epot(i)=1.d0
         endif
      enddo

      open(1,file="pot.out")
      do i=1,np
         write(1,100) vr(i),pot(i)
      enddo
      close(1)

C ---> Gamma's
      yr=0.d0
      if (istart.ne.0) then
         open(2,file="gamma.out")
         do i=1,np
            read(2,*) vrr,yr(i)
         enddo
         close(2)
      endif

C ==== MAIN LOOP ====

 1000 iter=iter+1

      cr=-1.d0-yr
      expy=dexp(yr)
      if(nclos.eq.0) then ! PY
        where(epot.gt.0.d0) cr=(epot-1.d0)*(yr+1.d0)
      elseif(nclos.eq.1) then ! HNC
        where(epot.gt.0.d0) cr=epot*expy-(yr+1.d0)
      elseif(nclos.eq.2) then ! MSA
        where(epot.gt.0.d0) cr=-beta*pot
      endif

      ck=cr
      call fft(ck,np,vr,vk,1)   
      yk=ck/(1.d0-rho*ck)-ck
      yn=yk
      call fft(yn,np,vr,vk,-1)   
      dy=sqrt(sum((yn-yr)**2)*dr)
      write(6,100) iter,dy
      yr=(1.d0-amix)*yn+amix*yr
      if (dy.gt.tol.and.iter.le.itermax) goto 1000

      sk=1.d0/(1.d0-rho*ck)
      gr=yr+cr+1.d0

      open(3,file="str.out")
      do i=1,np
         write(3,100) vr(i),gr(i)
      enddo
      close(3)

      open(4,file="sk.out")
      do i=1,np
         write(4,100) vk(i),sk(i)
      enddo
      close(4)

      open(5,file="gamma.out")
      do i=1,np
         write(5,100) vr(i),yr(i)
      enddo
      close(5)

c ---> Calculation of the internal energy
      alfa=pot*gr
      where(gr.le.0.0001d0) alfa=0.d0
      en=2.d0*pi*rho*dr*sum(alfa(:)*vr*vr)
      write(6,*) en
            
 100  format(1x,10(1x,g13.6))
 101  format('#',1x,10(1x,g13.6))
      end