############################################################
#                                                          #
#       Virtual Laboratory of Geomathematics in Python     #
#                                                          #
#             Analysis of Earthquake Energy                # 
#                      (27.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

import xlrd
import easygui as eg

# Ventana informacion
eg.msgbox("Within this script:\n 1) We selected a file with earthquakes frequency and magnitude information\n 2) The frequency for each magnitude is plotted\n 3) The limits for the magnitude intervale are selected using a cursor\n 4) The cumulative energy is plotted",title='Problem 3', 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 minimum magnitudes (1-100):',
                    title='Minimum Magnitud',
                    default=1,
                    lowerbound=1,
                    upperbound=100)
if num == None:
       print("No column selected")
       sys.exit()
magmin=read_excel(fichero,num)

numm = eg.integerbox(msg='Column number for the maximum magnitudes (1-100):',
                    title='Maximum Magnitud',
                    default=2,
                    lowerbound=1,
                    upperbound=100)
if numm == None:
       print("No column selected")
       sys.exit()
magmax=read_excel(fichero,numm)

numf = eg.integerbox(msg='Column number for the frequency (1-100):',
                    title='Frequency',
                    default=3,
                    lowerbound=1,
                    upperbound=100)
if numm == None:
       print("No column selected")
       sys.exit()

frequ=read_excel(fichero,numf)

magmed=[]
for i in range(len(magmax)):
      magmed.append((magmax[i]+magmin[i])*0.5)
     
# Limites en el eje x de la grafica de magnitudes y numero de intervalos
etlim=['Mag min', 'Mag max', 'Frequencys']
xmin=min(magmin)
xmax=max(magmax)
fmax=max(frequ)+100

# Frecuencia acumulada
fac=[]
faca=0
for i in range(len(frequ)):
    fac.append(faca+frequ[i])
    faca=fac[i]

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

# Histograma de las magnitudes
gs = GridSpec(2, 1, 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.semilogy(magmed,frequ,'bs')
ax1.set_title("    "+fichero)
ax1.set_ylabel("Frequency")
ax1.set_xlabel("Magnitude")

# Mensaje de  seleccion
textstring='Selected in the above frequency plot the two limits for the energy calculation'
axm=fig.add_subplot(gs[1])
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(ax1,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)
msum=[]
esum=[]
etot=[]
eant=0
j=0
for i in range(len(magmed)):
        if magmed[i]>amin and magmed[i]<amax:
              msum.append(magmed[i])
              esum.append(frequ[i]*(10**(5.24+1.44*magmed[i])))
              etot.append(eant+esum[j])
              eant=etot[j]
              j=j+1

# Dibuja los datos y el ajuste 
fig2=plt.figure('Energy')
ax2=fig2.add_subplot(111)
ax2.tick_params(direction='in')
ax2.yaxis.set_ticks_position('both')
ax2.xaxis.set_ticks_position('both')
# Datos
ax2.semilogy(msum,esum,'bs',label='Energy (mag)')
# Total
ax2.semilogy(msum,etot,'r-',label='Total Energy (<mag)')
ax2.legend(loc=2)
ax2.set_ylabel("Energy (J)")
ax2.set_xlabel("Magnitude")
ax2.set_title(fichero)

plt.show()