!____________________________________________________________________________
! simple-mtc.f90 
!
! numerically solve the randomly placed & oriented tile problem
! by generating many randomly rotated + shifted unit-cell square grids &
! checking how often none of the grid nodes fall within a central square tile 
!____________________________________________________________________________


!___________________________________________________________________________
! Main program.  Does nmax = 100 mil experiments and displays statistics	
!___________________________________________________________________________
	use ifport           ! intel fortran lib, contains rand() function
	implicit none
	integer :: n, idum = 1, overl, nonover = 0, nmax = 1e8
	real :: randnum, rot, dx, dy, err, res, res_theor, err_sig
	res_theor = 2.-6./(atan(1.)*4)
	do n = 1, nmax
!	use either one of the functions: randnum(idum), or rand() 
!		rot = 1.;	        ! uncomment this for fixed orientation (pi/4)
		rot = rand()*8.  	! full angle range, rot is in units of pi/4
		dx = rand() -.5;	dy = rand() -.5 
!		rot = randnum(idum)*8.  ! alternative random number generator
!		dx = randnum(idum)-.5;	dy = randnum(idum)-.5 
!		both rand() and randnum() give similar results 
!		yet another random number generator is encl. belowi, but not tested 

!   compute overl = # of gridpoints overlapping with tile 
		call  gridpts (rot,dx,dy, overl)
		if(mod(n,nmax/100)==0) print*,'rot,dr,overl',rot,dx,dy,overl
		if (overl == 0) nonover = nonover+1 ! count non-overlapping tiles
	end do
	res = nonover/real(nmax)
	err = 2./nmax**0.5   ! assume 2/N**0.5 relative error 
	print*,' stats: ',nonover,' nonoverl. among ',nmax,' trials'
	print*,' fraction = ', res,'+-',err*res
	err_sig = (res/res_theor -1.)/err
	print*,'rel err =',err_sig,' sigma'
	print*,'[expected 0.0901407 for random angles]' ! just quote the solution
	print*,'[expected 0.1715728 for fixed angle 45 deg]'   ! do not reveal it
	end



	subroutine gridpts (alpha_pi4, dx,dy, overlaps)
!   first three numbers are random orientation and shift (of a square grid)
!   function gridpts returns overlaps = # of gridpoints overlapping with tile 
	implicit none
	real, parameter :: pi_4 = atan(1.)
	integer :: i, j, overlaps
	real :: alpha_pi4, dx,dy, x, y, xx, yy, c, s 
	overlaps = 0
	c = cos(alpha_pi4 *pi_4);		s = sin(alpha_pi4 *pi_4)
!	instead of randomly rotating and shifting the tile, we do this with grid 
	do i = -1,2
	 x = -0.5+i       		! -1.5..+1.5 is the grid x-interval
	 do j = -1,2  
		y = -0.5+j       	! -1.5..+1.5 is the grid y-interval
		xx = x*c - y*s + dx ! (xx,yy) = rotated & shifted grid point
		yy = x*s + y*c + dy	
!	tile is fixed and resides in |x|,|y| < 0.5 area   
		if (abs(xx) < 0.5 .and. abs(yy) < 0.5) overlaps = overlaps +1
	 end do
	end do
	end  



!-----------optional random number generators

    function randnum (idum)
    ! serial replacement for ran1.f from F77 Num Recipes
    !
    ! “Minimal” random number generator of Park and Miller combined 
    ! with a Marsaglia shift sequence. 
    ! Returns a uniform random deviate between 0.0 and 1.0 (exclusive 
    ! of the endpoint values). This fully portable, scalar generator has
    ! the “traditional” (not Fortran 90) calling sequence with a random 
    ! deviate as the returned function value: call with idum a negative 
    ! integer*8 to initialize;   thereafter, do not alter idum except to 
    ! reinitialize. The period of this  generator is about 3.1e18 .
        implicit none
        integer, intent(INOUT) :: idum
        integer, parameter :: IA=16807,IM=2147483647,IQ=127773,IR=2836
        real, save :: am
        integer, save :: ix=-1, iy=-1, k
        real :: randnum 
        if (idum <= 0 .or. iy < 0) then
           am=nearest(1.0,-1.0)/IM
           iy=ior(ieor(888889999,abs(idum)),1)
           ix=ieor(777755555,abs(idum))
           idum=abs(idum)+1
        end if
        ix=ieor(ix,ishft(ix,13))
        ix=ieor(ix,ishft(ix,-17))
        ix=ieor(ix,ishft(ix,5))
        k=iy/IQ
        iy=IA*(iy-k*IQ)-IR*k
        if (iy < 0) iy=iy+IM
        randnum = am * ior(iand(IM,ieor(ix,iy)),1)        
     end function 



! another optional random number generator with two int seeds ix & iy
! they may need save attribute in the calling program declaration
! 1st call initializes the sequence, use ix=iy=0.
! Taken from Hoover's "Comutational Statistical Mechanics" 
! p.91 (pdf) (or 81 marked in book)
!
    function rand2 (ix,iy)
	i = 1029*ix + 1731
	j = i + 1029*iy + 507*ix - 1731
	ix = mod(i,2048)
	j = j +(i-ix)/2048
	iy = mod(j,2048)
	rand2 = (ix+2048*iy)/4194304.
    end