From 0ad528e6b86a77e0bac157f5e825ec3ded13d8fc Mon Sep 17 00:00:00 2001
From: Thorsten Liebig <Thorsten.Liebig@gmx.de>
Date: Wed, 25 Jul 2012 09:41:30 +0200
Subject: [PATCH] operator: use newton iteration to calculate numerical phase
 velocity

---
 FDTD/operator.cpp             | 46 ++++++++++++++++++++++++++++++++---
 matlab/Tutorials/RCS_Sphere.m |  2 +-
 2 files changed, 43 insertions(+), 5 deletions(-)

diff --git a/FDTD/operator.cpp b/FDTD/operator.cpp
index 1c8c31b..e917c22 100644
--- a/FDTD/operator.cpp
+++ b/FDTD/operator.cpp
@@ -1743,15 +1743,15 @@ void Operator::DeleteExtension(Operator_Extension* op_ext)
 
 double Operator::CalcNumericPhaseVelocity(unsigned int start[3], unsigned int stop[3], double propDir[3], float freq) const
 {
-	double phv=__C0__;
-
 	double average_mesh_disc[3];
-//	double k=2*PI*freq/__C0__;
+
+	//calculate average mesh deltas
 	for (int n=0;n<3;++n)
 	{
 		average_mesh_disc[n] = fabs(GetDiscLine(n,start[n])-GetDiscLine(n,stop[n]))*GetGridDelta() / (fabs(stop[n]-start[n]));
 	}
 
+	// if propagation is in a Cartesian direction, return analytic solution
 	for (int n=0;n<3;++n)
 	{
 		int nP = (n+1)%3;
@@ -1763,7 +1763,45 @@ double Operator::CalcNumericPhaseVelocity(unsigned int start[3], unsigned int st
 		}
 	}
 
-	cerr << "Operator::CalcNumericPhaseVelocity: Warning, propagation direction not in Cartesian direction, assuming phase velocity to be c0" << endl;
+	// else, do an newton iterative estimation
+	double k0=2*PI*freq/__C0__;
+	double k=k0;
+	double RHS = pow(sin(2*PI*freq*dT/2)/__C0__/dT,2);
+	double fk=1,fdk=0;
+	double old_phv=0;
+	double phv=__C0__;
+	double err_est = 1e-6;
+	int it_count=0;
+	while (fabs(old_phv-phv)>err_est)
+	{
+		++it_count;
+		old_phv=phv;
+		fk=0;
+		fdk=0;
+		for (int n=0;n<3;++n)
+		{
+			fk+= pow(sin(propDir[n]*k*average_mesh_disc[n]/2)/average_mesh_disc[n],2);
+			fdk+= propDir[n]*sin(propDir[n]*k*average_mesh_disc[n]/2)*cos(propDir[n]*k*average_mesh_disc[n]/2)/average_mesh_disc[n];
+		}
+		fk -= RHS;
+		k-=fk/fdk;
+
+		// do not allow a speed greater than c0 due to a numerical inaccuracy
+		if (k<k0)
+			k=k0;
+
+		phv=2*PI*freq/k;
+
+		//abort if iteration count is getting to high!
+		if (it_count>99)
+		{
+			cerr << "Operator::CalcNumericPhaseVelocity: Error, newton iteration estimation can't find a solution!!" << endl;
+			break;
+		}
+	}
+
+	if (g_settings.GetVerboseLevel()>1)
+		cerr << "Operator::CalcNumericPhaseVelocity: Newton iteration estimated solution: " << phv/__C0__ << "*c0 in " <<  it_count << " iterations." << endl;
 
 	return phv;
 }
diff --git a/matlab/Tutorials/RCS_Sphere.m b/matlab/Tutorials/RCS_Sphere.m
index 51a6555..0aeb130 100644
--- a/matlab/Tutorials/RCS_Sphere.m
+++ b/matlab/Tutorials/RCS_Sphere.m
@@ -54,7 +54,7 @@ CSX = AddSphere(CSX,'sphere',10,[0 0 0],sphere.rad);
 k_dir = [cos(inc_angle) sin(inc_angle) 0]; % plane wave direction
 E_dir = [0 0 1]; % plane wave polarization --> E_z
 
-CSX = AddPlaneWaveExcite(CSX, 'plane_wave', k_dir, E_dir);
+CSX = AddPlaneWaveExcite(CSX, 'plane_wave', k_dir, E_dir, f0);
 start = [-PW_Box/2 -PW_Box/2 -PW_Box/2];
 stop  = -start;
 CSX = AddBox(CSX, 'plane_wave', 0, start, stop);