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

namespace AMDA {
namespace Parameters {
	
	Tab2DData<float> mex_imah_calculate(const Tab2DData<float>& Hsp, const Tab2DData<float>& Gsp, const std::vector<double>& Energy, const short& Pacc, bool isFlux)
	{
		Tab2DData<float> lImaH(Hsp.getDim1Size(),96);
		float samplingTime = 0.1209;	
		float GeffH, GeffG;

		for(unsigned int i = 0; i < Hsp.getDim1Size(); ++i) { // Different elevations number                       			
			std::vector<float> imaHVector = Hsp[i];
			std::vector<float> imaGVector = Gsp[i];
			for (unsigned int iEn = 0; iEn < 96; iEn++) {// 96 Energy Steps
				
				(lImaH[i])[iEn] = 0.0;
				
				if ((float)Energy[iEn] > 0 &&  (imaHVector[iEn] + imaGVector[iEn]) > 0)
				{
					GeffH = gfl(Pacc, (float)Energy[iEn], false);
					GeffG = gfl(Pacc, (float)Energy[iEn], true);
					
					float coeffH = 0.0;
					if (GeffH > 0) {
						if (isFlux) // diff flux
							coeffH = 1.0/(GeffH * (float)Energy[iEn] * samplingTime);
						else // counts
							coeffH = 1.e-5/GeffH; 
					}
					
					float coeffG = 0.0;
					if (GeffG > 0) {
						if (isFlux) // diff flux
							coeffG = 1.0/(GeffG * (float)Energy[iEn] * samplingTime);
						else // counts
							coeffG = 1.e-5/GeffG; 
					}
					
					(lImaH[i])[iEn] = (imaHVector[iEn]*coeffH + imaGVector[iEn]*coeffG);
				}	 
			}
		} 
		
		return lImaH; 
	}
	
	/*
	 *   Calculate GeometricFactor for MEX / IMA / H+ special case (ghost!)
	 */
	float gfl(const short Pacc,  const float Energy, bool isGhost)
	{
		float Gl = -1.0;
		
		switch ( Pacc ) {
			case 0:
				Gl = exp(-0.378481*log(Energy) - 9.05248);
			break;
			case 4:
				if (isGhost) {
					if ((Energy < 800.0) && (Energy > 300.0))
						Gl = exp((14.0-13.1224)/(6.68461-6.21461)*(log(Energy) - 6.21461)-14.0);
					else if ((Energy >= 800.0) && (Energy < 1500.0))
						Gl = 2.0e-6;
					else if ((Energy >= 1500.0) && (Energy < 3000.0)) 
						Gl = exp((13.1224-14.4902)/(7.82405-7.31322)*(log(Energy) - 7.31322)-13.1224);
					else Gl = -4.0; 
					
					Gl = Gl /4.0;
				}
				else {
					if (Energy > 2000.0) 
						Gl = exp(-0.554935*log(Energy) - 6.40179);
					else if (Energy > 800.0)
						Gl = exp((15.6942-10.6198)/(7.60090-6.68461)*(log(Energy)-6.68461)-15.6942);
					else 
						Gl = -1.0;
				}
			break;
			case 7:
				if (isGhost) {
					if (Energy < 1500.0)
						Gl = 2.0e-5; 
					else Gl = -1.0;
				}
				else {
					if (Energy > 1500.0)
						Gl = exp(-0.611844*log(Energy) -5.31592);
					else 
						Gl = -1.0;
				}
			break;
			default:
				Gl = -1.0;
			break;
		}
		
		return Gl;
	}
}
}