HFSS Layered Impedance計算

這裡透過遞迴完成計算，公式可以參考HFSS Help Layered Impedance章節。

from cmath import pi, sinh, cosh, tan, sqrt

u0 = 4 * pi * 1e-7
e0 = 8.854187817e-12

x = [(1, 1, 0, 0, 8670000, 1e-6),
     (1, 1.00018, 0, 0, 1820000, 2e-6),
     (4.4, 1, 0.02, 0, 0, 3e-6),
     (1, 1, 0, 0, 0, 1e9)]

frequency = 1e9
w = 2 * pi * frequency
k0 = 2 * pi * frequency * sqrt(u0 * e0)

def calculate(x):
    erk1, urk1, lte, ltm, sigmak, dk = x[0]

    erk2 = -(sigmak) / (w * e0) - erk1 * lte
    urk2 = -urk1 * ltm

    erk = erk1 + 1j * erk2
    urk = urk1 + 1j * urk2

    Zwk = sqrt((u0 * urk) / (e0 * erk))

    if len(x) > 1:
        Zinput = calculate(x[1:])
        rk = 1j * k0 * sqrt(erk * urk)
        return Zwk * ((Zinput * cosh(rk * dk) + Zwk * sinh(rk * dk)) / (Zinput * sinh(rk * dk) + Zwk * cosh(rk * dk)))

    else:
        return Zwk
Z = calculate(x)

(圖一) 程式計算結果

(圖二) HFSS計算結果

如何讓IronPython呼叫Python輸出Base64字串並取得Base64字串

這裡的特點便是只透過記憶體傳遞資料，不須透過檔案讀寫便可以在IronPython及Python之間交換資料，可以字串形式傳遞資料或圖片。範例如下：

IronPython test1.py

import subprocess
import base64
import os
os.chdir(os.path.dirname(__file__))

os.environ['PATH'] = r'C:\Program Files\AnsysEM\AnsysEM20.2\Win64\commonfiles\CPython\3_7\winx64\Release\python'
proc = subprocess.Popen(['python', 'test2.py', '[1,2,3]', '[1,2,1]'], stdout=subprocess.PIPE, stderr=subprocess.STDOUT)
base64_img = proc.communicate()[0]

print(base64_img)
base64_img_bytes = base64_img.encode('utf-8')
with open('plot.png', 'wb') as file_to_save:
    decoded_image_data = base64.b64decode(base64_img_bytes)
    file_to_save.write(decoded_image_data)

Python test2.py

import matplotlib.pyplot as plt
import io
import sys
x, y = eval(sys.argv[1]), eval(sys.argv[2])

import base64
my_stringIObytes = io.BytesIO()
plt.plot(x, y)
plt.savefig(my_stringIObytes, format='png')
my_stringIObytes.seek(0)
my_base64_jpgData = base64.b64encode(my_stringIObytes.read())
result = str(my_base64_jpgData, encoding = "utf-8")

如何在新版AEDT使用使用matplotlib庫

在AEDT2020R2當中，我們的腳本可以連結AEDT安裝目錄當中的matplotlib庫來繪製圖表以完成後處理自動化。但是在2021R1之後執行腳本卻會產生問題。原因是2021R1之後的matplotlib必須使用tcl8.6及tk8.6。而AEDT的Python Lib目錄底下卻沒有提供這兩個目錄。因此我們可以在電腦當中試著搜尋這兩個目錄，如果存在，將這兩個目錄複製到Lib目錄底下即可。

將RLC網表轉換至矩陣做AC電路求解

上一篇我介紹了如何從電路矩陣計算出頻率響應，這裡的例子介紹如何將網表(netlist)轉成矩陣。這個範例程式碼只能處理由Vsin, R, L, C所構成的任意電路。從這個例子我們可以知道如何找到足夠多的方程式來構成足以求解的矩陣。這個範例的未知變數共14個。因此矩陣大小為14x14。

AEDT 原理圖與AC頻率響應

# -*- coding: utf-8 -*-
"""
Created on Thu Aug  5 07:42:35 2021

@author: mlin
"""
from cmath import rect
from math import radians
import numpy as np
from numpy import pi
import copy
import matplotlib.pyplot as plt

netlist = '''
L2 net_1 net_6 1e-08
V4 net_1 0 SIN(0 0 1000000000 0 0 0) AC 1 0
R8 0 net_3 10 
C11 net_4 net_6 1e-11
C14 net_6 net_1 1e-10
R17 0 net_4 10 
L20 net_6 net_3 1e-08
C23 net_7 net_3 1e-11
L26 0 net_7 1e-08
'''

current = {}
nodes = {}
for i in netlist.splitlines():
    if not i.strip():
        continue
    if i[0] == 'R':
        cmp_name, node1, node2, resistance = i.split()
        current[cmp_name] = (node1, node2, '-{}'.format(resistance))
    elif i[0] == 'L':
        cmp_name, node1, node2, inductance = i.split()
        current[cmp_name] = (node1, node2, '-1j*w*{}'.format(inductance))
    elif i[0] == 'C':
        cmp_name, node1, node2, capacitance = i.split()
        current[cmp_name] = (node1, node2, '-1/(1j*w*{})'.format(capacitance))
    elif i[0] == 'V':
        cmp_name, node1, node2, *_, mag, phase = i.split()
        current[cmp_name] = (node1, node2, str(rect(float(mag), radians(float(phase)))))
    else:
        continue

    for node, p in [(node1, '1'), (node2, '-1')]:
        try:
            nodes[node] += [(cmp_name, p)]
        except:
            nodes[node] = [(cmp_name, p)]

del (nodes['0'])
unknown_name = list(nodes.keys()) + list(current.keys())

unknown_array = []
for name in unknown_name:
    if name[0] == 'V':
        unknown_array.append(current[name][2])
    else:
        unknown_array.append('0')

matrix = []
for name1 in unknown_name:
    row = ['0'] * len(unknown_name)
    if name1 in nodes:
        for branch, polarization in nodes[name1]:
            row[unknown_name.index(branch)] = polarization
        matrix.append(row)
    else:
        node1, node2, value = current[name1]
        if name1[0] == 'V':
            value = '0'
        row[unknown_name.index(name1)] = value

        if node1 in unknown_name:
            row[unknown_name.index(node1)] = '1'
        if node2 in unknown_name:
            row[unknown_name.index(node2)] = '-1'
        matrix.append(row)

result = []
evaluated_matrix = copy.deepcopy(matrix)
freq = np.arange(0.1, 1e9 + 1e6, 1e6)
for f in freq:
    w = 2 * pi * f

    for i in range(len(unknown_name)):
        for j in range(len(unknown_name)):
            evaluated_matrix[i][j] = eval(matrix[i][j])

    A = np.array(evaluated_matrix)
    B = np.array([eval(i) for i in unknown_array])
    D = np.linalg.inv(A)
    E = np.dot(D, B)
    result.append(E)

result = list(zip(*result))
plt.plot(freq, [abs(i) for i in result[2]])
# plt.ylim(0, 1)
plt.grid()
plt.show()    

Python程式碼計算出的頻率響應

AC電路模擬的算法基礎

圖一是AEDT Designer當中一個簡單的RLC網路模擬及V2的電壓頻率響應。我手動將網表轉成矩陣形式A，邊界條件(1V電壓源)為B。將A做反矩陣運算得到D。D與B做相乘即可以得到[V1, V2, I1, I2, I3, I4]未知矩陣的解。將V2取出來畫出曲線與AEDT Designer所計算的結果相同。

當中運算量最大的是np.linalg.inv()反矩陣運算。這個例子可以說明為何反矩陣運算在頻域數值模擬當中扮演著非常重要的角色。關於Python反矩陣運算請參考：

https://dwightreid.com/blog/2015/09/21/python-how-to-solve-simultaneous-equations/

(圖一) 電路模擬結果

電壓電流方程組

電壓電流方程式矩陣表示法

import numpy as np
from numpy import pi

v2 = []
freq = np.arange(1e6, 1e9 + 1e7, 1e7)
for f in freq:
    w = 2 * pi * f
    R = 100
    Z1 = 1 / (1j * w * 1e-10)
    Z2 = 1j * w * 1e-7
    A = np.array([[1, 0, 0, 0, 0, 0],
                  [0, 0, 1, -1, -1, 0],
                  [0, 0, 0, 0, 1, -1],
                  [0, 1, 0, 0, 0, -R],
                  [1, 0, 0, -Z1, 0, 0],
                  [1, -1, 0, 0, -Z2, 0]])

    B = np.array([1, 0, 0, 0, 0, 0])
    D = np.linalg.inv(A)
    E = np.dot(D, B)
    v2.append(E[1])

import matplotlib.pyplot as plt

plt.plot(freq, [np.abs(i) for i in v2])
plt.ylim(0, 1)
plt.grid()
plt.show()

(圖二) Python計算結果

找出座標的角度(0-360)

 

from math import acos, sqrt, degrees, radians, cos, sin

pts = [(cos(radians(theta)), sin(radians(theta))) for theta in range(360)]

def getDegree(x, y):
    theta = acos(x/(sqrt(x**2 + y**2)))
    if y >= 0:
        return round(degrees(theta), 9)
    else:
        return round(360 - degrees(theta), 9)

for x, y in pts:
    print(getDegree(x, y))