############################################################
#                                                          #
#     Virtual Laboratory of Geomathematics in Python       #
#                                                          #
#   Analysis of Frequency Distribution of Earthquakes:     # 
#          Gutenbeg-Richter's Law   (26.09.2017)           #
#                                                          #                
#         Complutense University of Madrid, Spain          #
#                                                          #
#   THIS SCRIPT IS PROVIDED BY THE AUTHORS "AS IS" AND     #
#   CAN BE USED BY ANYONE FOR THE PURPOSES OF EDUCATION    #
#   AND RESEARCH.                                          #
#                                                          #
############################################################
import sys

from matplotlib.widgets import Cursor
from matplotlib.gridspec import GridSpec
import matplotlib.pylab as plt
import numpy as np
from math import log10

import xlrd
import easygui as eg

# Ventana informacion
eg.msgbox("In this script:\n 1) We selected a file with earthquakes magnitude information\n 2) An histogram of frequency is plotted\n 3) The cumulative number (for >mag) is also plotted in log(N) vs mag format\n 4) The limits for the Gutenberg-Richter fit are selected using a cursor\n 5) The cumulative number and the line of fit are plotted together",title='Problem 2', ok_button='Continue')


# Eleccion de fichero
seleccion = ["*.xls","*.xlsx","*.ods"]
fichero=eg.fileopenbox(msg="Please choose a file",title="Data file",default="*",filetypes=seleccion)
print(fichero)

# Funcion de apertura y lectura de la columna n de un fichero excel
def read_excel(fichero,n):
        book = xlrd.open_workbook(fichero)
        sh = book.sheet_by_index(0)

        # Metemos la columna n en una lista
        coln=[]
        for rx in range(1,sh.nrows):
               coln.append(float(sh.cell_value(colx=n-1,rowx=rx)))
        return coln

# Si no se ha seleccionado el fichero se cierra el programa
if fichero == None:
     print("No file selected")
     sys.exit()

# Lee las magnitudes
num = eg.integerbox(msg='Column number for the magnitudes (1-100):',
                    title='Magnitudes',
                    default=13,
                    lowerbound=1,
                    upperbound=100,
                    image=None)
if num == None:
       print("No column selected")
       sys.exit()
mag=read_excel(fichero,num)

# Limites en el eje x de la grafica de magnitudes y numero de intervalos
magmin=int(min(mag))
magmax=int(max(mag))+1
etlim=['Mag min', 'Mag max', 'Number of bins']
xlimi=eg.multenterbox(msg='Data for the histogram', title='Magnitude histogram',fields=etlim,values=(magmin,magmax,20))
xmin=float(xlimi[0])
xmax=float(xlimi[1])
nbin=int(xlimi[2])
tbin=(xmax-xmin)/float(nbin)

# Almacenamos los datos del histograma
histograma, bines=np.histogram(mag,bins=nbin,range=(xmin,xmax))
posb=[]
for i in range(1,len(bines)):
         posb.append((bines[i]+bines[i-1])*0.5)

# Frecuencia acumulada
fac=[]
freini=0
for i in range(len(posb)-1,-1,-1):
    freini=freini+histograma[i]
    fac.append(freini)
prev=posb[::-1]

# Dibujamos graficas
fig=plt.figure('Magnitude')
fig.subplots_adjust(wspace=0.3)

# Histograma de las magnitudes
gs = GridSpec(2, 2, height_ratios=[10,1])
ax1=fig.add_subplot(gs[:,0])
ax1.tick_params(direction='in')
ax1.yaxis.set_ticks_position('both')
ax1.xaxis.set_ticks_position('both')
ax1.hist(mag,bins=bines)
ax1.set_title("    "+fichero)
ax1.set_ylabel("Frequency")
ax1.set_xlabel("Magnitude")

# Grafica de log de la frecuencia acumulada
ax2=fig.add_subplot(gs[1])
ax2.tick_params(direction='in')
ax2.yaxis.set_ticks_position('both')
ax2.xaxis.set_ticks_position('both')
ax2.semilogy(prev,fac,'bs')
ax2.set_ylabel("Cumulative frequency (>mag)")
ax2.set_xlabel("Magnitude")

# Mensaje de  seleccion
textstring='Selected in the above cumulative frequency plot\n the two limits for the fit'
axm=fig.add_subplot(gs[3])
axm.text(0,0,textstring,fontsize=6,bbox={'facecolor':'blue', 'alpha':0.2, 'pad':1})
axm.axis('off')

# Seleccionamos rango de magnitudes para el ajuste
coords = []
cursor = Cursor(ax2,color='r',lw=1,horizOn=False,vertOn=True)

def onclick_x(event):
        if event.xdata != None and event.ydata != None:
              print(event.xdata)

        global coords
        coords.append(event.xdata)

        if len(coords)==2:
              fig.canvas.mpl_disconnect(cid)
              plt.close('Magnitude')
        return

cid = fig.canvas.mpl_connect('button_press_event', onclick_x)


plt.show()
plt.close()

## Limites para el ajuste
amin=min(coords)
amax=max(coords)
majus=[]
fajus=[]
for i in range(len(prev)):
        if prev[i]>amin and prev[i]<amax:
              majus.append(prev[i])
              fajus.append(log10(fac[i]))
              
fit=np.polyfit(majus,fajus,1)
m=fit[0]
b=fit[1]
ajuste=np.poly1d(fit)
la="log(y)="+str(m)+"x+"+str(b)
co=np.corrcoef(majus,fajus)[0,1]
print(la)
print('r=',co)

# Dibuja los datos y el ajuste
fig2=plt.figure('Fit')
ax3=fig2.add_subplot(111)
ax3.tick_params(direction='in')
ax3.yaxis.set_ticks_position('both')
ax3.xaxis.set_ticks_position('both')
# Datos
ax3.semilogy(prev,fac,'bs',label='data')
# Ajuste
ax3.semilogy(majus,10**ajuste(majus),'r-',label='fit: log(y)=%5.3fx+%5.3f\n    r=%6.4f' % (m, b, co))
ax3.legend(loc=1)
ax3.set_ylabel("Cumulative frequency (>mag)")
ax3.set_xlabel("Magnitude")
ax3.set_title(fichero)

plt.show()