diff --git a/FDTD/extensions/operator_ext_tfsf.cpp b/FDTD/extensions/operator_ext_tfsf.cpp
index 6bee532..8299430 100644
--- a/FDTD/extensions/operator_ext_tfsf.cpp
+++ b/FDTD/extensions/operator_ext_tfsf.cpp
@@ -43,6 +43,8 @@ void Operator_Ext_TFST::Init()
 				m_CurrAmp[n][l][c]=NULL;
 			}
 
+	m_Frequency = 0.0;
+	m_PhVel = __C0__;
 	Operator_Extension::Init();
 }
 
@@ -95,6 +97,7 @@ bool Operator_Ext_TFST::BuildExtension()
 		return false;
 	}
 
+	m_PhVel = __C0__;
 	CSPropExcitation* elec=NULL;
 	CSProperties* prop=NULL;
 	CSPrimitives* prim=NULL;
@@ -144,6 +147,20 @@ bool Operator_Ext_TFST::BuildExtension()
 			return false;
 		}
 
+		m_Frequency = elec->GetFrequency();
+//		if (m_Frequency<=0)
+//			m_Frequency = m_Op->GetExcitationSignal()->GetFrequencyOfInterest();
+		if (m_Frequency<=0)
+			m_PhVel=__C0__;
+		else
+			m_PhVel=m_Op->CalcNumericPhaseVelocity(m_Start,m_Stop,m_PropDir,m_Frequency);
+
+		if ((m_PhVel<0) || (m_PhVel>__C0__) || isnan(m_PhVel))
+		{
+			cerr << "Operator_Ext_TFST::BuildExtension: Warning, invalid phase velocity found, resetting to c0! " << endl;
+			m_PhVel = __C0__;
+		}
+
 		double origin[3];
 		unsigned int ui_origin[3];
 		for (int n=0;n<3;++n)
@@ -163,11 +180,11 @@ bool Operator_Ext_TFST::BuildExtension()
 
 			if (m_IncLow[n])
 			{
-				ui_origin[n] = m_Start[n];
+				ui_origin[n] = m_Start[n]-1;
 			}
 			else
 			{
-				ui_origin[n] = m_Stop[n];
+				ui_origin[n] = m_Stop[n]+1;
 			}
 			origin[n] = m_Op->GetDiscLine(n,ui_origin[n]);
 		}
@@ -244,16 +261,16 @@ bool Operator_Ext_TFST::BuildExtension()
 					if (m_ActiveDir[n][0])
 					{
 						m_Op->GetYeeCoords(nP,pos,coord,false);
-						dist = fabs(coord[0]-origin[0])*m_PropDir[0]+fabs(coord[1]-origin[1])*m_PropDir[1]+fabs(coord[2]-origin[2])*m_PropDir[2];
-						delay = dist*unit/__C0__/dT;
+						dist = fabs((coord[0]-origin[0])*m_PropDir[0])+fabs((coord[1]-origin[1])*m_PropDir[1])+fabs((coord[2]-origin[2])*m_PropDir[2]);
+						delay = dist*unit/m_PhVel/dT;
 						m_maxDelay = max((unsigned int)delay,m_maxDelay);
 						m_CurrDelay[n][0][1][ui_pos] = floor(delay);
 						m_CurrDelayDelta[n][0][1][ui_pos] = delay - floor(delay);
 						m_CurrAmp[n][0][1][ui_pos] = m_E_Amp[nP]*m_Op->GetEdgeLength(nP,pos);
 
 						m_Op->GetYeeCoords(nPP,pos,coord,false);
-						dist = fabs(coord[0]-origin[0])*m_PropDir[0]+fabs(coord[1]-origin[1])*m_PropDir[1]+fabs(coord[2]-origin[2])*m_PropDir[2];
-						delay = dist*unit/__C0__/dT;
+						dist = fabs((coord[0]-origin[0])*m_PropDir[0])+fabs((coord[1]-origin[1])*m_PropDir[1])+fabs((coord[2]-origin[2])*m_PropDir[2]);
+						delay = dist*unit/m_PhVel/dT;
 						m_maxDelay = max((unsigned int)delay,m_maxDelay);
 						m_CurrDelay[n][0][0][ui_pos] = floor(delay);
 						m_CurrDelayDelta[n][0][0][ui_pos] = delay - floor(delay);
@@ -268,16 +285,16 @@ bool Operator_Ext_TFST::BuildExtension()
 					{
 						pos[n] = m_Stop[n];
 						m_Op->GetYeeCoords(nP,pos,coord,false);
-						dist = fabs(coord[0]-origin[0])*m_PropDir[0]+fabs(coord[1]-origin[1])*m_PropDir[1]+fabs(coord[2]-origin[2])*m_PropDir[2];
-						delay = dist*unit/__C0__/dT;
+						dist = fabs((coord[0]-origin[0])*m_PropDir[0])+fabs((coord[1]-origin[1])*m_PropDir[1])+fabs((coord[2]-origin[2])*m_PropDir[2]);
+						delay = dist*unit/m_PhVel/dT;
 						m_maxDelay = max((unsigned int)delay,m_maxDelay);
 						m_CurrDelay[n][1][1][ui_pos] = floor(delay);
 						m_CurrDelayDelta[n][1][1][ui_pos] = delay - floor(delay);
 						m_CurrAmp[n][1][1][ui_pos] = m_E_Amp[nP]*m_Op->GetEdgeLength(nP,pos);
 
 						m_Op->GetYeeCoords(nPP,pos,coord,false);
-						dist = fabs(coord[0]-origin[0])*m_PropDir[0]+fabs(coord[1]-origin[1])*m_PropDir[1]+fabs(coord[2]-origin[2])*m_PropDir[2];
-						delay = dist*unit/__C0__/dT;
+						dist = fabs((coord[0]-origin[0])*m_PropDir[0])+fabs((coord[1]-origin[1])*m_PropDir[1])+fabs((coord[2]-origin[2])*m_PropDir[2]);
+						delay = dist*unit/m_PhVel/dT;
 						m_maxDelay = max((unsigned int)delay,m_maxDelay);
 						m_CurrDelay[n][1][0][ui_pos] = floor(delay);
 						m_CurrDelayDelta[n][1][0][ui_pos] = delay - floor(delay);
@@ -287,20 +304,10 @@ bool Operator_Ext_TFST::BuildExtension()
 						m_CurrAmp[n][1][1][ui_pos]*=m_Op->GetIV(nPP,pos);
 					}
 
-					if (m_IncLow[n])
-					{
-						if (m_ActiveDir[n][0])
-							m_CurrAmp[n][0][0][ui_pos]*=-1;
-						if (m_ActiveDir[n][1])
-							m_CurrAmp[n][1][1][ui_pos]*=-1;
-					}
-					else
-					{
-						if (m_ActiveDir[n][0])
-							m_CurrAmp[n][0][1][ui_pos]*=-1;
-						if (m_ActiveDir[n][1])
-							m_CurrAmp[n][1][0][ui_pos]*=-1;
-					}
+					if (m_ActiveDir[n][0])
+						m_CurrAmp[n][0][0][ui_pos]*=-1;
+					if (m_ActiveDir[n][1])
+						m_CurrAmp[n][1][1][ui_pos]*=-1;
 
 					if (pos[nP]==m_Stop[nP])
 					{
@@ -322,16 +329,16 @@ bool Operator_Ext_TFST::BuildExtension()
 					if (m_ActiveDir[n][0])
 					{
 						m_Op->GetYeeCoords(nP,pos,coord,true);
-						dist = fabs(coord[0]-origin[0])*m_PropDir[0]+fabs(coord[1]-origin[1])*m_PropDir[1]+fabs(coord[2]-origin[2])*m_PropDir[2];
-						delay = dist*unit/__C0__/dT;
+						dist = fabs((coord[0]-origin[0])*m_PropDir[0])+fabs((coord[1]-origin[1])*m_PropDir[1])+fabs((coord[2]-origin[2])*m_PropDir[2]);
+						delay = dist*unit/m_PhVel/dT + 1.0;
 						m_maxDelay = max((unsigned int)delay,m_maxDelay);
 						m_VoltDelay[n][0][1][ui_pos] = floor(delay);
 						m_VoltDelayDelta[n][0][1][ui_pos] = delay - floor(delay);
 						m_VoltAmp[n][0][1][ui_pos] = m_H_Amp[nP]*m_Op->GetEdgeLength(nP,pos,true);
 
 						m_Op->GetYeeCoords(nPP,pos,coord,true);
-						dist = fabs(coord[0]-origin[0])*m_PropDir[0]+fabs(coord[1]-origin[1])*m_PropDir[1]+fabs(coord[2]-origin[2])*m_PropDir[2];
-						delay = dist*unit/__C0__/dT;
+						dist = fabs((coord[0]-origin[0])*m_PropDir[0])+fabs((coord[1]-origin[1])*m_PropDir[1])+fabs((coord[2]-origin[2])*m_PropDir[2]);
+						delay = dist*unit/m_PhVel/dT + 1.0;
 						m_maxDelay = max((unsigned int)delay,m_maxDelay);
 						m_VoltDelay[n][0][0][ui_pos] = floor(delay);
 						m_VoltDelayDelta[n][0][0][ui_pos] = delay - floor(delay);
@@ -346,16 +353,16 @@ bool Operator_Ext_TFST::BuildExtension()
 					if (m_ActiveDir[n][1])
 					{
 						m_Op->GetYeeCoords(nP,pos,coord,true);
-						dist = fabs(coord[0]-origin[0])*m_PropDir[0]+fabs(coord[1]-origin[1])*m_PropDir[1]+fabs(coord[2]-origin[2])*m_PropDir[2];
-						delay = dist*unit/__C0__/dT;
+						dist = fabs((coord[0]-origin[0])*m_PropDir[0])+fabs((coord[1]-origin[1])*m_PropDir[1])+fabs((coord[2]-origin[2])*m_PropDir[2]);
+						delay = dist*unit/m_PhVel/dT + 1.0;
 						m_maxDelay = max((unsigned int)delay,m_maxDelay);
 						m_VoltDelay[n][1][1][ui_pos] = floor(delay);
 						m_VoltDelayDelta[n][1][1][ui_pos] = delay - floor(delay);
 						m_VoltAmp[n][1][1][ui_pos] = m_H_Amp[nP]*m_Op->GetEdgeLength(nP,pos,true);
 
 						m_Op->GetYeeCoords(nPP,pos,coord,true);
-						dist = fabs(coord[0]-origin[0])*m_PropDir[0]+fabs(coord[1]-origin[1])*m_PropDir[1]+fabs(coord[2]-origin[2])*m_PropDir[2];
-						delay = dist*unit/__C0__/dT;
+						dist = fabs((coord[0]-origin[0])*m_PropDir[0])+fabs((coord[1]-origin[1])*m_PropDir[1])+fabs((coord[2]-origin[2])*m_PropDir[2]);
+						delay = dist*unit/m_PhVel/dT + 1.0;
 						m_maxDelay = max((unsigned int)delay,m_maxDelay);
 						m_VoltDelay[n][1][0][ui_pos] = floor(delay);
 						m_VoltDelayDelta[n][1][0][ui_pos] = delay - floor(delay);
@@ -365,20 +372,10 @@ bool Operator_Ext_TFST::BuildExtension()
 						m_VoltAmp[n][1][1][ui_pos]*=m_Op->GetVI(nPP,pos);
 					}
 
-					if (m_IncLow[n])
-					{
-						if (m_ActiveDir[n][1])
-							m_VoltAmp[n][1][0][ui_pos]*=-1;
-						if (m_ActiveDir[n][0])
-							m_VoltAmp[n][0][1][ui_pos]*=-1;
-					}
-					else
-					{
-						if (m_ActiveDir[n][0])
-							m_VoltAmp[n][0][0][ui_pos]*=-1;
-						if (m_ActiveDir[n][1])
-							m_VoltAmp[n][1][1][ui_pos]*=-1;
-					}
+					if (m_ActiveDir[n][1])
+						m_VoltAmp[n][1][0][ui_pos]*=-1;
+					if (m_ActiveDir[n][0])
+						m_VoltAmp[n][0][1][ui_pos]*=-1;
 
 					if (pos[nP]==m_Stop[nP])
 					{
@@ -418,6 +415,7 @@ void Operator_Ext_TFST::ShowStat(ostream &ostr) const
 	Operator_Extension::ShowStat(ostr);
 	cout << "Active directions\t: " << "(" << m_ActiveDir[0][0] << "/" << m_ActiveDir[0][1] << ", " << m_ActiveDir[1][0] << "/" << m_ActiveDir[1][1]  << ", " << m_ActiveDir[2][0] << "/" << m_ActiveDir[2][1]  << ")" << endl;
 	cout << "Propagation direction\t: " << "(" << m_PropDir[0] << ", " << m_PropDir[1] << ", " << m_PropDir[2] << ")" << endl;
+	cout << "Rel. propagation speed\t: " << m_PhVel/__C0__ << "*c0  @ " << m_Frequency << " Hz" << endl;
 	cout << "E-field amplitude\t: " << "(" << m_E_Amp[0] << ", " << m_E_Amp[1] << ", " << m_E_Amp[2] << ")" << endl;
 	cout << "H-field amplitude\t: " << "(" << m_H_Amp[0]*__Z0__ << ", " << m_H_Amp[1]*__Z0__ << ", " << m_H_Amp[2]*__Z0__ << ")/Z0" << endl;
 	cout << "Box Dimensions\t\t: " << m_numLines[0] << " x " << m_numLines[1] << " x " << m_numLines[2] << endl;
diff --git a/FDTD/extensions/operator_ext_tfsf.h b/FDTD/extensions/operator_ext_tfsf.h
index c95f33b..61e1f03 100644
--- a/FDTD/extensions/operator_ext_tfsf.h
+++ b/FDTD/extensions/operator_ext_tfsf.h
@@ -60,6 +60,9 @@ protected:
 	double m_E_Amp[3];
 	double m_H_Amp[3];
 
+	double m_Frequency;
+	double m_PhVel;
+
 	unsigned int m_maxDelay;
 
 	// array setup [direction][low/high][component][ <mesh_position> ]
diff --git a/FDTD/operator.cpp b/FDTD/operator.cpp
index 170f4ab..1c8c31b 100644
--- a/FDTD/operator.cpp
+++ b/FDTD/operator.cpp
@@ -1740,3 +1740,30 @@ 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__;
+	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]));
+	}
+
+	for (int n=0;n<3;++n)
+	{
+		int nP = (n+1)%3;
+		int nPP = (n+2)%3;
+		if ((fabs(propDir[n])==1) && (propDir[nP]==0) && (propDir[nPP]==0))
+		{
+			double kx = asin(average_mesh_disc[0]/__C0__/dT*sin(2*PI*freq*dT/2))*2/average_mesh_disc[0];
+			return 2*PI*freq/kx;
+		}
+	}
+
+	cerr << "Operator::CalcNumericPhaseVelocity: Warning, propagation direction not in Cartesian direction, assuming phase velocity to be c0" << endl;
+
+	return phv;
+}
diff --git a/FDTD/operator.h b/FDTD/operator.h
index cb5a544..30b4ad9 100644
--- a/FDTD/operator.h
+++ b/FDTD/operator.h
@@ -147,6 +147,8 @@ public:
 
 	Operator_Ext_Excitation* GetExcitationExtension() const;
 
+	virtual double CalcNumericPhaseVelocity(unsigned int start[3], unsigned int stop[3], double propDir[3], float freq) const;
+
 protected:
 	//! use New() for creating a new Operator
 	Operator();
diff --git a/matlab/Tutorials/RCS_Sphere.m b/matlab/Tutorials/RCS_Sphere.m
index 9d27006..51a6555 100644
--- a/matlab/Tutorials/RCS_Sphere.m
+++ b/matlab/Tutorials/RCS_Sphere.m
@@ -20,13 +20,16 @@ unit = 1e-3; % all length in mm
 
 sphere.rad = 200;
 
+inc_angle = 0 /180*pi; %incident angle (to x-axis) in rad
+
 % size of the simulation box
 SimBox = 1000;
 PW_Box = 750;
 
 %% setup FDTD parameter & excitation function
-f_start = 200e6; % start frequency
+f_start =  50e6; % start frequency
 f_stop = 1000e6; % stop  frequency
+f0 = 500e6;
 
 FDTD = InitFDTD( 30000 );
 FDTD = SetGaussExcite( FDTD, 0.5*(f_start+f_stop), 0.5*(f_stop-f_start) );
@@ -34,12 +37,12 @@ BC = [1 1 1 1 1 1]*3;  % set boundary conditions
 FDTD = SetBoundaryCond( FDTD, BC );
 
 %% setup CSXCAD geometry & mesh
-% currently, openEMS cannot automatically generate a mesh
 max_res = c0 / f_stop / unit / 20; % cell size: lambda/20
 CSX = InitCSX();
 
-%create fixed lines for the simulation box, substrate and port
-mesh.x = SmoothMeshLines([-SimBox/2 0 SimBox/2], max_res);
+%create mesh
+smooth_mesh = SmoothMeshLines([0 SimBox/2], max_res);
+mesh.x = unique([-smooth_mesh smooth_mesh]);
 mesh.y = mesh.x;
 mesh.z = mesh.x;
 
@@ -48,7 +51,7 @@ CSX = AddMetal( CSX, 'sphere' ); % create a perfect electric conductor (PEC)
 CSX = AddSphere(CSX,'sphere',10,[0 0 0],sphere.rad);
 
 %% plane wave excitation
-k_dir = [1 0 0]; % plane wave direction --> x-direction
+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);
@@ -92,26 +95,44 @@ RunOpenEMS( Sim_Path, Sim_CSX);
 disp('Use Paraview to display the elctric fields dumped by openEMS');
 
 %%
-freq = 500e6;
-EF = ReadUI( 'et', Sim_Path, freq ); % time domain/freq domain voltage
+EF = ReadUI( 'et', Sim_Path, f0 ); % time domain/freq domain voltage
 Pin = 0.5*norm(E_dir)^2/Z0 .* abs(EF.FD{1}.val).^2;
 
 %%
-nf2ff = CalcNF2FF(nf2ff, Sim_Path, freq, [-180:2:180]*pi/180, [0 90]*pi/180,'Mode',1);
-RCS = 4*pi./Pin(1).*nf2ff.P_rad{1}(:,1);
-polar(nf2ff.theta,RCS);
-xlabel('z -->');
-ylabel('x -->');
+nf2ff = CalcNF2FF(nf2ff, Sim_Path, f0, pi/2, [-180:2:180]*pi/180, 'Mode',1);
+RCS = 4*pi./Pin(1).*nf2ff.P_rad{1}(:);
+polar(nf2ff.phi,RCS);
+xlabel('x -->');
+ylabel('y -->');
 hold on
 grid on
 
+drawnow
+
 %%
-disp( 'calculating 3D far field pattern and dumping to vtk (use Paraview to visualize)...' );
-thetaRange = (0:2:180);
-phiRange = (0:2:360) - 180;
-nf2ff = CalcNF2FF(nf2ff, Sim_Path, freq, thetaRange*pi/180, phiRange*pi/180,'Verbose',1,'Outfile','3D_Pattern.h5','Mode',1);
+freq = linspace(f_start,f_stop,100);
+EF = ReadUI( 'et', Sim_Path, freq ); % time domain/freq domain voltage
+Pin = 0.5*norm(E_dir)^2/Z0 .* abs(EF.FD{1}.val).^2;
 
-RCS = 4*pi./Pin(1).*nf2ff.P_rad{1};
-DumpFF2VTK([Sim_Path '/3D_Pattern.vtk'],RCS,thetaRange,phiRange,1e-3);
+nf2ff = CalcNF2FF(nf2ff, Sim_Path, freq, pi/2, pi+inc_angle, 'Mode',1);
+for fn=1:numel(freq)
+    back_scat(fn) = 4*pi./Pin(fn).*nf2ff.P_rad{fn}(1);
+end
 
-disp('Use Paraview to display the 3D scattering pattern');
+%%
+figure
+plot(freq/1e6,back_scat,'Linewidth',2);
+grid on;
+xlabel('frequency (MHz) \rightarrow');
+ylabel('RCS (m^2) \rightarrow');
+title('radar cross section');
+
+%%
+figure
+lambda = c0./freq;
+semilogy(sphere.rad*unit./lambda,back_scat/(pi*sphere.rad*unit*sphere.rad*unit),'Linewidth',2);
+ylim([10^-2 10^1])
+grid on;
+xlabel('sphere radius / wavelength \rightarrow');
+ylabel('RCS / (\pi a^2) \rightarrow');
+title('normalized radar cross section');