/* -*- Base: 10 ; Mode: C++ -*- */
/*------------------------------------------------------------------------
				     **	
  				FOST project
				     **
--------------------------------------------------------------------------
--------------------------------------------------------------------------
   				FILE LOG
		$Revision: 1.3 $    $Date: 2012-06-15 13:04:42 $
--------------------------------------------------------------------------
CREATION
	V.SAC

SUMMARY

DESCRIPTION

	The main function performs the following actions :
	<ul>
	<li>
	</ul>

------------------------------------------------------------------------*/

//=============================================================================
// 
// History of code
//
// creation 
//
// modification 
//=============================================================================

/**

*/
//=============================================================================
// Include section
//=============================================================================

// Standard libraries include files
//-----------------------------------------------------------------------------
#include <iostream>

// Include oriented definitions 
//-----------------------------------------------------------------------------

// Module Kernel include files
//-----------------------------------------------------------------------------

#include "DicError.hh"
#include "AMDA_exception.hh"
#include "ServicesServer.hh"
#include "PluginManager.hh"
#include "CurveFunctionWriter.hh"
#include "Boundaries.hh"

#include <math.h>

using namespace AMDA::Parameters;
using namespace plot;

// Other modules include files
//-----------------------------------------------------------------------------

// void venus_bowshock(double precision, CurveFunctionWriter::AttributeList& pAttributeList, CurveFunctionWriter::CurvePointList& curvePointList)
// {
// 	// param1 = 145.0 degrees
// 	// param2 = 1.35 Rv
// 	// param3 = 1.18
// 	// param4 = 0.55 Rv
// 	bowshock(precision, pAttributeList, curvePointList); 
// }
// 
// void mars_bowshock(double precision, CurveFunctionWriter::AttributeList& pAttributeList, CurveFunctionWriter::CurvePointList& curvePointList)
// {
// 	// param1 = 135.0 degrees
// 	// param2 = 2.04 Rm
// 	// param3 = 1.02
// 	// param4 = 0.55 Rm
// 	bowshock(precision, pAttributeList, curvePointList); 
// }
/**
 * @details This procedure return the curve definition of general bowshock
 *                r = L/(1+epsilon*cos(theta))
 */
void bowshock1(double precision, CurveFunctionWriter::AttributeList& pAttributeList, CurveFunctionWriter::CurvePointList& curvePointList)
{
	CurveFunctionWriter::AttributeList::iterator lparam1It = pAttributeList.find("param1");
	if ( lparam1It == pAttributeList.end()) 
	{
		BOOST_THROW_EXCEPTION(AMDA::AMDA_exception() << AMDA::errno_code(AMDA_INFORMATION_REQUEST_ERR) << AMDA::ex_msg(std::string("bowshock not found param1 attribute")));
	}
		
	CurveFunctionWriter::AttributeList::iterator lparam2It = pAttributeList.find("param2");
	if ( lparam2It == pAttributeList.end()) 
	{
		BOOST_THROW_EXCEPTION(AMDA::AMDA_exception() << AMDA::errno_code(AMDA_INFORMATION_REQUEST_ERR) << AMDA::ex_msg(std::string("bowshock not found param2 attribute")));
	}
	
	CurveFunctionWriter::AttributeList::iterator lparam3It = pAttributeList.find("param3");
	if ( lparam3It == pAttributeList.end()) 
	{
		BOOST_THROW_EXCEPTION(AMDA::AMDA_exception() << AMDA::errno_code(AMDA_INFORMATION_REQUEST_ERR) << AMDA::ex_msg(std::string("bowshock not found param3 attribute")));
	}
		
	
	for (int i = 0; i < precision; ++i)
	{
		double theta = (lparam1It->second * M_PI / 180.) * (double)(((double)i) / precision);
		double rr = lparam2It->second / (1.0 + lparam3It->second * cos(theta));
		CurveFunctionWriter::CurvePoint point;
		point.x = rr * cos(theta);
		point.y = rr * sin(theta);
		point.z = point.y;
		//std::cout << "BRE - " << point.x << " - " << point.y << " - " << point.z << std::endl;
		curvePointList.push_back(point);
	}
}

void bowshock(double precision, CurveFunctionWriter::AttributeList& pAttributeList, CurveFunctionWriter::CurvePointList& curvePointList)
{
	CurveFunctionWriter::AttributeList::iterator lparam1It = pAttributeList.find("param1");
	if ( lparam1It == pAttributeList.end()) 
	{
		BOOST_THROW_EXCEPTION(AMDA::AMDA_exception() << AMDA::errno_code(AMDA_INFORMATION_REQUEST_ERR) << AMDA::ex_msg(std::string("bowshock not found param1 attribute")));
	}
		
	CurveFunctionWriter::AttributeList::iterator lparam2It = pAttributeList.find("param2");
	if ( lparam2It == pAttributeList.end()) 
	{
		BOOST_THROW_EXCEPTION(AMDA::AMDA_exception() << AMDA::errno_code(AMDA_INFORMATION_REQUEST_ERR) << AMDA::ex_msg(std::string("bowshock not found param2 attribute")));
	}
	
	CurveFunctionWriter::AttributeList::iterator lparam3It = pAttributeList.find("param3");
	if ( lparam3It == pAttributeList.end()) 
	{
		BOOST_THROW_EXCEPTION(AMDA::AMDA_exception() << AMDA::errno_code(AMDA_INFORMATION_REQUEST_ERR) << AMDA::ex_msg(std::string("bowshock not found param3 attribute")));
	}
	
	CurveFunctionWriter::AttributeList::iterator lparam4It = pAttributeList.find("param4");
	if ( lparam4It == pAttributeList.end()) 
	{
		BOOST_THROW_EXCEPTION(AMDA::AMDA_exception() << AMDA::errno_code(AMDA_INFORMATION_REQUEST_ERR) << AMDA::ex_msg(std::string("bowshock not found param4 attribute")));
	}
	
	for (int i = 0; i < precision; ++i)
	{
		double theta = (lparam1It->second * M_PI / 180.) * (double)(((double)i) / precision);
		double rr = lparam2It->second / (1.0 + lparam3It->second * cos(theta));
		CurveFunctionWriter::CurvePoint point;
		point.x = lparam4It->second + rr * cos(theta);
		point.y = rr * sin(theta);
		point.z = point.y;
		//std::cout << "BRE - " << point.x << " - " << point.y << " - " << point.z << std::endl;
		curvePointList.push_back(point);
	}
}
/**
 * @details This procedure return the curve definition of earch  bowshock (Farris; pressure dependent )
 *                r = L/(1+epsilon*cos(theta))
 */
void bowshock_farris(double precision, CurveFunctionWriter::AttributeList& pAttributeList, CurveFunctionWriter::CurvePointList& curvePointList)
{
	CurveFunctionWriter::AttributeList::iterator lparam1It = pAttributeList.find("param1");
	if ( lparam1It == pAttributeList.end()) 
	{
		BOOST_THROW_EXCEPTION(AMDA::AMDA_exception() << AMDA::errno_code(AMDA_INFORMATION_REQUEST_ERR) << AMDA::ex_msg(std::string("bowshock not found param1 attribute")));
	}
		
	CurveFunctionWriter::AttributeList::iterator lparam2It = pAttributeList.find("ramPressure");
	if ( lparam2It == pAttributeList.end()) 
	{
		BOOST_THROW_EXCEPTION(AMDA::AMDA_exception() << AMDA::errno_code(AMDA_INFORMATION_REQUEST_ERR) << AMDA::ex_msg(std::string("bowshock not found param2 attribute")));
	}	
	 
       double epsilon = 0.81, L0 = 24.8, ramP = 1.8;
	
	double L, power = -0.166 ;
	L =  L0 * pow(lparam2It->second/ramP, power);
	
	for (int i = 0; i < precision; ++i)
	{
		double theta = (lparam1It->second * M_PI / 180.) * (double)(((double)i) / precision);
		double rr = L / (1.0 + epsilon * cos(theta));
		CurveFunctionWriter::CurvePoint point;
		point.x = rr * cos(theta);
		point.y = rr * sin(theta);
		point.z = point.y;		
		curvePointList.push_back(point);
	}
}
/**
 * @details This procedure return the curve definition of saturn  bowshock (Went; 2011; pressure dependent )
 *                 r = (1 + epsilon)*c1*Pdyn^(-1/c2) /(1 + epsilon*cos(theta))
 */
void saturn_bowshock_went(double precision, CurveFunctionWriter::AttributeList& pAttributeList, CurveFunctionWriter::CurvePointList& curvePointList)
{
	CurveFunctionWriter::AttributeList::iterator lparam1It = pAttributeList.find("param1");
	if ( lparam1It == pAttributeList.end()) 
	{
		BOOST_THROW_EXCEPTION(AMDA::AMDA_exception() << AMDA::errno_code(AMDA_INFORMATION_REQUEST_ERR) << AMDA::ex_msg(std::string("bowshock not found param1 attribute")));
	}
	
	CurveFunctionWriter::AttributeList::iterator lparam2It = pAttributeList.find("ramPressure");
	if ( lparam2It == pAttributeList.end()) 
	{
		BOOST_THROW_EXCEPTION(AMDA::AMDA_exception() << AMDA::errno_code(AMDA_INFORMATION_REQUEST_ERR) << AMDA::ex_msg(std::string("bowshock not found ramPressure attribute")));
	}	
	
	CurveFunctionWriter::AttributeList::iterator lparam3It = pAttributeList.find("epsilon");
	if ( lparam3It == pAttributeList.end()) 
	{
		BOOST_THROW_EXCEPTION(AMDA::AMDA_exception() << AMDA::errno_code(AMDA_INFORMATION_REQUEST_ERR) << AMDA::ex_msg(std::string("bowshock not found epsilon attribute")));
	}
	
	CurveFunctionWriter::AttributeList::iterator lparam4It = pAttributeList.find("c1");
	if ( lparam4It == pAttributeList.end()) 
	{
		BOOST_THROW_EXCEPTION(AMDA::AMDA_exception() << AMDA::errno_code(AMDA_INFORMATION_REQUEST_ERR) << AMDA::ex_msg(std::string("bowshock not found c1 attribute")));
	}
	
	CurveFunctionWriter::AttributeList::iterator lparam5It = pAttributeList.find("c2");
	if ( lparam5It == pAttributeList.end()) 
	{
		BOOST_THROW_EXCEPTION(AMDA::AMDA_exception() << AMDA::errno_code(AMDA_INFORMATION_REQUEST_ERR) << AMDA::ex_msg(std::string("bowshock not found c2 attribute")));
	}

	double power, L;
 
	power = -1.0 / lparam5It->second;
	L = (1.0 + lparam3It->second) * lparam4It->second * pow(lparam2It->second,power);
	 
	for (int i = 0; i < precision; ++i)
	{
		double theta = (lparam1It->second * M_PI / 180.) * (double)(((double)i) / precision);
		double rr = L / (1.0 + lparam3It->second * cos(theta));
		CurveFunctionWriter::CurvePoint point;
		point.x = rr * cos(theta);
		point.y = rr * sin(theta);
		point.z = point.y;		
		curvePointList.push_back(point);
	}
}

/**
 * @details This procedure return the curve definition of venus magnetopause.
 */
void venus_magnetopause(double precision, CurveFunctionWriter::AttributeList& pAttributeList, CurveFunctionWriter::CurvePointList& curvePointList)
{
	//param1 = 1.07
	CurveFunctionWriter::AttributeList::iterator lparam1It = pAttributeList.find("param1");
	if ( lparam1It == pAttributeList.end())
	{
		BOOST_THROW_EXCEPTION(AMDA::AMDA_exception() << AMDA::errno_code(AMDA_INFORMATION_REQUEST_ERR) << AMDA::ex_msg(std::string("venus_magnetopause not found param1 attribute")));
	}

	//param2 = -4
	CurveFunctionWriter::AttributeList::iterator lparam2It = pAttributeList.find("param2");
	if ( lparam2It == pAttributeList.end())
	{
		BOOST_THROW_EXCEPTION(AMDA::AMDA_exception() << AMDA::errno_code(AMDA_INFORMATION_REQUEST_ERR) << AMDA::ex_msg(std::string("venus_magnetopause not found param2 attribute")));
	}

	//param3 = 1.5
	CurveFunctionWriter::AttributeList::iterator lparam3It = pAttributeList.find("param3");
	if ( lparam3It == pAttributeList.end())
	{
		BOOST_THROW_EXCEPTION(AMDA::AMDA_exception() << AMDA::errno_code(AMDA_INFORMATION_REQUEST_ERR) << AMDA::ex_msg(std::string("venus_magnetopause not found param3 attribute")));
	}

	//param4 = 3
	CurveFunctionWriter::AttributeList::iterator lparam4It = pAttributeList.find("param4");
	if ( lparam4It == pAttributeList.end())
	{
		BOOST_THROW_EXCEPTION(AMDA::AMDA_exception() << AMDA::errno_code(AMDA_INFORMATION_REQUEST_ERR) << AMDA::ex_msg(std::string("venus_magnetopause not found param4 attribute")));
	}

	double b = lparam2It->second / (lparam3It->second - lparam1It->second);
	double a = -b * lparam1It->second;

	//double rr;
	for (int i = 0; i < precision; ++i)
	{
		CurveFunctionWriter::CurvePoint point;

		point.y = lparam4It->second * (double)(((double)i) / precision);

		if (point.y < lparam1It->second)
			point.x = sqrt(lparam1It->second * lparam1It->second - point.y * point.y);
		else
			point.x = a + b * point.y;

		point.z = point.y;

		//std::cout << "BRE - " << point.x << " - " << point.y << " - " << point.z << std::endl;
		curvePointList.push_back(point);
	}
}

/**
 * @details This procedure return the curve definition of mars magnetopause.
 */
void mars_magnetopause(double precision, CurveFunctionWriter::AttributeList& pAttributeList, CurveFunctionWriter::CurvePointList& curvePointList)
{
	//param1 = 1.2
	CurveFunctionWriter::AttributeList::iterator lparam1It = pAttributeList.find("param1");
	if ( lparam1It == pAttributeList.end())
	{
		BOOST_THROW_EXCEPTION(AMDA::AMDA_exception() << AMDA::errno_code(AMDA_INFORMATION_REQUEST_ERR) << AMDA::ex_msg(std::string("mars_magnetopause not found param1 attribute")));
	}

	//param2 = -0.0992 
	CurveFunctionWriter::AttributeList::iterator lparam2It = pAttributeList.find("param2");
	if ( lparam2It == pAttributeList.end())
	{
		BOOST_THROW_EXCEPTION(AMDA::AMDA_exception() << AMDA::errno_code(AMDA_INFORMATION_REQUEST_ERR) << AMDA::ex_msg(std::string("mars_magnetopause not found param2 attribute")));
	}

	//param3 = -0.2967
	CurveFunctionWriter::AttributeList::iterator lparam3It = pAttributeList.find("param3");
	if ( lparam3It == pAttributeList.end())
	{
		BOOST_THROW_EXCEPTION(AMDA::AMDA_exception() << AMDA::errno_code(AMDA_INFORMATION_REQUEST_ERR) << AMDA::ex_msg(std::string("mars_magnetopause not found param3 attribute")));
	}

	//param4 = 3
	CurveFunctionWriter::AttributeList::iterator lparam4It = pAttributeList.find("param4");
	if ( lparam4It == pAttributeList.end())
	{
		BOOST_THROW_EXCEPTION(AMDA::AMDA_exception() << AMDA::errno_code(AMDA_INFORMATION_REQUEST_ERR) << AMDA::ex_msg(std::string("mars_magnetopause not found param4 attribute")));
	}

	//double rr;
	for (int i = 0; i < precision; ++i)
	{
		CurveFunctionWriter::CurvePoint point;

		point.y = lparam4It->second * (double)(((double)i) / precision);
		 
		point.x = lparam2It->second * point.y * point.y * point.y * point.y + lparam3It->second * point.y * point.y + lparam1It->second;

		point.z = point.y;

		//std::cout << "BRE - " << point.x << " - " << point.y << " - " << point.z << std::endl;
		curvePointList.push_back(point);
	}
}
/**
 * @details This procedure return the curve definition of Shue magnetopause.
 */
void magnetopause_shue(double precision, CurveFunctionWriter::AttributeList& pAttributeList, CurveFunctionWriter::CurvePointList& curvePointList)
{
	 
	CurveFunctionWriter::AttributeList::iterator lparam1It = pAttributeList.find("param1");
	if ( lparam1It == pAttributeList.end())
	{
		BOOST_THROW_EXCEPTION(AMDA::AMDA_exception() << AMDA::errno_code(AMDA_INFORMATION_REQUEST_ERR) << AMDA::ex_msg(std::string("magnetopause not found param1 attribute")));
	}

	 
	CurveFunctionWriter::AttributeList::iterator lparam2It = pAttributeList.find("ramPressure");
	if ( lparam2It == pAttributeList.end())
	{
		BOOST_THROW_EXCEPTION(AMDA::AMDA_exception() << AMDA::errno_code(AMDA_INFORMATION_REQUEST_ERR) << AMDA::ex_msg(std::string("magnetopause not found param2 attribute")));
	}
	 
	CurveFunctionWriter::AttributeList::iterator lparam3It = pAttributeList.find("Bz");
	if ( lparam3It == pAttributeList.end())
	{
		BOOST_THROW_EXCEPTION(AMDA::AMDA_exception() << AMDA::errno_code(AMDA_INFORMATION_REQUEST_ERR) << AMDA::ex_msg(std::string("mars_magnetopause not found param3 attribute")));
	}
	
// 	   alpha=(.58-.007*Bz)*(1.+0.024*alog(RamPress))
//        r0=(10.22+1.29*tanh(0.184*(Bz+8.14)))*RamPress^(-1./6.6)
//        rr=r0*(2./(1.+COS(theta)))^alpha
       
	double alpha, r0, power = -0.152;
	
	alpha = (0.58 - 0.007 * lparam3It->second) * (1.0 + 0.024 * log(lparam2It->second));
	r0 = (10.22 + 1.29 * tanh(0.184 * (lparam3It->second + 8.14))) * pow(lparam2It->second,power);
		
	for (int i = 0; i < precision; ++i)
	{
		CurveFunctionWriter::CurvePoint point;
		
		double theta = (lparam1It->second * M_PI / 180.) * (double)(((double)i) / precision);
		double rr = r0 * pow((2.0 / (1.0 + cos(theta))),alpha);
		
		point.y = rr * sin(theta);
		 
		point.x = rr * cos(theta);

		point.z = point.y;

		//std::cout << "BRE - " << point.x << " - " << point.y << " - " << point.z << std::endl;
		curvePointList.push_back(point);
	}
}
/**
 * @details This procedure return the curve definition of Shue magnetopause for Mercure (Winslow; JGR;May 2013).
 */
void mercury_magnetopause(double precision, CurveFunctionWriter::AttributeList& pAttributeList, CurveFunctionWriter::CurvePointList& curvePointList)
{
	 
	CurveFunctionWriter::AttributeList::iterator lparam1It = pAttributeList.find("param1");
	if ( lparam1It == pAttributeList.end())
	{
		BOOST_THROW_EXCEPTION(AMDA::AMDA_exception() << AMDA::errno_code(AMDA_INFORMATION_REQUEST_ERR) << AMDA::ex_msg(std::string("magnetopause not found param1 attribute")));
	}

	 
	CurveFunctionWriter::AttributeList::iterator lparam2It = pAttributeList.find("rnose");
	if ( lparam2It == pAttributeList.end())
	{
		BOOST_THROW_EXCEPTION(AMDA::AMDA_exception() << AMDA::errno_code(AMDA_INFORMATION_REQUEST_ERR) << AMDA::ex_msg(std::string("magnetopause not found param2 attribute")));
	}
	 
	CurveFunctionWriter::AttributeList::iterator lparam3It = pAttributeList.find("alpha");
	if ( lparam3It == pAttributeList.end())
	{
		BOOST_THROW_EXCEPTION(AMDA::AMDA_exception() << AMDA::errno_code(AMDA_INFORMATION_REQUEST_ERR) << AMDA::ex_msg(std::string("mars_magnetopause not found param3 attribute")));
	}
	
//        rr=rnose*(2./(1.+COS(theta)))^alpha
       
	for (int i = 0; i < precision; ++i)
	{
		CurveFunctionWriter::CurvePoint point;
		
		double theta = (lparam1It->second * M_PI / 180.) * (double)(((double)i) / precision);
		double rr = lparam2It->second * pow((2.0 / (1.0 + cos(theta))),lparam3It->second);
		
		point.y = rr * sin(theta);
		 
		point.x = rr * cos(theta);

		point.z = point.y;		
		curvePointList.push_back(point);
	}
}
/**
 * @details This procedure return the curve definition of Shue magnetopause for Saturn (Arridge; JGR; 2006).
 */
void saturn_magnetopause(double precision, CurveFunctionWriter::AttributeList& pAttributeList, CurveFunctionWriter::CurvePointList& curvePointList)
{
	 
	CurveFunctionWriter::AttributeList::iterator lparam1It = pAttributeList.find("param1");
	if ( lparam1It == pAttributeList.end())
	{
		BOOST_THROW_EXCEPTION(AMDA::AMDA_exception() << AMDA::errno_code(AMDA_INFORMATION_REQUEST_ERR) << AMDA::ex_msg(std::string("magnetopause not found param1 attribute")));
	}
	 
	CurveFunctionWriter::AttributeList::iterator lparam2It = pAttributeList.find("ramPressure");
	if ( lparam2It == pAttributeList.end())
	{
		BOOST_THROW_EXCEPTION(AMDA::AMDA_exception() << AMDA::errno_code(AMDA_INFORMATION_REQUEST_ERR) << AMDA::ex_msg(std::string("magnetopause not found param2 attribute")));
	}
	 	 
	
//        rr=rnose*(2./(1.+COS(theta)))^alpha
//        rnose = a1*RamPress^(-a2)
//        alpha = a3 + a4*RamPress
       
	double a1 = 9.7, a2 = -0.24, a3 = 0.77, a4 = -1.5;
	double rnose, alpha;
	
	rnose = a1 * pow(lparam2It->second,a2);
	alpha = a3 + a4 * lparam2It->second;
	
	for (int i = 0; i < precision; ++i)
	{
		CurveFunctionWriter::CurvePoint point;
		
		double theta = (lparam1It->second * M_PI / 180.) * (double)(((double)i) / precision);
		double rr = rnose * pow((2.0 / (1.0 + cos(theta))),alpha);
		
		point.y = rr * sin(theta);
		 
		point.x = rr * cos(theta);

		point.z = point.y;		
		curvePointList.push_back(point);
	}
}