%matplotlib inline
from salishsea_tools import tidetools
import netCDF4 as NC
import matplotlib.pyplot as plt
import numpy as np

runnames = ['50s_15-21Sep','50s_22-25Sep','50s_26-29Sep','50s_30Sep-6Oct','50s_7-13Oct']
runlength = np.array([7.0,4.0,4.0,7.0,7.0])

mod_M2_eta_real1 = 0.0
mod_M2_eta_imag1 = 0.0
mod_K1_eta_real1 = 0.0
mod_K1_eta_imag1 = 0.0

for runnum in range(0,len(runnames)):
    harmT = NC.Dataset('/data/dlatorne/MEOPAR/SalishSea/results/'+runnames[runnum]+'/Tidal_Harmonics_eta.nc','r')
    #get imaginary and real components
    mod_M2_eta_real1 = mod_M2_eta_real1 + harmT.variables['M2_eta_real'][0,:,:]*runlength[runnum]
    mod_M2_eta_imag1 = mod_M2_eta_imag1 + harmT.variables['M2_eta_imag'][0,:,:]*runlength[runnum]
    mod_K1_eta_real1 = mod_K1_eta_real1 + harmT.variables['K1_eta_real'][0,:,:]*runlength[runnum]
    mod_K1_eta_imag1 = mod_K1_eta_imag1 + harmT.variables['K1_eta_imag'][0,:,:]*runlength[runnum]

totaldays = sum(runlength)
mod_M2_eta_real = mod_M2_eta_real1/totaldays
mod_M2_eta_imag = mod_M2_eta_imag1/totaldays
mod_K1_eta_real = mod_K1_eta_real1/totaldays
mod_K1_eta_imag = mod_K1_eta_imag1/totaldays

mod_M2_amp = np.sqrt(mod_M2_eta_real**2+mod_M2_eta_imag**2)
mod_M2_pha = -np.degrees(np.arctan2(mod_M2_eta_imag,mod_M2_eta_real))
mod_K1_amp = np.sqrt(mod_K1_eta_real**2+mod_K1_eta_imag**2)
mod_K1_pha = -np.degrees(np.arctan2(mod_K1_eta_imag,mod_K1_eta_real))

bathy, X, Y = tidetools.get_SS_bathy_data()
tidetools.plot_amp_map(X,Y,mod_M2_amp,'Combined 50s runs',0,'M2')
tidetools.plot_pha_map(X,Y,mod_M2_pha,'Combined 50s runs',0,'M2')
tidetools.plot_amp_map(X,Y,mod_K1_amp,'Combined 50s runs',0,'K1')
tidetools.plot_pha_map(X,Y,mod_K1_pha,'Combined 50s runs',0,'K1')