diff --git a/.pylint.d/Evaluate_dl1.stats b/.pylint.d/Evaluate_dl1.stats
new file mode 100644
index 0000000000000000000000000000000000000000..a73817f49a2ab6fd69ff4048e2e419c05be72a7c
Binary files /dev/null and b/.pylint.d/Evaluate_dl1.stats differ
diff --git a/All_US_errors_DoubleElbow_alpha_30.0_refl_1.npz b/All_US_errors_DoubleElbow_alpha_30.0_refl_1.npz
new file mode 100644
index 0000000000000000000000000000000000000000..3f16caa48f1ce929207f763d9881109784495b86
Binary files /dev/null and b/All_US_errors_DoubleElbow_alpha_30.0_refl_1.npz differ
diff --git a/All_US_errors_SingleElbow_alpha_30.0_refl_1.npz b/All_US_errors_SingleElbow_alpha_30.0_refl_1.npz
new file mode 100644
index 0000000000000000000000000000000000000000..707e2d7f8d0f2c0e481fee8bce25b8eb745dd17b
Binary files /dev/null and b/All_US_errors_SingleElbow_alpha_30.0_refl_1.npz differ
diff --git a/FlowMeters/DoubleElbow/Ultrasonic_diametral_1_30.npz b/FlowMeters/DoubleElbow/Ultrasonic_diametral_1_30.npz
index d8ca1ef8a08ec56a9ef856f1ca6e5d71a5d49149..673e4bd627c3451ecc1278c0815383ee697e039d 100644
Binary files a/FlowMeters/DoubleElbow/Ultrasonic_diametral_1_30.npz and b/FlowMeters/DoubleElbow/Ultrasonic_diametral_1_30.npz differ
diff --git a/FlowMeters/DoubleElbow/Ultrasonic_diametral_1_30_alt.npz b/FlowMeters/DoubleElbow/Ultrasonic_diametral_1_30_alt.npz
new file mode 100644
index 0000000000000000000000000000000000000000..d8ca1ef8a08ec56a9ef856f1ca6e5d71a5d49149
Binary files /dev/null and b/FlowMeters/DoubleElbow/Ultrasonic_diametral_1_30_alt.npz differ
diff --git a/FlowMeters/DoubleElbow/Ultrasonic_diametral_1_30_ref.npz b/FlowMeters/DoubleElbow/Ultrasonic_diametral_1_30_ref.npz
new file mode 100644
index 0000000000000000000000000000000000000000..088834358aed6ec10f1f4b9a901234145a783446
Binary files /dev/null and b/FlowMeters/DoubleElbow/Ultrasonic_diametral_1_30_ref.npz differ
diff --git a/FlowMeters/Double_S_Elbow/Ultrasonic_diametral_1_30.npz b/FlowMeters/Double_S_Elbow/Ultrasonic_diametral_1_30.npz
index 0c4525f9c60354e9193986ac8b25edacbee3845f..e27946b7809b3cc58e969d01b16ae76b652e21d2 100644
Binary files a/FlowMeters/Double_S_Elbow/Ultrasonic_diametral_1_30.npz and b/FlowMeters/Double_S_Elbow/Ultrasonic_diametral_1_30.npz differ
diff --git a/FlowMeters/SingleElbow/Ultrasonic_diametral_1_29.npz b/FlowMeters/SingleElbow/Ultrasonic_diametral_1_29.npz
deleted file mode 100644
index cdeabd042d4287bf2d4a664e29926e99cdeb9667..0000000000000000000000000000000000000000
Binary files a/FlowMeters/SingleElbow/Ultrasonic_diametral_1_29.npz and /dev/null differ
diff --git a/FlowMeters/SingleElbow/Ultrasonic_diametral_1_30.npz b/FlowMeters/SingleElbow/Ultrasonic_diametral_1_30.npz
index 25f9b7f403f65529ca2cce218d6fc510a4b00406..b59f4cc99f1cfeacf5c80fef1031a7bcf6377046 100644
Binary files a/FlowMeters/SingleElbow/Ultrasonic_diametral_1_30.npz and b/FlowMeters/SingleElbow/Ultrasonic_diametral_1_30.npz differ
diff --git a/FlowMeters/SingleElbow/transpose_Ultrasonic_diametral_1_30.npz b/FlowMeters/SingleElbow/transpose_Ultrasonic_diametral_1_30.npz
deleted file mode 100644
index b321c2077e6da909e7590c7dffb2aa5ee1162f0e..0000000000000000000000000000000000000000
Binary files a/FlowMeters/SingleElbow/transpose_Ultrasonic_diametral_1_30.npz and /dev/null differ
diff --git a/Flow_class.py b/Flow_class.py
index 646e2b073349543ea56f26ce82c34a3357ab754c..2e7092a09ebf2b80936750b35867e9c5b3ecc600 100644
--- a/Flow_class.py
+++ b/Flow_class.py
@@ -11,10 +11,11 @@ import scipy.interpolate as interpol
#import matplotlib.pyplot as plt
import json
from math import pi
-#import sys
+import sys
#from scipy.interpolate import Rbf
#sys.path.append('E:/NextCloud/FlowProfiles/')
-#sys.path.append('../../FlowProfiles')
+sys.path.append('./FlowProfiles/')
+sys.path.append('./HelpingFunctions/')
# import gh-class from global directory
from class_flowprofiles import flowprofiles
#sys.path.append('E:/NextCloud/GlobalPyFunc/')
@@ -26,8 +27,6 @@ import Integrate2D as int2d
#import pykrige as pk
import matplotlib.pyplot as plt
-import pdb
-
def do_reconstruct(A,phi,Nmodes=0,umean = 0):
# A are the POD-coefficients
# phi are the modes
@@ -163,6 +162,8 @@ class Elbow_profile():
#
coeffs, self.modes, self.umean = loadPODres(case)
self.Nmodes = self.modes.shape[-1]
+ # refine the resolution of the modes in radial direction by akima-spline interpolation
+ # self.refine_modes()
# split the modes in x, y, z ? better not because of the reconstruction?
#
if case=="SingleElbow":
@@ -189,7 +190,37 @@ class Elbow_profile():
self.make_fully()
#
# end Init
-
+ def refine_modes(self,refine_level = 2):
+ r = self.r[0]
+ for i in range(1,refine_level):
+ # the r-values that can be inserted
+ r_ins = r[0:-1] + np.diff(r)/2
+ r_ref = np.sort(np.concatenate((r,r_ins)))
+ r = r_ref
+ # prepare interpolation
+ modemat = self.modes.reshape(3,*self.r.shape,-1)
+ aki = interpol.Akima1DInterpolator(self.r[0],modemat,axis=2)
+ modematneu = aki(r_ref)
+ # overwrite the modes with finer grid
+ self.modes = modematneu.reshape(-1,self.Nmodes)
+ # refine Umean too
+ umeanmat = self.umean.reshape(3,*self.r.shape)
+ aki = interpol.Akima1DInterpolator(self.r[0],umeanmat,axis=2)
+ umeanmatneu = aki(r_ref)
+ self.umean = umeanmatneu.reshape(1,-1)
+ phi = self.phi[:,0]
+ # overwrite the private variables
+ self.phi,self.r = np.meshgrid(phi,r_ref,indexing='ij')
+ self.xm,self.ym = pol2cart(self.r,self.phi)
+ self.x = self.xm.ravel()
+ self.y = self.ym.ravel()
+ self.Nphi = self.r.shape[0]
+ self.Nr = self.r.shape[1]
+ #
+ self.int_weights = int2d.get_int_weights(x=self.x,y=self.y,method="tri")
+ #
+ self.make_fully()
+ #
def get_profile(self,Rk,dl,dist,x=None,y=None,addfully = False,asmat = True):
# returns the velocity profile at a certain distance and coordinates x,y, (if given)
# to evaluate the spline interpolation and get the coefficient matrix at the desired
@@ -230,10 +261,7 @@ class Elbow_profile():
# now reconstruct the velocities
# mx,my,mz = self.modes.reshape(3,*self.r.shape,-1)
modemat = self.modes.reshape(3,*self.r.shape,-1)
-
- # modes_ux_mat = self.modes_ux.reshape(*self.r.shape,-1)
- # modes_uy_mat = self.modes_uy.reshape(*self.r.shape,-1)
- # modes_uz_mat = self.modes_uz.reshape(*self.r.shape,-1)
+ #
if reverse:
#
modemat = modemat[:,:,::-1,:]
@@ -242,7 +270,6 @@ class Elbow_profile():
# calculate the weighting factors for the secondary flow components
wx = -D*np.cos(self.phi)
wy = -D*np.sin(self.phi)
- # print("mode.shape")
# if not reversed
else:
u = do_reconstruct_path(A, modemat, umean = self.umean.reshape(3,self.Nphi,self.Nr))
diff --git a/GUI_Elbow.py b/GUI_Elbow.py
index 4146b6b23fd71d1d4c99dde217f60cf9c4020fe6..8026169487ee997191bf9985ada11c8c99274207 100644
--- a/GUI_Elbow.py
+++ b/GUI_Elbow.py
@@ -245,18 +245,19 @@ class mainPanel(wx.Panel):
#name=CheckBoxNameStr)
self.checkbox_diff.Bind(wx.EVT_CHECKBOX, self.onChecked_diff)
#
- self.checkbox_uncert = wx.CheckBox(self, id=wx.ID_ANY, label="Show uncertainty") #, pos=DefaultPosition,
- #size=DefaultSize, style=0, validator=DefaultValidator,
- #name=CheckBoxNameStr)
- self.checkbox_uncert.Bind(wx.EVT_CHECKBOX, self.onChecked_uncert)
+ # self.checkbox_uncert = wx.CheckBox(self, id=wx.ID_ANY, label="Show uncertainty") #, pos=DefaultPosition,
+ # #size=DefaultSize, style=0, validator=DefaultValidator,
+ # #name=CheckBoxNameStr)
+ # self.checkbox_uncert.Bind(wx.EVT_CHECKBOX, self.onChecked_uncert)
# add a button to apply uncertainties
self.btn_uncertainty = wx.Button(self,-1,"Calc Uncert")
self.btn_uncertainty.Bind(wx.EVT_BUTTON,self.OnClicked_btn_uncertainty)
# add the elements to the boxsizer
# make the checkboxes and upper buttons
+ self.case_lbl = wx.StaticText(self, wx.ID_ANY, label =flow.case)
#hbox_check.AddSpacer(10)
hbox_check.Add(self.checkbox_diff,1, wx.ALIGN_TOP, border = b)
- hbox_check.Add(self.checkbox_uncert,1, wx.ALIGN_TOP, border = b)
+ hbox_check.Add(self.case_lbl,1, wx.ALIGN_TOP, border = b)
hbox_check.Add(self.btn_uncertainty,0,wx.ALIGN_TOP, border = b)
vbox_slider.Add(hbox_check,1, flag = wx.EXPAND | wx.TOP, border = 1)
#
@@ -450,7 +451,7 @@ class mainPanel(wx.Panel):
print("Flow meter successfully loaded from file")
return 1
except Exception as e:
- print("flow meter could not be loaded -> it will be calculated, this takes some time " + str(e))
+ print("flow meter could not be loaded -> it will be calculated for the whole data set, this takes some time " + str(e))
return 0
def OnClicked_btn_flowmeter(self,e):
global fm
@@ -798,7 +799,7 @@ class mainPanel(wx.Panel):
# the contour Plot in axis 0 (left side)
#ax.cla()
#self.__init__(self.parent)
-
+ self.case_lbl.SetLabel(flow.case)
self.ax0.cla()
self.ax1.cla()
@@ -1072,15 +1073,12 @@ class mainPanel(wx.Panel):
class Main(wx.Frame):
def __init__(self):
-
-
-
+ #
wx.Frame.__init__(self, parent = None,id = wx.ID_ANY, title = "Flowmeter behind Elbow", size = (1400,1000))
#
#splitter = wx.SplitterWindow(self)
self.mypanel = mainPanel(self)
-
-
+ #
# create a menu bar ------------------------------------------
menubar = wx.MenuBar()
fileMenu = wx.Menu()
@@ -1088,14 +1086,17 @@ class Main(wx.Frame):
fileItem_se = fileMenu.Append(wx.ID_ANY, 'Single Elbow', 'load the single 90 degree Elbow')
fileItem_de = fileMenu.Append(wx.ID_ANY, 'Double Elbow', 'load the double 90 degree Elbow Out-Of-Plane')
fileItem_dse = fileMenu.Append(wx.ID_ANY, 'Double_S_Elbow', 'load the double 90 degree Elbow In-Plane')
+ fileItem_ref = fileMenu.Append(wx.ID_ANY, 'Refine modes', 'refines the modes in radial direction for testing')
fileItem_q = fileMenu.Append(wx.ID_EXIT, 'Quit', 'Quit application')
- menubar.Append(fileMenu, '&File')
+ menubar.Append(fileMenu, '&Load Case')
self.SetMenuBar(menubar)
self.Bind(wx.EVT_MENU, self.OnQuit, fileItem_q)
self.Bind(wx.EVT_MENU, self.OnLoad_SE, fileItem_se)
self.Bind(wx.EVT_MENU, self.OnLoad_DE, fileItem_de)
self.Bind(wx.EVT_MENU, self.OnLoad_DSE, fileItem_dse)
+ self.Bind(wx.EVT_MENU, self.Onref, fileItem_ref)
+
# ----------------------------------------------------------
# bottom = BottomPanel(splitter, top)
# splitter.SplitHorizontally(top, bottom)
@@ -1112,6 +1113,9 @@ class Main(wx.Frame):
# self.mypanel = mainPanel(self)
# self.mypanel.__init__(self)
self.__init__()
+ #self.mypanel.case_lbl.SetLabel(flow.case)
+ self.mypanel.draw()
+ #
def OnLoad_DE(self,e):
global flow
flow = Elbow_profile('DoubleElbow')
@@ -1119,6 +1123,9 @@ class Main(wx.Frame):
wx.Frame.__init__(self, parent = None,id = wx.ID_ANY, title = "Flowmeter behind Double Elbow", size = (1400,1000))
self.mypanel = mainPanel(self)
self.mypanel.__init__(self)
+ #self.mypanel.case_lbl.SetLabel(flow.case)
+ self.mypanel.draw()
+ #
def OnLoad_DSE(self,e):
global flow
flow = Elbow_profile('Double_S_Elbow')
@@ -1126,7 +1133,13 @@ class Main(wx.Frame):
wx.Frame.__init__(self, parent = None,id = wx.ID_ANY, title = "Flowmeter behind Double S Elbow", size = (1400,1000))
self.mypanel = mainPanel(self)
self.mypanel.__init__(self)
-
+ #self.mypanel.case_lbl.SetLabel(flow.case)
+ self.mypanel.draw()
+ #
+ def Onref(self,e):
+ global flow
+ flow.refine_modes()
+ self.mypanel.draw()
# temp_dist = 10
# temp_Rc = 1.5
diff --git a/Get_Elbows_POD_int.py b/Get_Elbows_POD_int.py
deleted file mode 100644
index 6287299a69c39f3e41e13946e4b88d6af5fec279..0000000000000000000000000000000000000000
--- a/Get_Elbows_POD_int.py
+++ /dev/null
@@ -1,653 +0,0 @@
-# -*- coding: utf-8 -*-
-"""
-Created on Thu Apr 28 11:01:28 2022
-
-@author: weisse02
-"""
-
-# -*- coding: utf-8 -*-
-"""
-Created on Mon Dec 7 11:56:36 2020
-
-@author: Ichyc
-"""
-
-import os
-import numpy as np
-import scipy.interpolate as interpol
-#import matplotlib.pyplot as plt
-import json
-from math import pi
-import sys
-#from scipy.interpolate import Rbf
-#sys.path.append('E:/NextCloud/FlowProfiles/')
-#sys.path.append('../../FlowProfiles')
-# import gh-class from global directory
-from class_flowprofiles import flowprofiles
-#sys.path.append('E:/NextCloud/GlobalPyFunc/')
-from niceplot import *
-import Integrate2D as int2d
-sys.path.append('../Elbow_Variations/Auswertung/')
-from Fully_simu import get_fully_SA
-import pykrige.kriging_tools as kt
-import pykrige as pk
-import matplotlib.pyplot as plt
-
-#
-#
-#
-#
-
-
-def do_reconstruct(A,phi,Nmodes=0,umean = 0):
- if len(A.shape) == 1:
- A = A.reshape(1,A.shape[0])
- u_tilde = np.zeros((A.shape[0],phi.shape[0]))
- # if number of modes is not given take
- # all the modes given in phi and A
- if Nmodes == 0:
- Nmodes = phi.shape[1]
- for k in range(0,Nmodes):
- u_tilde = u_tilde + np.tensordot(A[:,k],phi[:,k], axes = 0)
- #
- return umean + np.real(u_tilde)
-
-def get_POD_coeffs(rk,dist,interpol_list,inttype):
- # A is the coefficient matrix that is recovert from interpolation A.shape[1] must be number of used modes
- rk = np.array(rk)
- dist = np.array(dist)
- A = np.zeros((*rk.shape,len(interpol_list) ))
- if len(A.shape) == 1:
- A = A.reshape(1,A.shape[0])
- A = A.T
- #
- for i, inter_i in enumerate(interpol_list):
- if inttype == "lin":
- A[i] = inter_i((rk,dist)).T
- elif inttype == "rbf":
- A[i] = inter_i(rk,dist).T
- elif inttype == "spline":
- A = np.zeros((rk.size,dist.size,len(interpol_list) ))
- if len(A.shape) == 1:
- A = A.reshape(1,A.shape[0])
- A = A.T
- # print(A.shape)
- A[i] = inter_i(rk,dist).T
- return A.T
-
-
-def do_reconstruct_add(A,phi,Nmodes,u_tilde,umean = 0):
- #u_tilde = np.zeros((A.shape[0],phi.shape[0]))
- for k in range(Nmodes-1,Nmodes):
- u_tilde = u_tilde + np.tensordot(A[:,k],phi[:,k], axes = 0)
- #
- return umean + np.real(u_tilde)
-
-
-def loadPODres():
- data_path = os.path.abspath("./POD_res/")
- #
- with open(os.path.join(data_path,"ux_PODres_Elbow.json")) as json_file:
- ux = json.load(json_file)
- with open(os.path.join(data_path,"uy_PODres_Elbow.json")) as json_file:
- uy = json.load(json_file)
- with open(os.path.join(data_path,"uz_PODres_Elbow.json")) as json_file:
- uz = json.load(json_file)
- # make numpy arrays from the lists
- keys = list(uz.keys())
- for key in keys:
- ux[key] = np.array(ux[key])
- uy[key] = np.array(uy[key])
- uz[key] = np.array(uz[key])
- #
- return ux, uy, uz
-
-def loadCoordinates():
- coord_file = os.path.abspath("./Dist_Coords.json")
- with open(coord_file) as json_file:
- coords = json.load(json_file)
- return coords
-
-
-
-#
-def find_match(strlist,s):
- matched_indexes = []
- i = 0
- length = len(strlist)
-
- while i < length:
- if s in strlist[i]:
- matched_indexes.append(i)
- i += 1
- # there is only one linear index make it a list
- # if not isinstance(matched_indexes, list):
- # lin_ind = [lin_ind]
- return matched_indexes
-
-
-
-class Elbow_profile():
- def __init__(self):
- # load the flow profiles behind the elbows
- delete_last=True
-
- # ------------------------------------------------------------
- coords = loadCoordinates()
- #
- self.D = coords["D"]
- self.dist = np.array(coords["dist"])/self.D
- self.dist = self.dist[:-1]
- # the courvature radius
- self.Rk = np.array(coords["Rk"])/self.D/1000
- self.r = np.array(coords["r"])
- self.phi = np.array(coords["phi"])
-
- # remove the last value of the matrix this is where r = 1
- if delete_last:
- self.r_wall = self.r[:,-1]
- self.phi_wall = self.phi[:,-1]
- self.x_wall,self.y_wal = pol2cart(self.r_wall,self.phi_wall)
- #
- self.r = self.r[:,0:-1]
- self.phi = self.phi[:,0:-1]
- # delete also the duplicated angle
- # r = r[0:-1,:]
- # phi_a = phi_a[0:-1,:]
- #
- self.xm,self.ym = pol2cart(self.r,self.phi)
- self.x = self.xm.ravel()
- self.y = self.ym.ravel()
- #
- # initialize all variables
- N0 = 3 # len(flownames)
- N1 = self.x.size
- N2 = self.dist.size
- N3 = self.Rk.size
- Nt = self.dist.size * self.Rk.size
- #
- pod_x,pod_y,pod_z = loadPODres()
- self.Nmodes = pod_x["modes"].shape[-1]
- self.modes_ux = pod_x["modes"]
- self.modes_uy = pod_y["modes"]
- self.modes_uz = pod_z["modes"]
- self.ux_mean = pod_x["umean"]
- self.uy_mean = pod_y["umean"]
- self.uz_mean = pod_z["umean"]
- #
- #
- pod_x["coeffs"] = pod_x["coeffs"].reshape((len(self.Rk),len(self.dist),self.Nmodes))
- pod_y["coeffs"] = pod_y["coeffs"].reshape((len(self.Rk),len(self.dist),self.Nmodes))
- pod_z["coeffs"] = pod_z["coeffs"].reshape((len(self.Rk),len(self.dist),self.Nmodes))
- #
- # plot for test
- # indrk = 7
- # inddist = 20
- # rktest = self.Rk[7]
- # disttest = self.dist[20]
- # A_test = pod_z["coeffs"][indrk,inddist,:]
- # print(A_test.shape)
- # print("--- Atest")
- # utest = do_reconstruct(A_test,pod_z["modes"],self.Nmodes,self.uz_mean ).ravel()
- # nicecontour(self.x,self.y,utest)
-
- polydeg = 3
- self.splineint_ux = []
- self.splineint_uy = []
- self.splineint_uz = []
- #
- for i in range(0,self.Nmodes):
- self.splineint_ux.append(interpol.RectBivariateSpline(self.dist,self.Rk,pod_x["coeffs"][:,:,i].T,kx= polydeg,ky=polydeg ))
- self.splineint_uy.append(interpol.RectBivariateSpline(self.dist,self.Rk,pod_y["coeffs"][:,:,i].T,kx= polydeg,ky=polydeg ))
- self.splineint_uz.append(interpol.RectBivariateSpline(self.dist,self.Rk,pod_z["coeffs"][:,:,i].T,kx= polydeg,ky=polydeg ))
-
- del pod_x
- del pod_y
- del pod_z
- #
- self.int_weights = int2d.get_int_weights(x=self.x,y=self.y,method="tri")
-
- self.make_fully()
- # The parametrized function to be plotted
- #self.get_profile = interpol.RegularGridInterpolator((self.Rk,self.dist), self.uall + self.u_fully)
- #
- def get_profile(self,dist,Rk,x=None,y=None,addfully = False):
- # returns the velocity profile at a certain distance and coordinates x,y, (if given)
- # to evaluate the spline interpolation and get the coefficient matrix at the desired
- # values of dist and Rk
- A_ux = np.array([splineint_i(dist,Rk) for splineint_i in self.splineint_ux] ).T
- A_uy = np.array([splineint_i(dist,Rk) for splineint_i in self.splineint_uy] ).T
- A_uz = np.array([splineint_i(dist,Rk) for splineint_i in self.splineint_uz] ).T
- #
- A_ux = A_ux.reshape(-1,A_ux.shape[-1])
- A_uy = A_uy.reshape(-1,A_uy.shape[-1])
- A_uz = A_uz.reshape(-1,A_uz.shape[-1])
- # now reconstruct the velocities
- print(A_ux.shape)
- print(A_uy.shape)
- print(A_uz.shape)
- ux = do_reconstruct(A_ux, self.modes_ux, umean = self.ux_mean)
- uy = do_reconstruct(A_uy, self.modes_uy, umean = self.uy_mean)
- uz = do_reconstruct(A_uz, self.modes_uz, umean = self.uz_mean)
- #
- if addfully:
- uz += self.u_fully.ravel()
- return ux,uy,uz
-
- def make_fully(self,Re = 5e4,ks=0.0,ret=False):
- # generating fully
- uvol = 1
- self.Re = Re
- visc = self.D*uvol/self.Re
- Q = uvol*pi*(self.D/2)**2 *3600
- #
- fully_class = flowprofiles(Q,self.D,visc,ks= ks*self.D)
- # fully profile
- self.u_fully = fully_class.get_u_u_vol(self.r)
- #plt.plot(r.ravel(),u_fully)
- #self.uall = self.uall + minus*self.u_fully
- #self.get_profile = interpol.RegularGridInterpolator((self.Rk,self.dist), self.uall)
- # return the profile if wanted
- if ret:
- return self.u_fully
- # def add_fully(self,add = True):
- # if add:
- # self.get_profile = interpol.RegularGridInterpolator((self.Rk,self.dist), self.uall + self.u_fully)
- # else:
- # self.get_profile = interpol.RegularGridInterpolator((self.Rk,self.dist), self.uall)
-#
-# #- ----------------------------------------------------- ---------------------
-
-
-
-
-if __name__ == "__main__":
- myflow = Elbow_profile()
- ux,uy,uz = myflow.get_profile(5,1.425,addfully=False)
- nicecontour(myflow.x,myflow.y,uz.ravel())
-
-
-
-# #- ----------------------------------------------------- ---------------------
-if 0: #__name__ == "__main__":
- #
- # load the coordinates
- delete_last = True
- #
- coords = loadCoordinates()
- #
- D = coords["D"]
- dist = np.array(coords["dist"][0:-1])/D
- # the courvature radius
- Rk = np.array(coords["Rk"])/D/1000
- r = np.array(coords["r"])
- # remove the last value of the matrix this is where r = 1
- if delete_last: r = r[:,0:-1]
- phi = np.array(coords["phi"])
- # remove the last value of the matrix this is where r = 1
- if delete_last: phi = phi[:,0:-1]
- xm, ym = pol2cart(r,phi)
- x = xm.ravel()
- y = ym.ravel()
- #
- # initialize all variables
- N0 = 3 # len(flownames)
- N1 = x.size
- N2 = dist.size
- N3 = Rk.size
- Nt = dist.size * Rk.size
- #
- pod_x,pod_y,pod_z = loadPODres()
- Nmodes = pod_x["modes"].shape[-1]
- #pod_x["modes"] = pod_x["modes"].reshape(*xm.shape,Nmodes)
- #
-
- #raise SystemExit(0) # stops the program to continue
-
- pod_x["coeffs"] = pod_x["coeffs"].reshape((len(Rk),len(dist),Nmodes))
- pod_y["coeffs"] = pod_y["coeffs"].reshape((len(Rk),len(dist),Nmodes))
- pod_z["coeffs"] = pod_z["coeffs"].reshape((len(Rk),len(dist),Nmodes))
- #
- # plot for test
- indrk = 7
- inddist = 20
- rktest = Rk[7]
- disttest = dist[20]
- A_test = pod_z["coeffs"][indrk,inddist,:]
- utest = do_reconstruct(A_test,pod_z["modes"],Nmodes,pod_z["umean"] ).ravel()
- nicecontour(x,y,utest)
- #
-
- # OK_x = pk.OrdinaryKriging(data[:, 0], data[:, 1], pod_x, variogram_model='linear',
- # verbose=False, enable_plotting=False)
- # # for uz
- # OK_y = pk.OrdinaryKriging(data[:, 0], data[:, 1], pod_y, variogram_model='linear',
- # verbose=False, enable_plotting=False)
- # # for uz
- # OK_z = pk.OrdinaryKriging(data[:, 0], data[:, 1], pod_z, variogram_model='linear',
- # verbose=False, enable_plotting=False)
-
-
- ## prepare the cross validation study
- save_path = "./POD_interpol/"
- # the reference grid:
- ddist,rrk = np.meshgrid(dist,Rk)
- # raise SystemExit(0)
- # tubsample the distances
- dist_red = []
- templen = 1000
- for n in range(2,len(dist)):
- tempdists = np.concatenate( (dist[0:-1:n], [dist[-1]]) )
- if tempdists.size < templen:
- dist_red.append(tempdists)
- templen = tempdists.size
- dist_red.reverse()
- dist_red = dist_red[0::2]
- # subsample the curvature radius
- Rk_red = []
- templen = 1000
- for n in range(2,len(Rk)):
- temprk = np.concatenate( (Rk[0:-1:n], [Rk[-1]]) )
- if temprk.size < templen:
- Rk_red.append(temprk)
- templen = temprk.size
- Rk_red.reverse()
- #
- #plt.figure()
- # for n in range(len(dist_red)-1,-1,-1):
- # dist_grid,rk_grid = np.meshgrid(dist_red[n],Rk_red[n])
- # plt.figure()
- # plt.plot(dist_grid.ravel(), rk_grid.ravel(),'*')
- # plt.xlabel("distance z/D")
- # plt.ylabel("Rc")
- # plt.grid()
- # plt.tight_layout()
- # Npoints = str(dist_grid.size)
- # print(Npoints)
- # plt.savefig(save_path + "Grid_red_N_"+ Npoints +".pdf")
- #
- # define a reference interpolator on the finest grid which gives the ground truth (nearest interpolation)
- ref_linint = []
- for i in range(0,Nmodes):
- # pure linear interpolation for the coefficients of each mode
- ref_linint.append(interpol.RegularGridInterpolator((dist,Rk), pod_z["coeffs"][:,:,i].T,bounds_error=False,fill_value=0))
- #
- for n in range(len(dist_red)):
- tempdist = dist_red[n]
- temprk = Rk_red[n]
- dist_grid,rk_grid = np.meshgrid(tempdist,temprk)
- # get the coefficient matrix as reference on the reduced grid
- A_temp = get_POD_coeffs(dist_grid,rk_grid,ref_linint, inttype="lin")
- Nsamples = dist_grid.size
- # Radial basis functions -----------------------------------------------------------
- # run through all possible rbf functions
- rbf_func = ['multiquadric', 'inverse',
- 'gaussian',
- 'linear',
- 'cubic',
- 'quintic',
- 'thin_plate']
- for j in range(0,len(rbf_func)):
- rbfint = []
- for i in range(0,Nmodes):
- rbfint.append(interpol.Rbf(dist_grid,rk_grid,A_temp[:,:,i].T,function = rbf_func[j]))
- A_rbf = get_POD_coeffs(ddist,rrk,rbfint, inttype="rbf")
- np.save(save_path + "rbf_N_" + str(Nsamples) + "_func_" + rbf_func[j] + ".npy" ,A_rbf)
- # do the same for kriging -----------------------------------------------------------
- krig_func = ["linear", "power", "gaussian", "spherical", "exponential"]
- for j in range(0,len(krig_func)):
- krigint = []
- for i in range(0,Nmodes):
- krigint.append(pk.OrdinaryKriging(dist_grid.ravel(), rk_grid.ravel(),A_temp[:,:,i].ravel(), variogram_model=krig_func[j],
- verbose=False, enable_plotting=False))
- A_krig = np.array([krigint_i.execute('grid', dist, Rk)[0].T for krigint_i in krigint]).T
- A_krig_std = np.array([krigint_i.execute('grid', dist, Rk)[1].T for krigint_i in krigint]).T
- np.save(save_path + "kriging_N_" + str(Nsamples) + "_func_" + krig_func[j] + ".npy" ,A_krig)
- np.save(save_path + "kriging_STD_N_" + str(Nsamples) + "_func_" + krig_func[j] + ".npy" ,A_krig_std )
- # for spline interpolation -----------------------------------------------------------
- # do spline only if Npoints >4
- if Nsamples > 12:
- polydeg = [1,2,3]
- for j in range(0,len(polydeg)):
- splineint = []
- for i in range(0,Nmodes):
- splineint.append(interpol.RectBivariateSpline(tempdist,temprk,A_temp[:,:,i].T,kx= polydeg[j],ky=polydeg[j] ))
- A_spline = np.array([splineint_i(dist,Rk) for splineint_i in splineint] ).T
- np.save(save_path + "Spline_N_" + str(Nsamples) + "_polydeg_" + str(polydeg[j]) + ".npy" ,A_spline)
- #
- # linear interpolation -----------------------------------------------------------
- linint = []
- for i in range(0,Nmodes):
- # pure linear interpolation for the coefficients of each mode
- linint.append(interpol.RegularGridInterpolator((tempdist,temprk), A_temp[:,:,i].T,bounds_error=False,fill_value=0))
- A_lin = get_POD_coeffs(ddist,rrk,linint, inttype="lin")
- np.save(save_path + "Linear_N_" + str(Nsamples) + ".npy" ,A_lin)
-
-
-
- # rktest += 0.1
- # disttest +=0.1
- # A_lin = get_POD_coeffs(disttest,rktest,linint,inttype="lin")
- # u_int = do_reconstruct(A_lin.ravel(),pod_z["modes"],A_lin.size,pod_z["umean"] ).ravel()
- # nicecontour(x,y,u_int.ravel())
- # plt.title("linear")
- # #
- # # test with rbf
- # A_rbf = get_POD_coeffs(disttest,rktest,rbfint,inttype="rbf")
- # u_int = do_reconstruct(A_rbf.ravel(),pod_z["modes"],A_rbf.size,pod_z["umean"] ).ravel()
- # nicecontour(x,y,u_int.ravel())
- # plt.title("RBF")
- # # test with spline
- # A_spline = get_POD_coeffs(disttest,rktest,splineint,inttype="rbf")
- # u_int = do_reconstruct(A_spline.ravel(),pod_z["modes"],A_spline.size,pod_z["umean"] ).ravel()
- # nicecontour(x,y,u_int.ravel())
- # plt.title("spline")
- # # test with kriging
- # #A_krig = get_POD_coeffs(disttest,rktest,krigint,inttype="rbf")
- # A_krig = np.array([krigint_i.execute('points', disttest, rktest)[0] for krigint_i in krigint])
- # u_int = do_reconstruct(A_krig.ravel(),pod_z["modes"],A_spline.size,pod_z["umean"] ).ravel()
- # nicecontour(x,y,u_int.ravel())
- # plt.title("kriging")
-
- # plt.figure()
- # plt.plot(A_test,'*-',label= "original")
- # plt.plot(A_lin.ravel(),'+:',label= "lin")
- # plt.plot(A_rbf.ravel(),'<:',label= "rbf")
- # plt.plot(A_spline.ravel(),'>--',label= "spline")
- # plt.grid()
- # plt.legend()
- # plt.tight_layout()
-
-
-
- raise SystemExit(0)
-
- # reconstruct a very fine A matrix
- Nfine = 100
- maxdist = 40
- rkfine = np.linspace(Rk.min(),Rk.max(),Nfine-1)
- distfine = np.linspace(dist.min(),np.sqrt(maxdist),Nfine+1)**2
- #
- distgrid, rkgrid = np.meshgrid(distfine,rkfine)
- # reconstruct the whole A matrix on the fine grid
- A_lin = get_POD_coeffs(distgrid,rkgrid,linint, inttype="lin")
- A_rbf = get_POD_coeffs(distgrid,rkgrid,rbfint, inttype="rbf")
- #A_spline = get_POD_coeffs(distfine,rkfine,splineint,inttype="spline")
- # this works:
- A_spline = np.array([splineint_i(distfine,rkfine) for splineint_i in splineint] ).T
- # kriging
- A_krig = np.array([krigint_i.execute('grid', distfine, rkfine)[0].T for krigint_i in krigint]).T
- A_ss = np.array([krigint_i.execute('grid', distfine, rkfine)[1].T for krigint_i in krigint]).T
- # prepare contour plots for coefficient matrix A
- #mode_i = 0 # the coefficients for the mode that is plotted
- #Atemp = pod_z["coeffs"][:,:,mode_i]
- maxdist_ind = np.where(dist<40)[0][-1] + 1
- # plot for the first 5 modes
- for mode_i in range(0,5):
- fig,ax = plt.subplots(1,5,figsize=(18,5))
- #
- nicecontour(ddist.ravel(),rrk.ravel(), pod_z["coeffs"][0:maxdist_ind,:,mode_i].ravel(), fig=ax[0],polar=False,equalsize=False)
- ax[0].set_title("original")
- #
- cs1 = nicecontour(distgrid.ravel(),rkgrid.ravel(), A_lin[:,:,mode_i].ravel(), fig=ax[1],polar=False,equalsize=False)
- ax[1].set_title("linear")
- #
- cs1 = nicecontour(distgrid.ravel(),rkgrid.ravel(), A_rbf[:,:,mode_i].ravel(), fig=ax[2],polar=False,equalsize=False)
- ax[2].set_title("RBF")
- #
- cs1 = nicecontour(distgrid.ravel(),rkgrid.ravel(), A_spline[:,:,mode_i].ravel(), fig=ax[3],polar=False,equalsize=False)
- ax[3].set_title("spline")
- #
- cs1 = nicecontour(distgrid.ravel(),rkgrid.ravel(), A_krig[:,:,mode_i].ravel(), fig=ax[4],polar=False,equalsize=False)
- ax[4].set_title("kriging")
-
- fig.tight_layout()
-
-
-
-
-# raise SystemExit(0)
-
-# for mode_i in range(0,40):
-# nicecontour(distgrid.ravel(),rkgrid.ravel(), A_spline[:,:,mode_i].ravel() - A_lin[:,:,mode_i].ravel(), polar=False,equalsize=False,withcolorbar=False)
-# plt.savefig("./nice_views/diff_spline_lin"+str(mode_i)+".png")
-# plt.close("all")
-
-
-
-
-# import numpy as np
-# from scipy import ndimage
-# import matplotlib.pyplot as plt
-
-# # Note that the output interpolated coords will be the same dtype as your input
-# # data. If we have an array of ints, and we want floating point precision in
-# # the output interpolated points, we need to cast the array as floats
-# data = np.arange(40).reshape((8,5)).astype(np.float)
-
-# # I'm writing these as row, column pairs for clarity...
-# coords = np.array([[1.2, 3.5], [6.7, 2.5], [7.9, 3.5], [3.5, 3.5]])
-# # However, map_coordinates expects the transpose of this
-# coords = coords.T
-
-# # The "mode" kwarg here just controls how the boundaries are treated
-# # mode='nearest' is _not_ nearest neighbor interpolation, it just uses the
-# # value of the nearest cell if the point lies outside the grid. The default is
-# # to treat the values outside the grid as zero, which can cause some edge
-# # effects if you're interpolating points near the edge
-# # The "order" kwarg controls the order of the splines used. The default is
-# # cubic splines, order=3
-# zi = ndimage.map_coordinates(data, coords, order=3, mode='nearest')
-
-# row, column = coords
-# nrows, ncols = data.shape
-# im = plt.imshow(data, interpolation='nearest', extent=[0, ncols, nrows, 0])
-# plt.colorbar(im)
-# plt.scatter(column, row, c=zi, vmin=data.min(), vmax=data.max())
-# for r, c, z in zip(row, column, zi):
-# plt.annotate('%0.3f' % z, (c,r), xytext=(-10,10), textcoords='offset points',
-# arrowprops=dict(arrowstyle='->'), ha='right')
-# plt.show()
-
-
-
-
-
-
-
-
-
-
-
-def get_corrfunc_const_2(dist):
- #return p[0]*x + p[1] + np.exp(-x**p[3] *p[4])
- #return p0*x + p1 + np.exp(-x**p3 *p4)*p5
- #p = [1.05819337, 1.77467745, 0.00353871, 1.00758248]
- #return p[0] * np.exp(-dist**p[1] *p[2])*p[3]
- p = [-1.67148330e+01, 1.21152710e+00, 2.76845394e+00, 3.88782903e-05]
- return p[1] * np.exp(-(dist-p[0])**p[2] *p[3])
-
-def get_corrfunc_const(dist,corrtype):
- #return p[0]*x + p[1] + np.exp(-x**p[3] *p[4])
- #return p0*x + p1 + np.exp(-x**p3 *p4)*p5
- if corrtype == 0:
- # constant factor
- p = [1.23934497, 0.04419126]
- elif corrtype == 1:
- # mean profile correction
- p = [1.21653678, 0.02861125]
- elif corrtype == 2:
- # deviation profile correction
- p = [0.86976819, 0.04046908]
- return p[0] * np.exp(-dist*p[1])
-
-def get_corrfunc_lin_meanconst(x):
- p = np.array([ 1.16627408, -0.02290368])
- x = x.clip(min=2.43, max=30.82)
- return p[0] + p[1]*x
-
-def get_corrfunc_sec(dist):
- p = np.array([1.02696228, -0.03255838])
- f = p[0] + p[1]*dist
- yind = np.where(f<0)[0] #
- if yind.size > 0:
- f[yind] = 0
- return f
-
-def get_corrfunc_mean_piecewiselin(x):
- fa_mean = np.array([1.0840439208753998, 0.9998665775415876, 0.936508879310944, 0.8761143416100763, 0.7280467671133543, 0.4028894012569246,0])
- dist = np.array([ 2.43, 5.56, 10.7 , 15.69, 20.82, 30.82,100])
- # peacewise interpolation
- pi = interpol.interp1d(dist,fa_mean)
- x = x.clip(min=2.43)
- return pi(x)
-
-
-def get_corrfunc_pathwise_const(dist,phi):
- def func_p_exp(x,p0,p1):
- #return p[0]*x + p[1] + np.exp(-x**p[3] *p[4])
- #return p0*x + p1 + np.exp(-x**p3 *p4)*p5
- return np.exp(-(x-p0)**2 /p1)
- #
- def func_lin(x,p0,p1):
- return p0 + p1*x
- #
- paras = np.array([[ 3.07451078e+00, 5.97263072e+02],
- [-4.64976255e+01, 4.16174476e+03],
- [ 4.49253567e-02, 3.57868450e-03]])
- p_coeffs = [func_p_exp(dist,*paras[0]) , func_p_exp(dist,*paras[1]) , func_lin(dist,*paras[2]) ]
- def func_exp(x,p0,p1,p2):
- return p0 - p1 * np.exp( (-(x/pi+0.5)**2)/(p2**2))
- #
- return func_exp(phi,*p_coeffs)
-
-def get_corrfunc_diff_const(dist):
- p = [ 6.71591469e-01, -1.43380241e+00, 7.15070860e+02]
- return p[0]*np.exp((-(dist + p[1])**2) /p[2])
-#
-def get_corrfunc_mean_const(dist):
- p = [ 1.04279875, -2.53149853, 798.25668037]
- return p[0]*np.exp((-(dist + p[1])**2) /p[2])
-
-
-
-
-def gaussfunc(r,phi,p0,p1,p2,p3,p4):
- p = (p0,p1,p2,p3,p4)
- return p[4] + p[2]*np.exp(-(r-p[3])**2 /p[0])*np.exp(-(phi + pi/2)**2 /p[1])
-
-# correction by multiplication of (1+u) with a exp function
-def get_corrfunc_gauss(r,phi,dist):
- dist = min(dist,31)
- p = [np.array([0.00570207, 0.05981588]),
- np.array([ 9.37053900e-07, -3.35839510e-05, 1.41653099e-03, -7.06934959e-03,
- 4.30957378e-02]),
- np.array([0.15368735]),
- np.array([-0.00838962, 0.83048114]),
- np.array([ -3.79461018e-07, 2.75128296e-05, -7.09316862e-04, 6.93034661e-03,
- 9.69319331e-01])]
- return gaussfunc(r,phi, np.polyval(p[0],dist), np.polyval(p[1],dist), np.polyval(p[2],dist), np.polyval(p[3],dist), np.polyval(p[4],dist))
-
-
-
-
diff --git a/HelpingFunctions/Ultrasonic_eval.py b/HelpingFunctions/Ultrasonic_eval.py
deleted file mode 100644
index 89d7aeda4e3408bcb785e0a13203dfffed76a7e6..0000000000000000000000000000000000000000
--- a/HelpingFunctions/Ultrasonic_eval.py
+++ /dev/null
@@ -1,453 +0,0 @@
-# -*- coding: utf-8 -*-
-"""
-Created on Wed Feb 3 15:17:28 2021
-
-@author: Ichyc
-"""
-
-import json
-import glob
-import os
-import numpy as np
-import matplotlib.pyplot as plt
-import matplotlib
-import scipy.interpolate as interpol
-import pandas as pd
-from matplotlib.ticker import MaxNLocator
-import sys
-# the mock-0.3.1 dir contains testcase.py, testutils.py & mock.py
-sys.path.append('E:/NextCloud/GlobalPyFunc/')
-from niceplot import *
-from Integrate2D import Tri_Int
-from math import pi
-from GetElbowProfile import GetElbowProfile, get_corrfunc_const
-from load_measurement import load_measurement
-
-# ------------------------------------------------------------
-# ------------------------------------------------------------
-# ------------------------------------------------------------
-
-# define the fonts
-# SMALL_SIZE = 8
-# MEDIUM_SIZE = 16
-# BIGGER_SIZE = 30
-
-# plt.rc('font', size=MEDIUM_SIZE) # controls default text sizes
-# plt.rc('axes', titlesize=MEDIUM_SIZE) # fontsize of the axes title
-# plt.rc('axes', labelsize=MEDIUM_SIZE) # fontsize of the x and y labels
-# plt.rc('xtick', labelsize=MEDIUM_SIZE) # fontsize of the tick labels
-# plt.rc('ytick', labelsize=MEDIUM_SIZE) # fontsize of the tick labels
-# plt.rc('legend', fontsize=MEDIUM_SIZE) # legend fontsize
-# plt.rc('figure', titlesize=BIGGER_SIZE) # fontsize of the figure title
-# # set the linewidth value globally
-# # show all the keys: plt.rcParams.keys()
-# #plt.rcParams['axes.linewidth'] = 1.5
-# plt.rcParams['lines.linewidth'] = 4
-# plt.rcParams['lines.markersize'] = 8
-# plt.rcParams['lines.markeredgewidth'] = 2
-# ------------------------------------------------------------
-# ------------------------------------------------------------
-# ------------------------------------------------------------
-
-
-
-new_rk = 1.425
-
-# ------------------------------------------------------------
-# load the coordinates of the simulation results
-rootdir = os.path.abspath("../")
-data_path= os.path.join(rootdir,"All_DistPlanes")
-with open(os.path.join(data_path,"Dist_Coords.json")) as json_file:
- coords = json.load(json_file)
-D = coords["D"]
-dist_s = np.array(coords["dist"])/D
-# the courvature radius
-Rk = np.array(coords["Rk"])/D/1000
-r_s = np.array(coords["r"])
-phi_s = np.array(coords["phi"])
-# remove the last value of the matrix this is where r = 1
-#r_s = r_s[:,0:-1]
-#phi = phi[:,0:-1]
-# calc the x and y values
-x_s,y_s = pol2cart(r_s.ravel(),phi_s.ravel())
-# -------------------------------------------------------------
-
-# the measurement distances
-dist_m = np.array([0.243 , 0.556, 1.07, 1.569, 2.082, 3.082])/D
-u_m = [0] * len(dist_m)
-# load all measurements
-for i in range(0,len(dist_m)):
- r_m, phi_m, z, u_m[i] = load_measurement(i, diff = True)
-x_m,y_m = pol2cart(r_m.ravel(), phi_m.ravel())
-minval = -0.3
-maxval = 0.3
-
-x,y, u_nocor = GetElbowProfile(new_rk, dist_m, apply_correction = 0 )
-x,y, u_cor_const = GetElbowProfile(new_rk, dist_m, apply_correction = 1 )
-x,y, u_cor_diff = GetElbowProfile(new_rk, dist_m, apply_correction = 3 )
-#
-u_nocor = u_nocor[0]
-u_cor_const = u_cor_const[0]
-u_cor_diff = u_cor_diff[0]
-#
-
-'''
-fig1, axs = plt.subplots(4, 6)
-fig1.set_size_inches(18, 12)
-for i,dist_i in enumerate(dist_m):
- new_dist = dist_i
- nicecontour(x_m, y_m, u_m[i].ravel(),minval = minval,maxval=maxval,fig = axs[0,i],withcolorbar=False )
- nicecontour(x, y, u_nocor[i,:],minval = minval,maxval=maxval,fig = axs[1,i],withcolorbar=False )
- nicecontour(x, y, u_cor_const[i,:],minval = minval,maxval=maxval,fig = axs[2,i],withcolorbar=False )
- nicecontour(x, y, u_cor_diff[i,:],minval = minval,maxval=maxval,fig = axs[3,i],withcolorbar=False )
- for j in range(0,4):
- axs[j,i].axis("equal")
- axs[j,i].set_xlim([-1, 1])
- axs[j,i].set_ylim([-1, 1])
- #
-#
-#fig1.colorbar(cs)
-#plt.tight_layout()
-plt.show()
-'''
-
-def make_whole_pathmat(r,phi,u):
- Npaths = u.shape[0]
- # check if the number is even
- if Npaths % 2:
- # its Odd, then delete the last
- u = u[0:-1,:]
- phi = phi[0:-1,:]
- r = r[0:-1,:]
- Npaths -=1
- M = int(Npaths/2)
- new_u = np.append( np.flip(u[0:M,1:],axis=1) , u[M:,:] ,axis=1)
- new_r = np.append( np.flip(-r[0:M,1:],axis=1) , r[M:,:] ,axis=1)
- new_phi = np.append( np.flip(phi[0:M,1:],axis=1), phi[M:,:],axis=1)
- #
- return new_r, new_phi, new_u
-#
-def path_int(r,u,N = int(1e4)):
- # integrate with akima spline interpolation
- aki = interpol.Akima1DInterpolator(r,u)
- rfine = np.linspace(min(r),max(r),N)
- return np.trapz(aki(rfine),rfine)/2
-#
-N_dist = len(dist_m)
-u_mg = [0] * N_dist
-u_sg_cor_const = [0] * N_dist
-u_sg_nocor = [0] * N_dist
-u_sg_cor_pathconst = [0] * N_dist
-#
-us_m_mean = np.zeros(N_dist)
-us_m_std = np.zeros(N_dist)
-#
-us_s_mean_nocor = np.zeros(N_dist)
-us_s_std_nocor = np.zeros(N_dist)
-#
-us_s_mean_cor_const = np.zeros(N_dist)
-us_s_std_cor_const = np.zeros(N_dist)
-#
-us_s_mean_cor_diff = np.zeros(N_dist)
-us_s_std_cor_diff = np.zeros(N_dist)
-
-
-for j,dist_i in enumerate(dist_m):
- r_mg , phi_mg, u_mg[j] = make_whole_pathmat(r_m,phi_m,u_m[j])
- r_sg , phi_sg, u_sg_nocor[j] = make_whole_pathmat(r_s,phi_s,u_nocor[j].reshape(r_s.shape))
- r_sg , phi_sg, u_sg_cor_const[j] = make_whole_pathmat(r_s,phi_s,u_cor_const[j].reshape(r_s.shape))
- r_sg , phi_sg, u_sg_cor_pathconst[j] = make_whole_pathmat(r_s,phi_s,u_cor_diff[j].reshape(r_s.shape))
- phi1_m = phi_mg[:,-1]
- phi1_s = phi_sg[:,-1]
- #
- #fig_lines, axl = axs = plt.subplots(1)
- us_m = np.zeros(phi1_m.shape)
- us_nocor = np.zeros(phi1_s.shape)
- us_cor_const = np.zeros(phi1_s.shape)
- us_cor_diff = np.zeros(phi1_s.shape)
- for i in range(0,r_mg.shape[0]):
- #axl.plot(r_mg[i,:],u_mg[i,:],'+:',label=str(np.round(phi_mg[i,-1]/pi ,decimals= 2)))
- us_m[i] = path_int(r_mg[i,:],u_mg[j][i,:])*100
- for i in range(0,r_sg.shape[0]):
- us_nocor[i] = path_int(r_sg[i,:], u_sg_nocor[j][i,:])*100
- us_cor_const[i] = path_int(r_sg[i,:], u_sg_cor_const[j][i,:])*100
- us_cor_diff[i] = path_int(r_sg[i,:], u_sg_cor_pathconst[j][i,:])*100
- #
- us_m_mean[j] = np.mean(us_m)
- us_m_std[j] = np.std(us_m)
- #
- us_s_mean_nocor[j] = np.mean(us_nocor)
- us_s_std_nocor[j] = np.std(us_nocor)
- #
- us_s_mean_cor_const[j] = np.mean(us_cor_const)
- us_s_std_cor_const[j] = np.std(us_cor_const)
- #
- us_s_mean_cor_diff[j] = np.mean(us_cor_diff)
- us_s_std_cor_diff[j] = np.std(us_cor_diff)
-
- #axl.grid("minor")
- #fig_lines.legend()
- #
- plt.figure()
- #us_m = us_m/2*100
- plt.plot(phi1_s,us_nocor,'g-',label="Simu")
- plt.plot(phi1_m,us_m,'r-',label="LDV")
- plt.plot(phi1_s,us_cor_const,'-',label="Simu cor const")
- plt.plot(phi1_s,us_cor_diff,'-',label="Simu cor mean_dev")
- plt.ylim(-19,0)
- plt.legend()
- plt.grid("minor")
- plt.title(str(np.round(dist_i)))
- plt.ylabel("Meter Error in %")
- plt.xlabel("Angular Position")
- plt.tight_layout()
- #plt.savefig("../Bilder_Ultrasonic/Error_over_angle_pos_" + str(j) + ".pdf")
- i_0 = np.where(phi1_m == 0)
-
-#
-plt.figure()
-plt.plot(dist_m,us_m_mean,'o:',label="LDV")
-plt.plot(dist_m,us_s_mean_nocor,'o:',label="Simu")
-plt.plot(dist_m,us_s_mean_cor_const,'o:',label="Simu cor const")
-plt.plot(dist_m,us_s_mean_cor_diff,'o:',label="cor diff const")
-plt.grid("minor")
-plt.legend()
-plt.tight_layout()
-
-# make the same with finer dist ------------------------------------------------
-
-
-def save_us_errors(new_rk,dist_s,r_s,phi_s):
- dist_s = dist_s[0:-1]
- x,y, u_nocor = GetElbowProfile(new_rk, dist_s, apply_correction = 0 )
- x,y, u_cor_const = GetElbowProfile(new_rk, dist_s, apply_correction = 1 )
- x,y, u_cor_diff = GetElbowProfile(new_rk, dist_s, apply_correction = 3 )
- #
- u_nocor = u_nocor[0]
- u_cor_const = u_cor_const[0]
- u_cor_diff =u_cor_diff[0]
- #
- N_dist = len(dist_s)
- u_mg = [0] * N_dist
- u_sg_cor_const = [0] * N_dist
- u_sg_nocor = [0] * N_dist
- u_sg_cor_pathconst = [0] * N_dist
- #
- us_s_mean_nocor = np.zeros(N_dist)
- us_s_std_nocor = np.zeros(N_dist)
- #
- us_s_mean_cor_const = np.zeros(N_dist)
- us_s_std_cor_const = np.zeros(N_dist)
- #
- us_s_mean_cor_diff = np.zeros(N_dist)
- us_s_std_cor_diff = np.zeros(N_dist)
-
-
- for j,dist_i in enumerate(dist_s):
- r_sg , phi_sg, u_sg_nocor[j] = make_whole_pathmat(r_s,phi_s,u_nocor[j].reshape(r_s.shape))
- r_sg , phi_sg, u_sg_cor_const[j] = make_whole_pathmat(r_s,phi_s,u_cor_const[j].reshape(r_s.shape))
- r_sg , phi_sg, u_sg_cor_pathconst[j] = make_whole_pathmat(r_s,phi_s,u_cor_diff[j].reshape(r_s.shape))
- phi1_s = phi_sg[:,-1]
- #
- #fig_lines, axl = axs = plt.subplots(1)
- us_nocor = np.zeros(phi1_s.shape)
- us_cor_const = np.zeros(phi1_s.shape)
- us_cor_diff = np.zeros(phi1_s.shape)
- # run through the angles
- for i in range(0,r_sg.shape[0]):
- us_nocor[i] = path_int(r_sg[i,:], u_sg_nocor[j][i,:])*100
- us_cor_const[i] = path_int(r_sg[i,:], u_sg_cor_const[j][i,:])*100
- us_cor_diff[i] = path_int(r_sg[i,:], u_sg_cor_pathconst[j][i,:])*100
- #
- '''
- mean_int = interpol.Akima1DInterpolator(dist_mod,us_mean)
- us_nocor[i] = np.trapz(u_sg_nocor[j][i,:],r_sg[i,:])/2*100
- us_cor_const[i] = np.trapz(u_sg_cor_const[j][i,:],r_sg[i,:])/2*100
- us_cor_diff[i] = np.trapz(u_sg_cor_pathconst[j][i,:],r_sg[i,:])/2*100
- '''
- #
- us_s_mean_nocor[j] = np.mean(us_nocor)
- us_s_std_nocor[j] = np.std(us_nocor)
- #
- us_s_mean_cor_const[j] = np.mean(us_cor_const)
- us_s_std_cor_const[j] = np.std(us_cor_const)
- #
- us_s_mean_cor_diff[j] = np.mean(us_cor_diff)
- us_s_std_cor_diff[j] = np.std(us_cor_diff)
-
- #axl.grid("minor")
- #fig_lines.legend()
- #
- '''
- plt.figure()
- #us_m = us_m/2*100
- plt.plot(phi1_s,us_nocor,'+:',label="Simu no cor")
- plt.plot(phi1_s,us_cor_const,'+:',label="Simu cor const")
- plt.plot(phi1_s,us_cor_diff,'+:',label="Simu cor pathwise")
- plt.ylim(-19,0)
- plt.legend()
- plt.grid("minor")
- plt.title(str(np.round(dist_i)))
- #i_0 = np.where(phi1_m == 0)
- '''
- return [us_s_mean_nocor,us_s_std_nocor] ,[us_s_mean_cor_const,us_s_std_cor_const],[us_s_mean_cor_diff,us_s_std_cor_diff]
-
-#Rk_mod = Rk[0:-1]
-
-# all the values are calculated if calc == True
-calc = False
-save_or_not = False
-if calc:
- all_errors_simu = [0] * len(Rk)
- all_errors_const = [0] * len(Rk)
- all_errors_const_diff = [0] * len(Rk)
-
-
- for i,rki in enumerate(Rk):
- all_errors_simu[i],all_errors_const[i],all_errors_const_diff[i] = save_us_errors(rki,dist_s,r_s,phi_s)
-
- np.savez('US_errors.npz',np.array(all_errors_simu))
- np.savez('US_errors_corrconst.npz',np.array(all_errors_const))
- np.savez('US_errors_const_diff.npz',np.array(all_errors_const_diff))
-else:
- # otherwise the values are loaded
- with np.load('US_errors.npz') as data:
- all_errors_simu = data["arr_0"] #.get("arr_0")
- with np.load('US_errors_corrconst.npz') as data:
- all_errors_const = data["arr_0"] #.get("arr_0")
- with np.load('US_errors_corrpathconst.npz') as data:
- all_errors_pathconst = data["arr_0"] #.get("arr_0")
- with np.load('US_errors_const_diff.npz') as data:
- all_errors_const_diff = data["arr_0"] #.get("arr_0")
-
-
-
-
-which_cor = 3
-# which correction is applied?
-if which_cor == 0:
- # no correction
- us_mean = all_errors_simu[:,0,:].T
- us_std = all_errors_simu[:,1,:].T
-elif which_cor == 1:
- # const correction
- us_mean = all_errors_const[:,0,:].T
- us_std = all_errors_const[:,1,:].T
-elif which_cor == 2:
- # const path correction
- us_mean = all_errors_pathconst[:,0,:].T
- us_std = all_errors_pathconst[:,1,:].T
-
-elif which_cor == 3:
- # const diff correction
- us_mean = all_errors_const_diff[:,0,:].T
- us_std = all_errors_const_diff[:,1,:].T
-
-dist_mod = dist_s[0:-1]
-mean_int = interpol.Akima1DInterpolator(dist_mod,us_mean)
-std_int = interpol.Akima1DInterpolator(dist_mod,us_std)
-
-#mean_int = interpol.interp1d( dist_mod,us_mean )
-#std_int = interpol.interp1d( dist_mod,us_std )
-
-fig_bars, ax_bars = plt.subplots(2, 3)
-fig_bars.set_size_inches(18, 9)
-fig_bars.tight_layout()
-ax_bars = ax_bars.ravel()
-
-#
-fig_fil, ax_fil = plt.subplots(2, 3)
-fig_fil.set_size_inches(18, 9)
-fig_fil.tight_layout()
-ax_fil = ax_fil.ravel()
-for i, disti in enumerate(dist_m):
-
-
- mean_errors = mean_int(disti)
- mean_std = std_int(disti)
- #mean_errors_10 = [bla[0][distN] for bla in all_errors_simu]
- #mean_std_10 = [bla[1][distN] for bla in all_errors_simu]
- #
- ax_bars[i].errorbar(Rk,mean_errors,yerr=mean_std,
- capsize=8,
- color='g',
- ecolor='y',
- #markerfacecolor='g',
- fmt='-',
- #markeredgewidth=2,
- label="simu")
- ax_bars[i].errorbar(new_rk,us_m_mean[i],yerr=us_m_std[i],fmt='ro',capsize=10,label="LDV")
- ax_bars[i].set_title("distance "+str(round(disti))+" D")
- ax_bars[i].set_xlim(0,10)
- #ax_bars[i].set_ylim(-17,0.5)
- ax_bars[i].grid("minor")
- if i>2: ax_bars[i].set_xlabel("Kurvature Radius in D")
- ax_bars[i].set_ylabel("Meter error in %")
-
-
-
- #with filled areas
- ax_fil[i].plot(Rk,mean_errors,'g-', label="mean_simu")
- ax_fil[i].fill_between(Rk, mean_errors-mean_std, mean_errors+mean_std,
- facecolor='green',
- alpha = 0.3,
- interpolate=True, label="std simu")
-
- ax_fil[i].errorbar(new_rk,us_m_mean[i],yerr=us_m_std[i],fmt='ro',capsize=10,label="LDV")
- ax_fil[i].set_title("distance "+str(round(disti))+" D")
- ax_fil[i].set_xlim(0,10)
- ax_fil[i].set_ylim(-20,0.5)
- ax_fil[i].grid("minor")
- if i>2: ax_fil[i].set_xlabel("Kurvature Radius in D")
- ax_fil[i].set_ylabel("Meter error in %")
-
-ax_bars[-1].legend(loc='lower right')
-ax_fil[-1].legend(loc='lower right')
-fig_fil.tight_layout()
-if save_or_not: fig_fil.savefig("../Bilder_Ultrasonic/Mean_Error_over_Rk_Piecewlincorr_" + str(which_cor) + ".pdf")
-fig_bars.tight_layout()
-if save_or_not: fig_bars.savefig("../Bilder_Ultrasonic/Mean_Error_over_Rk_Bar_Piecewlincorr_" + str(which_cor) + ".pdf")
-#
-plt.figure()
-
-plt.plot(dist_mod,us_mean[:,7],'g-', label="mean_simu")
-plt.fill_between(dist_mod, us_mean[:,7]-us_std[:,7], us_mean[:,7]+us_std[:,7],
- facecolor='green',
- alpha = 0.3,
- interpolate=True, label="std simu")
-
-# plt.errorbar(dist_mod,all_errors_simu[7,0,:],yerr=all_errors_simu[7,1,:],fmt='-', label="Simu")
-#plt.errorbar(dist_s,us_s_mean_cor_const,yerr=us_s_std_cor_const,fmt='-',label="Simu cor const")
-#plt.errorbar(dist_s,us_s_mean_cor_diff,yerr=us_s_std_cor_diff,fmt='-',label="Simu cor path const")
-plt.errorbar(dist_m,us_m_mean,yerr=us_m_std,fmt='ro:',capsize=10,label="LDV")
-plt.xlim(0,60)
-plt.ylim(-19,1)
-plt.grid("minor")
-plt.xlabel("Distance in Diameter")
-plt.ylabel("Meter error in %")
-plt.legend()
-plt.tight_layout()
-if save_or_not: plt.savefig("../Bilder_Ultrasonic/Mean_Error_over_dist_Piecewlincorr_" + str(which_cor) + ".pdf")
-
-
-# plot a 2-D contour plot for the error and the uncertainty
-end_ind = 64
-dd, rr = np.meshgrid(dist_mod[0:end_ind],Rk)
-nicecontour(dd.ravel(),rr.ravel(),(us_mean[0:end_ind,:].T).ravel(),maxval = 0, minval = -19 ,Nbins = 100,polar=False,withcontours = False, equalsize = False)
-plt.xlabel("Distance in D")
-plt.ylabel("Curvature Radius in D")
-plt.tight_layout()
-#plt.savefig(os.path.join('./', 'Numplot' + charnames[i] + '.png'))
-#
-# contour of the uncertainty
-end_ind = 64
-dd, rr = np.meshgrid(dist_mod[0:end_ind],Rk)
-nicecontour(dd.ravel(),rr.ravel(),(us_std[0:end_ind,:].T).ravel(),maxval = 10, minval = 0 ,Nbins = 100,polar=False,withcontours = False, equalsize = False)
-plt.xlabel("Distance in D")
-plt.ylabel("Curvature Radius in D")
-plt.tight_layout()
-
-
-
-
-
diff --git a/HelpingFunctions/getUSpath.py b/HelpingFunctions/getUSpath.py
deleted file mode 100644
index 4b41a847fd811e4e9723748d780514c67ce72ac1..0000000000000000000000000000000000000000
--- a/HelpingFunctions/getUSpath.py
+++ /dev/null
@@ -1,221 +0,0 @@
-# -*- coding: utf-8 -*-
-"""
-Created on Mon Sep 20 11:44:08 2021
-
-@author: weisse02
-"""
-import json
-#import glob
-#import os
-import numpy as np
-import matplotlib.pyplot as plt
-#import matplotlib
-import scipy.interpolate as interpol
-#import pandas as pd
-#from matplotlib.ticker import MaxNLocator
-#import sys
-# the mock-0.3.1 dir contains testcase.py, testutils.py & mock.py
-#sys.path.append('E:/NextCloud/GlobalPyFunc/')
-from niceplot import *
-#from Integrate2D import Tri_Int , get_L_norm
-from math import pi
-#sys.path.append('./Auswertung/')
-#from get_fully_SA import get_fully_SA
-#from GetElbowProfile import GetElbowProfile, get_corrfunc_const
-
-
-
-def set_axes_equal(ax: plt.Axes):
- """Set 3D plot axes to equal scale.
-
- Make axes of 3D plot have equal scale so that spheres appear as
- spheres and cubes as cubes. Required since `ax.axis('equal')`
- and `ax.set_aspect('equal')` don't work on 3D.
- """
- limits = np.array([
- ax.get_xlim3d(),
- ax.get_ylim3d(),
- ax.get_zlim3d(),
- ])
- origin = np.mean(limits, axis=1)
- radius = 0.5 * np.max(np.abs(limits[:, 1] - limits[:, 0]))
- _set_axes_radius(ax, origin, radius)
-
-def _set_axes_radius(ax, origin, radius):
- x, y, z = origin
- ax.set_xlim3d([x - radius, x + radius])
- ax.set_ylim3d([y - radius, y + radius])
- ax.set_zlim3d([z - radius, z + radius])
-
-
-
-# constructs the ultrasonic path for a given angle phi, alpha and type
-def get_US_points(pathtype = "diametrical", reflections = 0, alpha = pi/4, phi = 0,dz=0,R=1,Npoints = 100):
- # pathtype = "diametrical" means throuw the pipe center
- # or nondiametrical "non"
- # reflections = 0 for direct path, 1 for V-path, 2 for W-path
- # alpha is the axial angle
- # phi is the circumferential angle
- # dz is the axial position - optional
- # R is the pipe radius
- # Npoints = number of points that are samples along the path
- #
- D = 2*R
- # all points at the beginning, end reflection points of the paths
- p = np.zeros((reflections+2,3))
- # p0[3] is always dz
- p[0] = np.array([ R*np.cos(phi), R*np.sin(phi),dz])
- # direction vector
- #rv = np.array([-p[0,0],-p[0,1],np.sin(alpha) ])
- # the length of the whole path with all reflections
- L = 0
- # thepoints along the path on that it will be interpolated
- p_points = np.zeros((reflections+1,3,Npoints))
- # the direction vectors for each reflection
- rv = np.zeros((reflections+1,3))
- for i in range(1,reflections+2):
- Li = D/np.cos(alpha)
- p[i] = np.array([-p[i-1,0],-p[i-1,1], p[i-1,2] + Li* np.sin(alpha) ])
- L += Li
- # sample points along the line
- rvtemp = p[i]-p[i-1]
- #
- # p_points[i-1] = np.array([p[i-1,j] + np.linspace(0,1,Npoints)*rvtemp[j] for j in range(0,3) ])
- # without loop
- p_points[i-1] = (np.tensordot(np.linspace(0,1,Npoints),rvtemp,axes=0) + p[i-1]).T
- #np.array([p[i-1,j] + np.linspace(0,1,Npoints)*rvtemp[j] for j in range(0,3) ])
-
- rv[i-1] = rvtemp/np.linalg.norm(rvtemp)
-
- # normalize the z-values on diameter not on radius -scale with D
- p[:,2] /=D
- p_points[:,2,:] /=D
- # plot it nicely
- if 1:
- fig = plt.figure()
- ax = plt.axes(projection='3d')
- # plot circle
- fphi = np.linspace(0,2*pi,100)
- fx,fy = pol2cart(R, fphi)
- ax.plot(fx,fy,':')
- ax.plot3D(fx,fy,p[-1,2]*np.ones(len(fx)),':')
- ax.plot3D(p[:,0],p[:,1],p[:,2],'*-')
- #
- ax.scatter3D(p_points[:,0,:],p_points[:,1,:],p_points[:,2,:],c='r',marker='.',depthshade = False)
- #ax.plot3D(p[1,0],p[1,1],p[1,2],'*-')
- ax.set_xlabel("x")
- ax.set_ylabel("y")
- ax.set_zlabel("z")
- #ax.set_box_aspect((1,1,1))
- set_axes_equal(ax)
- #ax = fig.gca(projection='3d')
- #ax.set_aspect('equal')
- #ax.set_aspect('auto')
- #ax.axis("equal")
- fig.tight_layout()
- plt.show()
-
- return p_points, p , rvtemp
-
-
-if 0:
-
- # sample from elbow profiles
- pathtype = "diametrical";
- reflections = 3;
- alpha = 30/180*pi;
- phi = 0/180*pi;
- dz=0; # it is in D
- R=1;
- Npoints = 1000
-
- with open("../All_DistPlanes/Dist_Coords.json") as json_file:
- coords = json.load(json_file)
- #
- alldist = np.array(coords["dist"])/coords["D"]*2 # in Radius not in Diameters!!
- Rk = np.array(coords["Rk"])/1000/coords["D"]
- alldist = alldist[0:-8]
- #
- US_measure = [0]*len(alldist)
- for k,idist in enumerate(alldist):
- print(k)
- print(idist)
- p_points, p, rv = get_US_points(pathtype = pathtype, reflections = reflections, alpha = alpha, phi = phi,dz=idist,Npoints = 20)
- maxz = np.max(p[:,2])
- minz = np.min(p[:,2])
- #
- a = np.where(alldist<=minz)
- a = a[-1][-1]
- b = np.where(alldist>=maxz)
- b = b[0][0]
- #
- dist = alldist[a:b+1].copy()
- # GetElbowProfile(new_Rk, new_dist, new_x = None, new_y = None, crossFlow = False, delete_last = False , apply_correction = 0, allply_corr_sec = 0, Re = 5e4, with_fully = False ):
- #for i,rk in enumerate(Rk):
- #Rk = [1.425]
- US_measure[k] = [0]*len(Rk)
- for i,rk in enumerate(Rk):
- if i == 0: rk = np.round(rk,2)
- x,y, u = GetElbowProfile(rk,dist/2,crossFlow = True, apply_correction = 3, allply_corr_sec = 1 )
- # try RegularGridInterpolator .. with regular polar grid
- #
- # construct a regular grid:
- r0 = np.array(coords["r"][0])
- phi0 = np.array([phii[0] for phii in coords["phi"]])
- z0 = dist
- rr,pphi,zz = np.meshgrid(r0,phi0,z0) # np.meshgrid(x, y, z, indexing='ij', sparse=True)
- ux = u[0].T.reshape((len(phi0),len(r0),len(z0)))
- uy = u[1].T.reshape((len(phi0),len(r0),len(z0)))
- uz = u[2].T.reshape((len(phi0),len(r0),len(z0)))
- # the interpolators
- I_reg_x = interpol.RegularGridInterpolator((phi0,r0,z0),ux)
- I_reg_y = interpol.RegularGridInterpolator((phi0,r0,z0),uy)
- I_reg_z = interpol.RegularGridInterpolator((phi0,r0,z0),uz)
- #
- all_phi = phi0[0:-1].copy()
- US_measure[k][i] = [0]*len(all_phi)
- for j,phi_j in enumerate(all_phi):
- p_points, p, rv = get_US_points(pathtype = pathtype, reflections = reflections, alpha = alpha, phi = phi_j,dz=idist,Npoints = Npoints)
- points_r,points_phi = cart2pol(p_points[:,0,:], p_points[:,1,:])
- # perform interpolation
- u_path_x = I_reg_x((points_phi,points_r,p_points[:,2,:]))
- u_path_y = I_reg_y((points_phi,points_r,p_points[:,2,:]))
- u_path_z = I_reg_z((points_phi,points_r,p_points[:,2,:]))
-
- # calc the mean of the reflected paths
- if len(u_path_x.shape)>1:
- u_path_x = np.mean(u_path_x.T * rv[:,0],axis=1)
- u_path_y = np.mean(u_path_y.T * rv[:,1],axis=1)
- u_path_z = np.mean(u_path_z.T * rv[:,2],axis=1)
- ez = np.mean(rv[:,2])
- else:
- u_path_x *= rv[0]
- u_path_x *= rv[1]
- u_path_x *= rv[2]
- ez = rv[2]
- #
- u_path = u_path_x + u_path_y + u_path_z
- US_measure[k][i][j] = 1/ez * np.trapz(u_path,np.linspace(0,1,Npoints)) # np.mean(u_path)
- # end loop over phi
- # end loop over RK
- # end loop over distance
-
- #
- # plt.figure()
- # plt.plot(u_path.T,'*:')
- US_measure = np.array(US_measure)*100
-
-
- np.savez('All_US_errors_alpha_'+str(np.round(alpha/pi*180))+'_refl_'+str(int(reflections))+'_corr_3.npz',US_measure)
-
- US_measure_1 = US_measure[:,7,:]
- us_mean = np.mean(US_measure_1,axis=1)
- us_std = np.std(US_measure_1,axis=1)
- #alldist = alldist[0:-1]
- plt.figure()
- plt.plot(alldist/2,us_mean,'-')
- plt.fill_between(alldist/2, us_mean-us_std, us_mean+us_std,
- facecolor='green',
- alpha = 0.3,
- interpolate=True, label="std simu")
- plt.grid("minor")
\ No newline at end of file
diff --git a/HelpingFunctions/get_fully_SA.py b/HelpingFunctions/get_fully_SA.py
index 5318dcf9dea69eb1efef0bf379bc6533a3f739e4..9d41d1ebd53590640bf25f2bba0478330ac015b2 100644
--- a/HelpingFunctions/get_fully_SA.py
+++ b/HelpingFunctions/get_fully_SA.py
@@ -17,7 +17,7 @@ def get_fully_SA(r):
R = 0.1/2
uvol = 0.4
# get fully developed flow profile from Spalart-Allmaras Simulation
- fullyfile = os.path.abspath("E:/NextCloud/GlobalPyFunc/ProfileLine_SA.csv")
+ fullyfile = os.path.abspath("./HelpingFunctions/ProfileLine_SA.csv")
fullydata = pd.read_csv(fullyfile)
fullydata = fullydata.rename(columns={"U:0": "U_0", "U:1": "U_1" , "U:2": "U_2","Points:0": "x","Points:1": "y","Points:2": "z"})
u_fully_o = fullydata.U_2 / uvol
diff --git a/Screenshot_GUI.png b/Screenshot_GUI.png
new file mode 100644
index 0000000000000000000000000000000000000000..cc4013c0b9ef838f5174b9e5cf07e56adfcc5d2b
Binary files /dev/null and b/Screenshot_GUI.png differ
diff --git a/Ultrasonic_eval.py b/Ultrasonic_eval.py
new file mode 100644
index 0000000000000000000000000000000000000000..348a92e5af97fd7b15e92072f1bd5e13a73cadd2
--- /dev/null
+++ b/Ultrasonic_eval.py
@@ -0,0 +1,360 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Feb 3 15:17:28 2021
+
+@author: Ichyc
+"""
+
+import json
+import glob
+import os
+import numpy as np
+import matplotlib.pyplot as plt
+import matplotlib
+import scipy.interpolate as interpol
+import pandas as pd
+from matplotlib.ticker import MaxNLocator
+import sys
+# the mock-0.3.1 dir contains testcase.py, testutils.py & mock.py
+sys.path.append('E:/NextCloud/GlobalPyFunc/')
+from niceplot import *
+from Integrate2D import Tri_Int
+from math import pi
+from GetElbowProfile import GetElbowProfile, get_corrfunc_const
+from load_measurement import load_measurement
+import Flo_class as flowclass
+# ------------------------------------------------------------
+# ------------------------------------------------------------
+# ------------------------------------------------------------
+
+def make_whole_pathmat(r,phi,u):
+ Npaths = u.shape[0]
+ # check if the number is even
+ if Npaths % 2:
+ # its Odd, then delete the last
+ u = u[0:-1,:]
+ phi = phi[0:-1,:]
+ r = r[0:-1,:]
+ Npaths -=1
+ M = int(Npaths/2)
+ new_u = np.append( np.flip(u[0:M,1:],axis=1) , u[M:,:] ,axis=1)
+ new_r = np.append( np.flip(-r[0:M,1:],axis=1) , r[M:,:] ,axis=1)
+ new_phi = np.append( np.flip(phi[0:M,1:],axis=1), phi[M:,:],axis=1)
+ #
+ return new_r, new_phi, new_u
+#
+def path_int(r,u,N = int(1e4)):
+ # integrate with akima spline interpolation
+ aki = interpol.Akima1DInterpolator(r,u)
+ rfine = np.linspace(min(r),max(r),N)
+ return np.trapz(aki(rfine),rfine)/2
+#
+
+# ----------------------------------------------
+
+
+def save_us_errors(new_rk,dist_s,r_s,phi_s):
+ dist_s = dist_s[0:-1]
+ u_nocor = GetElbowProfile(new_rk, dist_s, apply_correction = 0 )
+ #
+ u_nocor = u_nocor[0]
+ #
+ N_dist = len(dist_s)
+ u_mg = [0] * N_dist
+ u_sg_nocor = [0] * N_dist
+ #
+ us_s_mean_nocor = np.zeros(N_dist)
+ us_s_std_nocor = np.zeros(N_dist)
+ #
+ for j,dist_i in enumerate(dist_s):
+ r_sg , phi_sg, u_sg_nocor[j] = make_whole_pathmat(r_s,phi_s,u_nocor[j].reshape(r_s.shape))
+ phi1_s = phi_sg[:,-1]
+ #
+ #fig_lines, axl = axs = plt.subplots(1)
+ us_nocor = np.zeros(phi1_s.shape)
+ us_cor_const = np.zeros(phi1_s.shape)
+ us_cor_diff = np.zeros(phi1_s.shape)
+ # run through the angles
+ for i in range(0,r_sg.shape[0]):
+ us_nocor[i] = path_int(r_sg[i,:], u_sg_nocor[j][i,:])*100
+ us_cor_const[i] = path_int(r_sg[i,:], u_sg_cor_const[j][i,:])*100
+ us_cor_diff[i] = path_int(r_sg[i,:], u_sg_cor_pathconst[j][i,:])*100
+ #
+ '''
+ mean_int = interpol.Akima1DInterpolator(dist_mod,us_mean)
+ us_nocor[i] = np.trapz(u_sg_nocor[j][i,:],r_sg[i,:])/2*100
+ us_cor_const[i] = np.trapz(u_sg_cor_const[j][i,:],r_sg[i,:])/2*100
+ us_cor_diff[i] = np.trapz(u_sg_cor_pathconst[j][i,:],r_sg[i,:])/2*100
+ '''
+ #
+ us_s_mean_nocor[j] = np.mean(us_nocor)
+ us_s_std_nocor[j] = np.std(us_nocor)
+ #
+ us_s_mean_cor_const[j] = np.mean(us_cor_const)
+ us_s_std_cor_const[j] = np.std(us_cor_const)
+ #
+ us_s_mean_cor_diff[j] = np.mean(us_cor_diff)
+ us_s_std_cor_diff[j] = np.std(us_cor_diff)
+
+ #axl.grid("minor")
+ #fig_lines.legend()
+ #
+ '''
+ plt.figure()
+ #us_m = us_m/2*100
+ plt.plot(phi1_s,us_nocor,'+:',label="Simu no cor")
+ plt.plot(phi1_s,us_cor_const,'+:',label="Simu cor const")
+ plt.plot(phi1_s,us_cor_diff,'+:',label="Simu cor pathwise")
+ plt.ylim(-19,0)
+ plt.legend()
+ plt.grid("minor")
+ plt.title(str(np.round(dist_i)))
+ #i_0 = np.where(phi1_m == 0)
+ '''
+ return [us_s_mean_nocor,us_s_std_nocor] ,[us_s_mean_cor_const,us_s_std_cor_const],[us_s_mean_cor_diff,us_s_std_cor_diff]
+#
+#
+
+if __name__=="main":
+
+ N_dist = len(dist_m)
+ u_mg = [0] * N_dist
+ u_sg_cor_const = [0] * N_dist
+ u_sg_nocor = [0] * N_dist
+ u_sg_cor_pathconst = [0] * N_dist
+ #
+ us_m_mean = np.zeros(N_dist)
+ us_m_std = np.zeros(N_dist)
+ #
+ us_s_mean_nocor = np.zeros(N_dist)
+ us_s_std_nocor = np.zeros(N_dist)
+ #
+ us_s_mean_cor_const = np.zeros(N_dist)
+ us_s_std_cor_const = np.zeros(N_dist)
+ #
+ us_s_mean_cor_diff = np.zeros(N_dist)
+ us_s_std_cor_diff = np.zeros(N_dist)
+
+
+ for j,dist_i in enumerate(dist_m):
+ r_mg , phi_mg, u_mg[j] = make_whole_pathmat(r_m,phi_m,u_m[j])
+ r_sg , phi_sg, u_sg_nocor[j] = make_whole_pathmat(r_s,phi_s,u_nocor[j].reshape(r_s.shape))
+ r_sg , phi_sg, u_sg_cor_const[j] = make_whole_pathmat(r_s,phi_s,u_cor_const[j].reshape(r_s.shape))
+ r_sg , phi_sg, u_sg_cor_pathconst[j] = make_whole_pathmat(r_s,phi_s,u_cor_diff[j].reshape(r_s.shape))
+ phi1_m = phi_mg[:,-1]
+ phi1_s = phi_sg[:,-1]
+ #
+ #fig_lines, axl = axs = plt.subplots(1)
+ us_m = np.zeros(phi1_m.shape)
+ us_nocor = np.zeros(phi1_s.shape)
+ us_cor_const = np.zeros(phi1_s.shape)
+ us_cor_diff = np.zeros(phi1_s.shape)
+ for i in range(0,r_mg.shape[0]):
+ #axl.plot(r_mg[i,:],u_mg[i,:],'+:',label=str(np.round(phi_mg[i,-1]/pi ,decimals= 2)))
+ us_m[i] = path_int(r_mg[i,:],u_mg[j][i,:])*100
+ for i in range(0,r_sg.shape[0]):
+ us_nocor[i] = path_int(r_sg[i,:], u_sg_nocor[j][i,:])*100
+ us_cor_const[i] = path_int(r_sg[i,:], u_sg_cor_const[j][i,:])*100
+ us_cor_diff[i] = path_int(r_sg[i,:], u_sg_cor_pathconst[j][i,:])*100
+ #
+ us_m_mean[j] = np.mean(us_m)
+ us_m_std[j] = np.std(us_m)
+ #
+ us_s_mean_nocor[j] = np.mean(us_nocor)
+ us_s_std_nocor[j] = np.std(us_nocor)
+ #
+ us_s_mean_cor_const[j] = np.mean(us_cor_const)
+ us_s_std_cor_const[j] = np.std(us_cor_const)
+ #
+ us_s_mean_cor_diff[j] = np.mean(us_cor_diff)
+ us_s_std_cor_diff[j] = np.std(us_cor_diff)
+
+ #axl.grid("minor")
+ #fig_lines.legend()
+ #
+ plt.figure()
+ #us_m = us_m/2*100
+ plt.plot(phi1_s,us_nocor,'g-',label="Simu")
+ plt.plot(phi1_m,us_m,'r-',label="LDV")
+ plt.plot(phi1_s,us_cor_const,'-',label="Simu cor const")
+ plt.plot(phi1_s,us_cor_diff,'-',label="Simu cor mean_dev")
+ plt.ylim(-19,0)
+ plt.legend()
+ plt.grid("minor")
+ plt.title(str(np.round(dist_i)))
+ plt.ylabel("Meter Error in %")
+ plt.xlabel("Angular Position")
+ plt.tight_layout()
+ #plt.savefig("../Bilder_Ultrasonic/Error_over_angle_pos_" + str(j) + ".pdf")
+ i_0 = np.where(phi1_m == 0)
+
+ #
+ plt.figure()
+ plt.plot(dist_m,us_m_mean,'o:',label="LDV")
+ plt.plot(dist_m,us_s_mean_nocor,'o:',label="Simu")
+ plt.plot(dist_m,us_s_mean_cor_const,'o:',label="Simu cor const")
+ plt.plot(dist_m,us_s_mean_cor_diff,'o:',label="cor diff const")
+ plt.grid("minor")
+ plt.legend()
+ plt.tight_layout()
+
+ # make the same with finer dist --
+
+
+ # all the values are calculated if calc == True
+ calc = False
+ save_or_not = False
+ if calc:
+ all_errors_simu = [0] * len(Rk)
+ all_errors_const = [0] * len(Rk)
+ all_errors_const_diff = [0] * len(Rk)
+
+
+ for i,rki in enumerate(Rk):
+ all_errors_simu[i],all_errors_const[i],all_errors_const_diff[i] = save_us_errors(rki,dist_s,r_s,phi_s)
+
+ np.savez('US_errors.npz',np.array(all_errors_simu))
+ np.savez('US_errors_corrconst.npz',np.array(all_errors_const))
+ np.savez('US_errors_const_diff.npz',np.array(all_errors_const_diff))
+ else:
+ # otherwise the values are loaded
+ with np.load('US_errors.npz') as data:
+ all_errors_simu = data["arr_0"] #.get("arr_0")
+ with np.load('US_errors_corrconst.npz') as data:
+ all_errors_const = data["arr_0"] #.get("arr_0")
+ with np.load('US_errors_corrpathconst.npz') as data:
+ all_errors_pathconst = data["arr_0"] #.get("arr_0")
+ with np.load('US_errors_const_diff.npz') as data:
+ all_errors_const_diff = data["arr_0"] #.get("arr_0")
+
+
+
+
+ which_cor = 3
+ # which correction is applied?
+ if which_cor == 0:
+ # no correction
+ us_mean = all_errors_simu[:,0,:].T
+ us_std = all_errors_simu[:,1,:].T
+ elif which_cor == 1:
+ # const correction
+ us_mean = all_errors_const[:,0,:].T
+ us_std = all_errors_const[:,1,:].T
+ elif which_cor == 2:
+ # const path correction
+ us_mean = all_errors_pathconst[:,0,:].T
+ us_std = all_errors_pathconst[:,1,:].T
+
+ elif which_cor == 3:
+ # const diff correction
+ us_mean = all_errors_const_diff[:,0,:].T
+ us_std = all_errors_const_diff[:,1,:].T
+
+ dist_mod = dist_s[0:-1]
+ mean_int = interpol.Akima1DInterpolator(dist_mod,us_mean)
+ std_int = interpol.Akima1DInterpolator(dist_mod,us_std)
+
+ #mean_int = interpol.interp1d( dist_mod,us_mean )
+ #std_int = interpol.interp1d( dist_mod,us_std )
+
+ fig_bars, ax_bars = plt.subplots(2, 3)
+ fig_bars.set_size_inches(18, 9)
+ fig_bars.tight_layout()
+ ax_bars = ax_bars.ravel()
+
+ #
+ fig_fil, ax_fil = plt.subplots(2, 3)
+ fig_fil.set_size_inches(18, 9)
+ fig_fil.tight_layout()
+ ax_fil = ax_fil.ravel()
+ for i, disti in enumerate(dist_m):
+
+
+ mean_errors = mean_int(disti)
+ mean_std = std_int(disti)
+ #mean_errors_10 = [bla[0][distN] for bla in all_errors_simu]
+ #mean_std_10 = [bla[1][distN] for bla in all_errors_simu]
+ #
+ ax_bars[i].errorbar(Rk,mean_errors,yerr=mean_std,
+ capsize=8,
+ color='g',
+ ecolor='y',
+ #markerfacecolor='g',
+ fmt='-',
+ #markeredgewidth=2,
+ label="simu")
+ ax_bars[i].errorbar(new_rk,us_m_mean[i],yerr=us_m_std[i],fmt='ro',capsize=10,label="LDV")
+ ax_bars[i].set_title("distance "+str(round(disti))+" D")
+ ax_bars[i].set_xlim(0,10)
+ #ax_bars[i].set_ylim(-17,0.5)
+ ax_bars[i].grid("minor")
+ if i>2: ax_bars[i].set_xlabel("Kurvature Radius in D")
+ ax_bars[i].set_ylabel("Meter error in %")
+
+
+
+ #with filled areas
+ ax_fil[i].plot(Rk,mean_errors,'g-', label="mean_simu")
+ ax_fil[i].fill_between(Rk, mean_errors-mean_std, mean_errors+mean_std,
+ facecolor='green',
+ alpha = 0.3,
+ interpolate=True, label="std simu")
+
+ ax_fil[i].errorbar(new_rk,us_m_mean[i],yerr=us_m_std[i],fmt='ro',capsize=10,label="LDV")
+ ax_fil[i].set_title("distance "+str(round(disti))+" D")
+ ax_fil[i].set_xlim(0,10)
+ ax_fil[i].set_ylim(-20,0.5)
+ ax_fil[i].grid("minor")
+ if i>2: ax_fil[i].set_xlabel("Kurvature Radius in D")
+ ax_fil[i].set_ylabel("Meter error in %")
+
+ ax_bars[-1].legend(loc='lower right')
+ ax_fil[-1].legend(loc='lower right')
+ fig_fil.tight_layout()
+ if save_or_not: fig_fil.savefig("../Bilder_Ultrasonic/Mean_Error_over_Rk_Piecewlincorr_" + str(which_cor) + ".pdf")
+ fig_bars.tight_layout()
+ if save_or_not: fig_bars.savefig("../Bilder_Ultrasonic/Mean_Error_over_Rk_Bar_Piecewlincorr_" + str(which_cor) + ".pdf")
+ #
+ plt.figure()
+
+ plt.plot(dist_mod,us_mean[:,7],'g-', label="mean_simu")
+ plt.fill_between(dist_mod, us_mean[:,7]-us_std[:,7], us_mean[:,7]+us_std[:,7],
+ facecolor='green',
+ alpha = 0.3,
+ interpolate=True, label="std simu")
+
+ # plt.errorbar(dist_mod,all_errors_simu[7,0,:],yerr=all_errors_simu[7,1,:],fmt='-', label="Simu")
+ #plt.errorbar(dist_s,us_s_mean_cor_const,yerr=us_s_std_cor_const,fmt='-',label="Simu cor const")
+ #plt.errorbar(dist_s,us_s_mean_cor_diff,yerr=us_s_std_cor_diff,fmt='-',label="Simu cor path const")
+ plt.errorbar(dist_m,us_m_mean,yerr=us_m_std,fmt='ro:',capsize=10,label="LDV")
+ plt.xlim(0,60)
+ plt.ylim(-19,1)
+ plt.grid("minor")
+ plt.xlabel("Distance in Diameter")
+ plt.ylabel("Meter error in %")
+ plt.legend()
+ plt.tight_layout()
+ if save_or_not: plt.savefig("../Bilder_Ultrasonic/Mean_Error_over_dist_Piecewlincorr_" + str(which_cor) + ".pdf")
+
+
+ # plot a 2-D contour plot for the error and the uncertainty
+ end_ind = 64
+ dd, rr = np.meshgrid(dist_mod[0:end_ind],Rk)
+ nicecontour(dd.ravel(),rr.ravel(),(us_mean[0:end_ind,:].T).ravel(),maxval = 0, minval = -19 ,Nbins = 100,polar=False,withcontours = False, equalsize = False)
+ plt.xlabel("Distance in D")
+ plt.ylabel("Curvature Radius in D")
+ plt.tight_layout()
+ #plt.savefig(os.path.join('./', 'Numplot' + charnames[i] + '.png'))
+ #
+ # contour of the uncertainty
+ end_ind = 64
+ dd, rr = np.meshgrid(dist_mod[0:end_ind],Rk)
+ nicecontour(dd.ravel(),rr.ravel(),(us_std[0:end_ind,:].T).ravel(),maxval = 10, minval = 0 ,Nbins = 100,polar=False,withcontours = False, equalsize = False)
+ plt.xlabel("Distance in D")
+ plt.ylabel("Curvature Radius in D")
+ plt.tight_layout()
+
+
+
+
+
diff --git a/getUSpath.py b/getUSpath.py
new file mode 100644
index 0000000000000000000000000000000000000000..ed83b529bcd9de2a81eb9bba7a01af32ca99d432
--- /dev/null
+++ b/getUSpath.py
@@ -0,0 +1,205 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Mon Sep 20 11:44:08 2021
+
+@author: weisse02
+"""
+import json
+#import glob
+#import os
+import numpy as np
+import matplotlib.pyplot as plt
+#import matplotlib
+import scipy.interpolate as interpol
+#import pandas as pd
+#from matplotlib.ticker import MaxNLocator
+import sys
+# the mock-0.3.1 dir contains testcase.py, testutils.py & mock.py
+sys.path.append('./HelpingFunctions/')
+from niceplot import *
+#from Integrate2D import Tri_Int , get_L_norm
+from math import pi
+#sys.path.append('./Auswertung/')
+#sys.path.append("../")
+import Flow_class as flowclass
+
+
+
+def set_axes_equal(ax: plt.Axes):
+ """Set 3D plot axes to equal scale.
+
+ Make axes of 3D plot have equal scale so that spheres appear as
+ spheres and cubes as cubes. Required since `ax.axis('equal')`
+ and `ax.set_aspect('equal')` don't work on 3D.
+ """
+ limits = np.array([
+ ax.get_xlim3d(),
+ ax.get_ylim3d(),
+ ax.get_zlim3d(),
+ ])
+ origin = np.mean(limits, axis=1)
+ radius = 0.5 * np.max(np.abs(limits[:, 1] - limits[:, 0]))
+ _set_axes_radius(ax, origin, radius)
+
+def _set_axes_radius(ax, origin, radius):
+ x, y, z = origin
+ ax.set_xlim3d([x - radius, x + radius])
+ ax.set_ylim3d([y - radius, y + radius])
+ ax.set_zlim3d([z - radius, z + radius])
+
+
+
+# constructs the ultrasonic path for a given angle phi, alpha and type
+def get_US_points(pathtype = "diametrical", reflections = 0, alpha = pi/4, phi = 0,dz=0,R=1,Npoints = 100,makeplot = False):
+ # pathtype = "diametrical" means through the pipe center
+ # or nondiametrical "non"
+ # reflections = 0 for direct path, 1 for V-path, 2 for W-path
+ # alpha is the axial angle
+ # phi is the circumferential angle
+ # dz is the axial position - optional
+ # R is the pipe radius
+ # Npoints = number of points that are samples along the path
+ #
+ D = 2*R
+ # all points at the beginning, end reflection points of the paths
+ p = np.zeros((reflections+2,3))
+ # p0[3] is always dz
+ p[0] = np.array([ R*np.cos(phi), R*np.sin(phi),dz])
+ # direction vector
+ #rv = np.array([-p[0,0],-p[0,1],np.sin(alpha) ])
+ # the length of the whole path with all reflections
+ L = 0
+ # thepoints along the path on that it will be interpolated
+ p_points = np.zeros((reflections+1,3,Npoints))
+ # the direction vectors for each reflection
+ rv = np.zeros((reflections+1,3))
+ for i in range(1,reflections+2):
+ Li = D/np.cos(alpha)
+ p[i] = np.array([-p[i-1,0],-p[i-1,1], p[i-1,2] + Li* np.sin(alpha) ])
+ L += Li
+ # sample points along the line
+ rvtemp = p[i]-p[i-1]
+ #
+ # p_points[i-1] = np.array([p[i-1,j] + np.linspace(0,1,Npoints)*rvtemp[j] for j in range(0,3) ])
+ # without loop
+ p_points[i-1] = (np.tensordot(np.linspace(0,1,Npoints),rvtemp,axes=0) + p[i-1]).T
+ #np.array([p[i-1,j] + np.linspace(0,1,Npoints)*rvtemp[j] for j in range(0,3) ])
+
+ rv[i-1] = rvtemp/np.linalg.norm(rvtemp)
+
+ # normalize the z-values on diameter not on radius -scale with D
+ # the points are already scaled on R here
+ # p[:,2] /=D
+ # p_points[:,2,:] /=D
+ # plot it nicely
+ if makeplot:
+ fig = plt.figure()
+ ax = plt.axes(projection='3d')
+ # plot circle
+ fphi = np.linspace(0,2*pi,100)
+ fx,fy = pol2cart(R, fphi)
+ ax.plot(fx,fy,':')
+ ax.plot3D(fx,fy,p[-1,2]*np.ones(len(fx)),':')
+ ax.plot3D(p[:,0],p[:,1],p[:,2],'*-')
+ #
+ ax.scatter3D(p_points[:,0,:],p_points[:,1,:],p_points[:,2,:],c='r',marker='.',depthshade = False)
+ #ax.plot3D(p[1,0],p[1,1],p[1,2],'*-')
+ ax.set_xlabel("x")
+ ax.set_ylabel("y")
+ ax.set_zlabel("z")
+ #ax.set_box_aspect((1,1,1))
+ set_axes_equal(ax)
+ #ax = fig.gca(projection='3d')
+ #ax.set_aspect('equal')
+ #ax.set_aspect('auto')
+ #ax.axis("equal")
+ fig.tight_layout()
+ plt.show()
+
+ return p_points, p , rv # rvtemp
+
+
+if __name__ == "__main__":
+
+ # sample from elbow profiles
+ pathtype = "diametrical";
+ reflections = 1;
+ alpha = 30/180*pi;
+ phi = 0/180*pi;
+ R=1;
+ Npoints = 1000
+
+ # load the flow case
+ case = "DoubleElbow"
+ flow = flowclass.Elbow_profile(case)
+ #
+ r0 = flow.r[0]
+ phi0 = flow.phi[:,0]
+ z0 = flow.dist * 2 # to transform in same unit as r, normalized by R not D
+ rr,pphi,zz = np.meshgrid(z0,phi0,r0,indexing='ij')
+ #
+ US_measure = np.zeros((len(flow.Rk),len(flow.dl),len(z0),len(phi0)))
+
+ for i,rk in enumerate(flow.Rk):
+ print(rk)
+ for l,dli in enumerate(flow.dl):
+ print(dli)
+ # u has the shape (Ndist, 3, Nphi,Nr)
+ # 3 are the three velocity components
+ u = flow.get_profile(Rk = rk, dl = dli, dist = flow.dist)
+ # reshape for better in interpolation make the 3 the last dimension
+ u = np.stack((u[:,0],u[:,1],u[:,2]),axis = 3)
+ # now u has dimension (Ndist, Nphi, Nr, 3 ) same as the regular grid
+ # the interpolator for all three velocities
+ I_reg = interpol.RegularGridInterpolator((z0,phi0,r0), u,bounds_error=False,fill_value=np.nan)
+ for k,idist in enumerate(z0):
+ # print(k)
+ # print(idist)
+ # this is only to get the min and max z-coordinates
+ # p_points, p, rv = get_US_points(pathtype = pathtype, reflections = reflections, alpha = alpha, phi = phi,dz=idist,Npoints = 20)
+ # maxz = np.max(p[:,2])
+ # minz = np.min(p[:,2])
+ #
+ #
+ for j,phi_j in enumerate(phi0):
+ p_points, p, rv = get_US_points(pathtype = pathtype, reflections = reflections, alpha = alpha, phi = phi_j,dz=idist,Npoints = Npoints)
+ points_r,points_phi = cart2pol(p_points[:,0,:], p_points[:,1,:])
+ # perform interpolation
+ u_path = I_reg((p_points[:,2,:],points_phi,points_r))
+ # calc the mean of the reflected paths
+ if len(u_path.shape)>2: # if reflected paths
+ u_path_x = np.mean(u_path[:,:,0].T * rv[:,0],axis=1)
+ u_path_y = np.mean(u_path[:,:,1].T * rv[:,1],axis=1)
+ u_path_z = np.mean(u_path[:,:,2].T * rv[:,2],axis=1)
+ ez = np.mean(rv[:,2])
+ else:
+ print("direct path")
+ u_path_x = u_path[:,0]*rv[0]
+ u_path_y = u_path[:,1]*rv[1]
+ u_path_z = u_path[:,2]*rv[2]
+ ez = rv[2]
+ #
+ u_path = u_path_x + u_path_y + u_path_z
+ US_measure[i,l,k,j] = 1/ez * np.trapz(u_path,np.linspace(0,1,Npoints)) # np.mean(u_path)
+ # end loop over phi
+ # end loop over RK
+ # end loop over distance
+
+ #
+ # plt.figure()
+ # plt.plot(u_path.T,'*:')
+ US_measure = US_measure
+ #
+ np.savez('All_US_errors_'+ case +'_alpha_'+str(np.round(alpha/pi*180))+'_refl_'+str(int(reflections))+'.npz',pathint = US_measure)
+ #
+ # US_measure_1 = US_measure[:,7,:]
+ # us_mean = np.mean(US_measure_1,axis=1)
+ # us_std = np.std(US_measure_1,axis=1)
+ # #alldist = alldist[0:-1]
+ # plt.figure()
+ # plt.plot(alldist/2,us_mean,'-')
+ # plt.fill_between(alldist/2, us_mean-us_std, us_mean+us_std,
+ # facecolor='green',
+ # alpha = 0.3,
+ # interpolate=True, label="std simu")
+ # plt.grid("minor")
\ No newline at end of file
diff --git a/plot_coeffs.py b/plot_coeffs.py
new file mode 100644
index 0000000000000000000000000000000000000000..763a90059cd81f3445600c8d6f7639a2e7f212a6
--- /dev/null
+++ b/plot_coeffs.py
@@ -0,0 +1,47 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Tue Feb 28 14:09:14 2023
+
+@author: weisse02
+"""
+
+from Flow_class import *
+
+
+MEDIUM_SIZE = 20
+plt.rc('font', size=MEDIUM_SIZE) # controls default text sizes
+plt.rc('axes', titlesize=MEDIUM_SIZE) # fontsize of the axes title
+plt.rc('axes', labelsize=MEDIUM_SIZE) # fontsize of the x and y labels
+plt.rc('xtick', labelsize=MEDIUM_SIZE) # fontsize of the tick labels
+plt.rc('ytick', labelsize=MEDIUM_SIZE) # fontsize of the tick labels
+plt.rc('legend', fontsize=MEDIUM_SIZE) # legend fontsize
+plt.rc('figure', titlesize=MEDIUM_SIZE)
+
+savepath = "C:/Users/weisse02/PTBbox/VirtMet/Auswertung_paper/Mode_Plots/Coeffs_plot/"
+case = "Double_S_Elbow"
+coeffs,modes,umean = loadPODres(case)
+
+flow = Elbow_profile(case)
+coeffs = coeffs.reshape(flow.Rk.size,flow.dl.size,flow.dist.size,coeffs.shape[-1])
+rkdist,distrk = np.meshgrid(flow.Rk,flow.dist,indexing="ij")
+rkdl,dlrk = np.meshgrid(flow.Rk,flow.dl,indexing="ij")
+dldist,distdl = np.meshgrid(flow.dl,flow.dist,indexing="ij")
+
+# plot over rk, dist
+nicecontour(rkdist.ravel(),distrk.ravel(),coeffs[:,0,:,0].ravel(),polar=False,equalsize =False)
+plt.xlabel("Rc in D")
+plt.ylabel("Distance z in D")
+plt.tight_layout()
+plt.savefig(savepath + case + "_M0_rcdist.pdf")
+
+nicecontour(rkdl.ravel(),dlrk.ravel(),coeffs[:,:,0,0].ravel(),polar=False,equalsize =False)
+plt.xlabel("Rc in D")
+plt.ylabel("dl in D")
+plt.tight_layout()
+plt.savefig(savepath + case + "_M0_rcdl_dist0.pdf")
+
+nicecontour(dldist.ravel(),distdl.ravel(),coeffs[0,:,:,0].ravel(),polar=False,equalsize =False)
+plt.xlabel("Rc in D")
+plt.ylabel("dl in D")
+plt.tight_layout()
+plt.savefig(savepath + case + "_M0_dldist_rk0.pdf")