diff --git a/matlab/private/ReadNF2FF.m b/matlab/private/ReadNF2FF.m
index 78bec86..eb0bcd9 100644
--- a/matlab/private/ReadNF2FF.m
+++ b/matlab/private/ReadNF2FF.m
@@ -27,6 +27,7 @@ if isOctave
         nf2ff.E_theta{n} = double(hdf.nf2ff.E_theta.FD.(['f' int2str(n-1) '_real']) +1i*hdf.nf2ff.E_theta.FD.(['f' int2str(n-1) '_imag']) );
         nf2ff.E_phi{n} = double(hdf.nf2ff.E_phi.FD.(['f' int2str(n-1) '_real']) +1i*hdf.nf2ff.E_phi.FD.(['f' int2str(n-1) '_imag']) );
         nf2ff.E_norm{n} = double(sqrt(abs(nf2ff.E_theta{n}).^2+abs(nf2ff.E_phi{n}).^2));
+        nf2ff.P_rad{n} = double(hdf.nf2ff.P_rad.FD.(['f' int2str(n-1)]));
     end
 else
     % matlab compatibility to older versions (will be removed soon)
@@ -39,6 +40,7 @@ else
             nf2ff.E_theta{n} = double(hdf5read(file,['/nf2ff/E_theta/FD/f' int2str(n-1) '_real']) + 1i*hdf5read(file,['/nf2ff/E_theta/FD/f' int2str(n-1) '_imag']));
             nf2ff.E_phi{n} = double(hdf5read(file,['/nf2ff/E_phi/FD/f' int2str(n-1) '_real']) + 1i*hdf5read(file,['/nf2ff/E_phi/FD/f' int2str(n-1) '_imag']));
             nf2ff.E_norm{n} = double(sqrt(abs(nf2ff.E_theta{n}).^2+abs(nf2ff.E_phi{n}).^2));
+            nf2ff.P_rad{n} = double(hdf5read(file,['/nf2ff/P_rad/FD/f' int2str(n-1)]));
         end
         return
     end
@@ -51,5 +53,6 @@ else
         nf2ff.E_theta{n} = double(h5read(file,['/nf2ff/E_theta/FD/f' int2str(n-1) '_real']) + 1i*h5read(file,['/nf2ff/E_theta/FD/f' int2str(n-1) '_imag']));
         nf2ff.E_phi{n} = double(h5read(file,['/nf2ff/E_phi/FD/f' int2str(n-1) '_real']) + 1i*h5read(file,['/nf2ff/E_phi/FD/f' int2str(n-1) '_imag']));
         nf2ff.E_norm{n} = double(sqrt(abs(nf2ff.E_theta{n}).^2+abs(nf2ff.E_phi{n}).^2));
+        nf2ff.P_rad{n} = double(h5read(file,['/nf2ff/P_rad/FD/f' int2str(n-1)]));
     end
 end
diff --git a/nf2ff/nf2ff.cpp b/nf2ff/nf2ff.cpp
index 9772865..4d449a8 100644
--- a/nf2ff/nf2ff.cpp
+++ b/nf2ff/nf2ff.cpp
@@ -349,6 +349,27 @@ bool nf2ff::Write2HDF5(string filename)
 			}
 		}
 	}
+
+	//dump radiated power
+	hdf_file.SetCurrentGroup("/nf2ff/P_rad/FD");
+	for (size_t fn=0;fn<m_freq.size();++fn)
+	{
+		stringstream ss;
+		ss << "f" << fn;
+		pos = 0;
+		float** field_data = GetRadPower(fn);
+		for (size_t j=0;j<m_numPhi;++j)
+			for (size_t i=0;i<m_numTheta;++i)
+			{
+				buffer[pos++]=field_data[i][j];
+			}
+		if (hdf_file.WriteData(ss.str(),buffer,dim,datasize)==false)
+		{
+			delete[] buffer;
+			cerr << "nf2ff::Write2HDF5: Error writing field data" << endl;
+			return false;
+		}
+	}
 	delete[] buffer;
 
 	//write frequency attribute
@@ -357,7 +378,7 @@ bool nf2ff::Write2HDF5(string filename)
 	buffer = new float[m_freq.size()];
 	//write radiated power attribute
 	for (size_t fn=0;fn<m_freq.size();++fn)
-		buffer[fn] = GetRadPower(fn);
+		buffer[fn] = GetTotalRadPower(fn);
 	hdf_file.WriteAtrribute("/nf2ff", "Prad",buffer,m_freq.size());
 
 	//write max directivity attribute
diff --git a/nf2ff/nf2ff.h b/nf2ff/nf2ff.h
index e59b8ed..bad28a4 100644
--- a/nf2ff/nf2ff.h
+++ b/nf2ff/nf2ff.h
@@ -37,11 +37,12 @@ public:
 
 	bool AnalyseFile(string E_Field_file, string H_Field_file);
 
-	float GetRadPower(size_t f_idx) const {return m_nf2ff.at(f_idx)->GetRadPower();}
+	float GetTotalRadPower(size_t f_idx) const {return m_nf2ff.at(f_idx)->GetTotalRadPower();}
 	float GetMaxDirectivity(size_t f_idx) const {return m_nf2ff.at(f_idx)->GetMaxDirectivity();}
 
 	complex<float>** GetETheta(size_t f_idx) const {return m_nf2ff.at(f_idx)->GetETheta();}
 	complex<float>** GetEPhi(size_t f_idx) const {return m_nf2ff.at(f_idx)->GetEPhi();}
+	float** GetRadPower(size_t f_idx) const {return m_nf2ff.at(f_idx)->GetRadPower();}
 
 	//! Write results to a hdf5 file
 	bool Write2HDF5(string filename);
diff --git a/nf2ff/nf2ff_calc.h b/nf2ff/nf2ff_calc.h
index c2cc046..7c0ccc0 100644
--- a/nf2ff/nf2ff_calc.h
+++ b/nf2ff/nf2ff_calc.h
@@ -75,11 +75,12 @@ public:
 	nf2ff_calc(float freq, vector<float> theta, vector<float> phi, vector<float> center);
 	~nf2ff_calc();
 
-	float GetRadPower() const {return m_radPower;}
+	float GetTotalRadPower() const {return m_radPower;}
 	float GetMaxDirectivity() const {return m_maxDir;}
 
 	complex<float>** GetETheta() const {return m_E_theta;}
 	complex<float>** GetEPhi() const {return m_E_phi;}
+	float** GetRadPower() const {return m_P_rad;}
 
 	unsigned int GetNumThreads() const {return m_numThreads;}
 	void SetNumThreads(unsigned int n) {m_numThreads=n;}