開発環境  ·  17日 12月 2024

    openEMSを再び(5)

     

    以前ギブアップしていたopenEMSにもう一度チャレンジする。


    実は前回のやつはAddLumpedElementってやつでキャパシタを入れていたので、キャパシタがない状態でシミュレーションしてみて、後からSパラに対してキャパシタを追加した場合と比較してみたい。
    しばらく無言で絵を貼り続ける
    で、Octaveのコードが、こう。
    1. % tutorial for hyp2mat - capacitor in a microstrip.
    2. %
    3. % run from openems matlab command prompt
    4. % See hyp2mat(1) - convert hyperlynx files to matlab scripts.
    5. % (C) 2011,2012 Thorsten Liebig <[email protected]>
    6. % Copyright 2012 Koen De Vleeschauwer.
    7. %
    8. % This file is part of hyp2mat.
    9. %
    10. % This program is free software: you can redistribute it and/or modify
    11. % it under the terms of the GNU General Public License as published by
    12. % the Free Software Foundation, either version 3 of the License, or
    13. % (at your option) any later version.
    14. %
    15. % This program is distributed in the hope that it will be useful,
    16. % but WITHOUT ANY WARRANTY; without even the implied warranty of
    17. % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
    18. % GNU General Public License for more details.
    19. %
    20. % You should have received a copy of the GNU General Public License
    21. % along with this program. If not, see <http://www.gnu.org/licenses/>.
    22. close all
    23. clear
    24. clc
    25. function port_return=Execute_main(Sim_Path,f_max,show_structure,port_active)
    26.   % initialize
    27.   physical_constants;
    28.   %f_max = 1e9;
    29.   % Importing printed circuit board
    30.   CSX = InitCSX();
    31.   CSX = ImportHyperLynx(CSX, 'small_loop_antenna2.hyp');
    32.   % Adding a component
    33. ## [pad1_material, pad1_start, pad1_stop] = GetHyperLynxPort(CSX, 'C1.1');
    34. ## [pad2_material, pad2_start, pad2_stop] = GetHyperLynxPort(CSX, 'C1.2');
    35. ##
    36. ## c1_height = pad1_stop - pad1_start;
    37. ## c1_start = [(pad1_start(1)+pad1_stop(1))/2, pad1_start(2), pad1_start(3)];
    38. ## c1_stop = [(pad2_start(1)+pad2_stop(1))/2, pad2_stop(2), pad2_stop(3)+c1_height];
    39. ##
    40. ## CSX = AddLumpedElement( CSX, 'Capacitor', 0, 'Caps', 1, 'C', 1e-12);
    41. ## CSX = AddBox( CSX, 'Capacitor', 0, c1_start, c1_stop );
    42.   % Adding excitation and load
    43.   [port1_material, port1_start, port1_stop] = GetHyperLynxPort(CSX, 'TP1.1');
    44.   [gnd1_material, gnd1_start, gnd1_stop] = GetHyperLynxPort(CSX, 'TP3.1');
    45.   [port2_material, port2_start, port2_stop] = GetHyperLynxPort(CSX, 'TP2.1');
    46.   [gnd2_material, gnd2_start, gnd2_stop] = GetHyperLynxPort(CSX, 'TP4.1');
    47.   if(port_active==1)
    48.     [CSX, port{1}] = AddLumpedPort( CSX, 999, 1, 50, gnd1_start, port1_stop, [0 0 -1], true);
    49.     [CSX, port{2}] = AddLumpedPort( CSX, 999, 2, 50, gnd2_start, port2_stop, [0 0 -1]);
    50.   else
    51.     [CSX, port{1}] = AddLumpedPort( CSX, 999, 1, 50, gnd1_start, port1_stop, [0 0 -1]);
    52.     [CSX, port{2}] = AddLumpedPort( CSX, 999, 2, 50, gnd2_start, port2_stop, [0 0 -1], true);
    53.   endif
    54.   % Setting up a mesh
    55.   unit = GetUnits(CSX);
    56.   substrate_epr = GetEpsilon(CSX);
    57.   resolution = c0 / f_max / sqrt(substrate_epr) / unit / 25;
    58.   AirBox = c0 / f_max / unit / 25;
    59.   z_top = port1_start(3);
    60.   z_bottom = gnd1_start(3);
    61.   z_middle = (z_top+z_bottom)/2;
    62.   mesh.x = [];
    63.   mesh.y = [];
    64.   mesh.z = [ z_middle ];
    65.   mesh = DetectEdges(CSX, mesh);
    66.   mesh.x = [min(mesh.x)-AirBox max(mesh.x)+AirBox mesh.x];
    67.   mesh.y = [min(mesh.y)-AirBox max(mesh.y)+AirBox mesh.y];
    68.   mesh.z = [min(mesh.z)-AirBox max(mesh.z)+2*AirBox mesh.z];
    69.   %mesh = SmoothMesh(mesh, resolution);
    70.   mesh = SmoothMesh(mesh, resolution, 'algorithm', [ 1 ]);
    71.   % Setting boundary conditions
    72.   FDTD = InitFDTD();
    73.   FDTD = SetGaussExcite(FDTD, f_max/2, f_max/2);
    74.   BC = {'PML_8' 'PML_8' 'PML_8' 'PML_8' 'PML_8' 'PML_8'};
    75.   FDTD = SetBoundaryCond(FDTD, BC );
    76.   mesh = AddPML(mesh, 8);
    77.   CSX = DefineRectGrid(CSX, unit, mesh);
    78.   % Simulation
    79.   %Sim_Path = 'tmp';
    80.   Sim_CSX = 'msl.xml';
    81.   disp([ 'Estimated simulation runtime: 6500 timesteps' ]); % inform user this may take a while...
    82.   WriteOpenEMS([Sim_Path '/' Sim_CSX], FDTD, CSX);
    83.   if(show_structure>0)
    84.     CSXGeomPlot([Sim_Path '/' Sim_CSX]);
    85.   endif
    86.     RunOpenEMS(Sim_Path, Sim_CSX);
    87.   port_return=port;
    88. endfunction
    89. f_max = 1e9;
    90. Sim_Path1='tmp1'
    91. [status, message, messageid] = rmdir(Sim_Path1, 's'); % clear previous directory
    92. [status, message, messageid] = mkdir(Sim_Path1 ); % create empty simulation folder
    93. Sim_Path2='tmp2'
    94. [status, message, messageid] = rmdir(Sim_Path2, 's'); % clear previous directory
    95. [status, message, messageid] = mkdir(Sim_Path2 ); % create empty simulation folder
    96. port=Execute_main(Sim_Path1,f_max,1,1)
    97. % Calculating the s-parameters
    98. f = linspace( 100e6, f_max, 1601 );
    99. port = calcPort( port, Sim_Path1, f, 'RefImpedance', 50);
    100. s11 = port{1}.uf.ref./ port{1}.uf.inc;
    101. s21 = port{2}.uf.ref./ port{1}.uf.inc;
    102. port=Execute_main(Sim_Path2,f_max,0,2)
    103. % Calculating the s-parameters
    104. port = calcPort( port, Sim_Path2, f, 'RefImpedance', 50);
    105. s12 = port{1}.uf.ref./ port{2}.uf.inc;
    106. s22 = port{2}.uf.ref./ port{2}.uf.inc;
    107. s_param=[];
    108. s_param(1,1,:) = s11;
    109. s_param(1,2,:) = s12;
    110. s_param(2,1,:) = s21;
    111. s_param(2,2,:) = s22;
    112. write_touchstone('s',f,s_param,'out.s2p');
    113. close all
    114. semilogx(f/1e6,20*log10(abs(s11)),'k-','LineWidth',2);
    115. hold on;
    116. grid on;
    117. semilogx(f/1e6,20*log10(abs(s21)),'r--','LineWidth',2);
    118. legend('S_{11}','S_{21}');
    119. ylabel('S-Parameter (dB)','FontSize',12);
    120. xlabel('frequency (GHz) \rightarrow','FontSize',12);
    121. ylim([-80 10]);
    122. print('sparam.png', '-dpng');
    123. % not truncated
    Hyperlynxファイルがこう。
    1. {VERSION=2.14}
    2. {UNITS=ENGLISH LENGTH}
    3. {BOARD "small_loop_antenna2.kicad_pcb"
    4.   (PERIMETER_SEGMENT X1=6.100000000 Y1=4.050000000 X2=4.900000000 Y2=4.050000000)
    5.   (PERIMETER_SEGMENT X1=4.900000000 Y1=4.050000000 X2=4.900000000 Y2=2.850000000)
    6.   (PERIMETER_SEGMENT X1=4.900000000 Y1=2.850000000 X2=6.100000000 Y2=2.850000000)
    7.   (PERIMETER_SEGMENT X1=6.100000000 Y1=2.850000000 X2=6.100000000 Y2=4.050000000)
    8. }
    9. {STACKUP
    10.   (SIGNAL T=0.00137795 P=0 C=1.724e-08 L="F.Cu" M=COPPER)
    11.   (DIELECTRIC T=0.03937 C=4.4 L="DE_F.Cu" M="FR4")
    12.   (SIGNAL T=0.00137795 P=0 C=1.724e-08 L="B.Cu" M=COPPER)
    13. }
    14. {DEVICES
    15.   (? REF="TP1" L="F.Cu")
    16.   (? REF="TP2" L="F.Cu")
    17.   (? REF="TP4" L="B.Cu")
    18.   (? REF="TP3" L="B.Cu")
    19. }
    20. {PADSTACK=0, 0.000000000
    21.   ("F.Cu", 1, 0.039370079, 0.039370079, 180.0, M)
    22. }
    23. {PADSTACK=1, 0.000000000
    24.   ("B.Cu", 1, 0.039370079, 0.039370079, 0.0, M)
    25. }
    26. {NET="GND"
    27.   (PIN X=5.4500000000 Y=4.0000000000 R="TP4.1" P=1)
    28.   (PIN X=5.5500000000 Y=4.0000000000 R="TP3.1" P=1)
    29.   (SEG X1=5.4500000000 Y1=4.0000000000 X2=5.5500000000 Y2=4.0000000000 W=0.0393700787 L="B.Cu")
    30. }
    31. {NET="Net-(TP1-Pad1)"
    32.   (PIN X=5.5500000000 Y=4.0000000000 R="TP1.1" P=0)
    33.   (PIN X=5.4500000000 Y=4.0000000000 R="TP2.1" P=0)
    34.   (SEG X1=5.4500000000 Y1=3.9000000000 X2=5.4000000000 Y2=3.8500000000 W=0.0393700787 L="F.Cu")
    35.   (SEG X1=5.4500000000 Y1=4.0000000000 X2=5.4500000000 Y2=3.9000000000 W=0.0393700787 L="F.Cu")
    36.   (SEG X1=5.1000000000 Y1=3.0000000000 X2=5.0000000000 Y2=3.1000000000 W=0.0393700787 L="F.Cu")
    37.   (SEG X1=5.9000000000 Y1=3.0000000000 X2=5.1000000000 Y2=3.0000000000 W=0.0393700787 L="F.Cu")
    38.   (SEG X1=6.0000000000 Y1=3.1000000000 X2=5.9000000000 Y2=3.0000000000 W=0.0393700787 L="F.Cu")
    39.   (SEG X1=6.0000000000 Y1=3.1000000000 X2=6.0000000000 Y2=3.6500000000 W=0.0393700787 L="F.Cu")
    40.   (SEG X1=5.4000000000 Y1=3.8000000000 X2=5.4000000000 Y2=3.8500000000 W=0.0393700787 L="F.Cu")
    41.   (SEG X1=5.6000000000 Y1=3.7500000000 X2=5.5500000000 Y2=3.8000000000 W=0.0393700787 L="F.Cu")
    42.   (SEG X1=5.0000000000 Y1=3.1000000000 X2=5.0000000000 Y2=3.6500000000 W=0.0393700787 L="F.Cu")
    43.   (SEG X1=5.0000000000 Y1=3.6500000000 X2=5.1000000000 Y2=3.7500000000 W=0.0393700787 L="F.Cu")
    44.   (SEG X1=5.9000000000 Y1=3.7500000000 X2=5.6000000000 Y2=3.7500000000 W=0.0393700787 L="F.Cu")
    45.   (SEG X1=5.3500000000 Y1=3.7500000000 X2=5.4000000000 Y2=3.8000000000 W=0.0393700787 L="F.Cu")
    46.   (SEG X1=5.1000000000 Y1=3.7500000000 X2=5.3500000000 Y2=3.7500000000 W=0.0393700787 L="F.Cu")
    47.   (SEG X1=6.0000000000 Y1=3.6500000000 X2=5.9000000000 Y2=3.7500000000 W=0.0393700787 L="F.Cu")
    48.   (SEG X1=5.5500000000 Y1=3.8000000000 X2=5.5500000000 Y2=4.0000000000 W=0.0393700787 L="F.Cu")
    49. }
    ボードサイズと層構成はこんな感じで認識されている。
    こんなパタンとして読み込まれている。
    そしてしばらく待つ、、、意外と早かった、、、6+6分。で、こうなる。
    で、Pythonで見てみる。前回と同じコード。
    1. import numpy as np
    2. import skrf as rf
    3. import matplotlib.pyplot as plt
    4. import scipy
    5.  
    6. def main():
    7.         
    8.     ANT1 = rf.Network('out.s2p')
    9.     ANT1.name='ANT1'
    10.     freq=ANT1.frequency
    11.     tl_media = rf.DefinedGammaZ0(freq, z0=50)
    12.     gnd = rf.Circuit.Ground(freq, name='gnd')
    13.     port1 = rf.Circuit.Port(freq, name='port1', z0=50)
    14.     port2 = rf.Circuit.Port(freq, name='port2', z0=50)
    15.     cnx=[[(port1,0),(ANT1,0)],
    16.          [(ANT1,1),(port2,0)],
    17.          [(gnd,0)]]
    18.     cir=rf.Circuit(cnx)
    19.     ntw=cir.network
    20.     ntw.frequency.unit='MHz'
    21.     fig=plt.figure()
    22.     ax1=fig.add_subplot(2,2,1)
    23.     ntw.plot_s_smith(ax=ax1,m=0, n=0, lw=2)
    24.     ax2_1=fig.add_subplot(2,2,2)
    25.     ax2_2=ax2_1.twinx()
    26.     ntw.plot_s_db(ax=ax2_1,m=1, n=0, lw=2, show_legend=False)
    27.     ntw.plot_s_vswr(ax=ax2_2,m=0, n=0, lw=2, show_legend=False, color='orangered')
    28.     ax2_2.set_ylim(1,6)
    29.     handler1, label1 = ax2_1.get_legend_handles_labels()
    30.     handler2, label2 = ax2_2.get_legend_handles_labels()
    31.     ax2_1.legend(handler1 + handler2, label1 + label2, loc=2, borderaxespad=0.)
    32.     ax3=fig.add_subplot(2,2,3)
    33.     ntw.plot_s_db(ax=ax3,m=0, n=1, lw=2)
    34.     ax4=fig.add_subplot(2,2,4)
    35.     ntw.plot_s_smith(ax=ax4,m=1, n=1, lw=2)
    36.     fig.tight_layout()
    37.     plt.show()
    38. if __name__ == "__main__":
    39.     main()
    で、こうなる。
    で、キャパシタを入れてみる。こんなかんじで
    1. import numpy as np
    2. import skrf as rf
    3. import matplotlib.pyplot as plt
    4. import scipy
    5.  
    6. def main():
    7.         
    8.     ANT1 = rf.Network('out.s2p')
    9.     ANT1.name='ANT1'
    10.     freq=ANT1.frequency
    11.     tl_media = rf.DefinedGammaZ0(freq, z0=50)
    12.     C11 = tl_media.capacitor(1e-12, name='C11')
    13.     gnd = rf.Circuit.Ground(freq, name='gnd')
    14.     port1 = rf.Circuit.Port(freq, name='port1', z0=50)
    15.     port2 = rf.Circuit.Port(freq, name='port2', z0=50)
    16.     cnx=[[(port1,0),(ANT1,0)],
    17.          [(ANT1,1),(C11,0)],
    18.          [(C11,1),(port2,0)],
    19.          [(gnd,0)]]
    20.     cir=rf.Circuit(cnx)
    21.     ntw=cir.network
    22.     ntw.frequency.unit='MHz'
    23.     fig=plt.figure()
    24.     ax1=fig.add_subplot(2,2,1)
    25.     ntw.plot_s_smith(ax=ax1,m=0, n=0, lw=2)
    26.     ax2_1=fig.add_subplot(2,2,2)
    27.     ax2_2=ax2_1.twinx()
    28.     ntw.plot_s_db(ax=ax2_1,m=1, n=0, lw=2, show_legend=False)
    29.     ntw.plot_s_vswr(ax=ax2_2,m=0, n=0, lw=2, show_legend=False, color='orangered')
    30.     ax2_2.set_ylim(1,6)
    31.     handler1, label1 = ax2_1.get_legend_handles_labels()
    32.     handler2, label2 = ax2_2.get_legend_handles_labels()
    33.     ax2_1.legend(handler1 + handler2, label1 + label2, loc=2, borderaxespad=0.)
    34.     ax3=fig.add_subplot(2,2,3)
    35.     ntw.plot_s_db(ax=ax3,m=0, n=1, lw=2)
    36.     ax4=fig.add_subplot(2,2,4)
    37.     ntw.plot_s_smith(ax=ax4,m=1, n=1, lw=2)
    38.     fig.tight_layout()
    39.     plt.show()
    40. if __name__ == "__main__":
    41.     main()
    実行すると、
    キャパシタまで一緒にOpenEMSで計算したのとほんのちょっとだけ違うけど、だいたいあってるよね。OpenEMSのAddLumpedElementがそれなりにちゃんと計算しているってことがわかる。ナイス！が、AddLumpedElementするために構造が複雑になってシミュレーションにめちゃ時間かかるので、使いどころはびみょー。実際のことろ実物とどのくらい一致するのかわからんけどOpenEMSなかなかおもしろい。
     

    カテゴリ： 非現実