name = 'GmO_TS13'
if name == 'GmO_Base':
    timeinc = 3659558400
elif name == 'GmO_TS13':
    timeinc = 3259872000
fT1 = NC.Dataset('../../myResults/NEMO36_Tides/'+name+'/PointAtkinson.nc','r')
ssh = fT1.variables['sossheig'][:]
time = (fT1.variables["time_counter"][:]-timeinc)/3600. 
print (ssh.shape)
plt.plot(time,ssh[:,0,0])

plt.plot(bathy_lon[415:418,281:285], bathy_lat[415:418, 281:285],'*')
plt.plot(bathy_lon[416, 283], bathy_lat[416, 283], 'ks')
plt.plot(-123.3176, 49.0422, 'ro')

east_lat = 49+2.531/60.
east_lon = -(123 + 19.059/60.)

print (east_lat, east_lon)

for i in [416,417]:
    for j in [282,283]:
        print (i, j, tidetools.haversine(east_lon, east_lat, bathy_lon[i,j], bathy_lat[i,j]))

print (lat[0,0], lon[0,0], bathy_lat[380, 240], bathy_lon[380, 240])

ieast = 35; jeast = 42
icent = 43; jcent = 25
cent_lat = 49.0401; cent_lon = -123.4261
plt.plot(v[:, 0, ieast, jeast])
plt.plot(v[:, 5, ieast, jeast])

print (np.mean(v[:, 0, ieast, jeast]), np.mean(v[:, 5, ieast, jeast]))

timeplus = np.zeros([5,960])
P1K1_amp = 0.31; P1K1_pha = -3
K2S2_amp = 0.27; K2S2_pha = 0
v_M2pha = np.zeros((sCO, sCO, 40)); v_M2amp = np.zeros((sCO, sCO, 40))
u_M2amp = np.zeros((sCO, sCO, 40)); u_M2pha = np.zeros((sCO, sCO, 40))
print (time[0:-1:2].shape)
timeplus[0] = time[0:-1:2]
timeplus[1] = P1K1_amp*np.ones_like(timeplus[0]); timeplus[2] = P1K1_pha*np.ones_like(timeplus[0])
timeplus[3] = K2S2_amp*np.ones_like(timeplus[0]); timeplus[4] = K2S2_pha*np.ones_like(timeplus[0])

vavg = np.zeros((960, sCO, sCO))
vavg_M2amp = np.zeros((sCO, sCO)); vavg_M2pha = np.zeros((sCO, sCO))
for ii in iichoice:
    for jj in jjchoice:
        m = v[:,:,ii,jj] == 0
        v_masked = np.ma.array(v[:,:,ii,jj], mask=m)
        vavg[:,ii,jj] = analyze.depth_average(v_masked, depths, depth_axis=1)
        fitted, cov = curve_fit(sixplustwo, timeplus, vavg[:,ii,jj]) 
        vavg_M2amp[ii,jj] = fitted[0]
        vavg_M2pha[ii,jj] = fitted[1]

uavg = np.zeros((960, sCO, sCO))
uavg_M2amp = np.zeros((sCO, sCO)); uavg_M2pha = np.zeros((sCO, sCO))
for ii in iichoice:
    for jj in jjchoice:
        m = u[:,:,ii,jj] == 0
        u_masked = np.ma.array(u[:,:,ii,jj], mask=m)
        uavg[:,ii,jj] = analyze.depth_average(u_masked, depths, depth_axis=1)
        fitted, cov = curve_fit(sixplustwo, timeplus, uavg[:,ii,jj]) 
        uavg_M2amp[ii,jj] = fitted[0]
        uavg_M2pha[ii,jj] = fitted[1]

for ii in iichoice:
    print (ii)
    for jj in jjchoice:
        for i in ddchoice:
            if np.abs(v[5,i,ii,jj]) > 0.00000001:
                fitted, cov = curve_fit(sixplustwo, timeplus, v[:,i,ii,jj]) 
                v_M2amp[ii,jj,i] = fitted[0]
                v_M2pha[ii,jj,i] = fitted[1]

plt.plot(v_M2amp[ieast, jeast], -depths[:],'o-')
plt.plot(v_M2amp[ieast, jeast-5], -depths[:],'o-')

for ii in iichoice:
    print (ii)
    for jj in jjchoice:
        for i in ddchoice:
            if np.abs(u[8,i,ii,jj]) > 0.00000001:
                fitted, cov = curve_fit(sixplustwo, timeplus, u[:,i,ii,jj]) 
                u_M2amp[ii,jj,i] = fitted[0]
                u_M2pha[ii,jj,i] = fitted[1]

plt.plot(u_M2amp[ieast, jeast], -depths[:],'o-')
plt.plot(u_M2amp[ieast, jeast-5], -depths[:],'o-')

scale = 0.4
print (major.shape)
grid_B = '../../NEMO-forcing/grid/bathy_downonegrid.nc'
fig, ax = plt.subplots(1, 1, figsize=(10,10))
viz_tools.set_aspect(ax) 
ii = ieast; jj = jeast; inode = 1
for ii, jj, inode in zip((ieast,icent), (jeast, jcent), (1,0)):
    ellsc = Ellipse(xy=(lon[ii,jj], lat[ii,jj]),
                width=scale * soon[inode][0],
                height=scale * soon[inode][1],
                angle= soon[inode][2],
                color='grey')
    ax.add_artist(ellsc)
    ellsc = Ellipse(xy=(lon[ii,jj], lat[ii,jj]),
                width=scale * node[inode][0],
                height=scale * node[inode][1],
                angle= node[inode][2],
                color='green')
    ax.add_artist(ellsc)
    ellsc = Ellipse(xy=(lon[ii,jj], lat[ii,jj]),
                width=scale * major[ii,jj],
                height=scale * minor[ii,jj],
                angle=theta[ii,jj]+29,
                color='red')
   
    ax.add_artist(ellsc)
viz_tools.plot_land_mask(ax, grid_B, coords='map')
viz_tools.plot_coastline(ax, grid_B, coords='map', isobath=5)
viz_tools.plot_coastline(ax, grid_B, coords='map', isobath=10)
viz_tools.plot_coastline(ax, grid_B, coords='map', isobath=20)

ax.set_xlim(-123.7, -123.)
ax.set_xticks((-123.6, -123.4, -123.2, -123.))
ax.set_xticklabels(('-123.6', '-123.4', '-123.2', '-123.'))
ax.set_yticks((48.9, 49, 49.1, 49.2, 49.3))
for tick in ax.xaxis.get_major_ticks():
                tick.label.set_fontsize(14)
for tick in ax.yaxis.get_major_ticks():
                tick.label.set_fontsize(14)
ax.set_xlabel("Longitude", fontsize=20)
ax.set_ylabel("Latitude", fontsize=20)

ax.set_ylim(48.89, 49.3)
ellsc = Ellipse(xy=(-123.13, 49.2),
                width=scale * 0.05,
                height=scale * 0.2,
                angle=29.,
                color='grey')
   
ax.add_artist(ellsc)
ax.text(-123.1, 49.2,'20 cm s$^{-1}$', fontsize=14, color='grey')
ax.text(-123.65, 49.28, 'New Model', fontsize=14, color='red')
ax.text(-123.65, 49.26, 'Observations', fontsize=14, color='green')
ax.text(-123.65, 49.24, 'Production Model', fontsize=14, color='grey')
ax.text(-123.3, 48.92, 'Bathymetry: 5 m, 10 m, 20 m')
ax.set_title('Barotropic Tidal Ellipses', fontsize=20)

scale = 0.2
print (major.shape)
grid_B = '../../NEMO-forcing/grid/bathy_downonegrid.nc'
fig, ax = plt.subplots(1, 1, figsize=(10,10))
viz_tools.set_aspect(ax) 
for ii in [20, 27, 35, 43, 50]:
    for jj in jjchoice:
        id = 0
        ellsc = Ellipse(
                xy=(lon[ii,jj], lat[ii,jj]),
                width=scale * major[ii,jj,id],
                height=scale * minor[ii,jj,id],
                angle=theta[ii,jj,id]+29,
                color='blue')
        ax.add_artist(ellsc)
        ellsc = Ellipse(
                xy=(lon[ii,jj], lat[ii,jj]),
                width=scale * major_avg[ii,jj],
                height=scale * minor_avg[ii,jj],
                angle=theta_avg[ii,jj]+29,
                color='red')
#        ax.add_artist(ellsc)
        
viz_tools.plot_land_mask(ax, grid_B, coords='map')
viz_tools.plot_coastline(ax, grid_B, coords='map', isobath=5)
viz_tools.plot_coastline(ax, grid_B, coords='map', isobath=10)
viz_tools.plot_coastline(ax, grid_B, coords='map', isobath=20)

ax.set_xlim(-123.7, -123.)
ax.set_xticks((-123.6, -123.4, -123.2, -123.))
ax.set_xticklabels(('-123.6', '-123.4', '-123.2', '-123.'))
ax.set_yticks((48.9, 49, 49.1, 49.2, 49.3))
for tick in ax.xaxis.get_major_ticks():
                tick.label.set_fontsize(14)
for tick in ax.yaxis.get_major_ticks():
                tick.label.set_fontsize(14)
ax.set_xlabel("Longitude", fontsize=20)
ax.set_ylabel("Latitude", fontsize=20)

ax.set_ylim(48.89, 49.3)
ellsc = Ellipse(xy=(-123.13, 49.2),
                width=scale * 0.05,
                height=scale * 0.2,
                angle=29.,
                color='grey')
   
ax.add_artist(ellsc)
ax.text(-123.1, 49.2,'20 cm s$^{-1}$', fontsize=14, color='grey')
#ax.text(-123.65, 49.28, 'Barotropic Tides', fontsize=14, color='red')
#ax.text(-123.65, 49.26, 'Surface Tides', fontsize=14, color='blue')
ax.text(-123.2, 48.91, 'Bathymetry: 5 m, 10 m, 20 m')
ax.set_title('30-m depth Tidal Ellipses', fontsize=20)

scale = 0.2
print (major.shape)
grid_B = '../../NEMO-forcing/grid/bathy_downonegrid.nc'
fig, ax = plt.subplots(1, 1, figsize=(15,15))
viz_tools.set_aspect(ax) 
id=0
cmap = plt.get_cmap('plasma')
cmap.set_bad('white')

mesh=ax.pcolormesh(lon[:-1,:-1], lat[:-1,:-1], CX[:,:,id], cmap=cmap, vmin=0, vmax=0.5)
        
viz_tools.plot_land_mask(ax, grid_B, coords='map')
viz_tools.plot_coastline(ax, grid_B, coords='map', isobath=5)
viz_tools.plot_coastline(ax, grid_B, coords='map', isobath=10)
viz_tools.plot_coastline(ax, grid_B, coords='map', isobath=20)

ax.set_xlim(-123.7, -123.)
ax.set_xticks((-123.6, -123.4, -123.2, -123.))
ax.set_xticklabels(('-123.6', '-123.4', '-123.2', '-123.'))
ax.set_yticks((48.9, 49, 49.1, 49.2, 49.3))
for tick in ax.xaxis.get_major_ticks():
                tick.label.set_fontsize(14)
for tick in ax.yaxis.get_major_ticks():
                tick.label.set_fontsize(14)
ax.set_xlabel("Longitude", fontsize=20)
ax.set_ylabel("Latitude", fontsize=20)

ax.set_ylim(48.89, 49.3)
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad="5%")

fig.colorbar(mesh, cax=cax)


#ax.text(-123.1, 49.2,'20 cm s$^{-1}$', fontsize=14, color='grey')
#ax.text(-123.65, 49.28, 'Barotropic Tides', fontsize=14, color='red')
#ax.text(-123.65, 49.26, 'Surface Tides', fontsize=14, color='blue')
ax.text(-123.16, 48.902, 'Bathymetry: 5 m, 10 m, 20 m')
ax.set_title('Cosine of Across Strait M2 Velocity', fontsize=25)

fig, axs = plt.subplots(1,2,figsize=(10,7))
ileft = 20
iright = 30
axs[0].plot(major[ileft, 20], -depths[:], 'o-', label="Farthest East", color='b', alpha=0.33)
for jj in [30, 40, 50]:
    print ((jj-20)/60.+0.33)
    axs[0].plot(major[ileft, jj], -depths[:], 'o-', color='b', alpha=(jj-20)/60.+0.33)
axs[0].plot(major[ileft, 60], -depths[:], 'o-', label="Farthest West", color='b', alpha=1)

axs[1].plot(major[iright, 10], -depths[:], 'o-', label="Farthest East", color='b', alpha=0.33)
for jj in [20, 30, 40, 50]:
    axs[1].plot(major[iright, jj], -depths[:], 'o-', color='b', alpha=(jj-20)/70.+0.33)
axs[1].plot(major[iright, 60], -depths[:], 'o-', label="Farthest West", color='b', alpha=1)
    

axs[0].set_ylim(-250,0)
axs[0].legend(loc='lower right')
axs[1].set_ylim(-80,0)
axs[1].legend(loc='lower right')
axs[0].set_title("Major Axis, Robert's Bank Line")
axs[1].set_title("Major Axis, Westham Island Line")

m = major == 0
major_masked = np.ma.array(major, mask=m)
mm = minor == 0
minor_masked = np.ma.array(minor, mask=mm)
mt = theta == 0
theta_masked = np.ma.array(theta, mask=mt)
fig, ax =plt.subplots(1,1,figsize=(10,10))
cmap = plt.get_cmap('plasma')
cmap.set_bad('burlywood')
mesh=ax.pcolormesh(lon[20], -depths[:], major_masked[20].transpose(), cmap=cmap, vmax=0.4, vmin=0.25)
ax.set_ylim((-20,0))
ax.set_xlim(-123.53, -123.1)
#ax.set_xlim(-123.45, -123.09)
#ax.set_xticks((-123.4, -123.3, -123.2, -123.1))
#ax.set_xticklabels(('-123.4','-123.3', '-123.2', '-123.1'))
for tick in ax.xaxis.get_major_ticks():
                tick.label.set_fontsize(14)
for tick in ax.yaxis.get_major_ticks():
                tick.label.set_fontsize(14)
ax.set_xlabel("Longitude", fontsize=20)
ax.set_ylabel("Depth (m)", fontsize=20)
cbar = fig.colorbar(mesh)
cbar.ax.tick_params(labelsize=14) 
cbar.set_label("Major Axis (m s$^{-1}$)", fontsize=20)
ax.set_title("Robert's Bank Line", fontsize=25)

m = major == 0
major_masked = np.ma.array(major,mask=m)
fig, ax =plt.subplots(1,1,figsize=(10,10))
cmap = plt.get_cmap('plasma')
cmap.set_bad('burlywood')
mesh=ax.pcolormesh(lon[27], -depths[:], major_masked[27].transpose(), cmap=cmap, vmax=0.4, vmin=0.25)
ax.set_ylim((-20,0))
ax.set_xlim(-123.53, -123.16)
#ax.set_xlim(-123.45, -123.09)
#ax.set_xticks((-123.4, -123.3, -123.2, -123.1))
#ax.set_xticklabels(('-123.4','-123.3', '-123.2', '-123.1'))
for tick in ax.xaxis.get_major_ticks():
                tick.label.set_fontsize(14)
for tick in ax.yaxis.get_major_ticks():
                tick.label.set_fontsize(14)
ax.set_xlabel("Longitude", fontsize=20)
ax.set_ylabel("Depth (m)", fontsize=20)
cbar = fig.colorbar(mesh)
cbar.ax.tick_params(labelsize=14) 
cbar.set_label("Major Axis (m s$^{-1}$)", fontsize=20)
ax.set_title("Middle Line", fontsize=25)

m = major == 0
major_masked = np.ma.array(major,mask=m)
fig, ax =plt.subplots(1,1,figsize=(10,10))
cmap = plt.get_cmap('plasma')
cmap.set_bad('burlywood')
mesh=ax.pcolormesh(lon[35], -depths[:]+0.5, CX[35].transpose(), cmap=cmap, vmax=0.35, vmin=0.)
ax.set_ylim((-50,0))
ax.set_xlim(-123.53, -123.16)
ax.set_xticks((-123.4, -123.3, -123.2))
ax.set_xticklabels(('-123.4','-123.3', '-123.2'))
for tick in ax.xaxis.get_major_ticks():
                tick.label.set_fontsize(14)
for tick in ax.yaxis.get_major_ticks():
                tick.label.set_fontsize(14)
ax.set_xlabel("Longitude", fontsize=20)
ax.set_ylabel("Depth (m)", fontsize=20)
cbar = fig.colorbar(mesh)
cbar.ax.tick_params(labelsize=14) 
cbar.set_label("Major Axis (m s$^{-1}$)", fontsize=20)
ax.set_title("Major Axis across Westham Island Line", fontsize=25)
ax.plot(east_lon,-5,'wo')

m = major == 0
major_masked = np.ma.array(major,mask=m)
p = phase == 0
phase_masked = np.ma.array(phase,mask=p)
fig, ax =plt.subplots(1,1,figsize=(10,10))
cmap = plt.get_cmap('plasma')
cmap.set_bad('burlywood')
mesh=ax.pcolormesh(lon[icent], -depths[:]+0.5, phase_masked[icent].transpose(), cmap=cmap)
ax.set_ylim((-20,0))
ax.set_xlim(-123.53, -123.16)
ax.set_xticks((-123.4, -123.3, -123.2))
ax.set_xticklabels(('-123.4','-123.3', '-123.2'))
for tick in ax.xaxis.get_major_ticks():
                tick.label.set_fontsize(14)
for tick in ax.yaxis.get_major_ticks():
                tick.label.set_fontsize(14)
ax.set_xlabel("Longitude", fontsize=20)
ax.set_ylabel("Depth (m)", fontsize=20)
cbar = fig.colorbar(mesh)
cbar.ax.tick_params(labelsize=14) 
cbar.set_label("Major Axis (m s$^{-1}$)", fontsize=20)
ax.set_title("Central Line", fontsize=25)

m = major == 0
major_masked = np.ma.array(major,mask=m)
p = phase == 0
phase_masked = np.ma.array(phase,mask=p)
fig, ax =plt.subplots(2,1,figsize=(10,15))
cmap = plt.get_cmap('plasma')
cmap.set_bad('burlywood')

for ia in [0,1]:
    if ia==0:
        mesh=ax[0].pcolormesh(lon[50], -depths[:]+0.5, CX[50].transpose(), cmap=cmap, vmax=0.35, vmin=0)
    else:
        mesh=ax[1].pcolormesh(lon[50], -depths[:]+0.5, phase_masked[50].transpose(), cmap=cmap, vmin=120)

    ax[ia].set_ylim((-20,0))
    ax[ia].set_xlim(-123.53, -123.16)
    ax[ia].set_xticks((-123.4, -123.3, -123.2))
    ax[ia].set_xticklabels(('-123.4','-123.3', '-123.2'))
    axe = ax[ia]
    for tick in axe.xaxis.get_major_ticks():
                tick.label.set_fontsize(14)
    for tick in axe.yaxis.get_major_ticks():
                tick.label.set_fontsize(14)
    ax[1].set_xlabel("Longitude", fontsize=20)
    ax[ia].set_ylabel("Depth (m)", fontsize=20)
    cbar = fig.colorbar(mesh, ax=ax[ia])
    cbar.ax.tick_params(labelsize=14) 
    cbar.set_label("Major Axis (m s$^{-1}$)", fontsize=20)
ax[0].set_title("Northern Line", fontsize=25)

m = major == 0
major_masked = np.ma.array(major,mask=m)
p = phase == 0
phase_masked = np.ma.array(phase,mask=p)
fig, ax =plt.subplots(2,1,figsize=(10,15))
cmap = plt.get_cmap('viridis')
cmap.set_bad('burlywood')

for ia in [0,1]:
    if ia==0:
        mesh=ax[0].pcolormesh(lon[57], -depths[:]+0.5, major_masked[57].transpose(), cmap=cmap, vmax=0.25)
    else:
        mesh=ax[1].pcolormesh(lon[57], -depths[:]+0.5, phase_masked[57].transpose(), cmap=cmap, vmin=120)

    ax[ia].set_ylim((-20,0))
    ax[ia].set_xlim(-123.53, -123.16)
    ax[ia].set_xticks((-123.4, -123.3, -123.2))
    ax[ia].set_xticklabels(('-123.4','-123.3', '-123.2'))
    axe = ax[ia]
    for tick in axe.xaxis.get_major_ticks():
                tick.label.set_fontsize(14)
    for tick in axe.yaxis.get_major_ticks():
                tick.label.set_fontsize(14)
    ax[1].set_xlabel("Longitude", fontsize=20)
    ax[ia].set_ylabel("Depth (m)", fontsize=20)
    cbar = fig.colorbar(mesh, ax=ax[ia])
    cbar.ax.tick_params(labelsize=14) 
    cbar.set_label("Major Axis (m s$^{-1}$)", fontsize=20)
ax[0].set_title("Further Northern Line", fontsize=25)

m = major == 0
major_masked = np.ma.array(major,mask=m)
p = phase == 0
phase_masked = np.ma.array(phase,mask=p)
fig, ax =plt.subplots(2,1,figsize=(10,15))
cmap = plt.get_cmap('viridis')
cmap.set_bad('burlywood')

for ia in [0,1]:
    if ia==0:
        mesh=ax[0].pcolormesh(lon[65], -depths[:]+0.5, major_masked[65].transpose(), cmap=cmap, vmax=0.25)
    else:
        mesh=ax[1].pcolormesh(lon[65], -depths[:]+0.5, phase_masked[65].transpose(), cmap=cmap, vmin=120)

    ax[ia].set_ylim((-20,0))
    ax[ia].set_xlim(-123.53, -123.16)
    ax[ia].set_xticks((-123.4, -123.3, -123.2))
    ax[ia].set_xticklabels(('-123.4','-123.3', '-123.2'))
    axe = ax[ia]
    for tick in axe.xaxis.get_major_ticks():
                tick.label.set_fontsize(14)
    for tick in axe.yaxis.get_major_ticks():
                tick.label.set_fontsize(14)
    ax[1].set_xlabel("Longitude", fontsize=20)
    ax[ia].set_ylabel("Depth (m)", fontsize=20)
    cbar = fig.colorbar(mesh, ax=ax[ia])
    cbar.ax.tick_params(labelsize=14) 
    cbar.set_label("Major Axis (m s$^{-1}$)", fontsize=20)
ax[0].set_title("Far Northern Line", fontsize=25)

print (depths[21])
base_CX, base_SX, base_CY, base_SY, base_ap, base_am, base_ep, base_em, base_major, base_minor, base_theta, base_phase, base_u_M2amp, base_u_M2pha, base_v_M2amp, base_v_M2pha = (
CX, SX, CY, SY, ap, am, ep, em, major, minor, theta, phase, u_M2amp, u_M2pha, v_M2amp, v_M2pha)

fig, axs = plt.subplots(1,2,figsize=(10,7))
for ii in [43]:
    for jj in jjchoice:
        axs[0].plot(major[ii, jj], -depths[:], 'o-', label=jj)
axs[0].set_ylim(-250,0)
axs[0].legend(loc='lower right')
for ii in [35]:
    for jj in jjchoice:
        axs[1].plot(major[ii, jj], -depths[:], 'o-', label=jj)
axs[1].set_ylim(-250,0)
axs[1].legend(loc='lower right')

fig, axs = plt.subplots(2,2,figsize=(15,15))
axs[0,0].plot(east[:-10,1],-east[:-10,0], color='grey', linewidth=3)
axs[0,0].plot(eastobs[1]/100., -eastobs[0], 'g', linewidth=3)
axs[0,0].plot(codar[1][0],0, 'go', markersize=10)
axs[0,0].plot(base_major[ieast, jeast, :30], -depths[:30], 'r', linewidth=3)
axs[0,0].plot(major[ieast, jeast, :30], -depths[:30], 'b', linewidth=3)
axs[0,0].set_xlim((0.1,0.4))
axs[0,0].set_title("M2: East: Major", fontsize=20)
axs[0,1].plot(east[:-10,2],-east[:-10,0], color='grey', label='Production', linewidth=3)
axs[0,1].plot(eastobs[2]/100., -eastobs[0], 'g', label='Observations', linewidth=3)
axs[0,1].plot(codar[1][1],0, 'go', markersize=10)
axs[0,1].plot(base_minor[ieast, jeast, :30], -depths[:30], 'r', label='Winter', linewidth=3)
axs[0,1].plot(minor[ieast, jeast, :30], -depths[:30], 'b', label='Freshet', linewidth=3)
axs[0,1].legend(loc='lower left', fontsize=14)
axs[0,1].set_xlim((-0.2,0.1))
axs[0,1].set_title("M2: East: Minor", fontsize=20)
axs[1,0].plot(east[:-10,3],-east[:-10,0], color='grey', linewidth=3)
axs[1,0].plot(eastobs[3], -eastobs[0], 'g', linewidth=3)
axs[1,0].plot(codar[1][2],0, 'go', markersize=10)
axs[1,0].plot(base_theta[ieast, jeast, :30]+29, -depths[:30], 'r', linewidth=3)
axs[1,0].plot(theta[ieast, jeast, :30]+29, -depths[:30],'b', linewidth=3)
axs[1,0].set_xlim((90,150))
axs[1,0].set_title("M2: East: Inclination", fontsize=20)
axs[1,1].plot(east[:-10,4],-east[:-10,0], color='grey', linewidth=3)
axs[1,1].plot(eastobs[4], -eastobs[0], 'g', linewidth=3)
axs[1,1].plot(codar[1][3],0, 'go', markersize=10)
axs[1,1].plot(base_phase[ieast, jeast, :30], -depths[:30], 'r', linewidth=3)
axs[1,1].plot(phase[ieast, jeast, :30], -depths[:30], 'b', linewidth=3)
axs[1,1].set_xlim((250,360))
axs[1,1].set_title("M2: East: Phase", fontsize=20)
axs[0,0].set_ylabel('Depth (m)', fontsize=14)
axs[1,0].set_ylabel("Depth (m)", fontsize=14)
for axe in (axs[0,0], axs[0,1], axs[1,0], axs[1,1]):
    for tick in axe.xaxis.get_major_ticks():
                tick.label.set_fontsize(14)
    for tick in axe.yaxis.get_major_ticks():
                tick.label.set_fontsize(14)
axs[0,0].set_xlabel("Major Axis (m/s)", fontsize=14)
axs[0,1].set_xlabel("Minor Axis (m/s)", fontsize=14)
axs[1,0].set_xlabel("Inclination (degrees)", fontsize=14)
axs[1,1].set_xlabel("Phase (degrees)", fontsize=14)

fig, axs = plt.subplots(2,2,figsize=(15,15))
axs[0,0].plot(cen[:-6,1],-cen[:-6,0], color='grey', linewidth=3)
axs[0,0].plot(central[1]/100., -central[0], 'g', linewidth=3)
axs[0,0].plot(codar[0][0],0, 'go', markersize=10)
axs[0,0].plot(base_major[icent, jcent, :35], -depths[:35], 'r', linewidth=3)
axs[0,0].plot(major[icent, jcent, :35], -depths[:35], 'b', linewidth=3)
axs[0,0].set_xlim((0.1,0.4))
axs[0,0].set_title("M2: Central: Major", fontsize=20)
axs[0,1].plot(cen[:-6,2],-cen[:-6,0], color='grey', label='Production', linewidth=3)
axs[0,1].plot(central[2]/100., -central[0], 'g', label='Observations', linewidth=3)
axs[0,1].plot(codar[0][1],0, 'go', markersize=10)
axs[0,1].plot(base_minor[icent, jcent, :35], -depths[:35], 'r', label='Winter', linewidth=3)
axs[0,1].plot(minor[icent, jcent, :35], -depths[:35], 'b', label='Freshet', linewidth=3)
axs[0,1].legend(loc='lower left', fontsize=14)
axs[0,1].set_xlim((-0.2,0.1))
axs[0,1].set_title("M2: Central: Minor", fontsize=20)
axs[1,0].plot(cen[:-6,3],-cen[:-6,0], color='grey', linewidth=3)
axs[1,0].plot(central[3], -central[0], 'g', linewidth=3)
axs[1,0].plot(codar[0][2],0, 'go', markersize=10)
axs[1,0].plot(base_theta[icent, jcent, :35]+29, -depths[:35], 'r', linewidth=3)
axs[1,0].plot(theta[icent, jcent, :35]+29, -depths[:35],'b', linewidth=3)
axs[1,0].set_xlim((90,150))
axs[1,0].set_title("M2: Central: Inclination", fontsize=20)
axs[1,1].plot(cen[:-6,4],-cen[:-6,0], color='grey', linewidth=3)
axs[1,1].plot(central[4], -central[0], 'g', linewidth=3)
axs[1,1].plot(codar[0][3],0, 'go', markersize=10)
axs[1,1].plot(base_phase[icent, jcent, :35], -depths[:35], 'r', linewidth=3)
axs[1,1].plot(phase[icent, jcent, :35], -depths[:35], 'b', linewidth=3)
axs[1,1].set_xlim((250,360))
axs[1,1].set_title("M2: Central: Phase", fontsize=20)
axs[0,0].set_ylabel('Depth (m)', fontsize=14)
axs[1,0].set_ylabel("Depth (m)", fontsize=14)
for axe in (axs[0,0], axs[0,1], axs[1,0], axs[1,1]):
    for tick in axe.xaxis.get_major_ticks():
                tick.label.set_fontsize(14)
    for tick in axe.yaxis.get_major_ticks():
                tick.label.set_fontsize(14)
axs[0,0].set_xlabel("Major Axis (m/s)", fontsize=14)
axs[0,1].set_xlabel("Minor Axis (m/s)", fontsize=14)
axs[1,0].set_xlabel("Inclination (degrees)", fontsize=14)
axs[1,1].set_xlabel("Phase (degrees)", fontsize=14)