2024年12月17日火曜日

openEMSを再び(4)

 

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


ちゃんとSパラを出力する。
TutorialのPatch_Antenna_Phased_Array.mでSパラを出力していたのでよくよく見てみると、AddLumpedPortをするときにexciteするポートを替えながら、都度シミュレーションしているみたい。よくよく考えたらネットアナでの測定もそうやる(ていうか装置が勝手にそうやってる)よね。
で、コードはこうなる(計算時間を短くするために周波数範囲を狭くした)。
  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.
  6. % (C) 2011,2012 Thorsten Liebig <thorsten.liebig@gmx.de>
  7. % Copyright 2012 Koen De Vleeschauwer.
  8. %
  9. % This file is part of hyp2mat.
  10. %
  11. % This program is free software: you can redistribute it and/or modify
  12. % it under the terms of the GNU General Public License as published by
  13. % the Free Software Foundation, either version 3 of the License, or
  14. % (at your option) any later version.
  15. %
  16. % This program is distributed in the hope that it will be useful,
  17. % but WITHOUT ANY WARRANTY; without even the implied warranty of
  18. % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  19. % GNU General Public License for more details.
  20. %
  21. % You should have received a copy of the GNU General Public License
  22. % along with this program. If not, see <http://www.gnu.org/licenses/>.
  23. close all
  24. clear
  25. clc
  27. function port_return=Execute_main(Sim_Path,f_max,show_structure,port_active)
  28.   % initialize
  29.   physical_constants;
  30.   %f_max = 1e9;
  32.   % Importing printed circuit board
  33.   CSX = InitCSX();
  34.   CSX = ImportHyperLynx(CSX, 'small_loop_antenna.hyp');
  36.   % Adding a component
  37.   [pad1_material, pad1_start, pad1_stop] = GetHyperLynxPort(CSX, 'C1.1');
  38.   [pad2_material, pad2_start, pad2_stop] = GetHyperLynxPort(CSX, 'C1.2');
  40.   c1_height = pad1_stop - pad1_start;
  41.   c1_start = [(pad1_start(1)+pad1_stop(1))/2, pad1_start(2), pad1_start(3)];
  42.   c1_stop = [(pad2_start(1)+pad2_stop(1))/2, pad2_stop(2), pad2_stop(3)+c1_height];
  44.   CSX = AddLumpedElement( CSX, 'Capacitor', 0, 'Caps', 1, 'C', 1e-12);
  45.   CSX = AddBox( CSX, 'Capacitor', 0, c1_start, c1_stop );
  47.   % Adding excitation and load
  48.   [port1_material, port1_start, port1_stop] = GetHyperLynxPort(CSX, 'TP1.1');
  49.   [gnd1_material, gnd1_start, gnd1_stop] = GetHyperLynxPort(CSX, 'TP3.1');
  50.   [port2_material, port2_start, port2_stop] = GetHyperLynxPort(CSX, 'TP2.1');
  51.   [gnd2_material, gnd2_start, gnd2_stop] = GetHyperLynxPort(CSX, 'TP4.1');
  53.   if(port_active==1)
  54.     [CSX, port{1}] = AddLumpedPort( CSX, 999, 1, 50, gnd1_start, port1_stop, [0 0 -1], true);
  55.     [CSX, port{2}] = AddLumpedPort( CSX, 999, 2, 50, gnd2_start, port2_stop, [0 0 -1]);
  56.   else
  57.     [CSX, port{1}] = AddLumpedPort( CSX, 999, 1, 50, gnd1_start, port1_stop, [0 0 -1]);
  58.     [CSX, port{2}] = AddLumpedPort( CSX, 999, 2, 50, gnd2_start, port2_stop, [0 0 -1], true);
  59.   endif
  61.   % Setting up a mesh
  62.   unit = GetUnits(CSX);
  63.   substrate_epr = GetEpsilon(CSX);
  64.   resolution = c0 / f_max / sqrt(substrate_epr) / unit / 25;
  65.   AirBox = c0 / f_max / unit / 25;
  67.   z_top = port1_start(3);
  68.   z_bottom = gnd1_start(3);
  69.   z_middle = (z_top+z_bottom)/2;
  71.   mesh.x = [];
  72.   mesh.y = [];
  73.   mesh.z = [ z_middle ];
  75.   mesh = DetectEdges(CSX, mesh);
  77.   mesh.x = [min(mesh.x)-AirBox max(mesh.x)+AirBox mesh.x];
  78.   mesh.y = [min(mesh.y)-AirBox max(mesh.y)+AirBox mesh.y];
  79.   mesh.z = [min(mesh.z)-AirBox max(mesh.z)+2*AirBox mesh.z];
  81.   %mesh = SmoothMesh(mesh, resolution);
  82.   mesh = SmoothMesh(mesh, resolution, 'algorithm', [ 1 ]);
  84.   % Setting boundary conditions
  86.   FDTD = InitFDTD();
  87.   FDTD = SetGaussExcite(FDTD, f_max/2, f_max/2);
  89.   BC = {'PML_8' 'PML_8' 'PML_8' 'PML_8' 'PML_8' 'PML_8'};
  90.   FDTD = SetBoundaryCond(FDTD, BC );
  92.   mesh = AddPML(mesh, 8);
  94.   CSX = DefineRectGrid(CSX, unit, mesh);
  96.   % Simulation
  98.   %Sim_Path = 'tmp';
  99.   Sim_CSX = 'msl.xml';
  101.   disp([ 'Estimated simulation runtime: 6500 timesteps' ]); % inform user this may take a while...
  103.   WriteOpenEMS([Sim_Path '/' Sim_CSX], FDTD, CSX);
  104.   if(show_structure>0)
  105.     CSXGeomPlot([Sim_Path '/' Sim_CSX]);
  106.   endif
  107.     RunOpenEMS(Sim_Path, Sim_CSX);
  108.   port_return=port;
  109. endfunction
  110. f_max = 1e9;
  111. Sim_Path1='tmp1'
  112. [status, message, messageid] = rmdir(Sim_Path1, 's'); % clear previous directory
  113. [status, message, messageid] = mkdir(Sim_Path1 ); % create empty simulation folder
  114. Sim_Path2='tmp2'
  115. [status, message, messageid] = rmdir(Sim_Path2, 's'); % clear previous directory
  116. [status, message, messageid] = mkdir(Sim_Path2 ); % create empty simulation folder
  118. port=Execute_main(Sim_Path1,f_max,1,1)
  119. % Calculating the s-parameters
  120. f = linspace( 100e6, f_max, 1601 );
  121. port = calcPort( port, Sim_Path1, f, 'RefImpedance', 50);
  122. s11 = port{1}.uf.ref./ port{1}.uf.inc;
  123. s21 = port{2}.uf.ref./ port{1}.uf.inc;
  125. port=Execute_main(Sim_Path2,f_max,0,2)
  126. % Calculating the s-parameters
  127. port = calcPort( port, Sim_Path2, f, 'RefImpedance', 50);
  128. s12 = port{1}.uf.ref./ port{2}.uf.inc;
  129. s22 = port{2}.uf.ref./ port{2}.uf.inc;
  131. s_param=[];
  132. s_param(1,1,:) = s11;
  133. s_param(1,2,:) = s12;
  134. s_param(2,1,:) = s21;
  135. s_param(2,2,:) = s22;
  136. write_touchstone('s',f,s_param,'out.s2p');
  138. close all
  140. semilogx(f/1e6,20*log10(abs(s11)),'k-','LineWidth',2);
  141. hold on;
  142. grid on;
  143. semilogx(f/1e6,20*log10(abs(s21)),'r--','LineWidth',2);
  144. legend('S_{11}','S_{21}');
  145. ylabel('S-Parameter (dB)','FontSize',12);
  146. xlabel('frequency (GHz) \rightarrow','FontSize',12);
  147. ylim([-80 10]);
  148. print('sparam.png', '-dpng');
  150. % not truncated
いちおう、今回の解析に使ったレイアウトファイルも載せておく。
  1. {VERSION=2.14}
  2. {UNITS=ENGLISH LENGTH}
  4. {BOARD "small_loop_antenna.kicad_pcb"
  5.   (PERIMETER_SEGMENT X1=6.100000000 Y1=4.050000000 X2=4.900000000 Y2=4.050000000)
  6.   (PERIMETER_SEGMENT X1=4.900000000 Y1=4.050000000 X2=4.900000000 Y2=2.850000000)
  7.   (PERIMETER_SEGMENT X1=4.900000000 Y1=2.850000000 X2=6.100000000 Y2=2.850000000)
  8.   (PERIMETER_SEGMENT X1=6.100000000 Y1=2.850000000 X2=6.100000000 Y2=4.050000000)
  9. }
  11. {STACKUP
  12.   (SIGNAL T=0.00137795 P=0 C=1.724e-08 L="F.Cu" M=COPPER)
  13.   (DIELECTRIC T=0.03937 C=4.4 L="DE_F.Cu" M="FR4")
  14.   (SIGNAL T=0.00137795 P=0 C=1.724e-08 L="B.Cu" M=COPPER)
  15. }
  17. {DEVICES
  18.   (? REF="TP4" L="B.Cu")
  19.   (? REF="TP3" L="B.Cu")
  20.   (? REF="TP2" L="F.Cu")
  21.   (? REF="TP1" L="F.Cu")
  22.   (? REF="C1" L="F.Cu")
  23. }
  25. {PADSTACK=0, 0.000000000
  26.   ("B.Cu", 1, 0.039370079, 0.039370079, 0.0, M)
  27. }
  29. {PADSTACK=1, 0.000000000
  30.   ("F.Cu", 1, 0.039370079, 0.039370079, 180.0, M)
  31. }
  33. {PADSTACK=2, 0.000000000
  34.   ("F.Cu", 2, 0.022047244, 0.024409449, 180.0, M)
  35. }
  37. {NET="GND"
  38.   (PIN X=5.4500000000 Y=4.0000000000 R="TP4.1" P=0)
  39.   (PIN X=5.5500000000 Y=4.0000000000 R="TP3.1" P=0)
  40.   (SEG X1=5.4500000000 Y1=4.0000000000 X2=5.5500000000 Y2=4.0000000000 W=0.0393700787 L="B.Cu")
  41. }
  43. {NET="Net-(C1-Pad1)"
  44.   (PIN X=5.5500000000 Y=4.0000000000 R="TP1.1" P=1)
  45.   (PIN X=5.4122047244 Y=3.9000000000 R="C1.1" P=2)
  46.   (SEG X1=5.4122047244 Y1=3.8622047244 X2=5.4000000000 Y2=3.8500000000 W=0.0236220472 L="F.Cu")
  47.   (SEG X1=5.4122047244 Y1=3.9000000000 X2=5.4122047244 Y2=3.8622047244 W=0.0236220472 L="F.Cu")
  48.   (SEG X1=5.9000000000 Y1=3.7500000000 X2=5.6000000000 Y2=3.7500000000 W=0.0393700787 L="F.Cu")
  49.   (SEG X1=5.6000000000 Y1=3.7500000000 X2=5.5500000000 Y2=3.8000000000 W=0.0393700787 L="F.Cu")
  50.   (SEG X1=6.0000000000 Y1=3.6500000000 X2=5.9000000000 Y2=3.7500000000 W=0.0393700787 L="F.Cu")
  51.   (SEG X1=6.0000000000 Y1=3.1000000000 X2=6.0000000000 Y2=3.6500000000 W=0.0393700787 L="F.Cu")
  52.   (SEG X1=5.5500000000 Y1=3.8000000000 X2=5.5500000000 Y2=4.0000000000 W=0.0393700787 L="F.Cu")
  53.   (SEG X1=5.4000000000 Y1=3.8000000000 X2=5.4000000000 Y2=3.8500000000 W=0.0393700787 L="F.Cu")
  54.   (SEG X1=5.3500000000 Y1=3.7500000000 X2=5.4000000000 Y2=3.8000000000 W=0.0393700787 L="F.Cu")
  55.   (SEG X1=5.1000000000 Y1=3.7500000000 X2=5.3500000000 Y2=3.7500000000 W=0.0393700787 L="F.Cu")
  56.   (SEG X1=5.0000000000 Y1=3.1000000000 X2=5.0000000000 Y2=3.6500000000 W=0.0393700787 L="F.Cu")
  57.   (SEG X1=5.0000000000 Y1=3.6500000000 X2=5.1000000000 Y2=3.7500000000 W=0.0393700787 L="F.Cu")
  58.   (SEG X1=5.1000000000 Y1=3.0000000000 X2=5.0000000000 Y2=3.1000000000 W=0.0393700787 L="F.Cu")
  59.   (SEG X1=5.9000000000 Y1=3.0000000000 X2=5.1000000000 Y2=3.0000000000 W=0.0393700787 L="F.Cu")
  60.   (SEG X1=6.0000000000 Y1=3.1000000000 X2=5.9000000000 Y2=3.0000000000 W=0.0393700787 L="F.Cu")
  61. }
  63. {NET="Net-(C1-Pad2)"
  64.   (PIN X=5.4500000000 Y=4.0000000000 R="TP2.1" P=1)
  65.   (PIN X=5.4500000000 Y=3.9000000000 R="C1.2" P=2)
  66.   (SEG X1=5.4500000000 Y1=3.9000000000 X2=5.4500000000 Y2=4.0000000000 W=0.0236220472 L="F.Cu")
  67. }
で、実行すると、1portあたり約1時間の計2時間で終了。いちおう、出力はこうなる。
で、出来上がったSパラを見てみる。が、Octaveでうまいこと表示させりゃいいんだけど、ここで、なぜかPythonに戻って、scikit-rfで表示する。コードはこんな感じ。
  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'
  11.     freq=ANT1.frequency
  12.     tl_media = rf.DefinedGammaZ0(freq, z0=50)
  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)
  17.     cnx=[[(port1,0),(ANT1,0)],
  18.          [(ANT1,1),(port2,0)],
  19.          [(gnd,0)]]
  20.     cir=rf.Circuit(cnx)
  21.     ntw=cir.network
  23.     ntw.frequency.unit='MHz'
  24.     fig=plt.figure()
  26.     ax1=fig.add_subplot(2,2,1)
  27.     ntw.plot_s_smith(ax=ax1,m=0, n=0, lw=2)
  29.     ax2_1=fig.add_subplot(2,2,2)
  30.     ax2_2=ax2_1.twinx()
  31.     ntw.plot_s_db(ax=ax2_1,m=1, n=0, lw=2, show_legend=False)
  32.     ntw.plot_s_vswr(ax=ax2_2,m=0, n=0, lw=2, show_legend=False, color='orangered')
  33.     ax2_2.set_ylim(1,6)
  34.     handler1, label1 = ax2_1.get_legend_handles_labels()
  35.     handler2, label2 = ax2_2.get_legend_handles_labels()
  36.     ax2_1.legend(handler1 + handler2, label1 + label2, loc=2, borderaxespad=0.)
  38.     ax3=fig.add_subplot(2,2,3)
  39.     ntw.plot_s_db(ax=ax3,m=0, n=1, lw=2)
  41.     ax4=fig.add_subplot(2,2,4)
  42.     ntw.plot_s_smith(ax=ax4,m=1, n=1, lw=2)
  43.     fig.tight_layout()
  44.     plt.show()
  46. if __name__ == "__main__":
  47.     main()
で、実行すると、
たぶん、いけてるんじゃないかな。
at
Labels: ,

0 件のコメント:

コメントを投稿

登録: コメントの投稿 (Atom)