From 0104f356e1a5ff40e9c896bdc7dbe6ab8fa9a422 Mon Sep 17 00:00:00 2001 From: Andrej Date: Wed, 26 May 2021 11:47:38 +0200 Subject: [PATCH] poskusil dodat iterativno metodo. nrw kljub poenostavitvi nestabilna --- s2eps.py | 102 ++++++++++++++++++++++++++++++++++++++++++++++--------- 1 file changed, 86 insertions(+), 16 deletions(-) diff --git a/s2eps.py b/s2eps.py index 51434f0..5093c8b 100644 --- a/s2eps.py +++ b/s2eps.py @@ -1,8 +1,10 @@ #!/usr/bin/python3 import numpy as np +from numpy import pi from operator import itemgetter import matplotlib.pyplot as plt +from scipy.optimize import fsolve def s2p_to_narray(file): """ prebere s2p file v narray. STOLPCI --> VRSTICE!!! """ @@ -20,28 +22,95 @@ def narray_to_s(narray): s['f'] = narray[0,:] return s -def s_to_eps(s, L): +def s_to_eps(s, L,fc): """ izracunam eps in tand iz s parametrov in dolzin """ s11, s21, s12, s22, f = itemgetter(*s.keys())(s) + c = 299792458 #hitrost svetlobe v vakuumu. - X = (s11**2 - s21**2 + 1)/2/s11 - G = X + np.sqrt(X**2 - 1) - Gm = X - np.sqrt(X**2 - 1) + K = (s11**2 - s21**2 + 1)/2/s11 + G = K + np.sqrt(K**2 - 1) + Gm = K - np.sqrt(K**2 - 1) G[np.abs(G) > 1] = Gm[np.abs(G) > 1] - P = (s11 + s21 - G)/(1-(s11+s21)*G) - Lambda2 = - (1/2/np.pi/L * np.log(1/P))**2 #izmisli resitev za korene kompleksnega logaritma - # argument korena mora biti 2*pi*n, kjer je n=L/lambda_g + T = (s11 + s21 - G)/(1-(s11+s21)*G) + Lambda2 = (1j/2/pi/L * np.log(T))**2 + l0g2 = ( (f/c)**2 - (fc/c)**2 )**-1 + print(l0g2*Lambda2) + test_plot(f, np.abs(l0g2*Lambda2)) + + return 0 + +#Nicolson-Ross-Weir metoda -> NUMERICNO NESTABILNA +#def s_to_eps(s, L): +# """ izracunam eps in tand iz s parametrov in dolzin """ +# s11, s21, s12, s22, f = itemgetter(*s.keys())(s) +# +# X = (s11**2 - s21**2 + 1)/2/s11 +# G = X + np.sqrt(X**2 - 1) +# Gm = X - np.sqrt(X**2 - 1) +# G[np.abs(G) > 1] = Gm[np.abs(G) > 1] +# P = (s11 + s21 - G)/(1-(s11+s21)*G) +# Lambda2 = - (1/2/np.pi/L * np.log(1/P))**2 #izmisli resitev za korene kompleksnega logaritma +# # argument korena mora biti 2*pi*n, kjer je n=L/lambda_g +# +# return measured_group_delay(f,P) +# +# +#naumi kako napisat sistem nelinearnih enacb in na kaj jih resevat! +#def iterative(s,fc,L,Lair): +# s11, s21, s12, s22, f = itemgetter(*s.keys())(s) +# c = 299792458 #hitrost svetlobe v vakuumu. +# +# p0 = 1j*np.sqrt( (2*np.pi*f/c)**2-(2*np.pi*fc/c)**2) +# z = lambda p,L: np.exp(-p*L) +# G = lambda p: (p0-p)/(p0+p) #predpostavim permeabilnost = 1 +# +# f1 = lambda p,L,Lair: np.abs(s21) - np.abs(z(p,L)*(1-G(p)**2)/(1-z(p,L)**2*G(p)**2)) +# f2 = lambda p,L,Lair: np.abs(s11) - np.abs(G(p)*(1-z(p,L)**2)/(1-z(p,L)**2*G(p)**2)) +# f3 = lambda p,L,Lair: (s21*s12-s11*s22)-np.exp(-p0*(Lair-L))*(z(p,L)**2-G(p)**2)/(1-z(p,L)**2*G(p)**2) +# +# p,l,l_air = fsolve([f1,f2,f3],[p0,L,Lair]) +# +# return p,l,l_air +# + +def sistem_enacb(s, fc, beta, L): + """ sistem enacb za iterativno iskanje dielektricnosti + predpostavim mu=1 + + X = [g, eps, T, G]""" + s11, s21, s12, s22, f = itemgetter(*s.keys())(s) + c = 299792458 #hitrost svetlobe v vakuumu. + + g0 = 2j*pi*f/c*np.sqrt(1-(fc/f)**2) + + # desne strani enacb + T = lambda g: np.exp(-g*L) + g = lambda eps: 2j*pi*f/c*np.sqrt(eps-(fc/f)**2) # g + G = lambda g: (g0-g)/(g0+g) + + # pomozne funkcije za zadnjo enacbo + h1 = lambda T,G: T*(1-G**2)/(1-T**2*G**2) + h2 = lambda T,G: G*(1-T**2)/(1-T**2*G**2) + + + f1 = lambda X: X[2]-T(X[0]) # T + f2 = lambda X: X[0]-g(X[1]) # g + f3 = lambda X: X[3]-G(X[0]) #G + f4 = lambda X: s21+beta*s11-h1(X[2],X[3])-beta*h2(X[2],X[3]) + + return lambda x: [f1(x),f2(x),f3(x),f4(x)] + +def iterative(s, L, fc): + c = 299792458 #hitrost svetlobe v vakuumu. + f = s['f'] + g0 = 2j*pi*f/c*np.sqrt(1-(fc/f)**2) + f = sistem_enacb(s, fc, 0.5, L) + return f - return measured_group_delay(f,P) def measured_group_delay(f, P): phase = np.unwrap(np.angle(P)) test_plot(f,phase) - # zgladim fazo s polinomsko aproksimacijo druge stopnje - #ph = np.polyfit(f,phase,2) - #print(ph) - #faza = ph[2]+ph[1]*f+ph[0]*f**2 - #test_plot(f,faza-phase) return -np.diff(phase)/np.diff(f)/2/np.pi def test_plot(x,*args): @@ -49,8 +118,9 @@ def test_plot(x,*args): plt.plot(x,y) plt.show() -a=s2p_to_narray('20.s2p') +### testni klici +a=s2p_to_narray('teflon.s2p') b = narray_to_s(a) locals().update(b) -c= s_to_eps(b,6e-2) -test_plot(f[:-1],c) +print(iterative(b,6e-2,6.557e9)) +#c = s_to_eps(b, 6e-2,6.557e9)