diff --git a/Analyse/PlotVoltage.m b/Analyse/PlotVoltage.m new file mode 100644 index 0000000..34e037b --- /dev/null +++ b/Analyse/PlotVoltage.m @@ -0,0 +1,21 @@ +%close all; +clear all; +clc + +tmp = load('../tmp/u1'); + +t = tmp(:,1); +u = tmp(:,2); + +L=numel(t); + +subplot(2,1,1); +plot(t,u); + +dt=t(2)-t(1); + +f = (1:L)/L/dt; +fu = fft(u)/L; +subplot(2,1,2); +plot(f(1:L/2),abs(fu(1:L/2))); + diff --git a/FDTD/engine.cpp b/FDTD/engine.cpp index 263af8d..e0fd0f4 100644 --- a/FDTD/engine.cpp +++ b/FDTD/engine.cpp @@ -36,11 +36,11 @@ bool Engine::IterateTS(unsigned int iterTS) for (unsigned int iter=0;iternumLines[2]-1;++pos[2]) + for (pos[0]=1;pos[0]numLines[0]-1;++pos[0]) { for (pos[1]=1;pos[1]numLines[1]-1;++pos[1]) { - for (pos[0]=1;pos[0]numLines[0]-1;++pos[0]) + for (pos[2]=1;pos[2]numLines[2]-1;++pos[2]) { //do the updates here //for x @@ -78,15 +78,15 @@ bool Engine::IterateTS(unsigned int iterTS) //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]+1][pos[2]] - volt[2][pos[0]][pos[1]][pos[2]] - volt[1][pos[0]][pos[1]][pos[2]+1] + volt[1][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]+1] - volt[0][pos[0]][pos[1]][pos[2]] - volt[2][pos[0]+1][pos[1]][pos[2]] + volt[2][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 x 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]+1][pos[1]][pos[2]] - volt[1][pos[0]][pos[1]][pos[2]] - volt[0][pos[0]][pos[1]+1][pos[2]] + volt[0][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]]); } } } diff --git a/FDTD/processvoltage.cpp b/FDTD/processvoltage.cpp index 5cb5a02..95ad474 100644 --- a/FDTD/processvoltage.cpp +++ b/FDTD/processvoltage.cpp @@ -38,5 +38,5 @@ void ProcessVoltage::Process() } // cerr << voltage << endl; voltages.push_back(voltage); - file << Eng->numTS << "\t" << voltage << endl; + file << (double)Eng->numTS*Op->GetTimestep() << "\t" << voltage << endl; } diff --git a/main.cpp b/main.cpp index 46c6eaa..1ea811a 100644 --- a/main.cpp +++ b/main.cpp @@ -12,6 +12,7 @@ using namespace std; void BuildMSL(ContinuousStructure &CSX); +void BuildDipol(ContinuousStructure &CSX); int main(int argc, char *argv[]) { @@ -19,8 +20,7 @@ int main(int argc, char *argv[]) //*************** setup/read geometry ************// ContinuousStructure CSX; -// CSX.ReadFromXML("csx-files/MSL.xml"); - BuildMSL(CSX); + BuildDipol(CSX); //*************** setup operator ************// CartOperator cop; @@ -38,7 +38,7 @@ int main(int argc, char *argv[]) // cop.DumpOperator2File("tmp/Operator"); cerr << "Nyquist number of timesteps: " << cop.GetNyquistNum(fmax) << endl; - unsigned int NrIter = cop.GetNyquistNum(fmax); + unsigned int NrIter = cop.GetNyquistNum(fmax)/3; cerr << "Time for operator: " << difftime(OpDoneTime,startTime) << endl; @@ -50,11 +50,13 @@ int main(int argc, char *argv[]) //*************** setup processing ************// ProcessVoltage PV(&cop,&eng); PV.OpenFile("tmp/u1"); - double start[]={0,0,0}; - double stop[]={0,50,0}; + double start[]={-500,-150,0}; + double stop[]={-500,150,0}; PV.DefineStartStopCoord(start,stop); unsigned int maxIter = 5000; + PV.Process(); +// NrIter=200; //*************** simulate ************// for (unsigned int i=0;iSetKappa(0.001); + CSX.AddProperty(mat); + + CSPrimBox* matbox = new CSPrimBox(CSX.GetParameterSet(),mat); + matbox->SetCoord(0,-1000.0);matbox->SetCoord(1,1000.0); + matbox->SetCoord(2,-1000.0);matbox->SetCoord(3,1000.0); + matbox->SetCoord(4,-1000.0);matbox->SetCoord(5,1000.0); + CSX.AddPrimitive(matbox); + + CSPropElectrode* elec = new CSPropElectrode(CSX.GetParameterSet()); + elec->SetExcitation(1,1); + elec->SetExcitType(0); +// elec->SetDelay(2.0e-9); + CSX.AddProperty(elec); + + CSPrimBox* box = new CSPrimBox(CSX.GetParameterSet(),elec); + box->SetCoord(0,-1.0);box->SetCoord(1,1.0); + box->SetCoord(2,-75.0);box->SetCoord(3,75.0); + box->SetCoord(4,-1.0);box->SetCoord(5,1.0); + CSX.AddPrimitive(box); + + CSRectGrid* grid = CSX.GetGrid(); + + for (int n=-1000;n<=1000;n+=20) + grid->AddDiscLine(2,(double)n); + for (int n=-1000;n<=1000;n+=20) + grid->AddDiscLine(0,(double)n); + for (int n=-1000;n<=1000;n+=20) + grid->AddDiscLine(1,(double)n); + + grid->SetDeltaUnit(1e-3); + + CSX.Write2XML("tmp/Dipol.xml"); +} + void BuildMSL(ContinuousStructure &CSX) { // CSPropMaterial* mat = new CSPropMaterial(CSX.GetParameterSet()); @@ -109,14 +149,14 @@ void BuildMSL(ContinuousStructure &CSX) CSRectGrid* grid = CSX.GetGrid(); - for (int n=-100;n<=100;n+=10) - grid->AddDiscLine(2,(double)n); - for (int n=-100;n<=100;n+=10) + for (int n=-300;n<=300;n+=20) grid->AddDiscLine(0,(double)n); - for (int n=0;n<=100;n+=10) + for (int n=0;n<=300;n+=20) grid->AddDiscLine(1,(double)n); + for (int n=-1000;n<=1000;n+=20) + grid->AddDiscLine(2,(double)n); grid->SetDeltaUnit(1e-3); - CSX.Write2XML("csx-files/MSL.xml"); + CSX.Write2XML("tmp/MSL.xml"); }