process fields: new SAR calculation

todo: needs much testing and evaluation

Signed-off-by: Thorsten Liebig <Thorsten.Liebig@gmx.de>

Common/processfields.cpp

@@ -73,10 +73,28 @@ string ProcessFields::GetFieldNameByType(DumpType type)
 		return "RotH-Field";
 	case SAR_LOCAL_DUMP:
 		return "SAR-local";
+	case SAR_1G_DUMP:
+		return "SAR_1g";
+	case SAR_10G_DUMP:
+		return "SAR_10g";
+	case SAR_RAW_DATA:
+		return "SAR_raw_data";
 	}
 	return "unknown field";
 }
 
+bool ProcessFields::NeedConductivity() const
+{
+	switch (m_DumpType)
+	{
+	case J_FIELD_DUMP:
+		return true;
+	default:
+		return false;
+	}
+	return false;
+}
+
 void ProcessFields::InitProcess()
 {
 	if (Enabled==false) return;

Common/processfields.h

@@ -40,7 +40,7 @@ public:
 	  Current dump types are electric field (E_FIELD_DUMP), magnetic field (H_FIELD_DUMP),
 	  (conduction) electric current density (kappa*E) (J_FIELD_DUMP) and total current density (rotH)
 	  */
-	enum DumpType { E_FIELD_DUMP, H_FIELD_DUMP, J_FIELD_DUMP, ROTH_FIELD_DUMP, SAR_LOCAL_DUMP};
+	enum DumpType { E_FIELD_DUMP=0, H_FIELD_DUMP=1, J_FIELD_DUMP=2, ROTH_FIELD_DUMP=3, SAR_LOCAL_DUMP=20, SAR_1G_DUMP=21, SAR_10G_DUMP=22, SAR_RAW_DATA=29};
 
 	virtual string GetProcessingName() const {return "common field processing";}
 
@@ -66,7 +66,7 @@ public:
 	void SetDumpMode2Cell() {SetDumpMode(Engine_Interface_Base::CELL_INTERPOLATE);}
 
 	//! Set dump type: 0 for E-fields, 1 for H-fields, 2 for D-fields, 3 for B-fields, 4 for J-fields, etc...
-	void SetDumpType(DumpType type) {m_DumpType=type;}
+	virtual void SetDumpType(DumpType type) {m_DumpType=type;}
 
 	double CalcTotalEnergyEstimate() const;
 
@@ -74,6 +74,8 @@ public:
 
 	static string GetFieldNameByType(DumpType type);
 
+	virtual bool NeedConductivity() const;
+
 protected:
 	DumpType m_DumpType;
 	FileType m_fileType;

Common/processfields_sar.cpp

@@ -19,37 +19,53 @@
 #include "operator_base.h"
 #include "tools/vtk_file_writer.h"
 #include "tools/hdf5_file_writer.h"
+#include "tools/sar_calculation.h"
 
 ProcessFieldsSAR::ProcessFieldsSAR(Engine_Interface_Base* eng_if) : ProcessFieldsFD(eng_if)
 {
+	m_UseCellKappa = false;
 }
 
 ProcessFieldsSAR::~ProcessFieldsSAR()
 {
-	for (size_t n = 0; n<m_FD_Fields.size(); ++n)
+	for (size_t n = 0; n<m_E_FD_Fields.size(); ++n)
-	{
 		Delete_N_3DArray(m_E_FD_Fields.at(n),numLines);
-		Delete_N_3DArray(m_J_FD_Fields.at(n),numLines);
-	}
 	m_E_FD_Fields.clear();
+
+	for (size_t n = 0; n<m_J_FD_Fields.size(); ++n)
+		Delete_N_3DArray(m_J_FD_Fields.at(n),numLines);
 	m_J_FD_Fields.clear();
 }
 
+void ProcessFieldsSAR::SetDumpType(DumpType type)
+{
+	if (type==SAR_RAW_DATA)
+		m_UseCellKappa = true;
+	ProcessFieldsFD::SetDumpType(type);
+}
+
+bool ProcessFieldsSAR::NeedConductivity() const
+{
+	if (m_UseCellKappa)
+		return false;
+	return	true;
+}
+
 void ProcessFieldsSAR::SetSubSampling(unsigned int subSampleRate, int dir)
 {
-	cerr << "ProcessFieldsSAR::SetSubSampling: Warning: Defining a sub-sampling for SAR calculation is not recommended!!! Expect false results!" << endl;
+	UNUSED(subSampleRate);UNUSED(dir);
-	ProcessFieldsFD::SetSubSampling(subSampleRate,dir);
+	cerr << "ProcessFieldsSAR::SetSubSampling: Warning: Defining a sub-sampling for SAR calculation is not allowed!!! Skipped!" << endl;
 }
 
 void ProcessFieldsSAR::SetOptResolution(double optRes, int dir)
 {
-	cerr << "ProcessFieldsSAR::SetOptResolution: Warning: Defining a sub-sampling for SAR calculation is not recommended!!! Expect false results!" << endl;
+	UNUSED(optRes);UNUSED(dir);
-	ProcessFieldsFD::SetOptResolution(optRes,dir);
+	cerr << "ProcessFieldsSAR::SetOptResolution: Warning: Defining a sub-sampling for SAR calculation is not allowed!!! Skipped!" << endl;
 }
 
 void ProcessFieldsSAR::InitProcess()
 {
-	if (m_DumpType!=SAR_LOCAL_DUMP)
+	if ((m_DumpType!=SAR_LOCAL_DUMP) && (m_DumpType!=SAR_1G_DUMP) && (m_DumpType!=SAR_10G_DUMP) && (m_DumpType!=SAR_RAW_DATA))
 	{
 		Enabled=false;
 		cerr << "ProcessFieldsSAR::InitProcess(): Error, wrong dump type... this should not happen... skipping!" << endl;
@@ -60,6 +76,15 @@ void ProcessFieldsSAR::InitProcess()
 		cerr << "ProcessFieldsSAR::InitProcess(): Warning, interpolation type is not supported, resetting to CELL!" << endl;
 		SetDumpMode2Cell();
 	}
+
+	if ((m_DumpType==SAR_RAW_DATA) && (m_fileType!=HDF5_FILETYPE))
+	{
+		Enabled=false;
+		cerr << "ProcessFieldsSAR::InitProcess(): Error, wrong file type for dumping raw SAR data! skipping" << endl;
+		return;
+
+	}
+
 	ProcessFieldsFD::InitProcess();
 
 	if (Enabled==false) return;
@@ -68,7 +93,8 @@ void ProcessFieldsSAR::InitProcess()
 	for (size_t n = 0; n<m_FD_Samples.size(); ++n)
 	{
 		m_E_FD_Fields.push_back(Create_N_3DArray<std::complex<float> >(numLines));
-		m_J_FD_Fields.push_back(Create_N_3DArray<std::complex<float> >(numLines));
+		if (!m_UseCellKappa)
+			m_J_FD_Fields.push_back(Create_N_3DArray<std::complex<float> >(numLines));
 	}
 }
 
@@ -85,6 +111,9 @@ int ProcessFieldsSAR::Process()
 	double T;
 	FDTD_FLOAT**** field_td=NULL;
 
+	//save dump type
+	DumpType save_dump_type = m_DumpType;
+
 	// calc E-field
 	m_DumpType = E_FIELD_DUMP;
 	field_td = CalcField();
@@ -111,32 +140,35 @@ int ProcessFieldsSAR::Process()
 	Delete_N_3DArray<FDTD_FLOAT>(field_td,numLines);
 
 	// calc J-field
-	m_DumpType = J_FIELD_DUMP;
+	if (!m_UseCellKappa)
-	field_td = CalcField();
-	T = m_Eng_Interface->GetTime(m_dualTime);
-	for (size_t n = 0; n<m_FD_Samples.size(); ++n)
 	{
-		std::complex<float> exp_jwt_2_dt = std::exp( (std::complex<float>)(-2.0 * _I * M_PI * m_FD_Samples.at(n) * T) );
+		m_DumpType = J_FIELD_DUMP;
-		exp_jwt_2_dt *= 2; // *2 for single-sided spectrum
+		field_td = CalcField();
-		exp_jwt_2_dt *= Op->GetTimestep() * m_FD_Interval; // multiply with timestep-interval
+		T = m_Eng_Interface->GetTime(m_dualTime);
-		field_fd = m_J_FD_Fields.at(n);
+		for (size_t n = 0; n<m_FD_Samples.size(); ++n)
-		for (pos[0]=0; pos[0]<numLines[0]; ++pos[0])
 		{
-			for (pos[1]=0; pos[1]<numLines[1]; ++pos[1])
+			std::complex<float> exp_jwt_2_dt = std::exp( (std::complex<float>)(-2.0 * _I * M_PI * m_FD_Samples.at(n) * T) );
+			exp_jwt_2_dt *= 2; // *2 for single-sided spectrum
+			exp_jwt_2_dt *= Op->GetTimestep() * m_FD_Interval; // multiply with timestep-interval
+			field_fd = m_J_FD_Fields.at(n);
+			for (pos[0]=0; pos[0]<numLines[0]; ++pos[0])
 			{
-				for (pos[2]=0; pos[2]<numLines[2]; ++pos[2])
+				for (pos[1]=0; pos[1]<numLines[1]; ++pos[1])
 				{
-					field_fd[0][pos[0]][pos[1]][pos[2]] += field_td[0][pos[0]][pos[1]][pos[2]] * exp_jwt_2_dt;
+					for (pos[2]=0; pos[2]<numLines[2]; ++pos[2])
-					field_fd[1][pos[0]][pos[1]][pos[2]] += field_td[1][pos[0]][pos[1]][pos[2]] * exp_jwt_2_dt;
+					{
-					field_fd[2][pos[0]][pos[1]][pos[2]] += field_td[2][pos[0]][pos[1]][pos[2]] * exp_jwt_2_dt;
+						field_fd[0][pos[0]][pos[1]][pos[2]] += field_td[0][pos[0]][pos[1]][pos[2]] * exp_jwt_2_dt;
+						field_fd[1][pos[0]][pos[1]][pos[2]] += field_td[1][pos[0]][pos[1]][pos[2]] * exp_jwt_2_dt;
+						field_fd[2][pos[0]][pos[1]][pos[2]] += field_td[2][pos[0]][pos[1]][pos[2]] * exp_jwt_2_dt;
+					}
 				}
 			}
 		}
+		Delete_N_3DArray<FDTD_FLOAT>(field_td,numLines);
 	}
-	Delete_N_3DArray<FDTD_FLOAT>(field_td,numLines);
 
 	//reset dump type
-	m_DumpType = SAR_LOCAL_DUMP;
+	m_DumpType = save_dump_type;
 
 	++m_FD_SampleCount;
 	return GetNextInterval();
@@ -147,85 +179,152 @@ void ProcessFieldsSAR::DumpFDData()
 	if (Enabled==false) return;
 	unsigned int pos[3];
 	unsigned int orig_pos[3];
-	FDTD_FLOAT*** SAR = Create3DArray<float>(numLines);
+	float*** SAR = Create3DArray<float>(numLines);
-	std::complex<float>**** E_field_fd = NULL;
-	std::complex<float>**** J_field_fd = NULL;
 	double coord[3];
-	double density;
 	ContinuousStructure* CSX = Op->GetGeometryCSX();
 	CSProperties* prop = NULL;
 	CSPropMaterial* matProp = NULL;
 
 	double power;
-	double l_pow;
 
-	for (size_t n = 0; n<m_FD_Samples.size(); ++n)
+	float*** cell_volume = Create3DArray<float>(numLines);
+	float*** cell_density = Create3DArray<float>(numLines);
+	float*** cell_kappa = NULL;
+	if (m_UseCellKappa)
+		cell_kappa = Create3DArray<float>(numLines);
+
+	bool found_UnIsotropic=false;
+
+	// calculate volumes and masses for all cells
+	for (pos[0]=0; pos[0]<numLines[0]; ++pos[0])
 	{
-		E_field_fd = m_E_FD_Fields.at(n);
+		orig_pos[0] = posLines[0][pos[0]];
-		J_field_fd = m_J_FD_Fields.at(n);
+		coord[0] = Op->GetDiscLine(0,orig_pos[0],true);
-		power = 0;
+		for (pos[1]=0; pos[1]<numLines[1]; ++pos[1])
-		for (pos[0]=0; pos[0]<numLines[0]; ++pos[0])
 		{
-			orig_pos[0] = posLines[0][pos[0]];
+			orig_pos[1] = posLines[1][pos[1]];
-			coord[0] = Op->GetDiscLine(0,orig_pos[0],true);
+			coord[1] = Op->GetDiscLine(1,orig_pos[1],true);
-			for (pos[1]=0; pos[1]<numLines[1]; ++pos[1])
+			for (pos[2]=0; pos[2]<numLines[2]; ++pos[2])
 			{
-				orig_pos[1] = posLines[1][pos[1]];
+				orig_pos[2] = posLines[2][pos[2]];
-				coord[1] = Op->GetDiscLine(1,orig_pos[1],true);
+				coord[2] = Op->GetDiscLine(2,orig_pos[2],true);
-				for (pos[2]=0; pos[2]<numLines[2]; ++pos[2])
+
+				cell_volume[pos[0]][pos[1]][pos[2]] = Op->GetCellVolume(orig_pos);
+				cell_density[pos[0]][pos[1]][pos[2]] = 0.0;
+
+				prop = CSX->GetPropertyByCoordPriority(coord,CSProperties::MATERIAL);
+				if (prop!=0)
 				{
-					orig_pos[2] = posLines[2][pos[2]];
+					matProp = dynamic_cast<CSPropMaterial*>(prop);
-					coord[2] = Op->GetDiscLine(2,orig_pos[2],true);
+					if (matProp)
-					l_pow  = abs(E_field_fd[0][pos[0]][pos[1]][pos[2]]) * abs(J_field_fd[0][pos[0]][pos[1]][pos[2]]);
-					l_pow += abs(E_field_fd[1][pos[0]][pos[1]][pos[2]]) * abs(J_field_fd[1][pos[0]][pos[1]][pos[2]]);
-					l_pow += abs(E_field_fd[2][pos[0]][pos[1]][pos[2]]) * abs(J_field_fd[2][pos[0]][pos[1]][pos[2]]);
-
-					power += 0.5*l_pow*Op->GetCellVolume(orig_pos);
-
-					prop = CSX->GetPropertyByCoordPriority(coord,CSProperties::MATERIAL);
-					SAR[pos[0]][pos[1]][pos[2]] = 0.0;
-					density=0.0;
-					if (prop!=0)
 					{
-						matProp = dynamic_cast<CSPropMaterial*>(prop);
+						found_UnIsotropic |= !matProp->GetIsotropy();
-						density = matProp->GetDensityWeighted(coord);
+						cell_density[pos[0]][pos[1]][pos[2]] = matProp->GetDensityWeighted(coord);
-						if (density>0)
+						if (m_UseCellKappa)
-						{
+							cell_kappa[pos[0]][pos[1]][pos[2]] = matProp->GetKappaWeighted(0,coord);
-							SAR[pos[0]][pos[1]][pos[2]] = l_pow*0.5/density;
-						}
 					}
 				}
 			}
 		}
+	}
+	if (found_UnIsotropic)
+		cerr << "ProcessFieldsSAR::DumpFDData(): Warning, found unisotropic material in SAR calculation... this is unsupported!" << endl;
 
-		if (m_fileType==VTK_FILETYPE)
+	float* cellWidth[3];
+	for (int n=0;n<3;++n)
+	{
+		cellWidth[n]=new float[numLines[n]];
+		for (unsigned int i=0;i<numLines[n];++i)
+			cellWidth[n][i]=Op->GetDiscDelta(n,posLines[n][i])*Op->GetGridDelta();
+	}
+
+	if (m_DumpType == SAR_RAW_DATA)
+	{
+		if (m_fileType!=HDF5_FILETYPE)
 		{
-			stringstream ss;
+			cerr << "ProcessFieldsSAR::DumpFDData(): Error, wrong file type, this should not happen!!! skipped" << endl;
-			ss << m_filename << fixed << "_f=" << m_FD_Samples.at(n);
+			return;
-
-			m_Vtk_Dump_File->SetFilename(ss.str());
-			m_Vtk_Dump_File->ClearAllFields();
-			m_Vtk_Dump_File->AddScalarField(GetFieldNameByType(m_DumpType),SAR);
-			if (m_Vtk_Dump_File->Write()==false)
-				cerr << "ProcessFieldsSAR::Process: can't dump to file...! " << endl;
 		}
-		else if (m_fileType==HDF5_FILETYPE)
+
+		size_t datasize[]={numLines[0],numLines[1],numLines[2]};
+		for (size_t n = 0; n<m_FD_Samples.size(); ++n)
 		{
 			stringstream ss;
 			ss << "f" << n;
-			size_t datasize[]={numLines[0],numLines[1],numLines[2]};
+			if (m_HDF5_Dump_File->WriteVectorField(ss.str(), m_E_FD_Fields.at(n), datasize)==false)
-			if (m_HDF5_Dump_File->WriteScalarField(ss.str(), SAR, datasize)==false)
+				cerr << "ProcessFieldsSAR::DumpFDData: can't dump to file...! " << endl;
-				cerr << "ProcessFieldsSAR::Process: can't dump to file...! " << endl;
-			float freq[1] = {(float)m_FD_Samples.at(n)};
-			if (m_HDF5_Dump_File->WriteAtrribute("/FieldData/FD/"+ss.str(),"frequency",freq,1)==false)
-				cerr << "ProcessFieldsSAR::Process: can't dump to file...! " << endl;
-			float pow[1] = {(float)power};
-			if (m_HDF5_Dump_File->WriteAtrribute("/FieldData/FD/"+ss.str(),"power",pow,1)==false)
-				cerr << "ProcessFieldsSAR::Process: can't dump to file...! " << endl;
 		}
-		else
-			cerr << "ProcessFieldsSAR::Process: unknown File-Type" << endl;
-	}
 
+		m_HDF5_Dump_File->SetCurrentGroup("/CellData");
+		if (m_UseCellKappa==false)
+			cerr <<  "ProcessFieldsSAR::DumpFDData: Error, cell conductivity data not available, this should not happen... skipping! " << endl;
+		else if (m_HDF5_Dump_File->WriteScalarField("Conductivity", cell_kappa, datasize)==false)
+			cerr << "ProcessFieldsSAR::DumpFDData: can't dump to file...! " << endl;
+		if (m_HDF5_Dump_File->WriteScalarField("Density", cell_density, datasize)==false)
+			cerr << "ProcessFieldsSAR::DumpFDData: can't dump to file...! " << endl;
+		if (m_HDF5_Dump_File->WriteScalarField("Volume", cell_volume, datasize)==false)
+			cerr << "ProcessFieldsSAR::DumpFDData: can't dump to file...! " << endl;
+	}
+	else
+	{
+		SAR_Calculation SAR_Calc;
+		SAR_Calc.SetNumLines(numLines);
+		if (m_DumpType == SAR_LOCAL_DUMP)
+			SAR_Calc.SetAveragingMass(0);
+		else if (m_DumpType == SAR_1G_DUMP)
+			SAR_Calc.SetAveragingMass(1e-3);
+		else if (m_DumpType == SAR_10G_DUMP)
+			SAR_Calc.SetAveragingMass(10e-3);
+		else
+		{
+			cerr << "ProcessFieldsSAR::DumpFDData: unknown SAR dump type...!" << endl;
+		}
+		SAR_Calc.SetCellDensities(cell_density);
+		SAR_Calc.SetCellWidth(cellWidth);
+		SAR_Calc.SetCellVolumes(cell_volume);
+		SAR_Calc.SetCellCondictivity(cell_kappa);
+
+		for (size_t n = 0; n<m_FD_Samples.size(); ++n)
+		{
+			SAR_Calc.SetEField(m_E_FD_Fields.at(n));
+			if (!m_UseCellKappa)
+				SAR_Calc.SetJField(m_J_FD_Fields.at(n));
+			power = SAR_Calc.CalcSARPower();
+			SAR_Calc.CalcSAR(SAR);
+
+			if (m_fileType==VTK_FILETYPE)
+			{
+				stringstream ss;
+				ss << m_filename << fixed << "_f=" << m_FD_Samples.at(n);
+
+				m_Vtk_Dump_File->SetFilename(ss.str());
+				m_Vtk_Dump_File->ClearAllFields();
+				m_Vtk_Dump_File->AddScalarField(GetFieldNameByType(m_DumpType),SAR);
+				if (m_Vtk_Dump_File->Write()==false)
+					cerr << "ProcessFieldsSAR::DumpFDData: can't dump to file...! " << endl;
+			}
+			else if (m_fileType==HDF5_FILETYPE)
+			{
+				stringstream ss;
+				ss << "f" << n;
+				size_t datasize[]={numLines[0],numLines[1],numLines[2]};
+				if (m_HDF5_Dump_File->WriteScalarField(ss.str(), SAR, datasize)==false)
+					cerr << "ProcessFieldsSAR::DumpFDData: can't dump to file...! " << endl;
+				float freq[1] = {(float)m_FD_Samples.at(n)};
+				if (m_HDF5_Dump_File->WriteAtrribute("/FieldData/FD/"+ss.str(),"frequency",freq,1)==false)
+					cerr << "ProcessFieldsSAR::DumpFDData: can't dump to file...! " << endl;
+				float pow[1] = {(float)power};
+				if (m_HDF5_Dump_File->WriteAtrribute("/FieldData/FD/"+ss.str(),"power",pow,1)==false)
+					cerr << "ProcessFieldsSAR::DumpFDData: can't dump to file...! " << endl;
+			}
+			else
+				cerr << "ProcessFieldsSAR::DumpFDData: unknown File-Type" << endl;
+		}
+	}
+	for (int n=0;n<3;++n)
+		delete[] cellWidth[n];
+	Delete3DArray(cell_volume,numLines);
+	Delete3DArray(cell_density,numLines);
+	Delete3DArray(cell_kappa,numLines);
 	Delete3DArray(SAR,numLines);
 }

Common/processfields_sar.h

@@ -26,6 +26,10 @@ public:
 	ProcessFieldsSAR(Engine_Interface_Base* eng_if);
 	virtual ~ProcessFieldsSAR();
 
+	virtual void SetDumpType(DumpType type);
+
+	virtual bool NeedConductivity() const;
+
 	virtual string GetProcessingName() const {return "SAR dump";}
 
 	virtual void InitProcess();
@@ -36,9 +40,14 @@ public:
 
 	virtual void SetOptResolution(double optRes, int dir=-1);
 
+	//! Set to true for using the conductivity found at the center of a cell, instead of default E*J
+	virtual void SetUseCellConductivity(bool val) {m_UseCellKappa=val;}
+
 protected:
 	virtual void DumpFDData();
 
+	bool m_UseCellKappa;
+
 	//! frequency domain electric field storage
 	vector<std::complex<float>****> m_E_FD_Fields;
 	//! frequency domain current density storage

openEMS.pro

@@ -155,6 +155,7 @@ SOURCES += tools/global.cpp \
     tools/array_ops.cpp \
     tools/ErrorMsg.cpp \
     tools/AdrOp.cpp \
+    tools/sar_calculation.cpp \
     tools/vtk_file_writer.cpp \
     tools/hdf5_file_writer.cpp
 
@@ -224,6 +225,7 @@ HEADERS += tools/ErrorMsg.h \
     tools/global.h \
     tools/useful.h \
     tools/aligned_allocator.h \
+    tools/sar_calculation.h \
     tools/vtk_file_writer.h \
     tools/hdf5_file_writer.h
 

openems.cpp

@@ -432,7 +432,7 @@ bool openEMS::SetupProcessing()
 						ProcField = new ProcessFieldsTD(NewEngineInterface());
 					else if ((db->GetDumpType()>=10) && (db->GetDumpType()<=13))
 						ProcField = new ProcessFieldsFD(NewEngineInterface());
-					else if (db->GetDumpType()==20)
+					else if ( ((db->GetDumpType()>=20) && (db->GetDumpType()<=22)) || (db->GetDumpType()==29) )
 						ProcField = new ProcessFieldsSAR(NewEngineInterface());
 					else
 						cerr << "openEMS::SetupFDTD: unknown dump box type... skipping!" << endl;
@@ -455,12 +455,16 @@ bool openEMS::SetupProcessing()
 							ProcField->SetDumpType((ProcessFields::DumpType)db->GetDumpType());
 
 						if (db->GetDumpType()==20)
-						{
 							ProcField->SetDumpType(ProcessFields::SAR_LOCAL_DUMP);
-						}
+						if (db->GetDumpType()==21)
+							ProcField->SetDumpType(ProcessFields::SAR_1G_DUMP);
+						if (db->GetDumpType()==22)
+							ProcField->SetDumpType(ProcessFields::SAR_10G_DUMP);
+						if (db->GetDumpType()==29)
+							ProcField->SetDumpType(ProcessFields::SAR_RAW_DATA);
 
 						//SetupMaterialStorages() has previewed storage needs... refresh here to prevent cleanup!!!
-						if ( ((db->GetDumpType()==2) || (db->GetDumpType()==12) || (db->GetDumpType()==20)) && Enable_Dumps )
+						if ( ProcField->NeedConductivity() && Enable_Dumps )
 							FDTD_Op->SetMaterialStoreFlags(1,true);
 
 						ProcField->SetDumpMode((Engine_Interface_Base::InterpolationType)db->GetDumpMode());

tools/sar_calculation.cpp

@@ -0,0 +1,632 @@
+/*
+*	Copyright (C) 2012 Thorsten Liebig (Thorsten.Liebig@gmx.de)
+*
+*	This program is free software: you can redistribute it and/or modify
+*	it under the terms of the GNU General Public License as published by
+*	the Free Software Foundation, either version 3 of the License, or
+*	(at your option) any later version.
+*
+*	This program is distributed in the hope that it will be useful,
+*	but WITHOUT ANY WARRANTY; without even the implied warranty of
+*	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*	GNU General Public License for more details.
+*
+*	You should have received a copy of the GNU General Public License
+*	along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+
+#include "sar_calculation.h"
+#include "cfloat"
+#include "array_ops.h"
+
+SAR_Calculation::SAR_Calculation()
+{
+	m_Vx_Used = NULL;
+	m_Vx_Valid = NULL;
+	m_DebugLevel = 1;
+	SetAveragingMethod(SIMPLE);
+	Reset();
+}
+
+void SAR_Calculation::Reset()
+{
+	Delete3DArray(m_Vx_Used,m_numLines);
+	m_Vx_Used = NULL;
+	Delete3DArray(m_Vx_Valid,m_numLines);
+	m_Vx_Valid = NULL;
+
+	m_avg_mass = 0;
+	m_numLines[0]=m_numLines[1]=m_numLines[2]=0;
+	m_cellWidth[0]=m_cellWidth[1]=m_cellWidth[2]=NULL;
+
+	m_cell_volume = NULL;
+	m_cell_density = NULL;
+	m_cell_conductivity = NULL;
+	m_E_field = NULL;
+	m_J_field = NULL;
+
+	Delete3DArray(m_Vx_Used,m_numLines);
+	m_Vx_Used = NULL;
+	Delete3DArray(m_Vx_Valid,m_numLines);
+	m_Vx_Valid = NULL;
+}
+
+void SAR_Calculation::SetAveragingMethod(SARAveragingMethod method)
+{
+	if (method==IEEE_C95_3)
+	{
+		m_massTolerance = 0.0001;
+		m_maxMassIterations = 100;
+		m_maxBGRatio = 0.1;
+		m_markPartialAsUsed = false;
+		m_UnusedRelativeVolLimit = 1.05;
+		m_IgnoreFaceValid = false;
+		return;
+	}
+	else if (method==IEEE_62704)
+	{
+		m_massTolerance = 0.05;
+		m_maxMassIterations = 100;
+		m_maxBGRatio = 1;
+		m_markPartialAsUsed = true;
+		m_UnusedRelativeVolLimit = 1;
+		m_IgnoreFaceValid = false;
+		return;
+	}
+	else if (method==SIMPLE)
+	{
+		m_massTolerance = 0.05;
+		m_maxMassIterations = 100;
+		m_maxBGRatio = 1;
+		m_markPartialAsUsed = true;
+		m_UnusedRelativeVolLimit = 1;
+		m_IgnoreFaceValid = true;
+		return;
+	}
+
+	// unknown method, fallback to simple
+	SetAveragingMethod(SIMPLE);
+}
+
+void SAR_Calculation::SetNumLines(unsigned int numLines[3])
+{
+	Delete3DArray(m_Vx_Used,m_numLines);
+	m_Vx_Used = NULL;
+	Delete3DArray(m_Vx_Valid,m_numLines);
+	m_Vx_Valid = NULL;
+
+	for (int n=0;n<3;++n)
+		m_numLines[n]=numLines[n];
+}
+
+void SAR_Calculation::SetCellWidth(float* cellWidth[3])
+{
+	for (int n=0;n<3;++n)
+		m_cellWidth[n]=cellWidth[n];
+}
+
+float*** SAR_Calculation::CalcSAR(float*** SAR)
+{
+	if (CheckValid()==false)
+	{
+		cerr << "SAR_Calculation::CalcSAR: SAR calculation is invalid due to missing values... Abort..." << endl;
+		return NULL;
+	}
+	if (m_avg_mass<=0)
+		return CalcLocalSAR(SAR);
+	return CalcAveragedSAR(SAR);
+}
+
+bool SAR_Calculation::CheckValid()
+{
+	for (int n=0;n<3;++n)
+		if (m_cellWidth[n]==NULL)
+			return false;
+	if (m_E_field==NULL)
+		return false;
+	if ((m_J_field==NULL) && (m_cell_conductivity==NULL))
+		return false;
+	if (m_cell_density==NULL)
+		return false;
+	if (m_avg_mass<0)
+		return false;
+	return true;
+}
+
+double SAR_Calculation::CalcSARPower()
+{
+	if (CheckValid()==false)
+	{
+		cerr << "SAR_Calculation::CalcSARPower: SAR calculation is invalid due to missing values... Abort..." << endl;
+		return 0;
+	}
+	double power=0;
+	unsigned int pos[3];
+	for (pos[0]=0; pos[0]<m_numLines[0]; ++pos[0])
+	{
+		for (pos[1]=0; pos[1]<m_numLines[1]; ++pos[1])
+		{
+			for (pos[2]=0; pos[2]<m_numLines[2]; ++pos[2])
+			{
+				power += CalcLocalPowerDensity(pos)*CellVolume(pos);
+			}
+		}
+	}
+	return power;
+}
+
+double SAR_Calculation::CalcLocalPowerDensity(unsigned int pos[3])
+{
+	double l_pow=0;
+	if (m_cell_conductivity==NULL)
+	{
+		l_pow  = abs(m_E_field[0][pos[0]][pos[1]][pos[2]]) * abs(m_J_field[0][pos[0]][pos[1]][pos[2]]);
+		l_pow += abs(m_E_field[1][pos[0]][pos[1]][pos[2]]) * abs(m_J_field[1][pos[0]][pos[1]][pos[2]]);
+		l_pow += abs(m_E_field[2][pos[0]][pos[1]][pos[2]]) * abs(m_J_field[2][pos[0]][pos[1]][pos[2]]);
+	}
+	else
+	{
+		l_pow  = m_cell_conductivity[pos[0]][pos[1]][pos[2]]*abs(m_E_field[0][pos[0]][pos[1]][pos[2]]) * abs(m_E_field[0][pos[0]][pos[1]][pos[2]]);
+		l_pow += m_cell_conductivity[pos[0]][pos[1]][pos[2]]*abs(m_E_field[1][pos[0]][pos[1]][pos[2]]) * abs(m_E_field[1][pos[0]][pos[1]][pos[2]]);
+		l_pow += m_cell_conductivity[pos[0]][pos[1]][pos[2]]*abs(m_E_field[2][pos[0]][pos[1]][pos[2]]) * abs(m_E_field[2][pos[0]][pos[1]][pos[2]]);
+	}
+	l_pow*=0.5;
+	return l_pow;
+}
+
+float*** SAR_Calculation::CalcLocalSAR(float*** SAR)
+{
+	unsigned int pos[3];
+	m_Valid=0;
+	m_Used=0;
+	m_Unused=0;
+	m_AirVoxel=0;
+	for (pos[0]=0; pos[0]<m_numLines[0]; ++pos[0])
+	{
+		for (pos[1]=0; pos[1]<m_numLines[1]; ++pos[1])
+		{
+			for (pos[2]=0; pos[2]<m_numLines[2]; ++pos[2])
+			{
+				if (m_cell_density[pos[0]][pos[1]][pos[2]]>0)
+				{
+					++m_Valid;
+					SAR[pos[0]][pos[1]][pos[2]]	= CalcLocalPowerDensity(pos)/m_cell_density[pos[0]][pos[1]][pos[2]];
+				}
+				else
+				{
+					++m_AirVoxel;
+					SAR[pos[0]][pos[1]][pos[2]] = 0;
+				}
+			}
+		}
+	}
+	return SAR;
+}
+
+int SAR_Calculation::FindFittingCubicalMass(unsigned int pos[3], float box_size, unsigned int start[3], unsigned int stop[3],
+					float partial_start[3], float partial_stop[3], float &mass, float &volume, float &bg_ratio, int disabledFace, bool ignoreFaceValid)
+{
+	unsigned int mass_iterations = 0;
+	float old_mass=0;
+	float old_box_size=0;
+	bool face_valid;
+	bool mass_valid;
+	bool voxel_valid;
+
+	//iterate over cubical sizes to find fitting volume to mass
+	while (mass_iterations<m_maxMassIterations)
+	{
+		// calculate cubical mass
+		face_valid = GetCubicalMass(pos, box_size/2, start, stop, partial_start, partial_stop, mass, volume, bg_ratio, disabledFace);
+
+		// check if found mass is valid
+		mass_valid = abs(mass-m_avg_mass)<=m_massTolerance*m_avg_mass;
+		voxel_valid = mass_valid && (face_valid==true) && (bg_ratio<m_maxBGRatio);
+
+		if ((face_valid==false) && (mass<m_avg_mass*(1.0-m_massTolerance)) && (ignoreFaceValid==false))
+		{
+			// this is an invalid cube with a too small total mass --> increasing the box will not yield a valid cube
+			return 1;
+		}
+		else if ((face_valid==false  || bg_ratio>=m_maxBGRatio) && (mass_valid))
+		{
+			// this is an invalid cube with a valid total mass
+			if (ignoreFaceValid)
+				return 0;
+			return 2;
+		}
+		if (voxel_valid)
+		{
+			// valid cube found
+			return 0;
+		}
+
+		// if no valid or finally invalid cube is found, calculate an alternaive cube size
+		if (mass_iterations==0)
+		{
+			// on first interation, try a relative resize
+			old_box_size=box_size;
+			box_size*=pow(m_avg_mass/mass,1.0/3.0);
+		}
+		else
+		{
+			// on later interations, try a newton approach
+			float new_box_size = box_size - (mass-m_avg_mass)/(mass-old_mass)*(box_size-old_box_size);
+			old_box_size = box_size;
+			box_size = new_box_size;
+		}
+		old_mass=mass;
+
+		++mass_iterations;
+	}
+
+	// m_maxMassIterations iterations are exhausted...
+	return -1;
+}
+
+bool SAR_Calculation::GetCubicalMass(unsigned int pos[3], float box_size, unsigned int start[3], unsigned int stop[3],
+									 float partial_start[3], float partial_stop[3], float &mass, float &volume, float &bg_ratio, int disabledFace)
+{
+	if ((box_size<=0) || isnan(box_size) || isinf(box_size))
+	{
+		cerr << "SAR_Calculation::GetCubicalMass: critical error: invalid averaging box size!! EXIT" << endl;
+		exit(-1);
+	}
+	bool face_valid=true;
+	for (int n=0;n<3;++n)
+	{
+		// check start position
+		start[n]=pos[n];
+		partial_start[n]=1;
+		float dist=m_cellWidth[n][pos[n]]/2;
+		if (disabledFace==2*n)
+			dist=box_size;
+		else
+		{
+			while (dist<box_size)
+			{
+				// box is outside of domain
+				if (start[n]==0)
+				{
+					partial_start[n]=-1;
+					break;
+				}
+				--start[n];
+				dist+=m_cellWidth[n][start[n]];
+			}
+
+			if (partial_start[n]!=-1)
+			{   // box ends inside stop[n] voxel
+				partial_start[n]=1-(dist-box_size)/m_cellWidth[n][start[n]];
+			}
+			if ((partial_start[n]!=-1) && (pos[n]==start[n]))
+				partial_start[n]=2*box_size/m_cellWidth[n][start[n]];
+		}
+
+		// check stop position
+		stop[n]=pos[n];
+		partial_stop[n]=1;
+		dist=m_cellWidth[n][pos[n]]/2;
+		if (disabledFace==2*n+1)
+			dist=box_size;
+		else
+		{
+			while (dist<box_size)
+			{
+				// box is outside of domain
+				if (stop[n]==m_numLines[n]-1)
+				{
+					partial_stop[n]=-1;
+					break;
+				}
+				++stop[n];
+				dist+=m_cellWidth[n][stop[n]];
+			}
+
+			if (partial_stop[n]!=-1)
+			{   // box ends inside stop[n] voxel
+				partial_stop[n]=1-(dist-box_size)/m_cellWidth[n][stop[n]];
+			}
+			if ((partial_stop[n]!=-1) && (pos[n]==stop[n]))
+				partial_stop[n]=2*box_size/m_cellWidth[n][stop[n]];
+		}
+	}
+
+	for (int n=0;n<3;++n)
+	{
+		if (partial_start[n]==-1)
+			face_valid=false;
+		if (partial_stop[n]==-1)
+			face_valid=false;
+	}
+
+	mass = 0;
+	volume = 0;
+	float bg_volume=0;
+	float weight[3];
+	unsigned int f_pos[3];
+	bool face_intersect[6] = {false,false,false,false,false,false};
+	for (f_pos[0]=start[0];f_pos[0]<=stop[0];++f_pos[0])
+	{
+		weight[0]=1;
+		if (f_pos[0]==start[0])
+			weight[0]*=abs(partial_start[0]);
+		if (f_pos[0]==stop[0])
+			weight[0]*=abs(partial_stop[0]);
+
+		for (f_pos[1]=start[1];f_pos[1]<=stop[1];++f_pos[1])
+		{
+			weight[1]=1;
+			if (f_pos[1]==start[1])
+				weight[1]*=abs(partial_start[1]);
+			if (f_pos[1]==stop[1])
+				weight[1]*=abs(partial_stop[1]);
+
+			for (f_pos[2]=start[2];f_pos[2]<=stop[2];++f_pos[2])
+			{
+				weight[2]=1;
+				if (f_pos[2]==start[2])
+					weight[2]*=abs(partial_start[2]);
+				if (f_pos[2]==stop[2])
+					weight[2]*=abs(partial_stop[2]);
+
+				mass += CellMass(f_pos)*weight[0]*weight[1]*weight[2];
+				volume += CellVolume(f_pos)*weight[0]*weight[1]*weight[2];
+
+				if (m_cell_density[f_pos[0]][f_pos[1]][f_pos[2]]==0)
+					bg_volume += CellVolume(f_pos)*weight[0]*weight[1]*weight[2];
+				else
+				{
+					for (int n=0;n<3;++n)
+					{
+						if (start[n]==f_pos[n])
+							face_intersect[2*n]=true;
+						if (stop[n]==f_pos[n])
+							face_intersect[2*n+1]=true;
+					}
+				}
+			}
+		}
+	}
+	//check if all bounds have intersected a material boundary
+	for (int n=0;n<6;++n)
+		face_valid *= face_intersect[n];
+
+	bg_ratio = bg_volume/volume;
+
+	return face_valid;
+}
+
+float SAR_Calculation::CalcCubicalSAR(float*** SAR, unsigned int pos[3], unsigned int start[3], unsigned int stop[3], float partial_start[3], float partial_stop[3], bool assignUsed)
+{
+	float power_mass=0;
+	float mass=0;
+	float weight[3];
+	unsigned int f_pos[3];
+	for (f_pos[0]=start[0];f_pos[0]<=stop[0];++f_pos[0])
+	{
+		weight[0]=1;
+		if (f_pos[0]==start[0])
+			weight[0]*=abs(partial_start[0]);
+		if (f_pos[0]==stop[0])
+			weight[0]*=abs(partial_stop[0]);
+
+		for (f_pos[1]=start[1];f_pos[1]<=stop[1];++f_pos[1])
+		{
+			weight[1]=1;
+			if (f_pos[1]==start[1])
+				weight[1]*=abs(partial_start[1]);
+			if (f_pos[1]==stop[1])
+				weight[1]*=abs(partial_stop[1]);
+
+			for (f_pos[2]=start[2];f_pos[2]<=stop[2];++f_pos[2])
+			{
+				weight[2]=1;
+				if (f_pos[2]==start[2])
+					weight[2]*=abs(partial_start[2]);
+				if (f_pos[2]==stop[2])
+					weight[2]*=abs(partial_stop[2]);
+
+				if (m_cell_density[f_pos[0]][f_pos[1]][f_pos[2]]>=0)
+				{
+					mass += CellMass(f_pos)*weight[0]*weight[1]*weight[2];
+					power_mass += CalcLocalPowerDensity(f_pos) * CellVolume(f_pos)*weight[0]*weight[1]*weight[2];
+				}
+			}
+		}
+	}
+	float vx_SAR = power_mass/mass;
+	if (SAR==NULL)
+		return vx_SAR;
+
+	SAR[pos[0]][pos[1]][pos[2]]=vx_SAR;
+
+	if (assignUsed==false)
+		return vx_SAR;
+
+	// assign SAR to all used voxel
+	bool is_partial[3];
+	for (f_pos[0]=start[0];f_pos[0]<=stop[0];++f_pos[0])
+	{
+		if ( ((f_pos[0]==start[0]) && (partial_start[0]!=1)) || ((f_pos[0]==stop[0]) && (partial_stop[0]!=1)) )
+			is_partial[0]=true;
+		else
+			is_partial[0]=false;
+
+		for (f_pos[1]=start[1];f_pos[1]<=stop[1];++f_pos[1])
+		{
+			if ( ((f_pos[1]==start[1]) && (partial_start[1]!=1)) || ((f_pos[1]==stop[1]) && (partial_stop[1]!=1)) )
+				is_partial[1]=true;
+			else
+				is_partial[1]=false;
+
+			for (f_pos[2]=start[2];f_pos[2]<=stop[2];++f_pos[2])
+			{
+				if ( ((f_pos[2]==start[2]) && (partial_start[2]!=1)) || ((f_pos[2]==stop[2]) && (partial_stop[2]!=1)) )
+					is_partial[2]=true;
+				else
+					is_partial[2]=false;
+
+				if ( (!is_partial[0] && !is_partial[1] && !is_partial[2]) || m_markPartialAsUsed)
+				{
+					if (!m_Vx_Valid[f_pos[0]][f_pos[1]][f_pos[2]] && (m_cell_density[f_pos[0]][f_pos[1]][f_pos[2]]>0))
+					{
+						m_Vx_Used[f_pos[0]][f_pos[1]][f_pos[2]]=true;
+						SAR[f_pos[0]][f_pos[1]][f_pos[2]]=max(SAR[f_pos[0]][f_pos[1]][f_pos[2]], vx_SAR);
+					}
+				}
+			}
+		}
+	}
+	return vx_SAR;
+}
+
+
+float*** SAR_Calculation::CalcAveragedSAR(float*** SAR)
+{
+	unsigned int pos[3];
+	m_Vx_Used = Create3DArray<bool>(m_numLines);
+	m_Vx_Valid = Create3DArray<bool>(m_numLines);
+
+	float voxel_volume;
+	float total_mass;
+	unsigned int start[3];
+	unsigned int stop[3];
+	float partial_start[3];
+	float partial_stop[3];
+	float bg_ratio;
+	int EC=0;
+
+	// debug counter
+	unsigned int cnt_case1=0;
+	unsigned int cnt_case2=0;
+	unsigned int cnt_NoConvergence=0;
+
+	m_Valid=0;
+	m_Used=0;
+	m_Unused=0;
+	m_AirVoxel=0;
+
+	for (pos[0]=0; pos[0]<m_numLines[0]; ++pos[0])
+	{
+		for (pos[1]=0; pos[1]<m_numLines[1]; ++pos[1])
+		{
+			for (pos[2]=0; pos[2]<m_numLines[2]; ++pos[2])
+			{
+				if (m_cell_density[pos[0]][pos[1]][pos[2]]==0)
+				{
+					SAR[pos[0]][pos[1]][pos[2]] = 0;
+					++m_AirVoxel;
+					continue;
+				}
+
+				// guess an initial box size and find a fitting cube
+				EC = FindFittingCubicalMass(pos, pow(m_avg_mass/m_cell_density[pos[0]][pos[1]][pos[2]],1.0/3.0)/2, start, stop,
+											partial_start, partial_stop, total_mass, voxel_volume, bg_ratio, -1, m_IgnoreFaceValid);
+
+				if (EC==0)
+				{
+					m_Vx_Valid[pos[0]][pos[1]][pos[2]] = true;
+					m_Vx_Used[pos[0]][pos[1]][pos[2]] = true;
+					++m_Valid;
+					CalcCubicalSAR(SAR, pos, start, stop, partial_start, partial_stop, true);
+				}
+				else if (EC==1)
+					++cnt_case1;
+				else if (EC==2)
+					++cnt_case2;
+				else if (EC==-1)
+					++cnt_NoConvergence;
+				else
+					cerr << "other EC" << EC << endl;
+			}
+		}
+	}
+	if (cnt_NoConvergence>0)
+	{
+		cerr << "SAR_Calculation::CalcAveragedSAR: Warning, for some voxel a valid averaging cube could not be found (no convergence)... " << endl;
+	}
+	if (m_DebugLevel>0)
+	{
+		cerr << "Number of invalid cubes (case 1): " << cnt_case1 << endl;
+		cerr << "Number of invalid cubes (case 2): " << cnt_case2 << endl;
+		cerr << "Number of invalid cubes (falied to converge): " << cnt_NoConvergence << endl;
+	}
+
+	// count all used and unused etc. + special handling of unused voxels!!
+	m_Used=0;
+	m_Unused=0;
+	for (pos[0]=0;pos[0]<m_numLines[0];++pos[0])
+	{
+		for (pos[1]=0;pos[1]<m_numLines[1];++pos[1])
+		{
+			for (pos[2]=0;pos[2]<m_numLines[2];++pos[2])
+			{
+				if (!m_Vx_Valid[pos[0]][pos[1]][pos[2]] && m_Vx_Used[pos[0]][pos[1]][pos[2]])
+					++m_Used;
+				if ((m_cell_density[pos[0]][pos[1]][pos[2]]>0) && !m_Vx_Valid[pos[0]][pos[1]][pos[2]] && !m_Vx_Used[pos[0]][pos[1]][pos[2]])
+				{
+					++m_Unused;
+
+					SAR[pos[0]][pos[1]][pos[2]] = 0;
+					float unused_volumes[6];
+					float unused_SAR[6];
+
+					float min_Vol=FLT_MAX;
+
+					// special handling of unused voxels:
+					for (int n=0;n<6;++n)
+					{
+						EC = FindFittingCubicalMass(pos, pow(m_avg_mass/m_cell_density[pos[0]][pos[1]][pos[2]],1.0/3.0)/2, start, stop,
+													partial_start, partial_stop, total_mass, unused_volumes[n], bg_ratio, n, true);
+						if ((EC!=0) && (EC!=2))
+						{
+							// this should not happen
+							cerr << "SAR_Calculation::CalcAveragedSAR: Error handling unused voxels, can't find fitting cubical averaging volume' " << endl;
+							unused_SAR[n]=0;
+						}
+						else
+						{
+							unused_SAR[n]=CalcCubicalSAR(NULL, pos, start, stop, partial_start, partial_stop, false);
+							min_Vol = min(min_Vol,unused_volumes[n]);
+						}
+					}
+					for (int n=0;n<6;++n)
+					{
+						if (unused_volumes[n]<=m_UnusedRelativeVolLimit*min_Vol)
+							SAR[pos[0]][pos[1]][pos[2]] = max(SAR[pos[0]][pos[1]][pos[2]],unused_SAR[n]);
+					}
+				}
+			}
+		}
+	}
+
+	if (m_Valid+m_Used+m_Unused+m_AirVoxel!=m_numLines[0]*m_numLines[1]*m_numLines[2])
+	{
+		cerr << "SAR_Calculation::CalcAveragedSAR: critical error, mismatch in voxel status count... EXIT" << endl;
+		exit(1);
+	}
+
+	if (m_DebugLevel>0)
+		cerr << "SAR_Calculation::CalcAveragedSAR: Stats: Valid=" << m_Valid << " Used=" << m_Used << " Unused=" << m_Unused << " Air-Voxel=" << m_AirVoxel << endl;
+
+	return SAR;
+}
+
+float SAR_Calculation::CellVolume(unsigned int pos[3])
+{
+	if (m_cell_volume)
+		return m_cell_volume[pos[0]][pos[1]][pos[2]];
+
+	float volume=1;
+	for (int n=0;n<3;++n)
+		volume*=m_cellWidth[n][pos[n]];
+	return volume;
+}
+
+float SAR_Calculation::CellMass(unsigned int pos[3])
+{
+	return m_cell_density[pos[0]][pos[1]][pos[2]]*CellVolume(pos);
+}
+

tools/sar_calculation.h

@@ -0,0 +1,119 @@
+/*
+*	Copyright (C) 2012 Thorsten Liebig (Thorsten.Liebig@gmx.de)
+*
+*	This program is free software: you can redistribute it and/or modify
+*	it under the terms of the GNU General Public License as published by
+*	the Free Software Foundation, either version 3 of the License, or
+*	(at your option) any later version.
+*
+*	This program is distributed in the hope that it will be useful,
+*	but WITHOUT ANY WARRANTY; without even the implied warranty of
+*	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+*	GNU General Public License for more details.
+*
+*	You should have received a copy of the GNU General Public License
+*	along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef SAR_CALCULATION_H
+#define SAR_CALCULATION_H
+
+#include <complex>
+
+class SAR_Calculation
+{
+public:
+	SAR_Calculation();
+
+	enum SARAveragingMethod { IEEE_C95_3, IEEE_62704, SIMPLE};
+
+	//! Reset and initialize all values (will keep all SAR settings)
+	void Reset();
+
+	//! Set the debug level
+	void SetDebugLevel(int level) {m_DebugLevel=level;}
+
+	//! Set the used averaging method
+	void SetAveragingMethod(SARAveragingMethod method);
+
+	//! Set number of lines in all direcitions. (mandatory information)
+	void SetNumLines(unsigned int numLines[3]);
+	//! Set cell width in all direcitions. (mandatory information for averaging)
+	void SetCellWidth(float* cellWidth[3]);
+
+	//! Set the averaging mash. (mandatory information for averaging)
+	void SetAveragingMass(float mass) {m_avg_mass=mass;}
+
+	//! Set the cell volumes (optional for speedup)
+	void SetCellVolumes(float*** cell_volume) {m_cell_volume=cell_volume;}
+
+	//! Set the cell densities (mandatory information)
+	void SetCellDensities(float*** cell_density) {m_cell_density=cell_density;}
+
+	//! Set the cell conductivities (mandatory if no current density field is given)
+	void SetCellCondictivity(float*** cell_conductivity) {m_cell_conductivity=cell_conductivity;}
+
+	//! Set the electric field (mandatory information)
+	void SetEField(std::complex<float>**** field) {m_E_field=field;}
+	//! Set the current density field (mandatory if no conductivity distribution is given)
+	void SetJField(std::complex<float>**** field) {m_J_field=field;}
+
+	//! Calculate the SAR, requires a preallocated 3D array
+	float*** CalcSAR(float*** SAR);
+
+	//! Calculate the total power dumped
+	double CalcSARPower();
+
+protected:
+	unsigned int m_numLines[3];
+	float* m_cellWidth[3];
+
+	float m_avg_mass;
+	float*** m_cell_volume;
+	float*** m_cell_density;
+	float*** m_cell_conductivity;
+	std::complex<float>**** m_E_field;
+	std::complex<float>**** m_J_field;
+
+	bool*** m_Vx_Used;
+	bool*** m_Vx_Valid;
+
+	unsigned int m_Valid;
+	unsigned int m_Used;
+	unsigned int m_Unused;
+	unsigned int m_AirVoxel;
+
+	int m_DebugLevel;
+
+	/*********** SAR calculation parameter and settings ***********/
+	float m_massTolerance;
+	unsigned int m_maxMassIterations;
+	float m_maxBGRatio;
+	bool m_markPartialAsUsed;
+	float m_UnusedRelativeVolLimit;
+	bool m_IgnoreFaceValid;
+
+	/*********** SAR calculations methods ********/
+	double CalcLocalPowerDensity(unsigned int pos[3]);
+
+	//! Calculate the local SAR
+	float*** CalcLocalSAR(float*** SAR);
+
+	/****** start SAR averaging and all necessary methods ********/
+	//! Calculate the averaged SAR
+	float*** CalcAveragedSAR(float*** SAR);
+
+	int FindFittingCubicalMass(unsigned int pos[3], float box_size, unsigned int start[3], unsigned int stop[3],
+						float partial_start[3], float partial_stop[3], float &mass, float &volume, float &bg_ratio, int disabledFace=-1, bool ignoreFaceValid=false);
+	bool GetCubicalMass(unsigned int pos[3], float box_size, unsigned int start[3], unsigned int stop[3],
+						float partial_start[3], float partial_stop[3], float &mass, float &volume, float &bg_ratio, int disabledFace=-1);
+
+	float CalcCubicalSAR(float*** SAR, unsigned int pos[3], unsigned int start[3], unsigned int stop[3], float partial_start[3], float partial_stop[3], bool assignUsed=false);
+	/****** end SAR averaging and all necessary methods ********/
+
+	bool CheckValid();
+	float CellVolume(unsigned int pos[3]);
+	float CellMass(unsigned int pos[3]);
+};
+
+#endif // SAR_CALCULATION_H