diff --git a/FDTD/engine_cylinder.cpp b/FDTD/engine_cylinder.cpp index 9f979f3..ea8a5dc 100644 --- a/FDTD/engine_cylinder.cpp +++ b/FDTD/engine_cylinder.cpp @@ -26,18 +26,191 @@ Engine_Cylinder* Engine_Cylinder::New(const Operator_Cylinder* op) Engine_Cylinder::Engine_Cylinder(const Operator_Cylinder* op) : Engine(op) { + cyl_Op = op; + if (cyl_Op->GetClosedAlpha()) + { + ++numLines[1]; //necessary for dobled voltage and current line in alpha-dir, operator will return one smaller for correct post-processing + } } Engine_Cylinder::~Engine_Cylinder() { + Reset(); } void Engine_Cylinder::Init() { Engine::Init(); + +// if (cyl_Op->GetClosedAlpha()) +// { +// unsigned int lastLine = Op->numLines[1]-1; //number of last alpha-line +// unsigned int pos[3]; +// for (pos[0]=0;pos[0]numLines[0];++pos[0]) +// { +// for (int n=0;n<3;++n) +// { +// delete[] volt[n][pos[0]][lastLine]; +// volt[n][pos[0]][lastLine] = volt[n][pos[0]][0]; +// delete[] curr[n][pos[0]][lastLine]; +// curr[n][pos[0]][lastLine] = curr[n][pos[0]][0]; +// } +// } +// } } void Engine_Cylinder::Reset() { - Engine::Reset(); +// if (cyl_Op->GetClosedAlpha()) +// { +// unsigned int lastLine = Op->numLines[1]-1; //number of last alpha-line +// unsigned int pos[3]; +// for (pos[0]=0;pos[0]numLines[0];++pos[0]) +// { +// for (int n=0;n<3;++n) +// { +// volt[n][pos[0]][lastLine] = NULL; +// curr[n][pos[0]][lastLine] = NULL; +// } +// } +// } } + +inline void Engine_Cylinder::CloseAlphaVoltages() +{ + unsigned int pos[3]; + // copy voltages from last alpha-plane to first + unsigned int last_A_Line = numLines[1]-1; + for (pos[0]=0;pos[0]GetClosedAlpha()==false) + return Engine::IterateTS(iterTS); + + for (unsigned int iter=0;iterGetClosedAlpha()==false) +// return Engine::UpdateVoltages(); +// +// //voltage updates +// for (pos[0]=0;pos[0]numLines[0];++pos[0]) +// { +// shift[0]=pos[0]; +// for (pos[1]=1;pos[1]numLines[1];++pos[1]) +// { +// shift[1]=pos[1]; +// for (pos[2]=0;pos[2]numLines[2];++pos[2]) +// { +// shift[2]=pos[2]; +// //do the updates here +// //for x +// volt[0][pos[0]][pos[1]][pos[2]] *= Op->vv[0][pos[0]][pos[1]][pos[2]]; +// volt[0][pos[0]][pos[1]][pos[2]] += Op->vi[0][pos[0]][pos[1]][pos[2]] * ( curr[2][pos[0]][pos[1]][pos[2]] - curr[2][pos[0]][pos[1]-shift[1]][pos[2]] - curr[1][pos[0]][pos[1]][pos[2]] + curr[1][pos[0]][pos[1]][pos[2]-shift[2]]); +// +// //for y +// volt[1][pos[0]][pos[1]][pos[2]] *= Op->vv[1][pos[0]][pos[1]][pos[2]]; +// volt[1][pos[0]][pos[1]][pos[2]] += Op->vi[1][pos[0]][pos[1]][pos[2]] * ( curr[0][pos[0]][pos[1]][pos[2]] - curr[0][pos[0]][pos[1]][pos[2]-shift[2]] - curr[2][pos[0]][pos[1]][pos[2]] + curr[2][pos[0]-shift[0]][pos[1]][pos[2]]); +// +// //for z +// volt[2][pos[0]][pos[1]][pos[2]] *= Op->vv[2][pos[0]][pos[1]][pos[2]]; +// volt[2][pos[0]][pos[1]][pos[2]] += Op->vi[2][pos[0]][pos[1]][pos[2]] * ( curr[1][pos[0]][pos[1]][pos[2]] - curr[1][pos[0]-shift[0]][pos[1]][pos[2]] - curr[0][pos[0]][pos[1]][pos[2]] + curr[0][pos[0]][pos[1]-shift[1]][pos[2]]); +// } +// } +// +// // copy voltages from last alpha-plane to first +// unsigned int last_A_Line = Op->numLines[1]-1; +// for (pos[2]=0;pos[2]numLines[2];++pos[2]) +// { +// volt[0][pos[0]][0][pos[2]] = volt[0][pos[0]][last_A_Line][pos[2]]; +// volt[1][pos[0]][0][pos[2]] = volt[1][pos[0]][last_A_Line][pos[2]]; +// volt[2][pos[0]][0][pos[2]] = volt[2][pos[0]][last_A_Line][pos[2]]; +// } +// +// } +//} +// +//inline void Engine_Cylinder::UpdateCurrents() +//{ +// if (cyl_Op->GetClosedAlpha()==false) +// return Engine::UpdateCurrents(); +// +// unsigned int pos[3]; +// for (pos[0]=0;pos[0]numLines[0]-1;++pos[0]) +// { +// for (pos[1]=0;pos[1]numLines[1]-1;++pos[1]) +// { +// for (pos[2]=0;pos[2]numLines[2]-1;++pos[2]) +// { +// //do the updates here +// //for x +// curr[0][pos[0]][pos[1]][pos[2]] *= Op->ii[0][pos[0]][pos[1]][pos[2]]; +// curr[0][pos[0]][pos[1]][pos[2]] += Op->iv[0][pos[0]][pos[1]][pos[2]] * ( volt[2][pos[0]][pos[1]][pos[2]] - volt[2][pos[0]][pos[1]+1][pos[2]] - volt[1][pos[0]][pos[1]][pos[2]] + volt[1][pos[0]][pos[1]][pos[2]+1]); +// +// //for y +// curr[1][pos[0]][pos[1]][pos[2]] *= Op->ii[1][pos[0]][pos[1]][pos[2]]; +// curr[1][pos[0]][pos[1]][pos[2]] += Op->iv[1][pos[0]][pos[1]][pos[2]] * ( volt[0][pos[0]][pos[1]][pos[2]] - volt[0][pos[0]][pos[1]][pos[2]+1] - volt[2][pos[0]][pos[1]][pos[2]] + volt[2][pos[0]+1][pos[1]][pos[2]]); +// +// //for z +// curr[2][pos[0]][pos[1]][pos[2]] *= Op->ii[2][pos[0]][pos[1]][pos[2]]; +// curr[2][pos[0]][pos[1]][pos[2]] += Op->iv[2][pos[0]][pos[1]][pos[2]] * ( volt[1][pos[0]][pos[1]][pos[2]] - volt[1][pos[0]+1][pos[1]][pos[2]] - volt[0][pos[0]][pos[1]][pos[2]] + volt[0][pos[0]][pos[1]+1][pos[2]]); +// } +// } +// // copy currents from first alpha-plane to last +// unsigned int last_A_Line = Op->numLines[1]-1; +// for (pos[2]=0;pos[2]numLines[2]-1;++pos[2]) +// { +// curr[0][pos[0]][last_A_Line][pos[2]] = curr[0][pos[0]][0][pos[2]]; +// curr[1][pos[0]][last_A_Line][pos[2]] = curr[1][pos[0]][0][pos[2]]; +// curr[2][pos[0]][last_A_Line][pos[2]] = curr[2][pos[0]][0][pos[2]]; +// } +// } +//} diff --git a/FDTD/engine_cylinder.h b/FDTD/engine_cylinder.h index af4d122..60ab5a3 100644 --- a/FDTD/engine_cylinder.h +++ b/FDTD/engine_cylinder.h @@ -25,13 +25,21 @@ class Engine_Cylinder : public Engine { public: static Engine_Cylinder* New(const Operator_Cylinder* op); + virtual ~Engine_Cylinder(); virtual void Init(); virtual void Reset(); + //!Iterate a number of timesteps + virtual bool IterateTS(unsigned int iterTS); + protected: Engine_Cylinder(const Operator_Cylinder* op); - ~Engine_Cylinder(); + + virtual inline void CloseAlphaVoltages(); + virtual inline void CloseAlphaCurrents(); + + const Operator_Cylinder* cyl_Op; }; #endif // ENGINE_CYLINDER_H diff --git a/FDTD/operator_cylinder.cpp b/FDTD/operator_cylinder.cpp index 51ad60b..c643162 100644 --- a/FDTD/operator_cylinder.cpp +++ b/FDTD/operator_cylinder.cpp @@ -45,6 +45,16 @@ void Operator_Cylinder::Reset() Operator::Reset(); } +inline unsigned int Operator_Cylinder::GetNumberOfLines(int ny) const +{ + //this is necessary for a correct field processing... cylindrical engine has to reset this by adding +1 + if (CC_closedAlpha && ny==1) + return numLines[1]-1; + + return numLines[ny]; +} + + bool Operator_Cylinder::SetGeometryCSX(ContinuousStructure* geo) { if (Operator::SetGeometryCSX(geo)==false) return false; @@ -53,10 +63,20 @@ bool Operator_Cylinder::SetGeometryCSX(ContinuousStructure* geo) // cerr << minmaxA -2*PI << " < " << (2*PI)/10/numLines[1] << endl; if (fabs(minmaxA-2*PI) < (2*PI)/10/numLines[1]) //check minmaxA smaller then a tenth of average alpha-width { - CC_closedAlpha = true; - --numLines[1]; cout << "Operator_Cylinder::SetGeometryCSX: Alpha is a full 2*PI => closed Cylinder..." << endl; - cerr << "Operator_Cylinder::SetGeometryCSX: closed cylinder not yet implemented... exit... " << endl; exit(1); + CC_closedAlpha = true; + discLines[1][numLines[1]-1] = discLines[1][0] + 2*PI; + cerr << "Operator_Cylinder::SetGeometryCSX: Warning, not handling the disc-line width and material averaging correctly yet for a closed cylinder..." << endl; + if (MainOp->GetIndexDelta(1,0)-MainOp->GetIndexDelta(1,numLines[1]-2) > (2*PI)/10/numLines[1]) + { + cerr << "Operator_Cylinder::SetGeometryCSX: first and last angle delta must be the same... deviation to large..." << MainOp->GetIndexDelta(1,0) - MainOp->GetIndexDelta(1,numLines[1]-2) << endl; + exit(1); + } + if (MainOp->GetIndexDelta(1,0)-MainOp->GetIndexDelta(1,numLines[1]-2) > 0) + { + cerr << "Operator_Cylinder::SetGeometryCSX: first and last angle delta must be the same... auto correction of deviation: " << MainOp->GetIndexDelta(1,0) - MainOp->GetIndexDelta(1,numLines[1]-2) << endl; + discLines[1][numLines[1]-2] = discLines[1][numLines[1]-1]-MainOp->GetIndexDelta(1,0); + } } else if (minmaxA>2*PI) {cerr << "Operator_Cylinder::SetGeometryCSX: Alpha Max-Min must not be larger than 2*PI!!!" << endl; Reset(); return false;} diff --git a/FDTD/operator_cylinder.h b/FDTD/operator_cylinder.h index de4c3ef..8a98a20 100644 --- a/FDTD/operator_cylinder.h +++ b/FDTD/operator_cylinder.h @@ -33,6 +33,11 @@ public: virtual void Reset(); + virtual unsigned int GetNumberOfLines(int ny) const; + + bool GetClosedAlpha() const {return CC_closedAlpha;} + bool GetR0Included() const {return CC_R0_included;} + protected: Operator_Cylinder(); virtual void Init(); diff --git a/matlab/examples/Coax_CylinderCoords.m b/matlab/examples/Coax_CylinderCoords.m new file mode 100644 index 0000000..94fca1b --- /dev/null +++ b/matlab/examples/Coax_CylinderCoords.m @@ -0,0 +1,129 @@ +close all; +clear all; +clc + +EPS0 = 8.85418781762e-12; +MUE0 = 1.256637062e-6; +C0 = 1/sqrt(EPS0*MUE0); +Z0 = sqrt(MUE0/EPS0); + +f0 = 0.5e9; +epsR = 1; + +abs_length = 500; +length = 3000; +port_dist = 1500; +rad_i = 100; +rad_a = 230; +max_mesh = 10; +max_alpha = max_mesh; +N_alpha = ceil(rad_a * 2*pi / max_alpha); +mesh_res = [max_mesh 2*pi/N_alpha max_mesh]; + +openEMS_Path = [pwd() '/../../']; +openEMS_opts = ''; +openEMS_opts = [openEMS_opts ' --disable-dumps']; +% openEMS_opts = [openEMS_opts ' --debug-material']; + +Sim_Path = 'tmp'; +Sim_CSX = 'coax.xml'; + +mkdir(Sim_Path); + +%setup FDTD parameter +FDTD = InitCylindricalFDTD(1e5,1e-5,'OverSampling',10); +FDTD = SetGaussExcite(FDTD,f0,f0); +BC = [0 0 1 1 0 0]; +FDTD = SetBoundaryCond(FDTD,BC); + +%setup CSXCAD geometry +CSX = InitCSX(); +mesh.x = rad_i : mesh_res(1) : rad_a; +mesh.y = linspace(0,2*pi,N_alpha); +% mesh.y = mesh.y + mesh_res(2)/2; +mesh.z = 0 : mesh_res(3) : length; +CSX = DefineRectGrid(CSX, 1e-3,mesh); + +%%%fake pml +finalKappa = 0.3/abs_length^4; +finalSigma = finalKappa*MUE0/EPS0/epsR; +CSX = AddMaterial(CSX,'pml'); +CSX = SetMaterialProperty(CSX,'pml','Kappa',finalKappa,'Epsilon',epsR); +CSX = SetMaterialProperty(CSX,'pml','Sigma',finalSigma); +CSX = SetMaterialWeight(CSX,'pml','Kappa',['pow(abs(z)-' num2str(length-abs_length) ',4)']); +CSX = SetMaterialWeight(CSX,'pml','Sigma',['pow(abs(z)-' num2str(length-abs_length) ',4)']); + +start = [rad_i mesh.y(1) length-abs_length]; +stop = [rad_a mesh.y(end) length]; +CSX = AddBox(CSX,'pml',0 ,start,stop); + + +CSX = AddMaterial(CSX,'fill'); +CSX = SetMaterialProperty(CSX,'fill','Epsilon',epsR); +start = [mesh.x(1) mesh.y(1) 0]; +stop = [mesh.x(end) mesh.y(end) length]; +CSX = AddBox(CSX,'fill',0 ,start,stop); + +start = [rad_i mesh.y(1) 0]; +stop = [rad_a mesh.y(end) 0]; + +CSX = AddExcitation(CSX,'excite',0,[1 0 0]); +weight{1} = '1/rho'; +weight{2} = 0; +weight{3} = 0; +CSX = SetExcitationWeight(CSX, 'excite', weight ); +CSX = AddBox(CSX,'excite',0 ,start,stop); + +%dump +CSX = AddDump(CSX,'Et_','DumpMode',0); +start = [mesh.x(1) , 0 , mesh.z(1)]; +stop = [mesh.x(end) , 0 , mesh.z(end)]; +CSX = AddBox(CSX,'Et_',0 , start,stop); + +CSX = AddDump(CSX,'Ht_','DumpType',1,'DumpMode',0); +CSX = AddBox(CSX,'Ht_',0,start,stop); + +% voltage calc (take a voltage average to be at the same spot as the +% current calculation) +CSX = AddProbe(CSX,'ut1_1',0); +start = [ rad_i 0 port_dist ];stop = [ rad_a 0 port_dist ]; +CSX = AddBox(CSX,'ut1_1', 0 ,start,stop); +CSX = AddProbe(CSX,'ut1_2',0); +start = [ rad_i 0 port_dist+mesh_res(3) ];stop = [ rad_a 0 port_dist+mesh_res(3) ]; +CSX = AddBox(CSX,'ut1_2', 0 ,start,stop); + +% current calc +CSX = AddProbe(CSX,'it1',1); +mid = 0.5*(rad_i+rad_a); +start = [ 0 mesh.y(1) port_dist ];stop = [ mid mesh.y(end) port_dist ]; +CSX = AddBox(CSX,'it1', 0 ,start,stop); + +%Write openEMS compatoble xml-file +WriteOpenEMS([Sim_Path '/' Sim_CSX],FDTD,CSX); + +%cd to working dir and run openEMS +savePath = pwd(); +cd(Sim_Path); %cd to working dir +command = [openEMS_Path 'openEMS.sh ' Sim_CSX ' ' openEMS_opts]; +disp(command); +system(command) +cd(savePath); + +UI = ReadUI({'ut1_1','ut1_2','it1'},'tmp/'); +u_f = (UI.FD{1}.val + UI.FD{2}.val)/2; %averaging voltages to fit current +i_f = UI.FD{3}.val; + +delta_t = UI.TD{3}.t(1) - UI.TD{1}.t(1); % half time-step (s) +i_f2 = i_f .* exp(-1i*2*pi*UI.FD{1}.f*delta_t); % compensate half time-step advance of H-field + +ZL = Z0/2/pi/sqrt(epsR)*log(rad_a/rad_i); %analytic line-impedance of a coax +plot(UI.FD{1}.f,ZL*ones(size(u_f)),'g'); +hold on; +grid on; +Z = u_f./i_f2; +plot(UI.FD{1}.f,real(Z),'Linewidth',2); +plot(UI.FD{1}.f,imag(Z),'r','Linewidth',2); +xlim([0 2*f0]); +legend('Z_L','\Re\{Z\}','\Im\{Z\}','Location','Best'); + +