g========+=========+=========+=========+=========+=========+=========+=$
C PROGRAM: levy
C TYPE   : main
C PURPOSE: perform Levy construction 
C          z_1 is the initial configuration
C          z_N is the final position ... and they are both fixed. 
C          The Levy construction allows to sample the intermediate
C          points z_2.....z_(N-1)  such that probability to sample them is
C          P(z_2, z_(N-1)) \sim
C          exp(-beta*(z_1-z_2)^2/2) x exp(-beta*(z_2-z_3)^2/2) ..
C                            ... x  exp(-beta*(z_(N-1)-z_N)^2/2)
C VERSION: 5-ma-2001 (bugs fixed in description, one in code) 
C========+=========+=========+=========+=========+=========+=========+=$
      parameter (N=50)
      dimension z(N),zmoy(N)
cc
cc    ...this is just an example
cc
      z(1)=0.
      z(N)=4.2
      beta=1.3
      Idum=-19199      

      do iter=1,10
cc
cc    ...end of example
cc
cc    to understand the following lines, you have to check the following:
cc
cc    suppose p(z_1,z_2) \sim exp(-(z_1 -z_2)^2/(2 tau_1)) and
cc    suppose p(z_2,z_3) \sim exp(-(z_2 -z_3)^2/(2 tau_2)) and
cc
cc    consider P(z1, z_2, z_3) =  p(z_1,z_2) x p(z_2,z_N)
cc
cc    first:  P(z_3) (fixing z_1, but integrating over z_2) is 
cc
cc            P(z_3) \sim exp(-(z_1 -z_3)^2/(2 (tau_1 + tau_2)))
cc
cc
cc    second: P(z_2) (fixing z_1 and z_3) is
cc
cc            P(z_2) \sim exp(-(z_2 -z_mean)^2/(2 tau_3)) 
cc
cc    with z_mean = (tau_2 z_1 + tau_1 z_3)/(tau_1 + tau_2)
cc    
cc    and 1/tau_3= 1/( 1/tau_1 + 1/tau_2)
cc
cc        ... this is easily done by writing out P(z_2) and P(z_3) in full
cc            detail and by completing squares.
cc
      zend=z(N)
      do j = 2, N-1
         zinit= z(j-1)
         tau_1=1./beta
         tau_2=(N-j)/beta
         zmean=(tau_2* zinit + tau_1 * zend)/(tau_1 + tau_2)
         sigmasq=1./( 1./tau_1 + 1./tau_2)
         sigma= sqrt(sigmasq)
         z(j)=zmean + gauss(Idum, sigma)
      end do
      do i=1,n
         write(*,*)i,z(i)
         zmoy(i)=zmoy(i)+z(i)
      end do
      write(*,*)
      end do
cc
cc    Here we choose one particular order for computing
cc    the z_i...Another popular choice consists in first computing
cc    z_(N/2)...then z_(N/4) and z_(3 N/4) ... then ... 
cc
      do i=1,n
         write(*,*)i,zmoy(i)/10.
      end do

      end
C========+=========+=========+=========+=========+=========+=========+=$
C PROGRAM: gauss
C TYPE   : subroutine
C PURPOSE: gaussian random numbers
C========+=========+=========+=========+=========+=========+=========+=$
      function gauss(Idum,sigma)
      data iset/0/
      save

      if (iset.eq.0) then
1       continue
        xx = ran3(Idum)
         x = 2. * xx - 1
        yy = ran3(Idum)
         y = 2. * yy - 1.
         r2tilde= x**2 + y**2
         if (r2tilde.ge.1.) go to 1
         fac = sqrt( -2. * log(r2tilde) / r2tilde)
         gset = x * fac
         gauss = y * fac
         iset = 1
      else
         gauss = gset
         iset = 0
      end if
      gauss = gauss * sigma
      end
C========+=========+=========+=========+=========+=========+=========+=$
C PROGRAM: ran3.f  
C TYPE   : function
C PURPOSE: generate random numbers
C I/O    :
C	
C VERSION: 14-NOV-95
C COMMENT: Initialize Idum with negative integer
C========+=========+=========+=========+=========+=========+=========+=$
      real function ran3(Idum)
      Parameter (mbig=1000000000,Mseed=161803398,Mz=0,fac=1./Mbig)
      Dimension MA(55)
      save
      if (Idum.lt.0.or.iff.eq.0) then
         iff=1
         mj=mseed-iabs(Idum)
         mj=mod(mj,mbig)
         ma(55)=mj
         mk=1
         do 11 i=1,54
            ii=mod(21*i,55)
            ma(ii)=mk
            mk=mj-mk
            if (mk.lt.mz) mk=mk+mbig
            mj=ma(ii)
 11      continue
         do 13 k=1,4
            do 12 i=1,55
               ma(i)=ma(i)-ma(1+mod(i+30,55))
               if (ma(i).lt.mz)ma(i)=ma(i)+mbig
 12         continue
 13      continue
         inext=0        
         inextp=31
         Idum=1
      end if
      inext=inext+1
      if (inext.eq.56) inext=1
      inextp=inextp+1
      if (inextp.eq.56) inextp=1
      mj=ma(inext)-ma(inextp)
      if (mj.lt.mz) mj=mj+mbig
      ma(inext)=mj
      ran3=mj*fac
      end