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

sdr_path = $!seeder_path!$
mus_path = $!musubi_path!$
omegas = $!musubi_omega!$
refineLevel = $!seeder_refineLevel!$
dx = $!seeder_dx_ref!$
height = $!seeder_height!$

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]
    nHeight = height/dx_loc
    fileID.write(str(1.0/2**refineLevel[ifile])+'\t'+str(diff)+'\t'+str(anal)+'\t'
    +str(diff/anal)+'\t'+str(diff/nHeight)+'\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(1.0/2**refineLevel[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, 'vel')
vel_slope_files = 'vel_slope.txt'

press_file = extractL2Norm(mus_path, 'press')
press_slope_files = 'press_slope.txt'

wss_file= extractL2Norm(mus_path, 'wss')
wss_slope_files = 'wss_slope.txt'


show_plot = False
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 ))
