program ising2
c
c ***********************************************************************
c       Simple Ising model simulation on square lattice.
c ***********************************************************************
c
	parameter (maxn=101,maxt=100001)
	implicit real*8 (a-h,o-z)
	common isp(0:maxn,0:maxn)
	dimension ord(maxt)
	character*12 of1,of2
	real etime,time(2)
10	print *,'Ising model on square lattice: lattice edge length?'
	read(5,*) ne
	if(ne.gt.maxn) go to 10
	print *,'initial spin config (1 up, -1 down), ranf seed ?'
	read(5,*) init,isd
	print *,'start and end temp (J/k), #steps ?'
	read(5,*) t1,t2,nt
	print *,'start and end fields (J/mu), #steps ?'
	read(5,*) h1,h2,nh
	print *,'ms0(init trans), ms(rec int), ns(#rec)?'
	read(5,*) ms0,ms,ns
	if(nt.gt.0.or.nh.gt.0) then
		print *,'reinitialize for each T/H ? [1 yes, 0 no]'
		read(5,*) f2
		print *,'transient per T/H ?'
		read(5,*) ms1
	endif
	if(ms0+(ns-1)*ms+1.gt.maxt.or.ms1+(ns-1)*ms+1.gt.maxt) then
		print *,'maxt is too small!'
		go to 10
	endif
	print *,'periodic (1) or free (2) boundary ?'
	read(5,*) ib
	print *,'output file for thermodynamic functions ?'
	read(5,15) of1
15	format(a12)
	print *,'output file for magnetization time series ?'
	read(5,15) of2
	time0=etime(time)
	if(nt.gt.0) then
		dt=(t2-t1)/nt
	else
		dt=0.d0
	endif
	if(nh.gt.0) then
		dh=(h2-h1)/nh
	else
		dh=0.d0
	endif
	write(6,35)
35	format(/2x,'Temperature',2x,'Magnetic fld',1x,'Magnetization',2x,'Energy',2x,
     c         'Susceptibility',1x,'Heat Capacity'/)
	m1=0
	call init_spin(ne,init,ord,m1,maxt,isd)
	call transient(t1,h1,ne,ms0,ib,ord,m1,maxt,isd)
	if(of1.ne.'') then
		open(1,file=of1)
		rewind(1)
	endif
	do 20 i=1,nt+1
		t=t1+(i-1)*dt
		do 30 j=1,nh+1
		h=h1+(j-1)*dh
		if(i.gt.1.or.j.gt.1) then
			if(f2.eq.1) call init_spin(ne,init,ord,m1,maxt,isd)
			call transient(t,h,ne,ms1,ib,ord,m1,maxt,isd)
		endif
	        call calc(t,h,or,en,chi,ch,ne,m1,ms,ns,ib,ord,maxt,isd)
		if(of1.ne.'') write(1,25) t,h,or,en,chi,ch
		write(6,25) t,h,or,en,chi,ch
25		format(1x,1p,6(e12.5,1x))
30		continue
20	continue
	if(of1.ne.'') close(1)
	if(of2.ne.'') then
		open(2,file=of2)
		rewind(2)
		do 50 m=1,m1
			o1=ord(m)/dfloat(ne)**2
			write(2,55) m,o1
55			format(1x,i6,1x,1p,e12.5)
50		continue
		close(2)
	endif
	time1=etime(time)
	write(6,45) time1-time0
45	format(/1x,'CPU time = ',1p,e12.5,' sec')
	stop
	end
c
c
c
	subroutine init_spin(ne,init,ord,m1,maxt,isd)
	parameter (maxn=101)
	implicit real*8 (a-h,o-z)
	real ranf
	common isp(0:maxn,0:maxn)
	dimension ord(1)
	do 70 j=0,ne+1
	isp(0,j)=0
70	isp(ne+1,j)=0
	do 80 i=1,ne
	isp(i,0)=0
80	isp(i,ne+1)=0
	m1=m1+1
	if(m1.gt.maxt) then
		ord(1)=ord(m1-1)
		m1=2
	endif
	if(init.eq.1.or.init.eq.-1) then
		do 10 i=1,ne
			do 15 j=1,ne
				isp(i,j)=init
15			continue
10		continue
		ord(m1)=dfloat(init*ne**2)
	else
		ord(m1)=0.d0
		do 20 i=1,ne
			do 25 j=1,ne
				if(ranf(isd).le.0.5) then
					isp(i,j)=1
				else
					isp(i,j)=-1
				endif
				ord(m1)=ord(m1)+isp(i,j)
25			continue
20		continue
	endif
	return
	end
c
c
c
	subroutine transient(t,h,ne,ms0,ib,ord,m1,maxt,isd)
	parameter (maxn=101)
	implicit real*8 (a-h,o-z)
	real ranf
	common isp(0:maxn,0:maxn)
	dimension ord(1)
	do 30 m=1,ms0
		m1=m1+1
		if(m1.gt.maxt) then
			ord(1)=ord(m1-1)
			m1=2
		endif
		ord(m1)=ord(m1-1)
		do 40 i=1,ne
			i1=i-1
			if(ib.eq.1.and.i.eq.1) i1=ne
			i2=i+1
			if(ib.eq.1.and.i.eq.ne) i2=1
			do 50 j=1,ne
				j1=j-1
				if(ib.eq.1.and.j.eq.1) j1=ne
				j2=j+1
				if(ib.eq.1.and.j.eq.ne) j2=1
				isum=isp(i1,j)+isp(i2,j)
     c					+isp(i,j1)+isp(i,j2)
				de=2*isp(i,j)*(isum+h)
				if(de.lt.0) then
					isp(i,j)=-isp(i,j)
                                        ord(m1)=ord(m1)+2*isp(i,j)
				elseif(ranf(isd).le.dexp(-de/t)) then
					isp(i,j)=-isp(i,j)
                                        ord(m1)=ord(m1)+2*isp(i,j)
				endif
50			continue
40		continue
30	continue
	return
	end
c
c
c
	subroutine calc(t,h,or,en,chi,ch,ne,m1,ms,ns,ib,ord,maxt,isd)
	parameter (maxn=101)
	implicit real*8 (a-h,o-z)
	real ranf
	common isp(0:maxn,0:maxn)
	dimension ord(1)
c
c initialize sum,sum2,sume,sume2 for the quantities to compute
c
	sum=0.d0
	sum2=0.d0
	sume=0.d0
	sume2=0.d0
	do 10 k=1,ns
	xmag=0.d0
	xe=0.d0
	if(k.eq.1) then
		m2=1
	else
		m2=ms+1
	endif
	do 30 m=1,m2
		if(m.lt.m2) then
			m1=m1+1
			if(m1.gt.maxt) then
				ord(1)=ord(m1-1)
				m1=2
			endif
			ord(m1)=ord(m1-1)
		endif
		do 40 i=1,ne
			i1=i-1
			if(ib.eq.1.and.i.eq.1) i1=ne
			i2=i+1
			if(ib.eq.1.and.i.eq.ne) i2=1
			do 50 j=1,ne
				j1=j-1
				if(ib.eq.1.and.j.eq.1) j1=ne
				j2=j+1
				if(ib.eq.1.and.j.eq.ne) j2=1
				isum=isp(i1,j)+isp(i2,j)
     c					+isp(i,j1)+isp(i,j2)			
				if(m.lt.m2) then
				   de=2*isp(i,j)*(isum+h)
				   if(de.lt.0) then
					isp(i,j)=-isp(i,j)
					ord(m1)=ord(m1)+2.d0*isp(i,j)
				   elseif(ranf(isd).le.dexp(-de/t)) then
					isp(i,j)=-isp(i,j)
					ord(m1)=ord(m1)+2.d0*isp(i,j)
				   endif
				else
				   xmag=xmag+isp(i,j)
				   xe=xe-0.5d0*isp(i,j)*isum-h*isp(i,j)
				endif
50			continue
40		continue
30	continue
c
c compute mag,energy add these to the sum,sum2
c
	sum=sum+xmag
	sum2=sum2+xmag**2
	sume=sume+xe
	sume2=sume2+xe**2
10	continue
c
c average and compute mag,energy,chi,ch and put them in the common
c
	or=sum/dfloat(ns)
	chi=(sum2/dfloat(ns)-or**2)/(t*dfloat(ne**2))
	or=or/dfloat(ne**2)
	en=sume/dfloat(ns)
	ch=(sume2/dfloat(ns)-en**2)/(t*ne)**2
	en=en/dfloat(ne**2)
	return
	end
c
c ***********************************************************************
c       Random number generator
c ***********************************************************************
c
      real function ranf(iseed)
      iseed=iseed*1566083941+1
      ranf=float(iseed)*2.328306e-10+0.5
      return
      end