#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Arbuckle injection around 1 point
"""

#Run with pandas v1.3.4, numpy v1.20.3, and matplotlib v3.4.3
from matplotlib import rcParams
rcParams['font.family'] = 'sans-serif'
rcParams['font.sans-serif'] = ['Times New Roman']
rcParams.update({'font.size': 12})
import matplotlib.pyplot as plt
import pandas as pd
#import fiona
import numpy as np

import datetime


def haversine(lat1, lon1, lat2, lon2):
    """
    Calculate the great circle distance between two points
    on the earth (specified in decimal degrees)

    All args must be of equal length.

    """
    lon1, lat1, lon2, lat2 = map(np.radians, [lon1, lat1, lon2, lat2])

    dlon = lon2 - lon1
    dlat = lat2 - lat1

    a = np.sin(dlat/2.0)**2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon/2.0)**2

    c = 2 * np.arcsin(np.sqrt(a))
    km = 6371 * c
    return km


injdat = pd.read_hdf('welldata.h5', '/injdat')

well_loc={
    1: (36.8711390000000,-98.3715000000000),
    2: (36.8539080000000,	-98.2899310000000),
    3:      (36.8124170000000,	-98.2954720000000),
    4:      (36.7839440000000,-98.3673060000000),
    5:      (36.7394170000000,	-98.0364440000000),
    6:      (36.6376670000000,	-97.9917500000000),
    7:    (36.0007380000000,	-96.8334643000000),
    8:    (35.9702401000000,	-96.8385227000000),
    9:    (36.0007130000000,	-97.1410830000000),
    10:    (35.6655126100000,	-96.7152025900000),
    11:    (36.2064570800000,	-96.4556269600000),
    12:    (35.7858900000000,	-97.3868690000000),
    13:    (36.5211111000000,	-97.0395555000000),
    14:    (36.5508480000000,	-97.0476170000000),
    15:    (36.2756923000000,	-97.5153323000000),
    }

for i in range(1,16):
    lat1, lon1 = well_loc[i]
    
    
    distvar = [200,5000,10000,15000] #radii from point, meters
    
    mydate = datetime.datetime.now()
    current = datetime.datetime(mydate.year, mydate.month, 1)
    next_month = datetime.datetime(mydate.year + int(mydate.month / 12), ((mydate.month % 12) + 1), 1)
    enddate = str(next_month.month)+'/1/'+str(next_month.year)
    
    fig=plt.figure(figsize=[10,3],linewidth=1.0)
    possibles = np.zeros((len(injdat['Longitude']),1))
    for dist in distvar:
        sums_up_all = []
        possibles = 1000*haversine(lat1,lon1,injdat['Latitude'].values,injdat['Longitude'].values)
        subinjdat = injdat[(possibles<dist)]
        subinjdat = subinjdat.reset_index(drop=True)
        print(len(subinjdat))
    
    
        timestep = pd.date_range(start='8/4/2016', end='3/4/2020', freq='2d')
        sums_up_all = pd.DataFrame()
        for f in np.arange(1,len(timestep)):
            liqp = np.nansum(subinjdat[(subinjdat['Report_Date']>timestep[f-1]) & (subinjdat['Report_Date']<timestep[f])]['Volume_BPD'])
            liqp2 = np.nansum(subinjdat[(subinjdat['Report_Date']>timestep[f-1]) & (subinjdat['Report_Date']<timestep[f])]['Pressure_PSI'])
            liqp2 = liqp2/len(subinjdat[(subinjdat['Report_Date']>timestep[f-1]) & (subinjdat['Report_Date']<timestep[f])])
            dfall = pd.DataFrame({'Date':timestep[f],'Volume':liqp,'PSI':liqp2},index=[f])
            #print(dfall)
            sums_up_all = sums_up_all.append(dfall)
        #plt.step(sums_up_all['Date'],sums_up_all['Volume']/1e6,label=str(dist)+' m')
        with open('W'+str(i)+'-'+str(dist)+' m-M.txt', 'w') as f:
            for strlin in range(1,len(sums_up_all)):
                f.write(str(sums_up_all.Date[strlin])+', '+str(sums_up_all.Volume[strlin])+', '+str(sums_up_all.PSI[strlin])+'\n')
    #plt.xlabel('Date')
    #plt.ylabel('SWD (MBa)')
    #plt.xlim([datetime.date(2016, 8, 1), datetime.date(2020, 4, 1)])
    #plt.legend()
    #plt.savefig('W'+str(i)+'smrange.png')
    #plt.show()
#plt.legend()
#plt.xlim([datetime.date(2016, 8, 4), datetime.date(2020, 3, 4)])