2024年12月13日金曜日

openEMSを再び(1)

 

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

以前、使ってみようとしてうまくいかんかったんだけど、もっかいチャレンジしてみる。この手のソフトって超絶高価すぎる売り物か、すごく限られたモデルしか解析できないフリーなものかのどっちかしかないと思う。openEMSは見た感じ、かなり自由に解析できるのにフリーってのが魅力なので、なんとか使えるようになりたい。
まずは、総本山(https://www.openems.de/)に行く。見た感じ、LinuxなやつはBuildせんかいってことになっている。"Linux Build instructions"をつつく。
そうすると、もうそのままやればいい感じに書いてあるが、ここで立ち止まる(そう考える人はまずいないと思う)。実はUbuntuはlibvtkのバージョンが違っているっぽくて、ここに書いているものそのままやればいいわけじゃない。もちろんSynapticで検索しながらやればlibvtk7じゃなくてlibvtk9いれればいいかな？ってところにたどり着く。が、じつはそういうのちゃんと書いてある。で、
で、
で、
ってなってるので、この通りに全部やる。といいたいが、liboctave-devはaptできないかもしれない。たぶんoctave-devで良いはず。
pythonライブラリのインストールはvenvに入れるってなると、
source ~/work_flash/venv1/bin/activate
ってして(activateがどこにいるかはそれぞれの環境による)から、
pip install numpy matplotlib cython3 python3-h5py
のほうがいいと思う。
で、やり切ったら、ページの最下部からNextを押す。
Clone, Build and Installってのが開く。ここはこの通りにやる。
で、pathを通すので~/.profileに
export PATH=”$HOME/opt/openEMS/bin:$PATH”
を追記して、
source ~/.profile
で反映。
で、octaveにもpath設定があるらしいので、
octave --gui
でoctaveを起動。
「編集」-->「Set Path」
で、
ってのを追加する。で、保存してreload。
で、tutorialをやってみるが、Documentationのページからだとpythonのtutorialしかなくて、今回はなぜかoctaveでやってみたいので、総本山のLegacy Wikiをつつく。
で、左のTutorialsをつつくと、そこそこ出てくるので、今回は、こちら
をやってみる。
コードはこちら
  1. close all
  2. clear
  3. clc
  5. %% Setup the physical constants and antenna parameters
  6. physical_constants;
  7. unit = 1e-3; % all length in mm
  9. % patch width in x-direction
  10. patch.width = 30; % resonant length
  11. % patch length in y-direction
  12. patch.length = 40;
  14. %substrate setup
  15. substrate.epsR = 3.38;
  16. substrate.kappa = 1e-3 * 2*pi*2.45e9 * EPS0*substrate.epsR;
  17. substrate.width = 60;
  18. substrate.length = 60;
  19. substrate.thickness = 1.524;
  20. substrate.cells = 4;
  22. %setup feeding
  23. feed.pos = -5.5; %feeding position in x-direction
  24. feed.width = 2; %feeding port width
  25. feed.R = 50; %feed resistance
  27. % size of the simulation box
  28. SimBox = [200 200 100];
  29. %% Setup FDTD parameters, excitation signals & boundary conditions
  30. f0 = 2e9; % center frequency
  31. fc = 1e9; % 20 dB corner frequency
  32. FDTD = InitFDTD('NrTS', 30000 );
  33. FDTD = SetGaussExcite( FDTD, f0, fc );
  34. BC = {'MUR' 'MUR' 'MUR' 'MUR' 'MUR' 'MUR'}; % boundary conditions
  35. FDTD = SetBoundaryCond( FDTD, BC );
  36. %% Setup the CSXCAD mesh
  37. max_res = c0 / (f0+fc) / unit / 20; % cell size: lambda/20
  38. CSX = InitCSX();
  40. %create fixed lines for the simulation box, substrate and port
  41. mesh.x = [-SimBox(1)/2 SimBox(1)/2 -substrate.width/2 substrate.width/2 -patch.width/2 patch.width/2 feed.pos];
  42. mesh.x = SmoothMeshLines( mesh.x, max_res, 1.4); % create a smooth mesh between specified fixed mesh lines
  44. mesh.y = [-SimBox(2)/2 SimBox(2)/2 -substrate.length/2 substrate.length/2 -feed.width/2 feed.width/2 -patch.length/2 patch.length/2];
  45. mesh.y = SmoothMeshLines( mesh.y, max_res, 1.4 );
  47. %create fixed lines for the simulation box and given number of lines inside the substrate
  48. mesh.z = [-SimBox(3)/2 linspace(0,substrate.thickness,substrate.cells) SimBox(3)/2 ];
  49. mesh.z = SmoothMeshLines( mesh.z, max_res, 1.4 );
  51. CSX = DefineRectGrid( CSX, unit, mesh );
  52. %% Setup the geometry
  53. %% create patch
  54. CSX = AddMetal( CSX, 'patch' ); % create a perfect electric conductor (PEC)
  55. start = [-patch.width/2 -patch.length/2 substrate.thickness];
  56. stop = [ patch.width/2 patch.length/2 substrate.thickness];
  57. CSX = AddBox(CSX,'patch',10,start,stop); % add a box-primitive to the metal property 'patch'
  59. %% create substrate
  60. CSX = AddMaterial( CSX, 'substrate' );
  61. CSX = SetMaterialProperty( CSX, 'substrate', 'Epsilon', substrate.epsR, 'Kappa', substrate.kappa );
  62. start = [-substrate.width/2 -substrate.length/2 0];
  63. stop = [ substrate.width/2 substrate.length/2 substrate.thickness];
  64. CSX = AddBox( CSX, 'substrate', 0, start, stop );
  66. %% create ground (same size as substrate)
  67. CSX = AddMetal( CSX, 'gnd' ); % create a perfect electric conductor (PEC)
  68. start(3)=0;
  69. stop(3) =0;
  70. CSX = AddBox(CSX,'gnd',10,start,stop);
  71. %% Setup the feeding port as a lumped port with 50 Ohms
  72. start = [feed.pos-.1 -feed.width/2 0];
  73. stop = [feed.pos+.1 +feed.width/2 substrate.thickness];
  74. [CSX port] = AddLumpedPort(CSX, 5 ,1 ,feed.R, start, stop, [0 0 1], true);
  75. %% Add a nf2ff box
  76. SimBox = SimBox - max_res * 4; %reduced SimBox size for nf2ff box
  77. [CSX nf2ff] = CreateNF2FFBox(CSX, 'nf2ff', -SimBox/2, SimBox/2);
  78. %% Create simulation folder
  79. %% Write xml simulation file
  80. %% Visualize the Geometry using AppCSXCAD
  81. %% Run openEMS
  82. %% prepare simulation folder
  83. Sim_Path = 'tmp';
  84. Sim_CSX = 'patch_ant.xml';
  86. [status, message, messageid] = rmdir( Sim_Path, 's' ); % clear previous directory
  87. [status, message, messageid] = mkdir( Sim_Path ); % create empty simulation folder
  89. %% write openEMS compatible xml-file
  90. WriteOpenEMS( [Sim_Path '/' Sim_CSX], FDTD, CSX );
  92. %% show the structure
  93. CSXGeomPlot( [Sim_Path '/' Sim_CSX] );
  95. %% run openEMS
  96. RunOpenEMS( Sim_Path, Sim_CSX );
  98. %% Post-Processing
  99. %% Read in port voltages and currents
  100. %% postprocessing & do the plots
  101. freq = linspace( max([1e9,f0-fc]), f0+fc, 501 );
  102. port = calcPort(port, Sim_Path, freq);
  103. %% Calculate & plot antenna input-impedance
  104. Zin = port.uf.tot ./ port.if.tot;
  105. s11 = port.uf.ref ./ port.uf.inc;
  106. P_in = 0.5 * port.uf.inc .* conj( port.if.inc ); % antenna feed power
  108. % plot feed point impedance
  109. figure
  110. plot( freq/1e6, real(Zin), 'k-', 'Linewidth', 2 );
  111. hold on
  112. grid on
  113. plot( freq/1e6, imag(Zin), 'r--', 'Linewidth', 2 );
  114. title( 'feed point impedance' );
  115. xlabel( 'frequency f / MHz' );
  116. ylabel( 'impedance Z_{in} / Ohm' );
  117. legend( 'real', 'imag' );
  118. %% Calculate & Plot S-Parameter and accepted power
  119. % plot reflection coefficient S11
  120. figure
  121. plot( freq/1e6, 20*log10(abs(s11)), 'k-', 'Linewidth', 2 );
  122. grid on
  123. title( 'reflection coefficient S_{11}' );
  124. xlabel( 'frequency f / MHz' );
  125. ylabel( 'reflection coefficient |S_{11}|' );
  127. drawnow
  128. %% Calculate & Plot antenna parameter & radiation pattern
  129. %% NFFF contour plots %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  130. %find resonance frequency from s11
  131. f_res_ind = find(s11==min(s11));
  132. f_res = freq(f_res_ind);
  134. % calculate the far field at phi=0 degrees and at phi=90 degrees
  135. disp( 'calculating far field at phi=[0 90] deg...' );
  137. nf2ff = CalcNF2FF(nf2ff, Sim_Path, f_res, [-180:2:180]*pi/180, [0 90]*pi/180);
  139. % display power and directivity
  140. disp( ['radiated power: Prad = ' num2str(nf2ff.Prad) ' Watt']);
  141. disp( ['directivity: Dmax = ' num2str(nf2ff.Dmax) ' (' num2str(10*log10(nf2ff.Dmax)) ' dBi)'] );
  142. disp( ['efficiency: nu_rad = ' num2str(100*nf2ff.Prad./real(P_in(f_res_ind))) ' %']);
  144. % normalized directivity as polar plot
  145. figure
  146. polarFF(nf2ff,'xaxis','theta','param',[1 2],'normalize',1)
  148. % log-scale directivity plot
  149. figure
  150. plotFFdB(nf2ff,'xaxis','theta','param',[1 2])
  151. % conventional plot approach
  152. % plot( nf2ff.theta*180/pi, 20*log10(nf2ff.E_norm{1}/max(nf2ff.E_norm{1}(:)))+10*log10(nf2ff.Dmax));
  154. drawnow
  156. %%
  157. disp( 'calculating 3D far field pattern and dumping to vtk (use Paraview to visualize)...' );
  158. thetaRange = (0:2:180);
  159. phiRange = (0:2:360) - 180;
  160. nf2ff = CalcNF2FF(nf2ff, Sim_Path, f_res, thetaRange*pi/180, phiRange*pi/180,'Verbose',1,'Outfile','3D_Pattern.h5');
  162. figure
  163. plotFF3D(nf2ff,'logscale',-20);
  166. E_far_normalized = nf2ff.E_norm{1} / max(nf2ff.E_norm{1}(:)) * nf2ff.Dmax;
  167. DumpFF2VTK([Sim_Path '/3D_Pattern.vtk'],E_far_normalized,thetaRange,phiRange,'scale',1e-3);
で、このまま実行すると、とまった？って思うけど、コマンドウィンドウでyes/noを入力しないと進まないようになっている。で、yesってして、
こういうの出るので、
ざっくり見て、閉じる(閉じないと先に進まない)
そすっと計算を進めている風になって、その後、こういうの
が出る。で、裏ではこのまま計算が進むんだけど、
まぁ、こんな感じで止まっちゃう。何かがないので何かを生成しようとしているっぽいけど失敗している感じ。
/home/hoge/opt/openEMS/share/openEMS/matlab/setup.m
に何やろうちしているか書いてある。
mkoctfile("h5readatt_octave.cc", ["-L" hdf5lib_dir], ["-I" hdf5inc_dir], "-lhdf5") これをやりたいらしい。ならば先にやって差し上げておこう。hdf5lib_dirとhdf5inc_dirはエラーメッセージから読み取ると、両方とも/usr/lib/x86_64-linux-gnu/hdf5/serialらしい。
なので、
cd /home/hoge/opt/openEMS/share/openEMS/matlab
mkoctfile h5readatt_octave.cc -L/usr/lib/x86_64-linux-gnu/hdf5/serial -I/usr/include/hdf5/serial -lhdf5
ってする。そすっと完了して、h5readatt_octave.octってのができる。で、再び実行してみると、最後までやりきって、こういう感じ
になる。たぶんうまくいった。
よし、何かに使ってみる、、、実物で検証できないテーマだな、、、家では狭すぎるし、田舎の誰もいない公園でも、法に触れることはやっちゃいかんし、、、

てかひさびさの更新！約２ヶ月！あれもこれもと思っているうちに時間がたってしまった。そして、考えていたうちのどれでもないテーマで書いてしまった、、、挙動不審の極み。認知症に足首くらいまで入ってるな。
at
Labels: ,

0 件のコメント:

コメントを投稿

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