#ifndef STATION_C
#define STATION_C
/**
*/
#include <iostream>
#include <fstream>
#include <stdlib.h>
#include <string.h>
#include <vector>
#include <cstdarg>
#include <string>
#include <math.h>

#include <stdio.h>
#include <unistd.h>

#include <limits.h>

#include <CCfits>

using namespace std ;
using namespace CCfits ;

#include "Station.h"

/**
* Constructor and default ctor. Defaults are the UVW platform zero 
* according to VLTSW20050050.
*/
Station::Station(const double lambda, const double phi, const string name, const double geoel)
	: Geocen(lambda,phi), geoRadius(rEarth+geoel), id(name), aperture(1.8)
{
	/* debugging
	* cout << setprecision(16) << lambda << " " << phi << endl ;
	*/
}

/**
* ctor with a standard station name
*/
Station::Station( const string name) : aperture(1.8)
{
	staIndx=init(name) ;
}

/**
* ctor with a standard station number
* @param idx 1 to 34 for a regular station or 0 for the UVW coordinate center
*/
Station::Station( const int idx) : aperture(1.8)
{
	if ( idx >=0 && idx <= 34)
	{
		init( Station::posN[idx] ) ;
		staIndx= idx ;
	}
	else
	{
		cerr << __FILE__ << " " << __LINE__ << " invalid station no " << idx << endl ;
		exit(EXIT_FAILURE) ;
	}
}


/**
* ctor
* @param[in] f the FITS file with the CONF STATION keywords in the PHDU
* @param[in] indx 0 for the T1 and 1 for the T2
*/
Station::Station(FITS &f, const int indx)
{
	PHDU & p = f.pHDU() ;
	string hkey = "HIERARCH ESO ISS CONF STATION" ;
	switch(indx)
	{
	case 0 :
		hkey += "1 " ;
		break ;
	default :
		hkey += "2 " ;
		break ;
	}
	/**
	* We read HIERARCH ESO ISS CONF STATION{1|2} in the primary FITS header.
	*/
	string sta_name ;
	p.readKey(hkey, sta_name) ;

	/**
	* We take also HIERARCH ESO ISS GEOELEV from the primary FITS header.
	*/
	double geoel ;
	p.readKey("HIERARCH ESO ISS GEOELEV", geoel) ;

	staIndx = init(sta_name,geoel) ;

	/**
	* We assume 1.8 m telescopes on rails and 8 meters on the U stations.
	* @todo take the actual diameters from the DIAMETER column of ARRAY_GEOMETRY. This would
	*   allow simulations with masked apertures.
	*/
	aperture = ( *id.data() == 'U' ) ? 8.1 : 1.8 ;

	/* debugging
	* cerr << id << " aperture " << aperture << endl ;
	*/
}

/**
* Compute the baseline vector to the other station.
* Eq on page 10 of
* \latexonly
* \cite{MatharArxiv04}.
* \endlatexonly
* \htmlonly
* <a href="http://arxiv.org/abs/astrop-ph/0411384">astro-ph/0411384</a>.
* \endhtmlonly
* See also http://www.faqs.org/faqs/geography/infosystems-faq/ Q5.1
* @param[in] oth the station that is somewhere else than this
* @return a vector with three components [m]
*    the component [0] is points N, the component [1] West, and the component [2] zenith
*/
vector<double> Station::baseVec( const Station & oth) const
{
	/* debugging
	* cout << "Phi " << geolat << " " << oth.geolat << endl ;
	* cout << "Lam " << geolon << " " << oth.geolon << " diff " << geolon-oth.geolon << endl ;
	* cout << "geoR " << geoRadius << endl ;
	*/
	vector<double> b ;
	const double deltalambda = geolon- oth.geolon ;	// difference in longitudes; sign doesn't matter

	/** Because the differences in the geographic coordinates are small,
	* a Taylor expansion is usually more accurate than the straight implementation
	* of the of the simple trigonometric functions. The angle differences deltalambda and
	* deltaphi can be up to 200 m/6378 km =3.e-5 rad. So a 2nd order Taylor expansion has a
	* relative accuracy of 9.e-10, a 3rd order Taylor expansion of 27.e-15. At baselines of up to 200m
	* the 3rd order is accurate to 6 pm and sufficient for all our purposes.
	*/
#if 0
	b.push_back( -sin(geolat)*cos(oth.geolat)*cos(deltalambda)+cos(geolat)*sin(oth.geolat) ) ;
	b.push_back( cos(oth.geolat)*sin(deltalambda) ) ;
	b.push_back( cos(geolat)*cos(oth.geolat)*cos(deltalambda)+sin(geolat)*sin(oth.geolat)-1. ) ;
#else
	/* Maple source code:
	* x := -sin(phi1)*cos(phi2)*cos(dl)+cos(phi1)*sin(phi2) ;
	* z := cos(phi1)*cos(phi2)*cos(dl)+sin(phi1)*sin(phi2)-1 ;
	* phi1 := phi2+dp ;
	* mtaylor(x,{dp,dl}) ;
	* mtaylor(z,{dp,dl}) ;
	*/
	const double deltaphi = geolat- oth.geolat ;	// difference in latitude
	b.push_back( deltaphi*(deltaphi*deltaphi/6.-1.)+ 0.5*deltalambda*deltalambda*cos(oth.geolat)
			*(sin(oth.geolat)+deltaphi*cos(oth.geolat)) ) ;
	b.push_back( cos(oth.geolat)*sin(deltalambda) ) ;
	b.push_back( -0.5*deltaphi*deltaphi
			+0.5*deltalambda*deltalambda*cos(oth.geolat)*(sin(oth.geolat)*deltaphi-cos(oth.geolat)) ) ;
#endif
	for(int c=0 ; c < 3; c++)
		b[c] *= geoRadius ;

	return b ;
}

/**
* The baseline angle from this Station to the other, equivalent to
* the list in <a href="http://www.eso.org/observing/etc/doc/vlti/baseline/vltistations.html">vltistations.html</a>
* up to the sign.
* as discussed in the note
* \latexonly
* \cite{MatharMIDIRep16}.
* \endlatexonly
* \htmlonly
* <a href="http://www.strw.leidenuniv.nl/~mathar/public/matharMIDI20051110.pdf">MIDI optical path differences and phases</a>.
* \endhtmlonly
* @param[in] oth the station representing T2
* @return the azimuth angle from \c this to \c oth, from North over West [rad]
* @remark the azimuth angle convention differs from the one used for most of the rest of the code.
* @since 2006-02-26
*/
double Station::baseAzi( const Station & oth) const
{
	/** We define a creat circle from \c this to \c oth. The normal of this plane is
	* the cross product of the two vectors fromthe Earth center to the stations. In Maple,
	* 	# T1
	* 	r1[0] := cos(lambda[0])*cos(Phi[0]) ;
	* 	r1[1] := sin(lambda[0])*cos(Phi[0]) ;
	* 	r1[2] := sin(Phi[0]) ;
	* 
	* 	# T2
	* 	r2[0] := cos(lambda[1])*cos(Phi[1]) ;
	* 	r2[1] := sin(lambda[1])*cos(Phi[1]) ;
	* 	r2[2] := sin(Phi[1]) ;
	* 
	* 	# normal of common plane
	* 	n[0] := r1[1]*r2[2]-r1[2]*r2[1] ;
	* 	n[1] := r1[2]*r2[0]-r1[0]*r2[2] ;
	* 	n[2] := r1[0]*r2[1]-r1[1]*r2[0] ;
	* 
	* 	# North direction at T1 (derivative w.r.t Phi)
	* 	N1[0] := -cos(lambda[0])*sin(Phi[0]) ;
	* 	N1[1] := -sin(lambda[0])*sin(Phi[0]) ;
	* 	N1[2] := cos(Phi[0]) ;
	* 
	* 	# North direction at T1 (derivative w.r.t lambda, then normalized)
	* 	E1[0] := -sin(lambda[0])*cos(Phi[0]) ;
	* 	E1[1] := cos(lambda[0])*cos(Phi[0]) ;
	* 	E1[2] := 0 ;
	* 	E1[0] := -sin(lambda[0]) ;
	* 	E1[1] := cos(lambda[0]) ;
	* 	E1[2] := 0 ;
	* 
	* 	# unit vector into the direction to turn T1 around the axis
	* 	rot[0] := n[1]*r1[2]-n[2]*r1[1] ;
	* 	rot[1] := n[2]*r1[0]-n[0]*r1[2] ;
	* 	rot[2] := n[0]*r1[1]-n[1]*r1[0] ;
	* 
	* 	# decompose rot into a vector along the E and the N direction
	* 	coscompN := rot[0]*N1[0]+rot[1]*N1[1]+rot[2]*N1[2] ;
	* 	coscompE := rot[0]*E1[0]+rot[1]*E1[1]+rot[2]*E1[2] ;
	* 
	* 	coscompN := simplify(coscompN,trig) ;
	* 	coscompE := simplify(coscompE,trig) ;
	*/
	double compN = cos(geolat) *sin(oth.geolat)-sin(geolat)*cos(oth.geolat)*cos(geolon-oth.geolon) ;
	double compE = cos(oth.geolat)*sin(oth.geolon-geolon) ;
	return -atan2(compE,compN) ;
}

/**
* decipher whether this is a station on the western or eastern hemisphere from the VLTI M16
* @return true if the station U coordinates is less than approximately 53 m
*/
bool Station::isWest() const
{
	switch ( id.c_str()[0] )
	{
	case 'A' :
	case 'B' :
	case 'C' :
	case 'D' :
	case 'E' :
	case 'G' :
		return true;
		break ;
	case 'H' :
	case 'I' :
	case 'J' :
	case 'K' :
	case 'L' :
	case 'M' :
		return false;
		break ;
	case 'U' :
		switch( id.c_str()[1])
		{
		case '1' :
		case '2' :
			return true ;
			break ;
		case '3' :
		case '4' :
			return false ;
			break ;
		default :
			cerr << __FILE__ << " " << __LINE__ << " invalid station id " << id << endl ;
			exit(EXIT_FAILURE) ;
		}
		break ;
	default :
		cerr << __FILE__ << " " << __LINE__ << " invalid station id " << id << endl ;
		exit(EXIT_FAILURE) ;
	}
}

/**
* decipher whether this is a station North or South of the DL tunnel.
* @return true if the station V coordinates are more North than approximately -40 m
* @since 2006-08-28
*/
bool Station::isNorth() const
{
	switch ( id.c_str()[0] )
	{
	case 'A' :
	case 'B' :
	case 'C' :
	case 'D' :
	case 'E' :
	case 'H' :
	case 'I' :
	case 'K' :
	case 'L' :
	case 'M' :
		return false;
		break ;
	case 'U' :
		return true;
		break ;
	case 'G' :
		switch( id.c_str()[1])
		{
		case '0' :
		case '1' :
			return false ;
			break ;
		case '2' :
			return true ;
			break ;
		default :
			cerr << __FILE__ << " " << __LINE__ << " invalid station id " << id << endl ;
			exit(EXIT_FAILURE) ;
		}
		break ;
	case 'J' :
		switch( id.c_str()[1])
		{
		case '1' :
		case '2' :
			return false ;
			break ;
		case '3' :
		case '4' :
		case '5' :
		case '6' :
			return false ;
			break ;
		default :
			cerr << __FILE__ << " " << __LINE__ << " invalid station id " << id << endl ;
			exit(EXIT_FAILURE) ;
		}
		break ;
	default :
		cerr << __FILE__ << " " << __LINE__ << " invalid station id " << id << endl ;
		exit(EXIT_FAILURE) ;
	}
}

/**
* @return true if the telescope is a UT, false if it is an AT.
* @since 2006-08-28
*/
bool Station::isUT() const
{
	return ( id.c_str()[0] == 'U' ) ;
}

/**
* @param oth the station to compare with
* @return true if the two stations are on the same rail
* @since 2006-02-26
*/
bool Station::sameRail(const Station & oth) const
{
	if ( id.c_str()[0] == 'U' )
		return false ;
	else
		return id.c_str()[0] == oth.id.c_str()[0] ;
}

int Station::init( const string name, const double geoel)
{
	for(int s=0 ; s < statNo ; s++)
		if ( name == posN[s] )
		{
			*this = Station(posLo[s],posLa[s],name,geoel) ;
			/* debugging
			* cout << name << " station " << s << endl ;
			*/
			aperture = ( s >= 31 ) ? 8.0 : 1.8 ;
			return (staIndx= s) ;
		}
	return (staIndx =-1) ;
}

/** pupil area of the telescope entrance
* @return the entrance aperture [m^2]
* This is simply a quarter of \f$\pi\f$ multiplied by the diameter squared.
*/
double Station::area() const
{
	return 0.25*M_PI*aperture*aperture ;
}

/** compute the cartesian coordinates in a geocentric frame
* @param indx 
* @param xyz
*/
void Station::spher2cart(double xyz[3]) const
{
	xyz[0] = geoRadius*cos(geolon)*cos(geolat) ;
	xyz[1] = geoRadius*sin(geolon)*cos(geolat) ;
	xyz[2] = geoRadius*sin(geolat) ;
}

/** compute the cartesian coordinates in a geocentric frame
* @param[in] indx station index to generate the spherical angles
* @param[in] r sphere radius [m] 
* @param[out] xyz Cartesian coordinates [m]
*/
void Station::spher2cart(const int indx, const double r, double xyz[3])
{
	xyz[0] = r*cos(posLo[indx])*cos(posLa[indx]) ;
	xyz[1] = r*sin(posLo[indx])*cos(posLa[indx]) ;
	xyz[2] = r*sin(posLa[indx]) ;
}

/** Generate the OIFITS format on standard output.
* This contains a \c OI_ARRAY table as described in
* \latexonly
* \cite{PaulsSPIE5491} and \url{http://www.mrao.cam.ac.uk/~jsy1001/exchange/}.
* \endlatexonly
* \htmlonly
* <a href="http://www.mrao.cam.ac.uk/~jsy1001/exchange/">OIFITS</a>.
* \endhtmlonly
* @since 2006-05-17
* @author Richard J. Mathar
*/
FITS * Station::oifitsInit(const double geoRref, double uvwRef[3])
{
	FITS *fits = new FITS("-",Write) ;
	fits->pHDU().makeThisCurrent() ;
	fits->pHDU().addKey("CREATOR",__FILE__,"Program to create this file (R. J. Mathar)") ;
	fits->pHDU().writeDate() ;

	vector<string> colNams ;
	vector<string> colUns ;
	vector<string> colFmt ;

	colNams.push_back("TEL_NAME") ; colUns.push_back("") ; colFmt.push_back("16A") ;
	colNams.push_back("STA_NAME") ; colUns.push_back("") ; colFmt.push_back("16A") ;
	colNams.push_back("STA_INDEX") ; colUns.push_back("1") ; colFmt.push_back("I") ;
	colNams.push_back("DIAMETER") ; colUns.push_back("m") ; colFmt.push_back("E") ;
	colNams.push_back("STAXYZ") ; colUns.push_back("m") ; colFmt.push_back("3D") ;

	fits->addTable("OI_ARRAY",0,colNams,colFmt,colUns) ;
	ExtHDU & ext = fits->extension("OI_ARRAY") ;
	ext.addKey("OI_REVN",1,"Revision number of the table definition") ;
	ext.addKey("ARRNAME","VLTI","Telescope array") ;
	ext.addKey("FRAME","GEOCENTRIC","earth-centered Cartesian reference frame") ;

	spher2cart(0,geoRref,uvwRef) ;
	ext.addKey("ARRAYX",uvwRef[0],"[m] array center coordinate") ;
	ext.addKey("ARRAYY",uvwRef[1],"[m] array center coordinate") ;
	ext.addKey("ARRAYZ",uvwRef[2],"[m] array center coordinate") ;

	ext.writeComment("OIFITS http://www.mrao.cam.ac.uk/~jsy1001/exchange/") ;
	ext.addKey("CREATOR",__FILE__,"Program to create this file (R. J. Mathar)") ;
	ext.writeDate() ;

	return fits ;
}

/** Generate the OIFITS format on standard output.
* This contains a \c OI_ARRAY table as described in
* \latexonly
* \cite{PaulsSPIE5491} and \url{http://www.mrao.cam.ac.uk/~jsy1001/exchange/}.
* \endlatexonly
* \htmlonly
* <a href="http://www.mrao.cam.ac.uk/~jsy1001/exchange/">OIFITS</a>.
* \endhtmlonly
* @since 2006-05-17
* @author Richard J. Mathar
*/
void Station::oifits(FITS * fits, const double uvwRef[3]) const
{
	ExtHDU & ext = fits->extension("OI_ARRAY") ;

	ext.makeThisCurrent() ;

	long nrows= ext.rows() ;
	vector<string> asc(1) ;

	// asc[0] = at(s).target ;
	asc[0] = ( id[0] == 'U') ? string("UT")+id[1] : "AT" ;
	ext.column("TEL_NAME").write(asc,nrows+1) ;

	asc[0] = id ;
	ext.column("STA_NAME").write(asc,nrows+1) ;

	short tmpI = staIndx ;
	ext.column("STA_INDEX").write(& tmpI,1L,nrows+1) ;

	double tmpD = aperture ;
	ext.column("DIAMETER").write(& tmpD,1L,nrows+1) ;

	double xyz[3] ;
	spher2cart(xyz) ;
	for(int d=0; d < 3; d++)
		xyz[d] -= uvwRef[d] ;
	ext.column("STAXYZ").write(xyz,3L,1L,nrows+1) ;
}

ostream & operator <<(ostream &os, const Station & some)
{
	os << some.id ;
	return os ;
}

const string Station::posN[] = {"O", "A0", "A1", "B0", "B1", "B2", "B3", "B4", "B5", "C0", "C1", "C2", "C3", "D0",
		"D1", "D2", "E0", "G0", "G1", "G2", "H0", "I1", "J1", "J2", "J3", "J4", "J5", "J6", "K0", "L0", "M0", "U1", "U2",
		"U3", "U4" } ;
const double Station::posLo[] = {-1.228798831, -1.228801200, -1.228800303, -1.228799896, -1.228798998, -1.228798549, -1.228798100,
                 -1.228797651, -1.228797202, -1.228798591, -1.228797693, -1.228797244, -1.228796795, -1.228795981,
                 -1.228794186, -1.228793288, -1.228793372, -1.228790762, -1.228787171, -1.228792109, -1.228785543,
                 -1.228781994, -1.228780282, -1.228778935, -1.228784771, -1.228785668, -1.228787015, -1.228788362,
                 -1.228780324, -1.228779019, -1.228777714, -1.228800386, -1.228796107, -1.228790929, -1.228780856 } ;
const double Station::posLa[] = {-0.429829904, -0.429838653, -0.429841025, -0.429838245, -0.429840617, -0.429841803, -0.429842989,
                 -0.429844175, -0.429845361, -0.429837837, -0.429840209, -0.429841395, -0.429842581, -0.429837021,
                 -0.429841765, -0.429844137, -0.429836204, -0.429835388, -0.429844877, -0.429831830, -0.429833756,
                 -0.429839279, -0.429836090, -0.429839649, -0.429824230, -0.429821857, -0.429818299, -0.429814741,
                 -0.429832124, -0.429831716, -0.429831308, -0.429833092, -0.429825122, -0.429819523, -0.429823005 } ;

// __oOo__
#endif // STATION_C