# -*- 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))