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

seeder_path = $!seeder_path!$
infiles = $!musubi_path!$
omegas = $!musubi_omega_diff!$
dx = $!seeder_dx!$
height = $!seeder_height!$
bc_type = $!seeder_bc_type!$
params = [bc_type,omegas]

nHeight =[]
for idx in dx:
  nHeight.append(height/idx)

labels = []
#pattern=[]
filedict=dict()
filenames = []
for ipar in product(*params):
  pattern = str(ipar[0])+'/mus_omega_diff'+str(ipar[1])
  filename = str(ipar[0])+'_mus_omega_diff'+str(ipar[1])
  #labels.append('bc\_type = '+str(ipar[0])+', $\omega$ = '+str(ipar[1]))
  #labels.append(str(ipar[0])+', $\omega\_diff$ = '+str(ipar[1]))
  labels.append('bc\_type='+str(ipar[0]))
  filenames.append(filename)
  for iFile in infiles:
    if pattern in iFile:
      if filename not in filedict:
        filedict[filename] = [iFile]
      else:
        filedict[filename].append(iFile)
      #  filenames.append(str(ipar[0])+'_mus_omega'+str(ipar[1]))

#print (filedict[filenames[0]])
#print (pattern)
#print (filenames)
#print (infiles[0])
#print ('moments/mus_omega0.8' in infiles[0])

def extractL2Norm(infiles, variable, fileprefix):
  filename = fileprefix+'_'+variable+"_l2Norm_omega.txt"
  fileID = open(filename, "w")
  fileID.write("#relative_dx \t difference_l2norm \t analy_l2norm \t rel_l2norm \t diff_l2norm_avg\n")

  slope_file = fileprefix+'_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]
    fileID.write(str(dx[ifile]/dx[0])+'\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[ifile]
      # constant for quadratic slope
      c2 = c1/dx[ifile]

      slopeID = open(slope_file, 'w')
      #print(dx[iFile], Air_error, c1, c2)
    slopeID.write( str(dx[ifile]/dx[0])+'\t'\
                  +str(c1*dx[ifile])+'\t'+str(c2*dx[ifile]*dx[ifile])+'\n')


  fileID.close()
  return filename

vel_file = []
for iFile in range(len(filenames)):
  print (filenames[iFile], filedict[filenames[iFile]])
  vel_file.append(extractL2Norm(filedict[filenames[iFile]], 'vel', filenames[iFile]))


show_plot = $!show_plot!$
markers = ['k-x','k-s','k-v','k-*','k-d','k-p']
plot  = []
font_size = 12
figsize=[5,4]
set_majorgrid = True
#font_family = 'sans_serif'
#font_type = 'sans_serif'

markers = ['r-o','b-s','k-*','g-d','k-v','kx','k-p','k-^','k+']

for iFile in range(len(filenames)):
  sPlot = False
  ePlot = False
  if iFile == 0:
    sPlot = True
  if iFile == (len(filenames)-1):
    ePlot = True

  plot.append(dict(kind = 'xy', 
              data=[vel_file[iFile]],
              row=[1,-1], col=[1,4],dtype='float', 
              startplot = sPlot, endplot = ePlot,
              #xlabel = 'log($\Delta x$)', 
              #ylabel = 'log(relative l2norm error)',
              label = labels[iFile], ls = markers[iFile], 
              #xmin=10, xmax=1e4,
              xlogscale = True, ylogscale = True,
              legend=dict(loc=2),
              #label = '', ls = markers[ifile],
              #ls = '-o', 
              format = 'pdf', figname = 'Vel_relativeError_differentBC', dpi=100 ))

