#include "mex_imao_calculate.hh"
#include <math.h>

namespace AMDA {
namespace Parameters {
	
	Tab2DData<float> mex_imao_calculate(const Tab2DData<float>& Hsp, const std::vector<double>& Energy, const short& Pacc)
	{
		Tab2DData<float> lImaH(Hsp.getDim1Size(),96);
		float samplingTime = 0.1209;	
		float Geff;

		for(int i = 0; i < Hsp.getDim1Size(); ++i) { // Different elevations number                       			
			std::vector<float> imaHVector = Hsp[i];
			
			for (unsigned int iEn = 0; iEn < 96; iEn++) {// 96 Energy Steps
				
				(lImaH[i])[iEn] = 0.0;
				
				if ((float)Energy[iEn] > 0.0 && imaHVector[iEn] > 0.0)
				{
					Geff = gfl(Pacc, (float)Energy[iEn]);
					
					float coeff = 0.0;
					if (Geff > 0.0) 
						coeff = 1.0/(Geff * (float)Energy[iEn] * samplingTime);

					(lImaH[i])[iEn] = (imaHVector[iEn]*coeff);
				}	 
			}
		} 
		
		return lImaH; 
	}
	
	/*
	 *   Calculate GeometricFactor for MEX / IMA / HEAVY IONS
	 */
	float gfl(const short Pacc,  const float Energy)
	{
		float Gl = -1.0;
		
		switch ( Pacc ) {
			case 0:
				Gl = exp(-0.378481*log(Energy) - 9.05248);
				
			break;
			case 4:
				if (Energy < 300.0)
					Gl = 7.0e-5;
				else
					Gl = exp(-0.69*log(Energy) - 5.63);
				 
			break;
			case 7:
				if (Energy < 600.0)
					Gl = 1.0e-4;
				else 
					Gl = exp(-0.611844*log(Energy) -5.31592);
				
			break;
			default:
				Gl = -1.0;
			break;
		}
		
		return Gl;
	}
}
}