!Program to simulate a solid at low temperature
!By Martin Eyles (9814937)
!Created 18/02/2000   -   Last Update 16/03/2000

!module containing my subroutines
MODULE subprogs
IMPLICIT NONE
!keep global the variables stated below
SAVE
	!grid3 is a pointer, that stores all the data for two solid blocks
	!grid1 and grid2, point to relavent sections of grid3, and represent
	!one block of solid each (NB - pointers (1,2) set to point to 3 in
	!routine dimgrid)
	!active set to point to which ever of the grid is being worked on
	!this is set along with control, in routine change_control
	!Finally - energy (which could alternatively be an allocatable array)
	!store the no. of sites with each quanta
	INTEGER, POINTER, DIMENSION(:,:) :: grid1,grid2,grid3,active
	INTEGER, POINTER, DIMENSION(:) :: energy
	INTEGER :: control

CONTAINS
!The actual subroutines

!Subroutine to find equilibrium in the active grid
SUBROUTINE eqm
	IMPLICIT NONE
	!i for itterations of loops, integrals store the sums of quanta
	!in the 2 loops, pen stores the no. of quanta at the peak.
	INTEGER :: i, integral1, integral2, pen
	
	!Outer loop repeats till eqm
	DO
		integral1=0
		integral2=0

		DO i=1,1000
			CALL swap
			!only add to integral every 20 iterations - speed up
			!the program
			IF (MODULO(i,20)==0) THEN
				CALL count
				!sum maximum, and zero energy
				integral1=integral1+energy(0)
				pen=MAXVAL(energy)
				integral1=integral1+pen
			END IF
		END DO
		DO i=1,1000
			CALL swap
			IF (MODULO(i,20)==0) THEN
				CALL count
				integral2=integral2+energy(0)
				pen=MAXVAL(energy)
				integral2=integral2+pen
			END IF
		END DO
		!compare two sums/integral then if close enough -> eqm reched
		!->exit do loop (no due to the large number of iterations summed
		!(1000/20=50), this is set to zero -> maximum accuracy, but slow.
		IF ((integral1<=(integral2+0)).AND.(integral1>=(integral2-0))) EXIT
	END DO

	!tell the user that eqm has been found
	PRINT*,'Equilibrium Found'
END SUBROUTINE eqm

SUBROUTINE dimgrid
!subroutine to set up block dimensions
	IMPLICIT NONE
	!user inputs of dimensions fro blocks
	INTEGER :: width1,width2,height
	
	!stop pointers pointing to grid3, so it can be dealocated
	NULLIFY (active,grid1,grid2)
	!deallocate grid 3 - to be re-allocated later
	DEALLOCATE (grid3)
	
	!prompt the user for values of grid sizes
	PRINT*,'Enter Grid 1 width'
	READ*,width1
	PRINT*,'Enter Grid 2 width'
	READ*,width2
	PRINT*,'Enter Grid Height'
	PRINT*,'(Same value used for 1 + 2)'
	READ*,height
	
	!allocate grid3, and point grids 1 and 2 to the appropriate
	!section thereof
	ALLOCATE (grid3((width1+width2),height))
	grid1 => grid3 (1:width1,1:height)
	grid2 => grid3 ((width1+1):(width1+width2),1:height)
	
	!point active to the grid that was being used before dealocation
	SELECT CASE(control)
		CASE(1)
			active=>grid1
		CASE(2)
			active=>grid2
		CASE(3)
			active=>grid3
	END SELECT
		
	!must set temperature values again
	CALL temp
END SUBROUTINE dimgrid

SUBROUTINE change_control
!subroutine to change the active grid
	IMPLICIT NONE
	!grid to change to
	INTEGER :: new
	
	PRINT*,'change to'
	DO
		READ*,new
		SELECT CASE(new)
			CASE(1)
				!they chose 1, so point to this, and set control
				!variable, for other routines
				control=new
				active=>grid1
				EXIT
			CASE(2)
				control=new
				active=>grid2
				EXIT
			CASE(3)
				control=new
				active=>grid3
				EXIT
			CASE DEFAULT
				!the input was invalid, so tell the user
				PRINT*,'grid must be 1,2 or 3 - (3 is 1&2 together)'
		END SELECT
	END DO
END SUBROUTINE change_control

SUBROUTINE temp
!routine to set grid temperature
	IMPLICIT NONE
	!inputs of temp and variable to store rand. no.
	REAL :: t1,t2,randno
	!no.s of quanta to place in grids, i as iteration, and
	!x and y are grid co-ords
	INTEGER :: quant1,quant2,i,x,y
	
	!erase grid data, ready for random quanta placement
	grid3=0
	
	!prompt for grid temp.
	PRINT*,'Enter grid1 temperature'
	READ*,t1
	PRINT*,'Enter grid2 temperature'
	READ*,t2
	
	!calculate no. of quanta to place in each grid
	quant1=INT(t1*SIZE(grid1))
	quant2=INT(t2*SIZE(grid2))
	
	!place quanta randomly on grid 1
	DO i=1,quant1
		CALL RANDOM_NUMBER(randno)
		x=(AINT(SIZE(grid1,1)*randno))+1
		CALL RANDOM_NUMBER(randno)
		y=(AINT(SIZE(grid1,2)*randno))+1
		grid1(x,y)=grid1(x,y)+1
	END DO

	!place quanta randomly on grid 2
	DO i=1,quant2
		CALL RANDOM_NUMBER(randno)
		x=(AINT(SIZE(grid2,1)*randno))+1
		CALL RANDOM_NUMBER(randno)
		y=(AINT(SIZE(grid2,2)*randno))+1
		grid2(x,y)=grid2(x,y)+1
	END DO
END SUBROUTINE temp

SUBROUTINE status
!print active grid, grid sizes, and grid temperatures
	IMPLICIT NONE
	!x,y dimensions of grid - used for all three, with values changed
	!i,j - iterate through grid     d1,d2,dy x size of grid 1&2
	!y size of grid (all same)
	INTEGER :: x,y,i,j,quanta,d1,d2,dy
	!temps of grids
	REAL :: temp1,temp2,temp3

	quanta=0
	x=SIZE(grid1,1)
	y=SIZE(grid1,2)
	!count quanta in grid
	DO i=1,x
		DO j=1,y
			quanta=quanta+grid1(i,j)
		END DO
	END DO	
	!calculate temperature, from quanta and grid area
	temp1=FLOAT(quanta)/(x*y)
	
	quanta=0
	x=SIZE(grid2,1)
	y=SIZE(grid2,2)
	DO i=1,x
		DO j=1,y
			quanta=quanta+grid2(i,j)
		END DO
	END DO	
	temp2=FLOAT(quanta)/(x*y)

	quanta=0
	x=SIZE(grid3,1)
	y=SIZE(grid3,2)
	DO i=1,x
		DO j=1,y
			quanta=quanta+grid3(i,j)
		END DO
	END DO	
	temp3=FLOAT(quanta)/(x*y)
	
	!find grid sizes
	d1=SIZE(grid1,1)
	d2=SIZE(grid2,1)
	dy=SIZE(grid3,2)	

	!output calculated info to user
	PRINT*,'Active Grid - ',control
	PRINT*
	PRINT*,'Grid Temperatures'
	PRINT*,'Grid1 - ',temp1
	PRINT*,'Grid2 - ',temp2
	PRINT*,'Grid3 - ',temp3
	PRINT*
	PRINT*,'Grid Dimensions'
	PRINT*,'Grid1 - ',d1,' * ',dy
	PRINT*,'Grid2 - ',d2,' * ',dy
	PRINT*,'Grid3 - ',d1+d2,' * ',dy
END SUBROUTINE status

SUBROUTINE advance
!routine to iterate the grid n times
	IMPLICIT NONE
	!n - no. of times to iterate - i, current iteration
	INTEGER :: n,i	

	!ask how many times to iterate
	PRINT*,'enter no. to advance'
	READ*,n
	
	!perform "swap" operation n times
	DO i=1,n
		CALL swap
	END DO
END SUBROUTINE advance

SUBROUTINE swap
!routine to randomly move quanta
	IMPLICIT NONE
	!x1,y1 - from co-ords   x2,y2 - to co-ords
	!xmax,ymax - grid dimensions
	INTEGER :: x1,y1,x2,y2,xmax,ymax
	!random no. variable
	REAL :: randno

	!find grid size
	xmax=SIZE(active,1)
	ymax=SIZE(active,2)

	!choose from position randomly
	DO
		CALL RANDOM_NUMBER(randno)
		x1=(AINT(xmax*randno))+1
		CALL RANDOM_NUMBER(randno)
		y1=(AINT(ymax*randno))+1
		!only valid if from position has quanta
		IF (active(x1,y1)>0) EXIT
	END DO

	!choose to position randomly
	DO
		CALL RANDOM_NUMBER(randno)
		x2=(AINT(xmax*randno))+1
		CALL RANDOM_NUMBER(randno)
		y2=(AINT(ymax*randno))+1
		!only valid if not same as from position
		IF ((x1/=x2).OR.(y1/=y2)) EXIT
	END DO

	!move one quanta from (x1,y1) to (x2,y2)
	active(x1,y1)=active(x1,y1)-1
	active(x2,y2)=active(x2,y2)+1
END SUBROUTINE swap

SUBROUTINE count
!routine to calculate quanta distribution
	IMPLICIT NONE
	INTEGER :: i,j,etemp,x,y,quanta

	DEALLOCATE (energy)
	x=SIZE(active,1)
	y=SIZE(active,2)

	!fin max no. of quanta in a position
	DO i=1,x
		DO j=1,y
			IF (active(i,j)>quanta) quanta=active(i,j)
		END DO
	END DO
	
	ALLOCATE (energy(0:quanta))

	!count quanta
	DO i=1,x
		DO j=1,y
			etemp=active(i,j)
			energy(etemp)=energy(etemp)+1
		END DO
	END DO

END SUBROUTINE count

SUBROUTINE graph
	IMPLICIT NONE
	!Define variables i for iterations, error for pgplot error status,
	!pgbeg to call pgplot, max_x and max_y for graph axis sizes
	INTEGER :: i, error, pgbeg, max_x, max_y
	!arrays for theoretical values, and x positions
	REAL, ALLOCATABLE, DIMENSION(:) :: theory
	INTEGER, ALLOCATABLE, DIMENSION(:) :: x
	
	CALL count
	
	max_y=0
	max_x=(SIZE(energy)-1)
	
	DO i=0,max_x
		IF (energy(i)>max_y) max_y=energy(i)
	END DO
	
	ALLOCATE (x(0:max_x))
	x=(/ (i, i=0,max_x) /)
	
	!Define device, and if an error in pgplot occurs, halt the program
	!NB /XW can be changed to /CPS to output a postscript file
	!(submitted with program)
	error = pgbeg(0,'/XW',1,1)
	IF (error /= 1) STOP

	!Define axes
	CALL pgenv(0.,FLOAT(max_x),0.,FLOAT(max_y),0,1)
	
	!Put title and axes labels on graph
	CALL pglab('Energy', 'No. of sites', 'Graph Of Sites per Energy')

	CALL pgsci(2)
	!Plot the graph
	CALL pgline(SIZE(energy),FLOAT(x),FLOAT(energy))
	
!-------------------------------------------------------------
!THEORETICAL GRAPH
	ALLOCATE (theory(0:max_x))
	!calculate theoretical
	DO i=0,max_x
		theory(i)=energy(0)*EXP(-LOG(2.0)*SIZE(active)*i/SUM(active))
	END DO
	CALL pgsci(3)
	!plot theoretcical
	CALL pgline(SIZE(energy),FLOAT(x),theory)
!-------------------------------------------------------------

	!Shutdown PGPlot properly
	CALL pgend
END SUBROUTINE graph

SUBROUTINE help
	!display a list of commands
	PRINT*,'-----------------------------------------------'
	PRINT*,'Solid - By Martin Eyles (9814937)'
	PRINT*
	PRINT*,'Command Summary'
	PRINT*,'X,x,Q,q   -   eXit/Quit'
	PRINT*,'H,h       -   display this Help'
	PRINT*,'A,a       -   change Active grid'
	PRINT*,'D,d       -   choose Dimensions of blocks'
	PRINT*,'E,e       -   find Equilibrium'
	PRINT*,'T,t       -   set Temperature of the grid'
	PRINT*,'G,g       -   plot Graph of quanta distribution'
	PRINT*,'N,n       -   perform random energy transfer'
	PRINT*,'              any Number of times'
	PRINT*,'S,s       -   show grid Status'
	PRINT*,'-----------------------------------------------'
END SUBROUTINE help
END MODULE subprogs

PROGRAM q1
!Main Program
	!allow acces to subroutines in module, and to global variables
	USE subprogs
	IMPLICIT NONE
	!Set up variable to store random seed
	INTEGER, DIMENSION(1) :: seed
	!Set up a keyboard input variable
	CHARACTER(LEN=1) :: user
	
	!set default grid size
	ALLOCATE (energy(0:100))
	ALLOCATE (grid3(100,50))
	grid1 => grid3(1:50,1:50)
	grid2 => grid3(51:100,1:50)
	active => grid1
	control=1

	!Initialize random number generator
	CALL SYSTEM_CLOCK(COUNT=seed(1))
	CALL RANDOM_SEED(PUT=seed)

	!Welcome The User
	PRINT*
	PRINT*,'***********************************'
	PRINT*,' Solid - By Martin Eyles (9814937)'
	PRINT*,' type h or H for Help'
	PRINT*,'***********************************'
	PRINT*

	DO
		!print cmd line style "solid", without new line
		WRITE(UNIT=*,FMT='(A8)',ADVANCE='NO')'Solid > '
		!take user input, and act acordingly
		READ*,user
		SELECT CASE(user)
			CASE('x','X','q','Q')
				!exit loop, which immediately exits th program
				EXIT
			CASE('h','H')
				CALL help
			CASE('a','A')
				CALL change_control
			!not documented - same as N or n then 1
			CASE('c','C')
				CALL swap
			CASE('d','D')
				CALL dimgrid
			CASE('e','E')
				CALL eqm
			CASE('g','G')
				CALL count
				CALL graph
			CASE('n','N')
				CALL advance
			CASE('s','S')
				CALL status
			CASE('t','T')
				CALL temp
			CASE('m','M')
				PRINT*,energy
			CASE DEFAULT
				!tell the user that the command isn't known
				!and tell them how to get help
				PRINT*,'Command not recognised'
				PRINT*,'type h or H for Help'
		END SELECT
	END DO

END PROGRAM q1