implicit none
      integer d1, L, MCS, i, j, k, FInew, dFI, D, index, m,o,t
      real ran1, P, E, dH, n, ne, no
      parameter (L = 20,D=70, MCS = 230000, P=3.14159265359, dFI = 10)
      parameter (ne = 1.7, no = 1.5)
      integer fi(1:D,1:L), nj(1:L), pj(1:L), ni(1:D), pi(1:D)
      real P2(0:360)
      open(8, file = 'L=20,D=70.txt', status = 'unknown')
      !-----------------------------------------------------------------

      d1=-1
      do i = 1, D
         do j = 1, L
              fi(i,j) = 0
         enddo
      enddo

      do i = 1, D
         ni(i) = i + 1
         pi(i) = i - 1
      enddo
      ni(D) = 1
      pi(1) = D

      do i = 1, L
         nj(i) = i + 1
         pj(i) = i - 1
      enddo

      do i = 0,360
         P2(i) = 0.5*(3*cos(((i-180)*P)/180.0)**2-1)
      enddo

      E=0.0
      do t=0,215
        write(*,*) E
        index = 0
        n = 0
        do k = 1, MCS
           do i = 1, D
              do j = 2, L-1
                 FInew = fi(i,j) + int((ran1(d1) - 0.5)*dFI)
                 if(FInew>90) Finew = Finew - 180
                 if(FInew<-90) Finew = Finew + 180
                 dH=-20*(P2(FInew-fi(pi(i),j)+180)+
     &          P2(FInew-fi(i,pj(j))+180)+P2(FInew-fi(ni(i),j)+180)+
     &          P2(FInew-fi(i,nj(j))+180))-E**2*P2(90-FInew +180)
     &          +20*(P2(fi(i,j)-fi(pi(i),j)+180)+
     &          P2(fi(i,j)-fi(i,pj(j))+180)+P2(fi(i,j)-fi(ni(i),j)+180)+
     &          P2(fi(i,j)-fi(i,nj(j))+180))+E**2*P2(90-fi(i,j)+180)
                 if(dH<0) then
                    fi(i,j) = FInew
                 else if(ran1(d1)<=exp(-dH)) then
                    fi(i,j) = FInew
                 endif
              enddo
           enddo
           if(k>30000.and.mod(k,100)==0) then
              index = index + 1
              do m = 1, D
                 do o = 1, L
                    n=n+(ne*no)/(no**2*cos(fi(m,o)/180.0*p)**2+
     &              ne**2*sin(fi(m,o)/180.0*p)**2)**(0.5)
                 enddo
              enddo
           endif
        enddo
        n = n/(index*L*D)
        write(8,*) E, n
        if(E>0.35.and.E<0.44) then
           E=E+0.01
        else
           E=E+0.02
        endif
      enddo
      end
      !.................................................................

      FUNCTION ran1(idum)
      INTEGER idum,IA,IM,IQ,IR,NTAB,NDIV
      REAL ran1,AM,EPS,RNMX
      PARAMETER (IA=16807,IM=2147483647,AM=1./IM,IQ=127773,IR=2836,
     *NTAB=32,NDIV=1+(IM-1)/NTAB,EPS=1.2e-7,RNMX=1.-EPS)
      INTEGER j,k,iv(NTAB),iy
      SAVE iv,iy
      DATA iv /NTAB*0/, iy /0/
      if (idum.le.0.or.iy.eq.0) then
        idum=max(-idum,1)
        do 11 j=NTAB+8,1,-1
          k=idum/IQ
          idum=IA*(idum-k*IQ)-IR*k
          if (idum.lt.0) idum=idum+IM
          if (j.le.NTAB) iv(j)=idum
11      continue
        iy=iv(1)
      endif
      k=idum/IQ
      idum=IA*(idum-k*IQ)-IR*k
      if (idum.lt.0) idum=idum+IM

      j=1+iy/NDIV
      iy=iv(j)
      iv(j)=idum
      ran1=min(AM*iy,RNMX)
      return
      END