import pyGlacier as pg
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation


fontSize = 7
plt.rcParams.update({'font.size': fontSize, 'legend.fontsize': fontSize-1,'legend.handlelength': 2, 'lines.linewidth': 0.5, 'text.usetex': True})
secondyear = 365.25*24*60*60
xRange = [0,12.5]

runname = 'RunFig5B_coupledSIA-SSA_kjetil'

a = pg.PostProcess(runname,file_format='mat',resolution=1)
a.t = a.t-500

fig = plt.figure(figsize = (18/2.54/2, 24/2.54/3.5)) #half page width in inches


indfind = np.where((a.t > 585.9) & (a.t < 586.1))

ind = [11475-13,11475,11475+26]
#ind = [5050,5056,5064]

for j in range(3):
	
	i = ind[j]

	ax1 = fig.add_axes([.13+.29*j, .72, .26, .27]); ax1.tick_params(labelbottom = False)
	ax2 = fig.add_axes([.13+.29*j, .42, .26, .27]); ax2.tick_params(labelbottom = False)
	ax3 = fig.add_axes([.13+.29*j, .12, .26, .27])

	ax1.cla()
	#ax1.set_title('t = '+str(round(a.t[i],1))+' years')
	ax1.text(.25,.82,'t = '+str(round(a.t[i],1))+' years', transform = ax1.transAxes, fontsize=fontSize-1)#, path_effects=[PathEffects.withStroke(linewidth=0.8,foreground='w')])



	glacier = ax1.fill(np.append(a.x,[a.x[-1],a.x[0]])/1000,np.append(a.b+a.H[i,:],[a.b[-1], a.b[-1]])/1000, facecolor = [170/255,211/255,233/255], edgecolor = [.7*170/255,.7*211/255,.7*233/255])
	base = ax1.fill(np.append(a.x,[a.x[-1],a.x[0]])/1000,np.append(a.b,[a.b[-1], a.b[-1]])/1000, facecolor = [.5,.5,.5], edgecolor = [.3,.3,.3])
	ax1.set_xlim(xRange[0],xRange[1]);
	ax1.set_ylim(.9,2.2);
	ax1.set_ylabel('y [km]')
	ax1.grid(linestyle=':',linewidth='0.5')

	ax2.cla()
	ax2.semilogy(a.x/1000,a.U[i,:])
	ax2.set_xlim(xRange[0],xRange[1]);
	ax2.set_ylim(1,5000)
	ax2.grid(linestyle=':',linewidth='0.5')

	effective_normal_stress = a.C*(a.rho[0]*a.g[0]*a.H[i,:]/1e6-(a.hydraulic_potential[i,:] - a.water_density[0]*a.g[0]*a.b/1e6) )
	tau_m = 2.0*a.H[i,:]/a.width[0] * (5.0/(a.A[0]*a.width[0]))**(1.0/3.0)*(np.abs(a.U_SSA[i,:])/secondyear)**(1.0/3.0) / 1e6
	tau_b = a.state_parameter[i,:]*(a.Sliding[i,:]/secondyear/a.As[0])**(1/3)/1e6
	ax3.cla()
	ax3.plot(a.x/1000,effective_normal_stress, label='$C \sigma_N$')
	ax3.plot(a.x/1000,tau_b, label = '$\\tau_b$')
	ax3.plot(a.x/1000,tau_m, label = '$\\tau_m$')
	ax3.set_xlim(xRange[0],xRange[1]);
	ax3.set_ylim(0,1);
	ax3.set_yticks([0.0,0.2,0.4,0.6,0.8])
	ax3.grid(linestyle=':',linewidth='0.5')
	
	if j==0:
		ax1.text(1,1,'a',fontsize='10')
		ax2.text(1,800,'b',fontsize='10') 		
		ax3.text(1,0.8,'c',fontsize='10')
		
	if j==1:
		ax1.text(1,1,'d',fontsize='10')
		ax2.text(1,800,'e',fontsize='10') 		
		ax3.text(1,0.8,'f',fontsize='10')
		
	if j==2:
		ax1.text(1,1,'g',fontsize='10')
		ax2.text(1,800,'h',fontsize='10') 		
		ax3.text(1,0.8,'i',fontsize='10')
	

	if j>0:
		ax1.set_ylabel('')
		ax1.tick_params(labelleft = False)
		ax2.set_xlabel('')
		ax2.tick_params(labelbottom = False)
		ax2.set_ylabel('')
		ax2.tick_params(labelleft = False)
		ax3.set_xlabel('x [km]',labelpad=1)
		ax3.set_ylabel('')
		ax3.tick_params(labelleft = False)

	else:
		ax1.set_ylabel('y [km]')
		ax2.set_ylabel('U [m/yr]')
		ax3.set_xlabel('x [km]',labelpad=1)
		ax3.set_ylabel('Stress [MPa]')
		ax3.legend(frameon=False)

plt.savefig('Figure_4.png',format='png',dpi=1500)