diff --git a/FDTD/engine.cpp b/FDTD/engine.cpp index b87d1b0..277234f 100644 --- a/FDTD/engine.cpp +++ b/FDTD/engine.cpp @@ -79,14 +79,12 @@ bool Engine::IterateTS(unsigned int iterTS) } //soft voltage excitation here (E-field excite) - for (unsigned int n=0;nE_Ex_Count;++n) + for (unsigned int n=0;nE_Exc_Count;++n) { - exc_pos = (int)numTS - (int)Op->E_Ex_delay[n]; + exc_pos = (int)numTS - (int)Op->E_Exc_delay[n]; exc_pos*= (exc_pos>0 && exc_pos<(int)Op->ExciteLength); // if (n==0) cerr << numTS << " => " << Op->ExciteSignal[exc_pos] << endl; - volt[0][Op->E_Ex_index[0][n]][Op->E_Ex_index[1][n]][Op->E_Ex_index[2][n]] += Op->E_Ex_amp[0][n]*Op->ExciteSignal[exc_pos]; - volt[1][Op->E_Ex_index[0][n]][Op->E_Ex_index[1][n]][Op->E_Ex_index[2][n]] += Op->E_Ex_amp[1][n]*Op->ExciteSignal[exc_pos]; - volt[2][Op->E_Ex_index[0][n]][Op->E_Ex_index[1][n]][Op->E_Ex_index[2][n]] += Op->E_Ex_amp[2][n]*Op->ExciteSignal[exc_pos]; + volt[Op->E_Exc_dir[n]][Op->E_Exc_index[0][n]][Op->E_Exc_index[1][n]][Op->E_Exc_index[2][n]] += Op->E_Exc_amp[n]*Op->ExciteSignal[exc_pos]; } //current updates diff --git a/FDTD/operator.cpp b/FDTD/operator.cpp index 91816ed..b232509 100644 --- a/FDTD/operator.cpp +++ b/FDTD/operator.cpp @@ -35,7 +35,9 @@ void Operator::Init() CSX = NULL; ExciteSignal = NULL; - E_Ex_delay = NULL; + E_Exc_delay = NULL; + E_Exc_amp=NULL; + E_Exc_dir=NULL; vv=NULL; vi=NULL; iv=NULL; @@ -43,8 +45,7 @@ void Operator::Init() for (int n=0;n<3;++n) { discLines[n]=NULL; - E_Ex_amp[n]=NULL; - E_Ex_index[n]=NULL; + E_Exc_index[n]=NULL; } MainOp=NULL; @@ -62,7 +63,9 @@ void Operator::Init() void Operator::Reset() { delete[] ExciteSignal; - delete[] E_Ex_delay; + delete[] E_Exc_delay; + delete[] E_Exc_dir; + delete[] E_Exc_amp; Delete_N_3DArray(vv,numLines); Delete_N_3DArray(vi,numLines); Delete_N_3DArray(iv,numLines); @@ -70,8 +73,7 @@ void Operator::Reset() for (int n=0;n<3;++n) { delete[] discLines[n]; - delete[] E_Ex_amp[n]; - delete[] E_Ex_index[n]; + delete[] E_Exc_index[n]; } delete MainOp; delete DualOp; @@ -487,8 +489,8 @@ bool Operator::Calc_ECPos(int n, unsigned int* pos, double* inEC) inEC[1] += 0; } - inEC[0]*=gridDelta/fabs(delta)/4*__EPS0__; - inEC[1]*=gridDelta/fabs(delta)/4; + inEC[0]*=gridDelta/fabs(delta)/4.0*__EPS0__; + inEC[1]*=gridDelta/fabs(delta)/4.0; //******************************* mu,sigma averaging *****************************// //shift down @@ -530,8 +532,8 @@ bool Operator::Calc_ECPos(int n, unsigned int* pos, double* inEC) inEC[3] = 0; } - inEC[2] = gridDelta * fabs(deltaP*deltaPP) * 2 * __MUE0__ / inEC[2]; - if (inEC[3]) inEC[3]=gridDelta*fabs(deltaP*deltaPP) * 2 / inEC[3]; + inEC[2] = gridDelta * fabs(deltaP*deltaPP) * 2.0 * __MUE0__ / inEC[2]; + if (inEC[3]) inEC[3]=gridDelta*fabs(deltaP*deltaPP) * 2.0 / inEC[3]; return true; } @@ -646,79 +648,86 @@ bool Operator::CalcEFieldExcitation() { if (dT==0) return false; vector vIndex[3]; - vector vExcit[3]; + vector vExcit; vector vDelay; + vector vDir; unsigned int ipos; - unsigned int pos[3]; + int pos[3]; double coord[3]; + double delta[3]; + double amp=0; for (pos[2]=0;pos[2]GetIndexDelta(2,pos[2])); for (pos[1]=0;pos[1]GetIndexDelta(1,pos[1])); for (pos[0]=0;pos[0]GetIndexDelta(0,pos[0])); coord[0] = discLines[0][pos[0]]; coord[1] = discLines[1][pos[1]]; coord[2] = discLines[2][pos[2]]; // CSProperties* prop = CSX->GetPropertyByCoordPriority(coord,(CSProperties::PropertyType)(CSProperties::ELECTRODE | CSProperties::METAL)); - CSProperties* prop = CSX->GetPropertyByCoordPriority(coord,(CSProperties::PropertyType)(CSProperties::ELECTRODE)); - if (prop) + CSProperties* prop = NULL; + for (int n=0;n<3;++n) { - CSPropElectrode* elec = prop->ToElectrode(); - if (elec!=NULL) - if ((elec->GetExcitType()==0) || (elec->GetExcitType()==1)) //soft or hard E-Field excite! + coord[n]+=delta[n]*0.5; + CSProperties* prop = CSX->GetPropertyByCoordPriority(coord,(CSProperties::PropertyType)(CSProperties::ELECTRODE)); + if (prop) + { + CSPropElectrode* elec = prop->ToElectrode(); + if (elec!=NULL) { - vDelay.push_back((unsigned int)(elec->GetDelay()/dT)); - for (int n=0;n<3;++n) + if ((elec->GetActiveDir(n)) && (pos[n]GetActiveDir(n)) && (pos[n]GetWeightedExcitation(n,coord)*delta*gridDelta); - else - vExcit[n].push_back(0); - if ((elec->GetExcitType()==1) && (elec->GetActiveDir(n))) //hard excite + amp = elec->GetWeightedExcitation(n,coord)*delta[n]*gridDelta; + if (amp!=0) { - vv[(n+1)%3][pos[0]][pos[1]][pos[2]] = 0; - vi[(n+1)%3][pos[0]][pos[1]][pos[2]] = 0; - vv[(n+2)%3][pos[0]][pos[1]][pos[2]] = 0; - vi[(n+2)%3][pos[0]][pos[1]][pos[2]] = 0; + vExcit.push_back(amp); + vDelay.push_back((unsigned int)(elec->GetDelay()/dT)); + vDir.push_back(n); + vIndex[0].push_back(pos[0]); + vIndex[1].push_back(pos[1]); + vIndex[2].push_back(pos[2]); + } + if (elec->GetExcitType()==1) //hard excite + { + vv[n][pos[0]][pos[1]][pos[2]] = 0; + vi[n][pos[0]][pos[1]][pos[2]] = 0; } } } + } + coord[n]-=delta[n]*0.5; } } } } - E_Ex_Count = vIndex[0].size(); - if (E_Ex_Count==0) + + E_Exc_Count = vExcit.size(); + cerr << "Operator::CalcEFieldExcitation: Found number of excitations points: " << E_Exc_Count << endl; + if (E_Exc_Count==0) cerr << "No E-Field excitation found!" << endl; for (int n=0;n<3;++n) { - delete[] E_Ex_index[n]; - E_Ex_index[n] = new unsigned int[E_Ex_Count]; - for (unsigned int i=0;i