program dipsys

*     Simulation of mixed magnetic dipolar system
*     -mu1 = 1 only concentration for step 1
*     -mu1 and mu2 for step 2

*     initilization of the system

      integer p,e,n,o,iter,itnum,sampl,t,ok,
     $                                          s,s1,chk0
      parameter (n = 24)

      real r,z,dy,dx,q,mu1,mu2,mind(n),dm(n,n),x(n),y(n),theta(n),ave,zerod,
     $                   dpl(n),alpha,rad(n),l,oldenergy,newenergy,errore

***      write(*,*) 'Input alpha parameter  : '
***      read(*,*) alpha
***      write(*,*) 'sampling frequency for the distances analisys : '
***      read(*,*) sampl


      alpha = .25
      sampl = 10000000

      o=1
      r=1
*      alpha = .25
      errore = 6
*           with r=1,n=20 l=alpha*0.2236068= .5*...=0.1118034
*                       n=40 0.3162277  alpha=2
      z=n
      l=alpha*sqrt(1/z)
      pi = 3.14159265
      mu1=1
      mu2=0
      open(1,file='os.dip')
************      open(9,file='mu.dip')
*      open(2,file='ns.dip')
      open(3,file='en.dip')
      do 2 e=1,n
            theta(e) = mu2
            if (e .lt. (n/2)) theta(e)=mu1
20          x(e) = (rand(0)-0.5)*2*r
	    y(e) = (rand(0)-0.5)*2*r
	    if ((x(e)**2+y(e)**2).ge.(r**2)) goto 20
       	    write(1,62) x(e),y(e),theta(e)
2     continue
      close(1)
*      COODINATES are CARTESIAN
*     system filled and random disposition of n particles done
*     CHOOSE RANDOMLY A PARTICLE AND MOVE IT RANDOM.
*     and iterate the process
*                ok=0 there isn't a defined separation in the lattice
      ok=0
      itnum=0
      DO WHILE (ok.eq.0)
	    ITNUM = ITNUM+1
	    p = rand(0)*n+1
*     calculate the energy of the choosen particle.
*     very IMPORTANT dx=0 dy=0 here !!
	    dx=0
	    dy=0
	    call energy(x,y,p,n,mu1,ene,dx,dy)
	    oldenergy = ene
*     n.b. here the movement depends on 2 random numbers
3           dx = (rand(0)-0.5)*2*r*l
	    dy= (rand(0)-0.5)*2*r*l
***check if the particle moves out of the system boundaries 
	    if (( (dx+x(p))**2+(dy+y(p))**2).ge.(r**2)) goto 3

***check if the particle tendes to superpose another when one of 
***                    two has mu=0
            chk0 = 0
            do 9 e=1,n
               if (e.ne.p) then
                  zerod = sqrt( (x(p)+dx-x(e))**2+(y(p)+dy-y(e))**2 ) 
                  if (zerod.lt.(1.0E-3)) chk0 = 1
               endif
 9             continue
            if (chk0 .eq. 1) goto 3   

*     new possible position is legal .... CHECK OF ENERGY
            if (theta(p).ne.0) then 
	       call energy(x,y,p,n,mu1,ene,dx,dy)
	       newenergy = ene
            else 
               newenergy = oldenergy
            endif

*     cfr the two values of energy new-old<0  make the move
	    if (newenergy.lt.oldenergy) then
**		  write(*,*) p,itnum,newenergy-oldenergy
		  x(p)=x(p)+dx
		  y(p)=y(p)+dy
		  write(3,*) itnum,(newenergy-oldenergy)
		  if (abs(newenergy-oldenergy).lt.(1.0E-12)) then
**                        write(*,*) newenergy-oldenergy
			stop
		  end if
	    else
*                 write(*,*) itnum,oldenergy,newenergy
	    endif
*     TRY TO BUILT THE DISTANCES MATRIX EVERY SAMPL-LOOP
	    if (itnum.eq.(sampl*o)) then
		  o=o+1
		  call dmatrix(x,y,dm,n)
		  open(4,file='dm.dip')
**		  open(5,file='av.dip')
**                  pause
                  open(2,file='ns.dip')
**                  open(6,file='status.dip')
                  write(4,62) ((dm(s,s1),s1=1,n),s=1,n)
                  write(2,62) (x(s),s=1,n), (y(s1),s1=1,n)
**     built the connectivity matrix assiciated with the lattice
                  
		  call dlattice(dm,mind,n,ok,ave,errore)
**                  write(6,*) p,itnum,newenergy-oldenergy
		  close(4)
**                  close(6)
**		  close(5)
                  close(2)
**                  write(*,*) p,itnum
      	    endif
      END DO
62    format (20f12.8)

*      do 200 e=1,n
*	    write(2,*) x(e),y(e)
*200   continue

*      close(6)
*      close(1)
*      close(2)
      close(3)
*      close(4)
*      close(5)
      end

      subroutine dlattice(dm,mind,num,ok,ave,err)
      real dm(num,num),mind(num),ave
      integer ok,m1,m2
	    ave=0
	    do m1=1,num
		  mind(m1)=dm(m1,1)
		  if (m1.eq.1) mind(m1)=dm(m1,2)
		  do m2=1,num

			if (m1.ne.m2) then
			   if (dm(m1,m2).lt.mind(m1)) mind(m1)=dm(m1,m2)
			end if
		  end do
		  ave=ave + mind(m1)
	    end do
	    ave=ave/num
	    ok=1
	    do m1=1,num
		  if (abs(((ave-mind(m1))/mind(m1))*100).gt.err) ok=0
***		  write(5,*) m1,ave,mind(m1)
	    end do
      end

*     COMPOSE THE DISTANCES MATRIX AND ORDER ITS ELEMENT !!!!
      subroutine dmatrix(x,y,dm,num)
      real x(num),y(num),dm(num,num),elem
      integer m1,m2,m,t,s1,s,nline
	    nline=num
	    do 50 il=1,num
		  dm(il,il)=0
		  do 60 jl=(il+1),num
			elem=sqrt((x(il)-x(jl))**2+(y(il)-y(jl))**2)
*     sorting ....
			m2=jl
			dm(il,m2)=elem
			dm(m2,il)=elem
60                continue
50          continue
      end

      subroutine energy(x,y,nump,num,mu1,ene,dx,dy)
      real r1,r2,ene,dist,x(num),y(num),mu1
      ene = 0
      do 10 i=1,num
	    if (i.ne.nump) then

	    dist=sqrt(( x(nump)+dx-x(i) )**2+( y(nump)+dy-y(i) )**2)
		  ene=ene+(mu1*mu1)/(dist*dist*dist)
	    end if
10      continue
      end