m2cce/internship
m2cce
/
internship
1
Fork 0
Code Issues Activity

Merge branch 'swash'

This commit is contained in:
Edgar P. Burkhart 2022-04-01 17:11:56 +02:00
parent 3d0e9cbca7 58ed9a4716
commit 5c20dc3db7
Signed by: edpibu
GPG Key ID: 9833D3C5A25BD227
4 changed files with 441 additions and 236 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

104
swash/processing/pa/PUV.py
View File

@ -16,17 +16,7 @@ from scipy.fftpack import fft, ifft
from scipy.signal import detrend
from numpy import angle
# fonction moyenne glissante
def lissage(Lx,Ly,p):
'''Fonction qui débruite une courbe par une moyenne glissante
sur 2P+1 points'''
Lxout=[]
Lyout=[]
Lxout = Lx[p: -p]
for index in range(p, len(Ly)-p):
average = np.mean(Ly[index-p : index+p+1])
Lyout.append(average)
return Lxout,Lyout
g = 9.81
# h = profondeur locale
# fs = fréquence echantillonage
@ -47,52 +37,78 @@ def PUV_dam(u,p,h,fs,zmesP,zmesU):
Ampu = np.abs(Yu[1 : N // 2]) / (N / 2.0)
Ampp = np.abs(Yp[1 : N // 2]) / (N / 2.0)
# pas de frequence deltaf
deltaf=1/(N/2)/delta*nyquist;
deltaf = 1 / (N / 2) / delta * nyquist
# phase
ThetaU=angle(Yu[1:N//2]);
ThetaP=angle(Yp[1:N//2]);
ThetaU = angle(Yu[1 : N // 2])
ThetaP = angle(Yp[1 : N // 2])
# calcul du coefficient de reflexion
nbf=len(xf);
nbf = len(xf)
if max(p) > 0:
# si pas de données de pression, jsute Sheremet
# attention : on commence par le deuxieme point, le premier correspondant a frequence nulle
iicutoff_xb = max(min(np.where(xf > 0.5))) # length(Fu_xb)) ?
# calcul de l'amplitude en deca de la frequence de coupure
k_xb = []
ii = 1
omega = []
ccc, ccs, cccp = [], [], []
aincident13, areflechi13 = [], []
aprog = []
r13 = []
ii = 0
while ii < iicutoff_xb:
k_xb[ii] = newtonpplus(xf[ii],h,0);
a1=Ampu[ii];
a3=Ampp[ii];
phi1=ThetaU[ii];
phi3=ThetaP[ii];
omega[ii]=2*pi*xf[ii];
cc=omega[ii]*cosh(xf[ii]*(zmesU+h))/sinh(xf[ii]*h);
ccp=omega[ii]*cosh(xf[ii]*(zmesP+h))/sinh(xf[ii]*h);
cs=omega[ii]*sinh(xf[ii]*(zmesU+h))/sinh(xf[ii]*h);
k_xb.append(newtonpplus(xf[ii], h, 0))
a1 = Ampu[ii]
a3 = Ampp[ii]
phi1 = ThetaU[ii]
phi3 = ThetaP[ii]
omega.append(2 * np.pi * xf[ii])
cc = omega[ii] * np.cosh(xf[ii] * (zmesU + h)) / np.sinh(xf[ii] * h)
ccp = omega[ii] * np.cosh(xf[ii] * (zmesP + h)) / np.sinh(xf[ii] * h)
cs = omega[ii] * np.sinh(xf[ii] * (zmesU + h)) / np.sinh(xf[ii] * h)
# Procedure de calcul des ai et ar sans courant
ccc[ii]=cc;
ccs[ii]=cs;
cccp[ii]=ccp;
ccc.append(cc)
ccs.append(cs)
cccp.append(ccp)
# calcul des amplitudes des ondes incidentes a partir de uhoriz et p
aincident13[ii]=0.5*sqrt(a1*a1/(cc*cc)+a3*a3*g*g*xf[ii]*xf[ii]/(omega[ii]*omega[ii]*ccp*ccp)+2*a1*a3*g*k_xb[ii]*cos(phi3-phi1)/(cc*ccp*omega[ii]))
areflechi13[ii]=0.5*sqrt(a1*a1/(cc*cc)+a3*a3*g*g*k_xb[ii]*k_xb[ii]/(omega[ii]*omega[ii]*ccp*ccp)-2*a1*a3*g*xf[ii]*cos(phi3-phi1)/(cc*ccp*omega[ii]))
aincident13.append(0.5 * np.sqrt(
a1 * a1 / (cc * cc)
+ a3
* a3
* g
* g
* xf[ii]
* xf[ii]
/ (omega[ii] * omega[ii] * ccp * ccp)
+ 2 * a1 * a3 * g * k_xb[ii] * np.cos(phi3 - phi1) / (cc * ccp * omega[ii])
))
areflechi13.append(0.5 * np.sqrt(
a1 * a1 / (cc * cc)
+ a3
* a3
* g
* g
* k_xb[ii]
* k_xb[ii]
/ (omega[ii] * omega[ii] * ccp * ccp)
- 2 * a1 * a3 * g * xf[ii] * np.cos(phi3 - phi1) / (cc * ccp * omega[ii])
))
cv = g * xf[ii] / (omega[ii] * ccp)
cp = 1 / cc
aprog[ii]= a3/(g*xf[ii]/(omega[ii]*ccp)); #hypothese progressive Drevard et al |u|/Cv
aprog.append(a3 / (g * xf[ii] / (omega[ii] * ccp)))
# hypothese progressive Drevard et al |u|/Cv
# aprog(ii)= a1/cp;
# |p|/Cp
if aincident13[ii] < 0.18 * ampliseuil:
r13[ii]=0
r13.append(0)
else:
r13[ii]=areflechi13[ii]/aincident13[ii]
r13.append(areflechi13[ii] / aincident13[ii])
# calcul des energies incidente et reflechie sans ponderation avec les voisins
Eprog=0.5*aprog**2/deltaf;
Eprog = 0.5 * aprog**2 / deltaf
Eixb = 0.5 * aincident13**2 / deltaf
Erxb = 0.5 * areflechi13**2 / deltaf
F = Fu_xb[1:iicutoff_xb]
return Eixb, Erxb, Eprog, F
'test'
"test"
# Calcul du vecteur d'onde en prÈsence d'un courant
@ -106,13 +122,23 @@ def newtonpplus(f,h,u) :
kh = 0.5
x = 0.001
u = -u
while (abs((kh - x)/x) > 0.00000001) :
i = 0
while abs((kh - x) / x) > 0.00001:
i += 1
if i > 10**5:
sys.exit(1)
x = kh
fx = x*math.tanh(x) - (h/g)*(2*math.pi*f-(u/h)*x)*(2*math.pi*f-(u/h)*x)
fprimx = math.tanh(x) + x*(1- (math.tanh(x))**2)+(2*u/g)*(2*math.pi*f-(u/h)*x)
fx = x * math.tanh(x) - (h / g) * (2 * math.pi * f - (u / h) * x) * (
2 * math.pi * f - (u / h) * x
)
fprimx = (
math.tanh(x)
+ x * (1 - (math.tanh(x)) ** 2)
+ (2 * u / g) * (2 * math.pi * f - (u / h) * x)
)
kh = x - (fx / fprimx)
k = kh / h
return k
# fin du calcul de k
# fin du calcul de k

183
swash/processing/pa/reflex3S.py
View File

@ -8,53 +8,67 @@ from matplotlib import animation
import cmath
from scipy.fft import fft
def newtonpplus(f, h, u):
# calcul de k:
g = 9.81
kh = 0.5
x = 0.001
u = -u
while (abs((kh - x)/x) > 0.00000001) :
while abs((kh - x) / x) > 0.00000001:
x = kh
fx = x*math.tanh(x) - (h/g)*(2*math.pi*f-(u/h)*x)*(2*math.pi*f-(u/h)*x)
fprimx = math.tanh(x) + x*(1- (math.tanh(x))**2)+(2*u/g)*(2*math.pi*f-(u/h)*x)
fx = x * math.tanh(x) - (h / g) * (2 * math.pi * f - (u / h) * x) * (
2 * math.pi * f - (u / h) * x
)
fprimx = (
math.tanh(x)
+ x * (1 - (math.tanh(x)) ** 2)
+ (2 * u / g) * (2 * math.pi * f - (u / h) * x)
)
kh = x - (fx / fprimx)
k = kh / h
return k
def newtonpmoins(f, h, u0):
# calcul de k:
g = 9.81;
kh = 0.5;
x = 0.01;
x = 6*h;
g = 9.81
kh = 0.5
x = 0.01
x = 6 * h
while (np.abs((kh - x)/x) > 0.00000001):
x = kh;
fx = x*math.tanh(x) - (h/g)*(2*math.pi*f-(u0/h)*x)*(2*math.pi*f-(u0/h)*x);
fprimx = math.tanh(x) + x*(1- (math.tanh(x))**2)+(2*u0/g)*(2*math.pi*f-(u0/h)*x);
kh = x - (fx/fprimx);
while np.abs((kh - x) / x) > 0.00000001:
x = kh
fx = x * math.tanh(x) - (h / g) * (2 * math.pi * f - (u0 / h) * x) * (
2 * math.pi * f - (u0 / h) * x
)
fprimx = (
math.tanh(x)
+ x * (1 - (math.tanh(x)) ** 2)
+ (2 * u0 / g) * (2 * math.pi * f - (u0 / h) * x)
)
kh = x - (fx / fprimx)
k = kh / h
return k
# Calcul du vecteur d'onde a partir de la frÈquence
# kh : vecteur d'onde * profondeur d'eau
def newtonpropa(hi, f):
# calcul de k:
g=9.81;
si = (2*math.pi*f)**2/g;
kh = 0.5;
x = 0.001;
while (np.abs((kh - x)/x) > 0.00000001) :
x = kh;
fx = x*math.tanh(x) - si*hi;
fprimx = math.tanh(x) + x*(1- (math.tanh(x))**2);
kh = x - (fx/fprimx);
kpropa = kh/hi;
g = 9.81
si = (2 * math.pi * f) ** 2 / g
kh = 0.5
x = 0.001
while np.abs((kh - x) / x) > 0.00000001:
x = kh
fx = x * math.tanh(x) - si * hi
fprimx = math.tanh(x) + x * (1 - (math.tanh(x)) ** 2)
kh = x - (fx / fprimx)
kpropa = kh / hi
return kpropa
def reflex3S(x1, x2, x3, xs1, xs2, xs3, h, mean_freq, fmin, fmax):
# Analyse avec transformee de fourier d un signal en sinus
# calcul du coefficient de reflexion en presence d un courant u
@ -98,47 +112,46 @@ def reflex3S(x1,x2,x3,xs1,xs2,xs3,h,mean_freq,fmin,fmax) :
# xs1=0;xs2=0.80;xs3=1.30;
# ENTREES DONNEES DES 3 SONDES AMONT et des 2 SONDES AVAL
ampliseuil=0.005;
ampliseuil = 0.005
#'check'
# pause
# PAS DE TEMPS
deltat1=1/mean_freq;
deltat2=1/mean_freq;
deltat3=1/mean_freq;
deltat1 = 1 / mean_freq
deltat2 = 1 / mean_freq
deltat3 = 1 / mean_freq
# transformees de Fourier
Y1 = fft(x1,len(x1));
N1 = len(Y1);
Y2 = fft(x2,len(x2));
N2 = len(Y2);
Y3 = fft(x3,len(x3));
N3 = len(Y3);
Y1 = fft(x1, len(x1))
N1 = len(Y1)
Y2 = fft(x2, len(x2))
N2 = len(Y2)
Y3 = fft(x3, len(x3))
N3 = len(Y3)
# amplitudes normalisees, soit coef de fourier
amplitude1=np.abs(Y1[1:N1//2])/(N1//2);
nyquist = 1/2;
freq1 = (np.arange(1, (N1//2)+1, 1)-1)/(N1//2)/deltat1*nyquist;
amplitude2=np.abs(Y2[1:N2//2])/(N2//2);
nyquist = 1/2;
freq2 = (np.arange(1, (N2//2)+1, 1)-1)/(N2//2)/deltat2*nyquist;
amplitude3=np.abs(Y3[1:N3//2])/(N3//2);
nyquist = 1/2;
freq3 = (np.arange(1, (N3//2)+1, 1)-1)/(N3//2)/deltat3*nyquist;
amplitude1 = np.abs(Y1[1 : N1 // 2]) / (N1 // 2)
nyquist = 1 / 2
freq1 = (np.arange(1, (N1 // 2) + 1, 1) - 1) / (N1 // 2) / deltat1 * nyquist
amplitude2 = np.abs(Y2[1 : N2 // 2]) / (N2 // 2)
nyquist = 1 / 2
freq2 = (np.arange(1, (N2 // 2) + 1, 1) - 1) / (N2 // 2) / deltat2 * nyquist
amplitude3 = np.abs(Y3[1 : N3 // 2]) / (N3 // 2)
nyquist = 1 / 2
freq3 = (np.arange(1, (N3 // 2) + 1, 1) - 1) / (N3 // 2) / deltat3 * nyquist
# recherche de la phase
theta1=np.angle(Y1[1:N1//2]);
theta2=np.angle(Y2[1:N2//2]);
theta3=np.angle(Y3[1:N3//2]);
theta1 = np.angle(Y1[1 : N1 // 2])
theta2 = np.angle(Y2[1 : N2 // 2])
theta3 = np.angle(Y3[1 : N3 // 2])
# pas de frequence deltaf
deltaf=1/(N1//2)/deltat1*nyquist;
nbrefreq=len(freq1);
deltaf = 1 / (N1 // 2) / deltat1 * nyquist
nbrefreq = len(freq1)
# Caracteristiques fondamentaux,sondes canaux 1 et 3
# distances entre les sondes
x12=xs2-xs1;
x13=xs3-xs1;
x23=xs3-xs2;
x12 = xs2 - xs1
x13 = xs3 - xs1
x23 = xs3 - xs2
# Debut calcul des coefficients de reflexion
indmin=np.min(np.where(freq1>0.02));
indfmin=np.min(np.where(freq1>fmin));
indfmax=np.max(np.where(freq1<fmax));
indmin = np.min(np.where(freq1 > 0.02))
indfmin = np.min(np.where(freq1 > fmin))
indfmax = np.max(np.where(freq1 < fmax))
T = []
fre = []
@ -183,16 +196,34 @@ def reflex3S(x1,x2,x3,xs1,xs2,xs3,h,mean_freq,fmin,fmax) :
delta13m = -kmoins * x13
delta23m = -kmoins * x23
# calcul du coefficient de reflexion a partir des sondes 1 et 2
aincident12.append(math.sqrt(a1*a1+a2*a2-2*a1*a2*math.cos(phi12+delta12p))/(2*np.abs(math.sin((delta12p+delta12m)/2))))
areflechi12.append(math.sqrt(a1*a1+a2*a2-2*a1*a2*math.cos(phi12-delta12m))/(2*np.abs(math.sin((delta12p+delta12m)/2))))
aincident12.append(
math.sqrt(a1 * a1 + a2 * a2 - 2 * a1 * a2 * math.cos(phi12 + delta12p))
/ (2 * np.abs(math.sin((delta12p + delta12m) / 2)))
)
areflechi12.append(
math.sqrt(a1 * a1 + a2 * a2 - 2 * a1 * a2 * math.cos(phi12 - delta12m))
/ (2 * np.abs(math.sin((delta12p + delta12m) / 2)))
)
# r12(jj)=areflechi12(jj)/aincident12(jj);
# calcul du coefficient de reflexion a partir des sondes 2 et 3
aincident23.append(math.sqrt(a2*a2+a3*a3-2*a2*a3*math.cos(phi23+delta23p))/(2*np.abs(math.sin((delta23p+delta23m)/2))))
areflechi23.append(math.sqrt(a2*a2+a3*a3-2*a2*a3*math.cos(phi23-delta23m))/(2*np.abs(math.sin((delta23p+delta23m)/2))))
aincident23.append(
math.sqrt(a2 * a2 + a3 * a3 - 2 * a2 * a3 * math.cos(phi23 + delta23p))
/ (2 * np.abs(math.sin((delta23p + delta23m) / 2)))
)
areflechi23.append(
math.sqrt(a2 * a2 + a3 * a3 - 2 * a2 * a3 * math.cos(phi23 - delta23m))
/ (2 * np.abs(math.sin((delta23p + delta23m) / 2)))
)
# r23(jj)=areflechi23(jj)/aincident23(jj);
# calcul du coefficient de reflexion a partir des sondes 1 et 3
aincident13.append(math.sqrt(a1*a1+a3*a3-2*a1*a3*math.cos(phi13+delta13p))/(2*np.abs(math.sin((delta13p+delta13m)/2))))
areflechi13.append(math.sqrt(a1*a1+a3*a3-2*a1*a3*math.cos(phi13-delta13m))/(2*np.abs(math.sin((delta13p+delta13m)/2))))
aincident13.append(
math.sqrt(a1 * a1 + a3 * a3 - 2 * a1 * a3 * math.cos(phi13 + delta13p))
/ (2 * np.abs(math.sin((delta13p + delta13m) / 2)))
)
areflechi13.append(
math.sqrt(a1 * a1 + a3 * a3 - 2 * a1 * a3 * math.cos(phi13 - delta13m))
/ (2 * np.abs(math.sin((delta13p + delta13m) / 2)))
)
# r13.append(areflechi13[jj]/aincident13[jj])
# calcul du coefficient de reflexion par methode des 3 sondesavec moindres carres
delta1m = 0
@ -201,12 +232,33 @@ def reflex3S(x1,x2,x3,xs1,xs2,xs3,h,mean_freq,fmin,fmax) :
delta1p = 0
delta2p = delta12p
delta3p = delta13p
s1=cmath.exp(-1j*2*delta1m)+cmath.exp(-1j*2*delta2m)+cmath.exp(-1j*2*delta3m)
s2=cmath.exp(+1j*2*delta1p)+cmath.exp(+1j*2*delta2p)+cmath.exp(+1j*2*delta3p)
s12=cmath.exp(1j*(delta1p-delta1m))+cmath.exp(1j*(delta2p-delta2m))+cmath.exp(1j*(delta3p-delta3m))
s3=a1*cmath.exp(-1j*(phi1+delta1m))+a2*cmath.exp(-1j*(phi2+delta2m))+a3*cmath.exp(-1j*(phi3+delta3m))
s4=a1*cmath.exp(-1j*(phi1-delta1p))+a2*cmath.exp(-1j*(phi2-delta2p))+a3*cmath.exp(-1j*(phi3-delta3p))
s1 = (
cmath.exp(-1j * 2 * delta1m)
+ cmath.exp(-1j * 2 * delta2m)
+ cmath.exp(-1j * 2 * delta3m)
)
s2 = (
cmath.exp(+1j * 2 * delta1p)
+ cmath.exp(+1j * 2 * delta2p)
+ cmath.exp(+1j * 2 * delta3p)
)
s12 = (
cmath.exp(1j * (delta1p - delta1m))
+ cmath.exp(1j * (delta2p - delta2m))
+ cmath.exp(1j * (delta3p - delta3m))
)
s3 = (
a1 * cmath.exp(-1j * (phi1 + delta1m))
+ a2 * cmath.exp(-1j * (phi2 + delta2m))
+ a3 * cmath.exp(-1j * (phi3 + delta3m))
)
s4 = (
a1 * cmath.exp(-1j * (phi1 - delta1p))
+ a2 * cmath.exp(-1j * (phi2 - delta2p))
+ a3 * cmath.exp(-1j * (phi3 - delta3p))
)
s5 = s1 * s2 - s12 * s12
ai.append(abs((s2 * s3 - s12 * s4) / s5))
ar.append(abs((s1 * s4 - s12 * s3) / s5))
# refl[jj]=ar[jj]/ai[jj];
@ -216,5 +268,4 @@ def reflex3S(x1,x2,x3,xs1,xs2,xs3,h,mean_freq,fmin,fmax) :
Ereflechi123.append(0.5 * ar[count] * ar[count] / deltaf)
count += 1
return ai, ar, Eincident123, Ereflechi123, indfmin, indfmax, fre

83
swash/processing/post_reflex.py Normal file
View File

@ -0,0 +1,83 @@
import argparse
import configparser
import logging
import pathlib
import matplotlib.pyplot as plt
import numpy as np
import scipy.signal as sgl
from .reflex3s import reflex3s
from .pa.reflex3S import reflex3S
from .pa.PUV import PUV_dam
parser = argparse.ArgumentParser(description="Post-process swash output")
parser.add_argument("-v", "--verbose", action="count", default=0)
parser.add_argument("-c", "--config", default="config.ini")
args = parser.parse_args()
logging.basicConfig(level=max((10, 20 - 10 * args.verbose)))
log = logging.getLogger("post")
log.info("Starting post-processing")
config = configparser.ConfigParser()
config.read(args.config)
inp = pathlib.Path(config.get("post", "inp"))
root = pathlib.Path(config.get("swash", "out"))
log.info(f"Reading data from '{inp}'")
x = np.load(inp.joinpath("xp.npy"))
t = np.load(inp.joinpath("tsec.npy"))
botl = np.load(inp.joinpath("botl.npy"))
watl = np.load(inp.joinpath("watl.npy"))
vel = np.load(inp.joinpath("vel.npy"))
x0 = config.getint("post", "x0")
arg_x0 = np.abs(x - x0).argmin()
t0 = config.getfloat("post", "t0")
arg_t0 = np.abs(t - t0).argmin()
dt = np.diff(t).mean()
f = 1 / dt
idx = [arg_x0 - 5, arg_x0, arg_x0 + 7]
wl = watl[arg_t0:, idx]
# Reflex3S
res = reflex3S(*wl.T, *x[idx], botl[idx].mean(), f, 0.02, 0.2)
f_, ai_, ar_ = reflex3s(wl.T, x[idx], botl[idx].mean(), f)
ai, ar, _, _, _, _, fre = res
# Custom PUV
eta = watl[t > t0, arg_x0]
u = vel[t > t0, 0, arg_x0]
phi_eta = sgl.welch(eta, f)
phi_u = sgl.welch(u, f)
phi_eta_u = sgl.csd(eta, u, f)
R = np.sqrt(
(np.abs(phi_eta[1]) + np.abs(phi_u[1]) - 2 * phi_eta_u[1].real)
/ (np.abs(phi_eta[1]) + np.abs(phi_u[1]) + 2 * phi_eta_u[1].real)
)
out = pathlib.Path(config.get("post", "out")).joinpath(f"t{t0}x{x0}")
log.info(f"Saving plots in '{out}'")
out.mkdir(parents=True, exist_ok=True)
plt.plot(fre, np.asarray(ar) / np.asarray(ai))
plt.xlim((0, 0.3))
plt.ylim((0, 1))
plt.savefig(out.joinpath("reflex3s_pa.pdf"))
plt.clf()
plt.plot(f_, ar_ / ai_)
plt.plot(phi_eta[0], R)
plt.xlim((0, 0.3))
plt.ylim((0, 1))
plt.savefig(out.joinpath("reflex3s.pdf"))
# pressk = np.load(inp.joinpath("pressk.npy"))
# press = pressk[:, 1, :]
# res = PUV_dam(vel[arg_t0:, 0, 500], press[arg_t0:, 500], botl[500], f, botl[500]/2, botl[500]/2)
# print(res)

45
swash/processing/reflex3s.py Normal file
View File

@ -0,0 +1,45 @@
import numpy as np
from scipy import fft
from scipy import signal as sgl
from scipy import optimize as opti
def reflex3s(eta, x, h, f):
eta_psd = np.stack([fft.rfft(z) for z in eta])
f_psd = fft.rfftfreq(eta.shape[1], 1 / f)
eta_amp = np.abs(eta_psd) / (f_psd.size - 1)
eta_phase = np.angle(eta_psd)
g = 9.81
k = np.asarray(
[
opti.root_scalar(
f=lambda k: k * np.tanh(k) - (2 * np.pi * f) ** 2 / g * h,
fprime=lambda k: np.tanh(k) + k * (1 - np.tanh(k) ** 2),
x0=0.5,
).root
/ h
for f in f_psd
]
)
s1 = 1 + np.exp(2j * k * (x[1] - x[0])) + np.exp(2j * k * (x[2] - x[0]))
s2 = 1 + np.exp(-2j * k * (x[1] - x[0])) + np.exp(-2j * k * (x[2] - x[0]))
s12 = 3
s3 = (
eta_amp[0] * np.exp(-1j * (eta_phase[0]))
+ eta_amp[1] * np.exp(-1j * (eta_phase[1] - k * (x[1] - x[0])))
+ eta_amp[2] * np.exp(-1j * (eta_phase[2] - k * (x[2] - x[0])))
)
s4 = (
eta_amp[0] * np.exp(-1j * (eta_phase[0]))
+ eta_amp[1] * np.exp(-1j * (eta_phase[1] + k * (x[1] - x[0])))
+ eta_amp[2] * np.exp(-1j * (eta_phase[2] + k * (x[2] - x[0])))
)
s5 = s1 * s2 - s12**2
ai = np.abs((s2 * s3 - s12 * s4) / s5)
ar = np.abs((s1 * s4 - s12 * s3) / s5)
return f_psd, ai, ar