from visual import *
from visual.graph import *
from random import *

SITES=20
MAXQUAKE=SITES*SITES
FCRIT=4.0
TOTAL_TIME=4000
DELTA_F=0.01

pic1=gdisplay(x=0,y=0,width=400,height=400,title="Earthquake size distribution",
             xtitle="log(size)",ytitle="-log(N(size))",xmax=4.0,xmin=0,ymin=0,
             ymax=15.0,
             foreground=color.black, background=color.white)
plot1=gdots(gdisplay=pic1,color=color.blue)


hdist2 = gdisplay(x=400, y=0, ymax = 0.5, xmin=0,xmax= MAXQUAKE/2,
             width=400, height=400, title="Relative Frequency of Quakes",
             xtitle="quake_size", ytitle="relative frequency")
hdist2.display.userspin=1

observation2 = gvbars(color=color.red)
for i in arange(0,MAXQUAKE+1):
    observation2.plot(pos=(i,1e-1))

force=zeros((SITES+1,SITES+1),Float)
move=zeros((SITES,SITES))
quake_size=zeros( (MAXQUAKE+1) )


def check(height,move):
    unstable=false
    for i in range (0,SITES):
        for j in range (0,SITES):
             if (force[i,j] > FCRIT):
                move[i,j]=true
                unstable=true
             else:
                move[i,j]=false
    return(unstable)

def relax(force,move):
    for i in range (0,SITES):
        for j in range (0,SITES):
            if (move[i,j]==true):
                force[i,j] -= FCRIT
                if(i<SITES-1):
                    force[i+1,j] += 1
                if (i>0):
                    force[i-1,j] += 1
                if (j<SITES-1):
                    force[i,j+1] += 1
                if (j>0):
                    force[i,j-1] += 1
                    
## free boundary conditions: force at boundary is always zero

def initial(force,quake_size):
    for i in range (0,SITES):  
         for j in range (0,SITES):
            force[i,j]=2.0*random()
            
    for i in range (0,SITES):
        force[i,SITES]=0
        force[SITES,i]=0
        
    for i in range (0,MAXQUAKE):
        quake_size[i]=0
    

initial(force,quake_size)
time=0
total_quakes=0
max_quake_size=1

started=0
while(time<TOTAL_TIME):
    quake=0
    for i in range (0,SITES):
        for j in range (0, SITES):
            force[i,j] += DELTA_F
            
    time += 1
    unstable=true
    while (unstable):
        unstable=check(force,move)
        if(unstable):
            relax(force,move)
            quake += 1
            started=1

    

    if (started):
    #emulate hidden functions for gvbars
        for s in arange(MAXQUAKE+1):
            val=float(quake_size[s])/float(max(max(quake_size),1))
            observation2.vbars[s].pos[1]=val/2.
            observation2.vbars[s].height=val
            observation2.gdisplay.resize(observation2.vbars[s].pos[0],observation2.vbars[s].pos[1])    

        if (quake>=max_quake_size):
           max_quake_size=quake
        
        quake_size[quake] += 1

        total_quakes += quake
    

print "total number of quakes",total_quakes
print "biggest quake",max_quake_size

for i in range(1,max_quake_size):
    tmp=float(quake_size[i])
    if(tmp>0.0):
      plot1.plot(pos=(log(float(i)),-log(tmp/float(time))))
        

print quake_size