############################################################
#                                                          #
#      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.                                          #
#                                                          #
############################################################
#!/usr/bin/env python
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

# Ventana informacion
print("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")


# Eleccion de fichero
fichero=raw_input('Please choose an excel file (including path)\n')
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 = int(input('Column number for the magnitudes (1-100):'))
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
xmin=int(min(mag))
xmax=int(max(mag))+1
nbin=20
ans1=raw_input('Data for the histogram:\n Min. mag.='+str(xmin)+'? (y/n) ')
if ans1 != 'y':
     xmin=int(input('New value for the minimum magnitude?\n'))
print('Min. mag.='+str(xmin))
ans2=raw_input('Max. mag.='+str(xmax)+'? (y/n) ')
if ans2 != 'y':
     xmax=int(input('New value for the maximum magnitude?\n'))
print('Max. mag.='+str(xmax))
ans3=raw_input('Number of bins=20? (y/n) ')
if ans3 != 'y':
     nbin=int(input('New value for the number of histogram bins?\n'))
print('Number of histogram bins='+str(nbin))

tbin=(xmax-xmin)/float(nbin)

# Almacenamos los datos del histograman
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()