From 074779e53f20e8868bf86020e2633aef0d4c0c18 Mon Sep 17 00:00:00 2001
From: weisse02 <andreas.weissenbrunner@ptb.de>
Date: Tue, 25 Oct 2022 15:53:36 +0200
Subject: [PATCH] gui works again for single elbows but is relatively slow,
 double Elbow still not working

---
 Flow_class.py | 203 ++++++++++++++++++++++++++++++++++++++++++--------
 GUI_Elbow.py  | 118 ++++++++++++++++++-----------
 2 files changed, 247 insertions(+), 74 deletions(-)

diff --git a/Flow_class.py b/Flow_class.py
index f21fd7c..d5dbc9c 100644
--- a/Flow_class.py
+++ b/Flow_class.py
@@ -143,6 +143,7 @@ class Elbow_profile():
         self.dist    = self.dist[:-1]
         # the courvature radius
         self.Rk      = np.array(coords["Rk"])/self.D #/1000
+        self.Rk[-1]  = np.round(self.Rk[-1],3)
         self.r       = np.array(coords["r"])
         self.phi     = np.array(coords["phi"])
         # remove the last value of the matrix this is where r = 1 
@@ -169,7 +170,7 @@ class Elbow_profile():
         if "dl" in coords.keys():
             self.dl = np.array(coords["dl"])/self.D
         else:
-            self.dl = np.array([0])
+            self.dl = np.array([0,1])
         
         return 1
     
@@ -186,6 +187,7 @@ class Elbow_profile():
         #
         pod_x,pod_y,pod_z = loadPODres(case)
         self.Nmodes    = pod_x["modes"].shape[-1]
+        self.case = case
         #
         if case=="SingleElbow":
             # save the modes in local variables and add zeros at the wall coordinates
@@ -198,9 +200,9 @@ class Elbow_profile():
             self.uz_mean  =  np.concatenate((pod_z["umean"].reshape(self.Nphi,self.Nr-1),np.zeros((self.Nphi,1))),axis=1).ravel()
             #
             #
-            # 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))
+            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))
         else: 
             self.modes_ux =  pod_x["modes"]
             self.modes_uy =  pod_y["modes"]
@@ -210,9 +212,9 @@ class Elbow_profile():
             self.uy_mean  =  pod_y["umean"]
             self.uz_mean  =  pod_z["umean"]
             #
-        coeffs_x = pod_x["coeffs"].reshape((len(self.Rk),len(self.dl),len(self.dist),self.Nmodes))
-        coeffs_y = pod_y["coeffs"].reshape((len(self.Rk),len(self.dl),len(self.dist),self.Nmodes))
-        coeffs_z = pod_z["coeffs"].reshape((len(self.Rk),len(self.dl),len(self.dist),self.Nmodes))
+            coeffs_x = pod_x["coeffs"].reshape((len(self.Rk),len(self.dl),len(self.dist),self.Nmodes))
+            coeffs_y = pod_y["coeffs"].reshape((len(self.Rk),len(self.dl),len(self.dist),self.Nmodes))
+            coeffs_z = pod_z["coeffs"].reshape((len(self.Rk),len(self.dl),len(self.dist),self.Nmodes))
         #
         # plot for test 
         # indrk    = 7
@@ -230,21 +232,31 @@ class Elbow_profile():
         self.regint_uy   = []
         self.regint_uz   = []
         # the spline interpolation attention it can only evaluate ordered input arrays 
-        for i in range(0,self.Nmodes):
-            # self.splineint_ux.append(interpol.RectBivariateSpline(self.dist,self.Rk,self.dl,pod_x["coeffs"][:,:,i].T,kx= polydeg,ky=polydeg ))
-            # self.splineint_uy.append(interpol.RectBivariateSpline(self.dist,self.Rk,self.dl,pod_y["coeffs"][:,:,i].T,kx= polydeg,ky=polydeg ))
-            # self.splineint_uz.append(interpol.RectBivariateSpline(self.dist,self.Rk,self.dl,pod_z["coeffs"][:,:,i].T,kx= polydeg,ky=polydeg ))
-            self.regint_ux.append(interpol.RegularGridInterpolator((self.Rk,self.dl,self.dist), coeffs_x[:,:,:,i] ))
-            self.regint_uy.append(interpol.RegularGridInterpolator((self.Rk,self.dl,self.dist), coeffs_y[:,:,:,i] ))
-            self.regint_uz.append(interpol.RegularGridInterpolator((self.Rk,self.dl,self.dist), coeffs_z[:,:,:,i] ))
+        if case=="SingleElbow":
+            for i in range(0,self.Nmodes):
+                # self.splineint_ux.append(interpol.RectBivariateSpline(self.dist,self.Rk,self.dl,pod_x["coeffs"][:,:,i].T,kx= polydeg,ky=polydeg ))
+                # self.splineint_uy.append(interpol.RectBivariateSpline(self.dist,self.Rk,self.dl,pod_y["coeffs"][:,:,i].T,kx= polydeg,ky=polydeg ))
+                # self.splineint_uz.append(interpol.RectBivariateSpline(self.dist,self.Rk,self.dl,pod_z["coeffs"][:,:,i].T,kx= polydeg,ky=polydeg ))
+                self.regint_ux.append(interpol.RegularGridInterpolator((self.Rk,self.dist), pod_x["coeffs"][:,:,i],method = 'linear',bounds_error = False,fill_value=None))#,bounds_error = False,fill_value=None ))
+                self.regint_uy.append(interpol.RegularGridInterpolator((self.Rk,self.dist), pod_y["coeffs"][:,:,i],method = 'linear',bounds_error = False,fill_value=None))#,bounds_error = True,fill_value=None ))
+                self.regint_uz.append(interpol.RegularGridInterpolator((self.Rk,self.dist), pod_z["coeffs"][:,:,i],method = 'linear',bounds_error = False,fill_value=None))#,bounds_error = True,fill_value=None ))
+        else: 
+            for i in range(0,self.Nmodes):
+                # self.splineint_ux.append(interpol.RectBivariateSpline(self.dist,self.Rk,self.dl,pod_x["coeffs"][:,:,i].T,kx= polydeg,ky=polydeg ))
+                # self.splineint_uy.append(interpol.RectBivariateSpline(self.dist,self.Rk,self.dl,pod_y["coeffs"][:,:,i].T,kx= polydeg,ky=polydeg ))
+                # self.splineint_uz.append(interpol.RectBivariateSpline(self.dist,self.Rk,self.dl,pod_z["coeffs"][:,:,i].T,kx= polydeg,ky=polydeg ))
+                self.regint_ux.append(interpol.RegularGridInterpolator((self.Rk,self.dl,self.dist), coeffs_x[:,:,:,i],bounds_error = False,fill_value=None ))
+                self.regint_uy.append(interpol.RegularGridInterpolator((self.Rk,self.dl,self.dist), coeffs_y[:,:,:,i],bounds_error = False,fill_value=None ))
+                self.regint_uz.append(interpol.RegularGridInterpolator((self.Rk,self.dl,self.dist), coeffs_z[:,:,:,i],bounds_error = False,fill_value=None ))
+            del coeffs_x
+            del coeffs_y
+            del coeffs_z 
         #
         #
         del pod_x 
         del pod_y 
         del pod_z
-        del coeffs_x
-        del coeffs_y
-        del coeffs_z 
+        
         #
         self.int_weights = int2d.get_int_weights(x=self.x,y=self.y,method="tri")
         
@@ -258,9 +270,14 @@ class Elbow_profile():
         # 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([int_i((Rk,dl,dist)) for int_i in self.regint_ux] ).T
-        A_uy = np.array([int_i((Rk,dl,dist)) for int_i in self.regint_uy] ).T
-        A_uz = np.array([int_i((Rk,dl,dist)) for int_i in self.regint_uz] ).T
+        if self.case=="SingleElbow":
+            A_ux = np.array([int_i((Rk,dist)) for int_i in self.regint_ux] ).T
+            A_uy = np.array([int_i((Rk,dist)) for int_i in self.regint_uy] ).T
+            A_uz = np.array([int_i((Rk,dist)) for int_i in self.regint_uz] ).T
+        else:
+            A_ux = np.array([int_i((Rk,dl,dist)) for int_i in self.regint_ux] ).T
+            A_uy = np.array([int_i((Rk,dl,dist)) for int_i in self.regint_uy] ).T
+            A_uz = np.array([int_i((Rk,dl,dist)) for int_i in self.regint_uz] ).T
         #
         # print("this is A for dist = " + str(dist) + " and Rk = " + str(Rk) + ":")
         # print(A_uz)
@@ -287,7 +304,67 @@ class Elbow_profile():
         #
         return ux,uy,uz
     
-    def get_path_profile_diam(self,Rk,dl,dist,reverse =False, addfully = False,R = 1):
+    def get_path_profile_diam_se(self,Rk,dist,reverse =False, addfully = False,R = 1):
+        # 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
+        #
+        D = 2*R
+        # for interpolation dist must be strictly increasing
+        #if reverse: dist = dist[::-1] 
+        A_ux = np.array([int_i((Rk,dist)) for int_i in self.regint_ux] ).T
+        A_uy = np.array([int_i((Rk,dist)) for int_i in self.regint_uy] ).T
+        A_uz = np.array([int_i((Rk,dist)) for int_i in self.regint_uz] ).T
+        #
+        # print("this is A for path-- dist = " + str(dist) + " and Rk = " + str(Rk) + ":")
+        # print(A_uz)
+        # 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])
+        # print(A_ux.shape)
+        #print("this is shape of A: --")
+        #print(A_uz.shape)
+        # now reconstruct the velocities
+        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: 
+            # A_ux = A_ux[::-1]
+            # A_uy = A_uy[::-1]
+            # A_uz = A_uz[::-1]
+            # change the modes to be reverse to fit to the coeffMatrix A
+            modes_ux_mat = modes_ux_mat[:,::-1,:]
+            modes_uy_mat = modes_uy_mat[:,::-1,:]
+            modes_uz_mat = modes_uz_mat[:,::-1,:]
+            #
+            ux = do_reconstruct_path(A_ux, modes_ux_mat, umean = self.ux_mean.reshape(self.Nphi,self.Nr)[:,::-1])
+            uy = do_reconstruct_path(A_uy, modes_uy_mat, umean = self.uy_mean.reshape(self.Nphi,self.Nr)[:,::-1])
+            uz = do_reconstruct_path(A_uz, modes_uz_mat, umean = self.uz_mean.reshape(self.Nphi,self.Nr)[:,::-1])
+            # 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:
+            ux = do_reconstruct_path(A_ux, modes_ux_mat, umean = self.ux_mean.reshape(self.Nphi,self.Nr))#[::-1]
+            uy = do_reconstruct_path(A_uy, modes_uy_mat, umean = self.uy_mean.reshape(self.Nphi,self.Nr))#[::-1]
+            uz = do_reconstruct_path(A_uz, modes_uz_mat, umean = self.uz_mean.reshape(self.Nphi,self.Nr))#[::-1]
+            wx = D*np.cos(self.phi)
+            wy = D*np.sin(self.phi)
+        # if the dimensions of a is larger:
+        # b1 = bla.reshape(*bla.shape[0:-1],self.xm.shape)
+        if addfully:
+            uz += self.u_fully
+        # if dist[0] == 0:
+        #     print("Reverse is " + str(reverse))
+        #     print("wx")
+        #     print(wx[:,0])
+        #     print("wy")
+        #     print(wy[:,0])
+        # return the radial velocity part and the axial velocity
+        return (wx*ux + wy*uy,uz)
+    
+    def get_path_profile_diam_de(self,Rk,dl,dist,reverse =False, addfully = False,R = 1):
         # 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
@@ -439,13 +516,65 @@ class Elbow_profile():
         return (dz0,dz1), L, wz, (dl0,dl1)
     
     
-    def get_pathint_allphi(self,dz,dL,wz,Rc,dl,L,makeplot=False):
+    def get_pathint_allphi_se(self,dz,dL,wz,Rc,L,makeplot=False):
+        # only for singleElbow case
+        Int = 0 #np.zeros(dz[0].shape[1]) 
+        ind0  = np.where(self.phi[:,0] == 0)[0][0] # 40 # index where phi==0
+        for ipath in range(dz[0].shape[0]):
+            # now reconstruct the profiles from the POD coeffs and integrate them
+            upaths_r0, upaths_z0 = self.get_path_profile_diam_se(Rk = Rc, dist = dz[0][ipath],reverse=True) 
+            upaths_r1, upaths_z1 = self.get_path_profile_diam_se(Rk = Rc, dist = dz[1][ipath],)
+            # integrate the velocities along the paths 
+            # the radial velocity is already weighted
+            # the integrals of the half-paths until the center point
+            #upaths_r0 = 0
+            #upaths_r1 = 0
+            Int0 = np.trapz(upaths_r0 + upaths_z0*wz  , x = dL[0][ipath]) 
+            Int1 = np.trapz(upaths_r1 + upaths_z1*wz , x = dL[1][ipath])
+            # the integral over the whole path
+            Int_path = Int0 + np.concatenate((Int1[ind0:-1],Int1[0:ind0+1]))
+            
+            # if the pathnumber is even:
+            if np.mod(ipath,2): 
+                # rotation of the paths, means an index shift
+                Int += np.concatenate((Int_path[ind0:-1],Int_path[0:ind0+1]))
+            else:
+                # add the pathintegrals
+                Int += Int_path
+            #
+            if makeplot and dz[0][0,0]==0:
+                # the velocities rearanged to fit for plotting
+                uz1temp = np.concatenate((upaths_z1[ind0:-1,:],upaths_z1[0:ind0+1,:]))
+                # the whole path
+                uzpath = np.concatenate((upaths_z0,uz1temp),axis=1)
+                if np.mod(ipath,2): 
+                    uzpath = np.concatenate((uzpath[ind0:-1,:],uzpath[0:ind0+1,:]),axis=0)
+                plt.figure()     
+                plt.plot(np.concatenate((dL[0][ipath],dL[1][ipath])), uzpath[60],':*')
+                # plt.plot(dL[0][ipath],upaths_z0[60],'+:')
+                # upaths_z1 = np.concatenate((upaths_z1[ind0:-1,:],upaths_z1[0:ind0+1,:]))
+                # plt.plot(dL[1][ipath],upaths_z1[60],'o:')
+                # plot over the radius coordinates
+                # self.ax_meter.plot(1-flow.r[0][::-1],upaths_z0[0],'*:')
+                # self.ax_meter.plot(1+flow.r[0],upaths_z1[0],'<:')
+                # or plot over the axial distance
+                #plt.plot(dz[0][ipath],upaths_z0[0],'*:')
+                #plt.plot(dz[1][ipath],upaths_z1[0],'<:')
+                plt.ylim([-0.35,0.35])
+                plt.grid()
+                plt.tight_layout()
+            # end if
+        # end loop over path-reflections
+        Int /= L*wz
+        return Int
+    def get_pathint_allphi_de(self,dz,dL,wz,Rc,dl,L,makeplot=False):
+        # only for double elbow case
         Int = 0 #np.zeros(dz[0].shape[1]) 
         ind0  = np.where(self.phi[:,0] == 0)[0][0] # 40 # index where phi==0
         for ipath in range(dz[0].shape[0]):
             # now reconstruct the profiles from the POD coeffs and integrate them
-            upaths_r0, upaths_z0 = self.get_path_profile_diam(Rk = Rc,dl = dl , dist = dz[0][ipath],reverse=True) 
-            upaths_r1, upaths_z1 = self.get_path_profile_diam(Rk = Rc,dl = dl , dist = dz[1][ipath],)
+            upaths_r0, upaths_z0 = self.get_path_profile_diam_de(Rk = Rc,dl = dl , dist = dz[0][ipath],reverse=True) 
+            upaths_r1, upaths_z1 = self.get_path_profile_diam_de(Rk = Rc,dl = dl , dist = dz[1][ipath])
             # integrate the velocities along the paths 
             # the radial velocity is already weighted
             # the integrals of the half-paths until the center point
@@ -490,6 +619,7 @@ class Elbow_profile():
         Int /= L*wz
         return Int
     
+    
     def get_pathint(self, fm,R=1,Rcs = None,dls = None,dists=None,makeplot = False):
         # calculates the diametrical pathintegrals 
         # for all angles, distances, curvature radii and inbetween dists dl
@@ -505,12 +635,23 @@ class Elbow_profile():
         print("this is L " + str(L))
         print("this is dL " + str(dL[-1][-1,-1]))
         # run through all distances
-        Int = np.zeros((len(Rcs),len(dls),len(dists),self.Nphi))
-        for i,rci in enumerate(Rcs):
-            for k,dlk in enumerate(dls):
-                for j,distj in enumerate(dists):
-                    #run through all the path reflections and integrate them for all angles
-                    Int[i,k,j] = self.get_pathint_allphi(dz + distj,dL,weights,rci,dlk,L,makeplot=makeplot)
-        return Int
+        if self.case == "SingleElbow":
+            Int = np.zeros((len(Rcs),len(dists),self.Nphi))
+            for i,rci in enumerate(Rcs):
+                    for j,distj in enumerate(dists):
+                        #run through all the path reflections and integrate them for all angles
+                        Int[i,j] = self.get_pathint_allphi_se(dz + distj,dL,weights,rci,L,makeplot=makeplot)
+                        # print(rci)
+                        # print(distj)
+            regInt = interpol.RegularGridInterpolator((self.Rk,self.dist,self.phi[:,0]),Int,bounds_error = False,fill_value=None)
+        else:
+            Int = np.zeros((len(Rcs),len(dls),len(dists),self.Nphi))
+            for i,rci in enumerate(Rcs):
+                for k,dlk in enumerate(dls):
+                    for j,distj in enumerate(dists):
+                        #run through all the path reflections and integrate them for all angles
+                        Int[i,k,j] = self.get_pathint_allphi_de(dz + distj,dL,weights,rci,dlk,L,makeplot=makeplot)
+            regInt = interpol.RegularGridInterpolator((self.Rk,self.dl,self.dist,self.phi[:,0]),Int,bounds_error = False,fill_value=None)
+        return Int , regInt
 # end flow class
 # 
diff --git a/GUI_Elbow.py b/GUI_Elbow.py
index ac972d1..806e5dd 100644
--- a/GUI_Elbow.py
+++ b/GUI_Elbow.py
@@ -209,7 +209,7 @@ class mainPanel(wx.Panel):
         self.bounds_ks   = [0.0,1e-1]
         #
         # make the float slider for the parameters 
-        sliderbox_dl,   self.sld_dl,   self.slider_label_dl_value,  self.txtbox_dl_u   = self.float_slider("inbetween distance in D",self.dl,flow.dl.min(),flow.dl.max(),factor=100)
+        sliderbox_dl,   self.sld_dl,   self.slider_label_dl_value,  self.txtbox_dl_u   = self.float_slider("inbetween distance in D",self.dl,flow.dl.min(),flow.dl.max(),factor=1000)
         sliderbox_Rc,   self.sld_Rc,   self.slider_label_Rc_value,  self.txtbox_Rc_u   = self.float_slider("Curvature Radius in D",self.Rc,flow.Rk.min(),flow.Rk.max(),factor=1000)
         sliderbox_dist, self.sld_dist, self.slider_label_dist_value,self.txtbox_dist_u = self.float_slider("Distance in D",self.dist,flow.dist.min(),flow.dist.max(),factor=1)
         sliderbox_phi,  self.sld_phi,  self.slider_label_phi_value, self.txtbox_phi_u  = self.float_slider("Angle in grad",0,-180,180,factor= 1)
@@ -261,6 +261,7 @@ class mainPanel(wx.Panel):
         vbox_slider.Add(hbox_check,1, flag = wx.EXPAND | wx.TOP, border = 1)
         #
         # add all sliderboxes to the vbox 
+        # if flow.case != "SingleElbow":
         vbox_slider.Add(sliderbox_dl  ,1, wx.EXPAND | wx.ALL, b)
         vbox_slider.Add(sliderbox_Rc  ,1, wx.EXPAND | wx.ALL, b)
         vbox_slider.Add(sliderbox_dist,1, wx.EXPAND | wx.ALL, b)
@@ -442,11 +443,11 @@ class mainPanel(wx.Panel):
         print(fm.options[0][fm.selectedOptions[0]] + ', ' + fm.options[1][fm.selectedOptions[1]] + ', ' + fm.selectedOptions[2] + ', ' + fm.selectedOptions[3])
         #
         # calculate the flow meter integrals for all paths distances and curvature radii
-        self.pathint = flow.get_pathint(fm, R=1, Rcs = None,dls= None, dists=None,makeplot = False)
+        self.pathint , self.regint_fm = flow.get_pathint(fm, R=1, Rcs = None,dls= None, dists=None,makeplot = False)
         print("pathint shape")
         print(self.pathint.shape)
         # a regular grid interpolator for the flow meter values
-        self.regint_fm = interpol.RegularGridInterpolator((flow.Rk,flow.dl,flow.dist,flow.phi[:,0]),self.pathint)
+        #self.regint_fm = interpol.RegularGridInterpolator((flow.Rk,flow.dl,flow.dist,flow.phi[:,0]),self.pathint)
         self.draw_meter()
         # now integrate the paths with the weights ... tomorrow
         #pdb.set_trace()
@@ -515,7 +516,6 @@ class mainPanel(wx.Panel):
                 # self.slider_label_dist_value.SetValue(self.dist)
                 # self.slider_label_phi_value.SetValue(flow.phi[self.phi,0])
                 # self.slider_label_Re_value.SetValue(self.Re)
-                
                 self.draw()
                 if self.pathint.any():
                     self.draw_meter()
@@ -557,7 +557,7 @@ class mainPanel(wx.Panel):
         
     def OnSliderScroll_dl(self,e):
         obj = e.GetEventObject() 
-        self.dl = obj.GetValue() 
+        self.dl = obj.GetValue()/1000
         self.slider_label_dl_value.SetValue(str(self.dl))
         #self.draw()
     def OnSliderScroll_dist(self, e):
@@ -599,19 +599,50 @@ class mainPanel(wx.Panel):
         if self.pathint.any():
             self.fm_mean = np.zeros_like(self.pathint)
             self.fm_std  = np.zeros_like(self.pathint)
-            for i,rc in enumerate(flow.Rk):
-                for j,dl in enumerate(flow.dl):
-                    for k,dist in enumerate(flow.dist):
-                        for p,phi in enumerate(flow.phi[:,0]):
-                            self.fm_mean[i,j,k,p], self.fm_std [i,j,k] = self.get_uncertainty(rc,dl,dist,phi)
+            
+            if flow.case == "SingleElbow":
+                for i,rc in enumerate(flow.Rk):
+                    for j,dist in enumerate(flow.dist):
+                        for k,phi in enumerate(flow.phi[:,0]):
+                            self.fm_mean[i,j,k], self.fm_std[i,j,k] = self.get_uncertainty(rc,0,dist,phi)
+                        #
                     #
+                    print(rc)
                 #
-                print(rc)
+                self.regint_mean = interpol.RegularGridInterpolator((flow.Rk,flow.dist,flow.phi[:,0]),self.fm_mean,bounds_error = False,fill_value=None)
+                self.regint_std  = interpol.RegularGridInterpolator((flow.Rk,flow.dist,flow.phi[:,0]),self.fm_std,bounds_error = False,fill_value=None)
+            else:
+                for i,rc in enumerate(flow.Rk):
+                    for j,dl in enumerate(flow.dl):
+                        for k,dist in enumerate(flow.dist):
+                            for p,phi in enumerate(flow.phi[:,0]):
+                                self.fm_mean[i,j,k,p], self.fm_std[i,j,k,p] = self.get_uncertainty(rc,dl,dist,phi)
+                            #
+                        #
+                    #
+                    print(rc)
+                #
+                self.regint_mean = interpol.RegularGridInterpolator((flow.Rk,flow.dl,flow.dist,flow.phi[:,0]),self.fm_mean,bounds_error = False,fill_value=None)
+                self.regint_std  = interpol.RegularGridInterpolator((flow.Rk,flow.dl,flow.dist,flow.phi[:,0]),self.fm_std,bounds_error = False,fill_value=None)
             #
-            self.regint_mean = interpol.RegularGridInterpolator((flow.Rk,flow.dl,flow.dist,flow.phi[:,0]),self.fm_mean)
-            self.regint_std  = interpol.RegularGridInterpolator((flow.Rk,flow.dl,flow.dist,flow.phi[:,0]),self.fm_std)
             self.draw_meter()
         #
+    def interpol_meter(self,Rc,dl,dist,phi):
+        y_u = np.zeros(0)
+        if flow.case == "SingleElbow":
+            if self.fm_mean.any():
+                y   = self.regint_mean((Rc,dist,phi))
+                y_u = self.regint_std((Rc,dist,phi))
+            else:
+                y   = self.regint_fm((Rc,dist,phi))
+        else:
+            if self.fm_mean.any():
+                y   = self.regint_mean((Rc,dl,dist,phi))
+                y_u = self.regint_std((Rc,dl,dist,phi))
+            else:
+                y   = self.regint_fm((Rc,dl,dist,phi))
+        return y, y_u
+        
     def draw_meter(self):
         #get_path = interpol.Akima1DInterpolator(flow.phi[:,0],self.pathint,axis=2)
         #pathintphi = get_path(self.phi)
@@ -621,33 +652,21 @@ class mainPanel(wx.Panel):
         y_u = np.zeros(0)
         if plot_over == 0: # this is "distance"
             x = flow.dist
-            if self.fm_mean.any():
-                y   = self.regint_mean((self.Rc,self.dl,flow.dist,self.phi))
-                y_u = self.regint_std((self.Rc,self.dl,flow.dist,self.phi))
-            else:
-                y = self.regint_fm((self.Rc,flow.dist,self.phi))
+            y, y_u = self.interpol_meter(self.Rc,self.dl,flow.dist,self.phi)
             xlabel = "distance in z/D"
             # the selected x-value from the sliders
             x_selected = self.dist
             #
         elif plot_over == 1: # this is "Rc"
             x = flow.Rk
-            if self.fm_mean.any():
-                y   = self.regint_mean((flow.Rk,self.dl,self.dist,self.phi))
-                y_u = self.regint_std((flow.Rk,self.dl,self.dist,self.phi))
-            else:    
-                y = self.regint_fm((flow.Rk,self.dl,self.dist,self.phi))
+            y, y_u = self.interpol_meter(flow.Rk,self.dl,self.dist,self.phi)
             xlabel = "curvature radius in Rc/D"
             # the selected x-value from the sliders
             x_selected = self.Rc
             #
         elif plot_over == 2: # this is "phi"
             x = flow.phi[:,0]/pi*180
-            if self.fm_mean.any():
-                y   = self.regint_mean((self.Rc,self.dl,self.dist,flow.phi[:,0]))
-                y_u = self.regint_std((self.Rc,self.dl,self.dist,flow.phi[:,0]))
-            else:
-                y = self.regint_fm((self.Rc,self.dl,self.dist,flow.phi[:,0]))
+            y,y_u   = self.interpol_meter(self.Rc,self.dl,self.dist,flow.phi[:,0])
             xlabel = "angle in grad"
             # the selected x-value from the sliders
             x_selected = self.phi/pi*180
@@ -885,13 +904,14 @@ class mainPanel(wx.Panel):
         self.fully_std = np.concatenate((np.flip(s),s))
         # np.tensordot(a,ii,axes=0).T
     #
-    def get_uncertainty(self,Rc,dist,phi):
+    def get_uncertainty(self,Rc,dl,dist,phi):
         # Rc = flow.Rk[-1]
         # test 
         # collect integration values for dist and rc and phi
         Ndist = 10 #2* self.bounds_phi[]/(flow.Nphi-1)
         Nrc   = 10
         Nphi  = 81
+        Ndl   = 10
         #np.where(flow.phi[:,0] < self.phi)
         #dphi  = 3* np.min(self.phi-  self.bounds_phi[0]) /(flow.Nphi-1)
         if self.dist_u > 0:
@@ -903,7 +923,7 @@ class mainPanel(wx.Panel):
         #
         if self.Rc_u > 0: 
             lim0   = max(self.bounds_Rc[0], Rc - self.Rc_u)
-            lim1   = min(self.bounds_Rc[1], Rc + self.Rc_u) - np.finfo(float).eps
+            lim1   = min(self.bounds_Rc[1], Rc + self.Rc_u)  - np.finfo(float).eps
             all_Rc = np.linspace(lim0,lim1 , Nrc)
             # print("Rc_lim1")
             # print(lim1)
@@ -921,18 +941,26 @@ class mainPanel(wx.Panel):
         else: 
             all_phi = phi
         #
-        # print("all_rc:")
-        # print(all_Rc)
-        # print("all_dist:")
-        # print(all_dist)
-        # print("all_phi:")
-        # print(all_phi)
-        # gather all 
-        rcgrid, distgrid, phigrid = np.meshgrid(all_Rc,all_dist,all_phi)
-        # ugrid = flow.get_profile(distgrid,rcgrid)
-        all_fm = self.regint_fm((rcgrid,distgrid,phigrid))
-        # flatten the first two dimensions
-        # ugrid   = ugrid.reshape(-1,*ugrid.shape[2:])
+        if flow.case ==  "SingleElbow":
+            # gather all 
+            rcgrid, distgrid, phigrid = np.meshgrid(all_Rc,all_dist,all_phi)
+            # ugrid = flow.get_profile(distgrid,rcgrid)
+            all_fm = self.regint_fm((rcgrid,distgrid,phigrid))
+            # flatten the first two dimensions
+            # ugrid   = ugrid.reshape(-1,*ugrid.shape[2:])
+        else: 
+            # also variate the inbetween distance dl
+            if self.dl_u > 0:
+                lim0   = max(self.bounds_dl[0], dl - self.dl_u)
+                lim1   = min(self.bounds_dl[1], dl + self.dl_u) # - np.finfo(float).eps
+                all_dl = np.linspace(lim0,lim1 , Ndl)
+            else: 
+                all_dl = dl
+            # and ghater all
+            
+            rcgrid,dlgrid, distgrid, phigrid = np.meshgrid(all_Rc,all_dl,all_dist,all_phi)
+            # ugrid = flow.get_profile(distgrid,rcgrid)
+            all_fm = self.regint_fm((rcgrid,dlgrid,distgrid,phigrid))
         return np.mean(all_fm), np.std(all_fm)
         
     #
@@ -977,11 +1005,15 @@ class Main(wx.Frame):
         self.Close()
     #
     def OnLoad_SE(self,e):
+        global flow
         flow = Elbow_profile('SingleElbow')
         wx.Frame.__init__(self, parent = None,id = wx.ID_ANY, title = "Flowmeter behind Elbow", size = (1400,1000))
+        self.mypanel = mainPanel(self)
     def OnLoad_DE(self,e):
+        global flow
         flow = Elbow_profile('DoubleElbow')
         wx.Frame.__init__(self, parent = None,id = wx.ID_ANY, title = "Flowmeter behind Elbow", size = (1400,1000))
+        self.mypanel = mainPanel(self)
 
 
 # temp_dist = 10
@@ -991,7 +1023,7 @@ class Main(wx.Frame):
 
 if __name__ == "__main__":
     #
-    flow = Elbow_profile('DoubleElbow')
+    flow = Elbow_profile('SingleElbow')
     fm   = FlowMeter()
     globdz = 0
     #
-- 
GitLab