import os
import math
import glob
import numpy as np
from itertools import *

sdr_path = $!seeder_path!$ 
mus_path = $!musubi_path!$
nHeight = $!seeder_nHeight!$
height = 0.41
dx_ref = height/nHeight[0]
omega = $!musubi_omega!$
mus_path_1.8 = []
#mus_path_1.9 = []
mus_path_new = []
for ipath in range(len(mus_path)):
	if omega[0] in mus_path[ipath]:
		mus_path_1.8.append(mus_path[ipath])
#	elif omega[1] in mus_path[ipath]:
#		mus_path_1.9.append(mus_path[ipath]) 

mus_path_new.append(mus_path_1.8)
#mus_path_new.append(mus_path_1.9)

for iomega in range(len(omega)):
	
	os.system('mkdir omega'+omega[iomega])
		def extractL2Norm(infiles, variable):
		         filename = variable+"_l2Norm.txt"
		         fileID = open(filename, "w")
		         fileID.write("#dx/dx_ref \t difference_l2norm \t analy_l2norm \t rel_l2norm, diff_l2norm_avg\n")
		
		         slope_file = variable+'_slope.txt'
		         slopeID = open(slope_file, 'w')
		         slopeID.write('#relative dx 1st-order 2nd-order')
		         for ifile in range(len(infiles)):  
		           dataFile = glob.glob(infiles[ifile]+'/tracking/*'+variable+'_l2norm*.res')
		           data = np.genfromtxt(dataFile[0], dtype = None)
		           if data.size == data.shape[0]:
		             data.resize(1,data.size)
		           diff = data[-1,1]
		           anal = data[-1,2]
		           #dx_loc = dx/2**refineLevel[ifile]
		           dx_loc = height/nHeight[ifile]
		           #nHeight = height/dx_loc
		           fileID.write(str(dx_ref/dx[ifile])+'\t'+str(diff)+'\t'+str(anal)+'\t'
		           +str(diff/anal)+'\t'+str(diff/nHeight[ifile])+'\n')  
		           #print(data[-1,1],data[-1,2])
		
		           if ifile == 0:
		             # constant for linear slope
		             c1 = (diff/anal)/dx_loc
		             # constant for quadratic slope
		             c2 = c1/dx_loc
		
		             slopeID = open(slope_file, 'w')
		             #print(dx[iFile], Air_error, c1, c2)
		           slopeID.write( str(dx_ref/(height/nHeight[ifile]))+'\t'\
		                         +str(c1*dx_loc)+'\t'+str(c2*dx_loc*dx_loc)+'\n')
		
		           fileID.close()
		           slopeID.close()
		           return filename
		
	        vel_file = extractL2Norm(mus_path_new[iomega], 'vel')
		vel_slope_files = 'vel_slope.txt'
		
		press_file = extractL2Norm(mus_path_new[iomega], 'press')
		press_slope_files = 'press_slope.txt'
		
		wss_file= extractL2Norm(mus_path_new[iomega], 'wss')
		wss_slope_files = 'wss_slope.txt'
		
		
		show_plot = True
		markers = ['k-x','k-s','k-v','k-*','k-d','k-p']
		plot  = []
		font_size = 12
		figsize=[5,4]
		set_majorgrid = True
		set_minorgrid = True
		#font_family = 'sans_serif'
		#font_type = 'sans_serif'
		
		markers = ['r-o','b-s','k-*','g-d','k-v','kx','k-p','k-^','k+']
		
		plot.append(dict(kind = 'xy', 
		            data=[vel_slope_files],
		            col=[1,2],dtype='float',ls='k--', 
		            label = '1st-order',
		            xlogscale = True, 
		            ylogscale = True,
		            startplot = True))
		plot.append(dict(kind = 'xy', 
		            data=[vel_slope_files],
		            col=[1,3],dtype='float',ls='k-.',
		            xlogscale = True, 
		            ylogscale = True,
		            label='2nd-order'
		            ))
		
		plot.append(dict(kind = 'xy', 
		            data=[vel_file],
		            col=[1,4],dtype='float', 
		            startplot = False, endplot = True,
		            xlabel = '$\Delta x/\Delta x_{ref}$', 
		            ylabel = 'Velocity relative l2norm error',
		            ls = markers[0], 
		            xlogscale = True, ylogscale = True,
		            legend=dict(loc=2),
		            format = 'png', figname = 'Vel_relativeError', dpi=100 ))
		
		
		plot.append(dict(kind = 'xy', 
		            data=[press_slope_files],
		            col=[1,2],dtype='float',ls='k--', 
		            label = '1st-order',
		            xlogscale = True, 
		            ylogscale = True,
		            startplot = True))
		plot.append(dict(kind = 'xy', 
		            data=[press_slope_files],
		            col=[1,3],dtype='float',ls='k-.',
		            xlogscale = True, 
		            ylogscale = True,
		            label='2nd-order'
		            ))
		
		plot.append(dict(kind = 'xy', 
		            data=[press_file],
		            col=[1,4],dtype='float', 
		            startplot = False, endplot = True,
		            xlabel = '$\Delta x/\Delta x_{ref}$', 
		            ylabel = 'Pressure relative l2norm error',
		            ls = markers[0], 
		            legend=dict(loc=2),
		            xlogscale = True, ylogscale = True,
		            format = 'png', figname = 'press_relativeError', dpi=100 ))
		
		plot.append(dict(kind = 'xy', 
		            data=[wss_slope_files],
		            col=[1,2],dtype='float',ls='k--', 
		            label = '1st-order',
		            xlogscale = True, 
		            ylogscale = True,
		            startplot = True))
		plot.append(dict(kind = 'xy', 
		            data=[wss_slope_files],
		            col=[1,3],dtype='float',ls='k-.',
		            xlogscale = True, 
		            ylogscale = True,
		            label='2nd-order'
		            ))
		
		plot.append(dict(kind = 'xy', 
		            data=[wss_file],
		            col=[1,4],dtype='float', 
		            startplot = False, endplot = True,
		            xlabel = '$\Delta x/\Delta x_{ref}$', 
		            ylabel = 'WSS relative l2norm error',
		            ls = markers[0], 
		            xlogscale = True, ylogscale = True,
		            legend=dict(loc=2),
		            format = 'png', figname = 'wss_relativeError', dpi=100 ))
	
		os.system('mv vel_slope.txt omega'+omega[iomega], 'mv press_slope.txt omega'+omega[iomega], 'mv wss_slope.txt omega'+omega[iomega]) 		
	
		# show_plot = True
		# marker_slope = ['k-x','k-s','k-v','k-*','k-d','k-p']
		# plot  = []
		# font_size = 12
		# figsize=[5,4]
		# set_majorgrid = True
		# set_minorgrid = True
		# markers = ['r-o','b-s','g-*','r-d','b-v','gx','r-p','b-^','g+']
		# vel_slope_file = 'vel_simple_slope.txt'
		# plot.append(dict(kind = 'xy', 
		#             data=[vel_slope_file],
		#             col=[1,2],dtype='float',ls=marker_slope[0], 
		#             label = '1st-order\_',#+str(names[iplot]),
		#             xlogscale = True, 
		#             ylogscale = True,
		#             startplot = True,
		#             endplot = False))
		# plot.append(dict(kind = 'xy', 
		#             data=[vel_slope_file],
		#             col=[1,3],dtype='float',ls=marker_slope[1],
		#             label='2nd-order\_',#+str(names[iplot]),
		#             xlogscale = True, 
		#             ylogscale = True,
		#             #startplot = stplot,
		#             endplot = False,
		#             ))
		# for iplot in range(len(functions)):
		#     #vel_slope_files.append('vel_'+str(functions[iplot])+'_slope.txt')
		#     stplot = False
		#     eplot = False
		#     if iplot == len(functions)-1:
		#         eplot = True
		# 
		#     if iplot != len(functions)-1:
		#         plot.append(dict(kind = 'xy', 
		#                     data=['vel_'+str(functions[iplot])+'_l2norm.txt'],
		#                     col=[1,4],dtype='float', label = str(names[iplot])+'\_error',
		#                     startplot = stplot, endplot = eplot,
		#                     ls = markers[iplot], ))
		#     else:
		#         plot.append(dict(kind = 'xy', 
		#                     data=['vel_'+str(functions[iplot])+'_l2norm.txt'],
		#                     col=[1,4],dtype='float', label = str(names[iplot])+'\_error',
		#                     startplot = stplot, endplot = eplot,
		#                     ls = markers[iplot], 
		#                     xlabel = '$\Delta x/\Delta x_{ref}$', 
		#                     ylabel = 'Velocity relative l2norm error',
		#                     xlogscale = True, ylogscale = True,
		#                     legend=dict(loc=4),
		#                     format = 'png', figname = 'Vel_relativeError', dpi=100 ))
