/////////////////
// CPP HEADERS //
/////////////////
#include <iostream>
#include <fstream>

#include "math.h"

////////////
/// BOOST //
////////////
#include <boost/tokenizer.hpp>
#include <boost/lexical_cast.hpp>
#include <boost/regex.hpp>
#include "boost/date_time/gregorian/gregorian.hpp"
#include "boost/date_time/posix_time/posix_time.hpp"

///////////////
// TINYXML-2 //
///////////////
#include "tinyxml2.h"

////////////
/// DDLIB //
////////////
#include "DD.h"

//////////
/// GSL //
//////////
#include "gsl_spline.h"

//////////////////////////
// FUNCTIONS PROTOTYPES //
//////////////////////////
#include "tools.h"

////////////////
// NAMESPACES //
///////////////
using namespace std;
using namespace tinyxml2;



////////////////////////////
/// FUNCTIONS DEFINITIONS //
////////////////////////////

/**
 * @brief      Print mission features
 *
 * @param      satelliteName   Name of the satellite
 * @param      satelliteVi     Virtual instrument's name of the satellite
 * @param      satelliteOrbit  Name of the netCDF variable for the satellite
 * @param 	   satelliteRes    Satellite orbit resolution
 * @param      plasmaName      Name of the plasma dataset
 * @param      plasmaVi        Virtual instrument's name of the plasma dataset
 * @param      plasmaVNorm     Name of the netCDF variable for plasma dataset
 * @param      plasmaRamP      Ram pressure parameter into the netCDF file
 * @param      plasmaVel       Velocity parameter into the netCDF file
 * @param      plasmaMag       Magnetic field parameter into the netCDF file
 * @param      indexName       Name of the Sym H index
 * @param      indexVi         Virtual instrument's name of the Sym H index
 * @param      indexSym        Sym H index parameter into the netCDF file
 * @param      start           Start time
 * @param      stop            Stop time
 * @param      ddStart         DD style start time
 * @param      ddStop          DD style stop time
 * @param      ddInterval      DD style time interval
 */
void printMission(string& satelliteName, string& satelliteVi, 
				string& satelliteOrbit, string& satelliteRes,
				string& plasmaName, string& plasmaVi, string& plasmaVNorm, 
				string& plasmaRamP, string& plasmaVel, string& plasmaMag,
				string& indexName, string& indexVi, string& indexSym, 
				string& start, string& stop, string& ddStart, string& ddStop, string& ddInterval)
{
	cout << "********** MISSION ************" << endl;
	cout << "Satellite name : " << satelliteName << endl;
	cout << "Satellite VI : " << satelliteVi << endl;
	cout << "Satellite orbit param : " << satelliteOrbit << endl;
	cout << "Satellite orbit resolution : " << satelliteRes << endl;
	cout << "Plasma name : " << plasmaName << endl;
	cout << "Plasma VI : " << plasmaVi << endl;
	cout << "Plasma V norm param : " << plasmaVNorm << endl;
	cout << "Plasma ram pressure : " << plasmaRamP << endl;
	cout << "Plasma velocity : " << plasmaVel << endl;
	cout << "Plasma B mag : " << plasmaMag << endl;
	cout << "Sym index name : " << indexName << endl;
	cout << "Sym index VI : " << indexVi << endl;
	cout << "Sym index param : " << indexSym << endl;
	cout << "Start time : " << start << endl;
	cout << "Stop time : " << stop << endl;
	cout << "DD Start time : " << ddStart << endl;
	cout << "DD Stop time : " << ddStop << endl;
	cout << "DD Interval : " << ddInterval << endl;
	cout << "********** END OF MISSION ************" << endl << endl;
}



/**
 * @brief      Get mission parameters
 *
 * @param      filename        Name of the XML config file
 * @param      satelliteName   Name of the satellite
 * @param      satelliteVi     Virtual instrument's name for the satellite
 * @param      satelliteOrbit  Name of the parameter into the netCDF file for the satellite's orbit
 * @param      satelliteRes	   Time resolution for satellite orbit
 * @param      plasmaName      Name of plasma dataset
 * @param      plasmaVi        Virtual instrument's name for the plasma dataset
 * @param      plasmaVNorm     Name of the parameter into the netCDF file for the plasma's velocity norm
 * @param      plasmaRamP      Name of the parameter into the netCDF file for the plasma's ram pressure
 * @param      plasmaVel       Name of the parameter into the netCDF file for the plasma's velocity vector
 * @param      plasmaMag       Name of the parameter into the netCDF file for the plasma's magnetic field vector
 * @param      indexName       Name of the Sym index
 * @param      indexVi         Virtual instrument's name for the Sym index
 * @param      indexSym        Name of the parameter into the netCDF file for the Sym index
 * @param      start           Start time
 * @param      stop            Stop time
 */
void loadMission(string& filename, 
	string& satelliteName, string& satelliteVi, 
	string& satelliteOrbit, string& satelliteRes,
	string& plasmaName, string& plasmaVi, string& plasmaVNorm,
	string& plasmaRamP, string& plasmaVel, string& plasmaMag,
	string& indexName, string& indexVi, string& indexSym,
	string& start, string& stop)
{

	// LOADING THE DOM
	XMLDocument doc;
	
	if ( doc.LoadFile( filename.c_str() ) != 0 )
	{
		cout << "[ERROR] loadMission : Failed to load " << filename << endl;
		
		exit(EXIT_FAILURE);
	}
	

	// GET PARAMETERS VALUES
	XMLElement* pRoot = doc.FirstChildElement( "MISSION" );
	
	satelliteName              = pRoot->FirstChildElement( "SATELLITE" )->FirstChildElement( "NAME" )->GetText();
	satelliteVi                = pRoot->FirstChildElement( "SATELLITE" )->FirstChildElement( "VI" )->GetText();
	satelliteOrbit = pRoot->FirstChildElement( "SATELLITE" )->FirstChildElement( "ORBIT" )->GetText();
	satelliteRes = pRoot->FirstChildElement( "SATELLITE" )->FirstChildElement( "RES" )->GetText();

	plasmaName              = pRoot->FirstChildElement( "PLASMA" )->FirstChildElement( "NAME" )->GetText();
	plasmaVi                = pRoot->FirstChildElement( "PLASMA" )->FirstChildElement( "VI" )->GetText();
	plasmaVNorm = pRoot->FirstChildElement( "PLASMA" )->FirstChildElement( "VNORM" )->GetText();
	plasmaRamP = pRoot->FirstChildElement( "PLASMA" )->FirstChildElement( "RAMP" )->GetText();
	plasmaVel = pRoot->FirstChildElement( "PLASMA" )->FirstChildElement( "VEL" )->GetText();
	plasmaMag = pRoot->FirstChildElement( "PLASMA" )->FirstChildElement( "MAG" )->GetText();

	indexName = pRoot->FirstChildElement( "INDEX" )->FirstChildElement( "NAME" )->GetText();
	indexVi = pRoot->FirstChildElement( "INDEX" )->FirstChildElement( "VI" )->GetText();
	indexSym = pRoot->FirstChildElement( "INDEX" )->FirstChildElement( "SYM" )->GetText();

	start             = pRoot->FirstChildElement( "START" )->GetText();
	stop              = pRoot->FirstChildElement( "STOP" )->GetText();

}


/**
 * @brief      Split date time into 2 parts : date and clock time
 *
 * @param[in]  strDate  Date to split
 *
 * @return     Array containing the splitted date : [0] -> Date / [1] -> Clock time
 */
vector<string> tokenizedDate(string strDate)
{
	vector<string> tDate;

  	boost::char_separator<char> sep("T");
	boost::tokenizer< boost::char_separator<char> > tok(strDate, sep); 
  	boost::tokenizer< boost::char_separator<char> >::iterator tok_iter;
  
	for(tok_iter = tok.begin(); tok_iter != tok.end(); ++tok_iter)
	{
		tDate.push_back(*tok_iter);
	}

	// tDate[1] = tDate[1].substr(0, tDate[1].size()-1);

  	return tDate;
}


/**
 * @brief      Check if date is in good format : yyyy-mm-ddThh:mm:ss.mls (modify if a pattern is matched)
 *
 * @param      date  Date to check for
 *
 * @return     True : good format / False : bad format
 */
bool checkDateFormat(string& date)
{
	boost::regex regex0("^[0-9]{4}[-]{1}[0-9]{2}[-]{1}[0-9]{2}[T]{1}[0-9]{2}[:]{1}[0-9]{2}[:]{1}[0-9]{2}$");
	boost::regex regex1("^[0-9]{4}[-]{1}[0-9]{2}[-]{1}[0-9]{2}[T]{1}[0-9]{2}[:]{1}[0-9]{2}[:]{1}[0-9]{2}[.]{1}[0-9]{3}$");
    boost::regex regex2("^[0-9]{4}[-]{1}[0-9]{2}[-]{1}[0-9]{2}[T]{1}[0-9]{2}[:]{1}[0-9]{2}[:]{1}[0-9]{2}[.]{1}[0-9]{6}$");
    boost::regex regex3("^[0-9]{4}[-]{1}[0-9]{2}[-]{1}[0-9]{2}[T]{1}[0-9]{2}[:]{1}[0-9]{2}[:]{1}[0-9]{2}[.]{1}[0-9]{9}$");
    boost::regex regex4("^[0-9]{4}[-]{1}[0-9]{2}[-]{1}[0-9]{2}[T]{1}[0-9]{2}[:]{1}[0-9]{2}[:]{1}[0-9]{2}[.]{1}[0-9]{3}[Z]{1}$");
    
    bool match0 = boost::regex_match(date, regex0);
    bool match1 = boost::regex_match(date, regex1);
    bool match2 = boost::regex_match(date, regex2);
    bool match3 = boost::regex_match(date, regex3);
    bool match4 = boost::regex_match(date, regex4);

    if (match0)
    {
    	date += ".000";
    	return true; 
    }
    else if (match1)
    {
    	return true; 
    }
    else if (match2 || match3)
    {
    	date = date.substr(0, 23);
    	return true; 
    }
    else if (match4)
    {
    	date = date.substr(0, date.size()-1);
    	return true; 
    }
    else
    {
    	return false; 
    }
    
}


/**
 * @brief      Date to DD date : yyyydoyhhmmssmls
 *
 * @param[in]  date  Date to transform
 *
 * @return     DD date
 */
string toDDTime(string date)
{
	if ( !checkDateFormat(date) )
	{
		string badDateFormatErr = "[ERROR] Tools::toDDTime -> Start time is bad formatted " + date + " : cannot build DD start time";
      	throw std::string(badDateFormatErr);
	} 

	string dateTime  = (tokenizedDate(date))[0];
	string clockTime = (tokenizedDate(date))[1];
	string gregoDate = dateTime.substr(0, 4)+dateTime.substr(5, 2)+dateTime.substr(8, 2);
    
    boost::gregorian::date d( boost::gregorian::from_undelimited_string(gregoDate) );
    string doy = boost::lexical_cast<std::string>( d.day_of_year()-1 ); 
    
    if ( doy.size() < 1 )
    {
    	throw std::string("[ERROR] Tools::toDDTime -> Failed to build day of year");
    }
    else if ( doy.size() < 2 )
    {
      doy = "00" + doy;
    }
    else if( doy.size() < 3 )
    {
      doy = "0" + doy;
    }
    else if ( doy.size() > 3 )
    {
    	throw std::string("[ERROR] Tools::toDDTime -> Failed to build day of year");
    }
    
    string ddTime = dateTime.substr(0, 4)+doy+clockTime.substr(0, 2)+clockTime.substr(3, 2)+clockTime.substr(6, 2)+clockTime.substr(9);
    
	return ddTime;
}


/**
 * @brief      Make DD style time intervall : yyyydoyhhmmssmls
 *
 * @param[in]  startTime  Start time
 * @param[in]  stopTime   Stop time
 *
 * @return     DD style time intervall
 */
string makeDDIntervall(string startTime, string stopTime)
{
	if ( !(checkDateFormat(startTime) && checkDateFormat(stopTime)) )
	{
		string badBoundsErr = "[ERROR] makeDDIntervall -> Failed to build DD intervall because of bad time boundaries : " + startTime + " / " + stopTime;
		throw std::string(badBoundsErr);
	}

	vector<string> startVec = tokenizedDate(startTime);
	vector<string> stopVec = tokenizedDate(stopTime);

	boost::posix_time::ptime ptStart(boost::posix_time::time_from_string(startVec[0]+ " " + startVec[1]));
	boost::posix_time::ptime ptStop(boost::posix_time::time_from_string(stopVec[0]+ " " + stopVec[1]));
	boost::posix_time::time_period tp(ptStart, ptStop);
	
	if ( tp.is_null() )
	{
		string badBoundsErr = "[ERROR] makeDDIntervall -> Failed to build DD intervall because of bad time boundaries : " + startTime + " / " + stopTime;
		throw std::string(badBoundsErr);
	}
	
	string days = boost::lexical_cast<string>( tp.length().hours()/24 );
	string hours = boost::lexical_cast<string>( tp.length().hours()%24 );
	string minutes = boost::lexical_cast<string>( tp.length().minutes() );
	string seconds = boost::lexical_cast<string>( tp.length().seconds() );
	string mls = boost::lexical_cast<string>( tp.length().fractional_seconds() );
	
	if ( days.size() < 2 )
	{
		days = "00" + days;
	}
	else if ( days.size() < 3 )
	{
		days = "0" + days;
	}

	if ( hours.size() < 2 )
	{
		hours = "0" + hours;
	}

	if ( minutes.size() < 2 )
	{
		minutes = "0" + minutes;
	}

	if ( seconds.size() < 2 )
	{
		seconds = "0" + seconds;
	}

	if ( mls.size() < 2 )
	{
		mls = "00" + mls;
	}
	else if ( mls.size() == 6 )
	{
		mls = mls.substr(0, 3);
	}
	else if ( mls.size() == 9 )
	{
		mls = mls.substr(0, 3);
	}

	string ddIntervall = "0000" + days + hours + minutes + seconds + mls;

	return ddIntervall;
}



/**
 * @brief      Get plasma data from DD base
 *
 * @param[in]  plasmaVi        Name of the VI into DD base
 * @param[in]  ddStart         Start time DD time format
 * @param[in]  ddInterval      Time interval DD time format
 * @param[in]  timeParameter   Name of time parameter into NC file
 * @param[in]  vnormParameter  Name of velocity norm parameter into NC file
 * @param[in]  rampParameter   Name of ram pressure parameter into NC file
 * @param[in]  velParameter    Name of velocity parameter into NC file
 * @param[in]  magParameter    Name of magnetic field parameter into NC file
 * @param      plasmaTimeArray Time array
 * @param      vnormArray      Velocity norm array
 * @param      rampArray       Ram pressure array
 * @param      vxArray         Vx component of the velocity array
 * @param      vyArray         Vy component of the velocity array
 * @param      vzArray         Vz component of the velocity array
 * @param      bxArray         Bx component of the velocity array
 * @param      byArray         By component of the velocity array
 * @param      bzArray         Bz component of the velocity array
 */
void getPlasmaData(string plasmaVi, string ddStart, string ddInterval,
	string timeParameter, string vnormParameter, string rampParameter,
	string velParameter, string magParameter,
	vector<string>& plasmaTimeArray, vector<float>& vnormArray,
	vector<float>& rampArray, 
	vector<float>& vxArray, vector<float>& vyArray, vector<float>& vzArray,
	vector<float>& bxArray, vector<float>& byArray, vector<float>& bzArray)
{
	int id = DD_SetVariable( const_cast<char*>( plasmaVi.c_str() ) );

	char *st = const_cast<char*>( ddStart.c_str() );

	int error = DD_SetTime(id, st);

	if (error < 0)
	{
		std::string ddSetTimeErr = "[ERROR] getPlasmaData : bad time pointer init in DD_SetTime -> err value : " +  error;
		throw std::string(ddSetTimeErr);
	}

	// GET DATA
	char *params[5];
	
	params[0] = const_cast<char*>( timeParameter.c_str() );
	params[1] = const_cast<char*>( vnormParameter.c_str() );
	params[2] = const_cast<char*>( rampParameter.c_str() );
	params[3] = const_cast<char*>( velParameter.c_str() );
	params[4] = const_cast<char*>( magParameter.c_str() );
	
	char *timeIntervall = const_cast<char*>( ddInterval.c_str() );
	
	DD_data_t *data;

	int status = 0;

	do {

		status = DD_GetMultiData(id, 5, static_cast<char**>(params), timeIntervall, &data, 1);
		
		if (status < 0)
		{
			std::string ddGetDataErr = "[ERROR] getPlasmaData : failed to get data -> status : " + status;
			throw std::string(ddGetDataErr);
		}

		for (int j = 0; j < data->VarNumber; ++j)
		{
			plasmaTimeArray.push_back( static_cast<char*>(data[0].Variables[j]) );
		  	vnormArray.push_back( *(static_cast<float*>(data[1].Variables[j])) );
			rampArray.push_back( *(static_cast<float*>(data[2].Variables[j])) );
			vxArray.push_back( static_cast<float*>(data[3].Variables[j])[0] );
		  	vyArray.push_back( static_cast<float*>(data[3].Variables[j])[1] );
		  	vzArray.push_back( static_cast<float*>(data[3].Variables[j])[2] );
		  	bxArray.push_back( static_cast<float*>(data[4].Variables[j])[0] );
		  	byArray.push_back( static_cast<float*>(data[4].Variables[j])[1] );
		  	bzArray.push_back( static_cast<float*>(data[4].Variables[j])[2] );

		}


	} while (status == MOREDATA);
	
    DD_Close(id);
}


/**
 * @brief      Get orbit data from DD base
 *
 * @param[in]  orbitVi         Name of the VI into DD base
 * @param[in]  orbitRes        Resolution of the satellite orbit
 * @param[in]  ddStart         Start time DD time format
 * @param[in]  ddInterval      Time interval DD time format
 * @param[in]  timeParameter   Name of time parameter into NC file
 * @param[in]  orbitParameter  Name of orbit parameter into NC file
 * @param      orbitArray      Time array
 * @param      x               X component array
 * @param      y               Y component array
 * @param      z               Z component array
 */
void getOrbitData(string satelliteVi, string orbitRes,
	string ddStart, string ddInterval,
	string timeParameter, string orbitParameter, 
	vector<string>& satelliteTimeArray,
	vector<float>& x, vector<float>& y, vector<float>& z)
{
	int id = DD_SetVariable( const_cast<char*>( satelliteVi.c_str() ) );

    char *st = const_cast<char*>( ddStart.c_str() );

	int error = DD_SetTime(id, st);

	if (error < 0)
	{
		std::string ddSetTimeErr = "[ERROR] getOrbitData : bad time pointer init in DD_SetTime -> err value : " +  error;
		throw std::string(ddSetTimeErr);
	}

	char *params[2];
	
	params[0] = const_cast<char*>( timeParameter.c_str() );
	params[1] = const_cast<char*>( orbitParameter.c_str() );
	
	char *timeIntervall = const_cast<char*>( ddInterval.c_str() );

	DD_data_t *data;

	int status = 0;

	do {

		status = DD_GetMultiData(id, 2, static_cast<char**>(params), timeIntervall, &data, 1);
		
		if (status < 0)
		{
			std::string ddGetDataErr = "[ERROR] getOrbitData : failed to get data -> status : " + status;
			throw std::string(ddGetDataErr);
		}

		int intOrbitStep = static_cast<int>( boost::lexical_cast<int>(orbitRes) );
		int samplingStep;

		if ( intOrbitStep < 300 )
		{	
			// Resampling to 5min resolution for resolution strictly superior to 300s
			// samplingStep = intOrbitStep/60;
			samplingStep = static_cast<int>( 300/intOrbitStep );

		} else {
			// Else : interpolation to do ...
			samplingStep = 1;
		}

		int dataType = data[1].type;
		
		for (int j = 0; j < data->VarNumber; j+=samplingStep)	
		{
			satelliteTimeArray.push_back( static_cast<char*>(data[0].Variables[j]) );
		  	
			if ( dataType == 3 ) // DOUBLE
			{
				x.push_back( static_cast<float>(static_cast<double*>(data[1].Variables[j])[0]) );
		  		y.push_back( static_cast<float>(static_cast<double*>(data[1].Variables[j])[1]) );
		  		z.push_back( static_cast<float>(static_cast<double*>(data[1].Variables[j])[2]) );

			}
			else if ( dataType == 2 ) // FLOAT
			{
				
				x.push_back( static_cast<float*>(data[1].Variables[j])[0] );
		  		y.push_back( static_cast<float*>(data[1].Variables[j])[1] );
		  		z.push_back( static_cast<float*>(data[1].Variables[j])[2] );

			}
			else
			{
				std::string ddGetDataTypeErr = "[ERROR] getOrbitData : failed to get data -> Unknown data type (not double nor float)";
				throw std::string(ddGetDataTypeErr);	
			}
		}


	} while (status == MOREDATA);
	
    DD_Close(id);
}


/**
 * @brief      Get Sym H index data from DD base
 *
 * @param[in]  symIndexVi         Name of the VI into DD base
 * @param[in]  ddStart            Start time DD time format
 * @param[in]  ddInterval         Time interval DD time format
 * @param[in]  timeParameter      Name of time parameter into NC file
 * @param[in]  symIndexParameter  Name of Sym H index parameter into NC file
 * @param      symIndexTimeArray  Time array
 * @param      symIndexArray      Sym H index array
 */
void getSymIndexData(string symIndexVi, string ddStart, string ddInterval,
	string timeParameter, string symIndexParameter, 
	vector<string>& symIndexTimeArray, vector<float>& symIndexArray)
{
	int id = DD_SetVariable( const_cast<char*>( symIndexVi.c_str() ) );

    char *st = const_cast<char*>( ddStart.c_str() );

	int error = DD_SetTime(id, st);

	if (error < 0)
	{
		std::string ddSetTimeErr = "[ERROR] getSymIndexData : bad time pointer init in DD_SetTime -> err value : " +  error;
		throw std::string(ddSetTimeErr);
	}

	char *params[2];
	
	params[0] = const_cast<char*>( timeParameter.c_str() );
	params[1] = const_cast<char*>( symIndexParameter.c_str() );
	
	char *timeIntervall = const_cast<char*>( ddInterval.c_str() );

	DD_data_t *data;

	int status = 0;

	do {

		status = DD_GetMultiData(id, 2, static_cast<char**>(params), timeIntervall, &data, 1);
		
		if (status < 0)
		{
			std::string ddGetDataErr = "[ERROR] getSymIndexData : failed to get data -> status : " + status;
			throw std::string(ddGetDataErr);
		}

		// 5 min to 1 min resampling
		for (int j = 0; j < data->VarNumber; j+=5)
		{
			symIndexTimeArray.push_back( static_cast<char*>(data[0].Variables[j]) );
		  	symIndexArray.push_back( static_cast<float*>(data[1].Variables[j])[0] );
		}


	} while (status == MOREDATA);
	
    DD_Close(id);
}


/**
 * @brief      Remove NAN values
 *
 * @param      n              Size of arrays
 * @param      rampArray      Ram pressure array
 * @param      vnormArray     Velocity norm array
 * @param      vx             Vx component array
 * @param      vy             Vy component array
 * @param      vz             Vz component array
 * @param      bx             Bx component array
 * @param      by             By component array
 * @param      bz             Bz component array
 * @param      symIndexArray  Sym H index array
 * @param      flag           Quality flag array
 */
void cleanData(int& n, vector<float>& rampArray, vector<float>& vnormArray,
				vector<float>& vx, vector<float>& vy, vector<float>& vz,
				vector<float>& bx, vector<float>& by, vector<float>& bz,
				vector<float>& symIndexArray, vector<int>& flag)		
{
	for (int i = 0; i < n; ++i)
	{
		// No value for plasma data
		if ( isnan(rampArray[i]) )
		{
			rampArray[i] = 1.5;
			vnormArray[i] = 500.0;
			vx[i] = -500.0;
			vy[i] = -29.78;
			// vy[i] = 0.0;
			vz[i] = 0.0;

			flag[i] = -1;
		}

		// No value for magnetic field data
		if ( isnan(by[i]) )
		{
			bx[i] = 0.0;
	  		by[i] = 0.0;
	  		bz[i] = 1.0;

	  		flag[i] = -1;
		}

		// No value for Sym H index data
		if ( isnan(symIndexArray[i]) )
		{
			symIndexArray[i] = -10.0;

			flag[i] = -1;
		}
	}
}

void writeInputsFile(string& filename, int& size, 
					vector<string>& plasmaTimeArray,
					vector<float>& rampArray, vector<float>& vnormArray,
					vector<float>& vx, vector<float>& vy, vector<float>& vz,
					vector<float>& bx, vector<float>& by, vector<float>& bz,
					vector<float>& x, vector<float>& y, vector<float>& z,
					vector<float>& symIndexArray, 
					vector<int>& flag)
{
	ofstream inFile(filename.c_str(), ios::out | ios::trunc);

	string year, day, hour, minute;

    if(inFile)
    {
    	for (int i = 0; i < size; ++i)
		{
			year = plasmaTimeArray[i].substr(0, 4);
			day = plasmaTimeArray[i].substr(4, 3);
			hour = plasmaTimeArray[i].substr(7, 2);
			minute = plasmaTimeArray[i].substr(9, 2);
			
			int yearInt = boost::lexical_cast<int>(year);

			int dayInt = boost::lexical_cast<int>(day);
			dayInt++;

			int hourInt = boost::lexical_cast<int>(hour);

			int minuteInt = boost::lexical_cast<int>(minute);

			inFile << yearInt << " " 
					<< dayInt << " " 
					<< hourInt << " " 
					<< minuteInt << " "
					<< x[i] << " " << y[i] << " " << z[i] << " "
					<< bx[i] << " " << by[i] << " " << bz[i] << " "
					<< vx[i] << " " << vy[i] + 29.78 << " " << vz[i] << " "
					// << vx[i] << " " << vy[i] << " " << vz[i] << " "
					<< rampArray[i] << " "
					<< symIndexArray[i] << " "
					<< flag[i]
					<< endl;
		}

		inFile.close();
    }
}




/**
 * @brief      Linear interpolation for satellite orbit data
 *
 * @param[in]  satTime  Time array of the satellite
 * @param      x        X component array
 * @param      y        Y component array
 * @param      z        Z component array
 */
void linearInterpolation(vector<string> satTime,  
						vector<float>& x, vector<float>& y, vector<float>& z)
{
	int n = satTime.size();
	double *p_t = new double[n];
	double *p_x = new double[n];
	double *p_y = new double[n];
	double *p_z = new double[n];

	for (int i = 0; i < n; ++i)
	{
		p_t[i] = boost::lexical_cast<double>(satTime[i]);
		p_x[i] = boost::lexical_cast<double>(x[i]);
		p_y[i] = boost::lexical_cast<double>(y[i]);
		p_z[i] = boost::lexical_cast<double>(z[i]);
	}

    gsl_interp *interpolation = gsl_interp_alloc (gsl_interp_linear, n);

   	gsl_interp_init(interpolation, p_t, p_x, n);
    gsl_interp_init(interpolation, p_t, p_y, n);
    gsl_interp_init(interpolation, p_t, p_z, n);

    gsl_interp_accel * accelerator =  gsl_interp_accel_alloc();

	double startTime = p_t[0];
	double stopTime  = p_t[n-1];
	double step      = 300;

	vector<float> newX;
	vector<float> newY;
	vector<float> newZ;
	double nx, ny, nz;

    for (double ti = startTime; ti < stopTime; ti += step)
    {
    	nx = gsl_interp_eval(interpolation, p_t, p_x, ti, accelerator);
    	ny = gsl_interp_eval(interpolation, p_t, p_y, ti, accelerator);
    	nz = gsl_interp_eval(interpolation, p_t, p_z, ti, accelerator);
    	
    	newX.push_back( boost::lexical_cast<float>(nx) );
    	newY.push_back( boost::lexical_cast<float>(ny) );
    	newZ.push_back( boost::lexical_cast<float>(nz) );
    }
    
	x = newX;
	y = newY;
	z = newZ;

    delete p_t;
    delete p_x;
    delete p_y;
    delete p_z;

    gsl_interp_free(interpolation);
    gsl_interp_accel_free(accelerator);
    
}


/**
 * @brief      Bad looking tool to add new data for Geotail orbit (x2)
 *
 * @param      x     X component array
 * @param      y     Y component array
 * @param      z     Z component array
 */
void GTL_interp(vector<float>& x, vector<float>& y, vector<float>& z)
{
	vector<float> newX, newY, newZ;
	int size = x.size();

	for (int i = 0; i < size; ++i)
	{
		newX.push_back(x[i]);
		newX.push_back(x[i]);

		newY.push_back(y[i]);
		newY.push_back(y[i]);

		newZ.push_back(z[i]);
		newZ.push_back(z[i]);
	}

	x = newX;
	y = newY;
	z = newZ;
}