/** @addtogroup tools
*/

#ifndef PRWDS_C
#define PRWDS_C

/**@{*/

/**
*
* @file
*
* @author Richard J. Mathar
* @version $Id: prWds.C,v 1.2 2006/03/06 17:47:29 mathar Exp $
*
*/

#include <math.h>
#include <iostream>
#include <iomanip>
#include <sstream>
#include <libgen.h>
#include <unistd.h>
#include <sys/types.h>
#include <regex.h>
#include <fstream>

// STL headers
#include <string>

#include "units.h"

using namespace std ;

/**
 * a right ascension
*/
class RA {
public:
	string wdsra ;	// 9 characters of the WDS entry
	double rad ;	// the value converted to radians

	RA(const char *wdsradenom ) : wdsra(wdsradenom)
	{
		initrad() ;
	}

	RA(const string & wdsradenom ) : wdsra(wdsradenom)
	{
		initrad() ;
	}

	RA(const double ra) : rad(ra)
	{
		if( rad < 0.)
			rad += 2.*M_PI ;
		initstr() ;
	}

	operator double() const
	{
		return rad ;
	}

protected:
private:
	/**
	* convert the string to radians
	*/
	void initrad()
	{
		int h= atoi(string(wdsra,0,2).c_str()) ;
		int m= atoi(string(wdsra,2,2).c_str()) ; 
		double s= atof(string(wdsra,4).c_str()) ; 
#ifdef DEBUG
		cout << wdsra << " is " << h << " " << m << " " << s << endl ;
#endif
		rad = (h+m/60. + s/3600.)*M_PI/12. ;
	}

	/**
	* and convert the radians to string
	*/
	void initstr()
	{
		char resul[11] ;	// including \n
		double hrs = rad ;
		/* radians to hours, where 24 hrs = 2 pi */
		hrs *= 12./M_PI ;
		int hh = (int)hrs ;
		hrs -= hh; 
		hrs *= 60 ;	// residual minutes
		int mm = (int)hrs ;
		hrs -= mm; 
		hrs *= 60 ;	// residual seconds
		sprintf(resul,"%02d%02d%06.3f",hh,mm,hrs) ;
		wdsra = resul ;
	}
} ;

/**
 * a declination
*/
class DEC {
public:
	string wdsdec ;	// the 9 characters of the WDS entry
	double rad ;	// the value converted to radians

	/**
	* ctor with a string
	* @param wdsdecdenom a string in the +-ddmmss.sss format
	*/
	DEC(const char *wdsdecdenom ) : wdsdec(wdsdecdenom)
	{
		initrad() ;
	}

	/**
	* ctor with a string
	* @param wdsdecdenom a string in the +-ddmmss.sss format
	*/
	DEC(const string & wdsdecdenom ) : wdsdec(wdsdecdenom)
	{
		initrad() ;
	}

	/**
	* ctor with a string
	* @param dec the declination [rad]
	*/
	DEC(const double dec) : rad(dec)
	{
		initstr() ;
	}

	/**
	* return an APES type of +-??d ??m ??.?? representation
	* @return the value [rad]
	*/
	operator double() const
	{
		return rad ;
	}
protected:
private:
	/**
	* convert the string to the representation of the number (rad)
	*/
	void initrad()
	{
		int d= atoi(string(wdsdec,1,2).c_str()) ; 	// excluding the sign
		int m= atoi(string(wdsdec,3,2).c_str()) ; 
		double s= atof(string(wdsdec,5).c_str()) ; 
#ifdef DEBUG
		cout << wdsdec << " is " << d << " " << m << " " << s << endl ;
#endif
		rad = DEG2RAD(d+m/60. + s/3600.) ;
		if ( (wdsdec.c_str())[0] == '-' )
			rad *= -1. ;
	}

	/**
	* and convert the radians to the WDS type of string
	*/
	void initstr()
	{
		char resul[12] ;	// sign, including plus \n
		/* radians to degrees, where 180 deg = pi */
		double deg = RAD2DEG(fabs(rad)) ;
		int dd = (int)deg ;
		deg -= dd; 
		deg *= 60 ;	// residual minutes
		int mm = (int)deg ;
		deg -= mm; 
		deg *= 60 ;	// residual seconds
		sprintf(resul,"%c%02d%02d%06.3f",(rad>=0.)? '+': '-', dd,mm,deg) ;
		wdsdec = resul ;
	}
} ;

/**
 * a position
*/
class RaDec {
public:
	RA ra;	// the right ascension component
	DEC dec;	// the declination component

	/**
	* ctor from a 18-character WDS 2000 coordinate
	* @param[in] coo 18 letters in a hhmmss.ss +-ddmmss.s format of the WDS
	*/
	RaDec(const string & coo) : ra( string(coo,0,9) ), dec( string(coo,9,9) )
	{
	}

	/**
	* ctor from two values
	* @param ra RA in radians
	* @param dec DEC in radians
	*/
	RaDec(const double rasc, const double decli) : ra(rasc), dec(decli)
	{
	}

	/**
	* move the coordinate by position angle theta and distance rho
	* @param theta position angle (deg)
	* @param rho position distance (arcsec)
	*/
	RaDec move(const string theta, const string rho)
	{
#ifdef DEBUG
		cout << "Separation " << rho << " as\n" ;
		cout << "Separation " << theta << " deg\n" ;
#endif
		double th = atof(theta.c_str()) ;	// degrees
		double r = atof(rho.c_str()) ;	// as
		// convert to radian
		th = DEG2RAD(th) ;
		r = ARCSEC2RAD(r) ;

		return move(th,r) ;
	}

protected:
private:
	/* Maple code to derive the final formula:
	with(linalg) ;
	assume(ra,real) ;
	assume(dec,real) ;
	assume(raa,real) ;
	assume(deca,real) ;
	assume(the,real) ;

	# primary star on (ra,dec) spherical coordinates
	ps := vector([cos(dec)*cos(ra),cos(dec)*sin(ra),sin(dec)]) ;
	# tangential vector to the North is derivative w.r.t. dec
	psN := vector([-sin(dec)*cos(ra),-sin(dec)*sin(ra),cos(dec)]) ;
	# tangential vector to the East pinned at the PS location
	psE := crossprod(psN,ps) ;

	# build a linear combination of the two components in the
	# tangential plane, rotated by the positional angle 'the=theta'
	# relative to the North. This is the direction 'ssder' where
	# the secondary star is found after a differential rotation.
	tmp1 := scalarmul(psN,cos(the)) ;
	tmp2 := scalarmul(psE,sin(the)) ;
	ssder := matadd(tmp1,tmp2) ;

	# the rotation axis is the cross product between the 
	# outwards PS direction vector and the direction towards the SS.
	ax := crossprod(ps,ssder) ;
	ax := map(simplify,ax) ;

	# test that the axis is indeed perpendicular to the ps.
	# dotprod(ps,ax) ;
	# simplify(%) ;

	# The plane in which the rotation takes place is spanned by
	# the direction to the PS and another vector perpendicular to
	# the rotation axis.
	psperp := crossprod(ax,ps) ;

	# In this plane we build a linear combination with the star
	# distance 'rho' to find the SS position.
	tmp1 := scalarmul(ps,cos(rho)) ;
	tmp2 := scalarmul(psperp,sin(rho)) ;
	ss := matadd(tmp1,tmp2) ;
	ss := map(simplify,ss) ;

	# and this is
	ss := [cos(rho) cos(dec~) cos(ra~) - sin(rho) cos(the~) sin(dec~) cos(ra~) - sin(rho) sin(the~) sin(ra~),

	    cos(rho) cos(dec~) sin(ra~) - sin(rho) cos(the~) sin(dec~) sin(ra~) + sin(rho) sin(the~) cos(ra~),

	    cos(rho) sin(dec~) + sin(rho) cos(dec~) cos(the~)]
	*/
	RaDec move(double th, double r)
	{
		double xyz[3] ;
		xyz[0] = cos(r)*cos(dec)*cos(ra)-sin(r)*cos(th)*sin(dec)*cos(ra)-sin(r)*sin(th)*sin(ra);
		xyz[1] = cos(r)*cos(dec)*sin(ra)-sin(r)*cos(th)*sin(dec)*sin(ra)+sin(r)*sin(th)*cos(ra);
		xyz[2] = cos(r)*sin(dec)+sin(r)*cos(dec)*cos(th);
		double newra = atan2(xyz[1],xyz[0]) ;
		double newdec = asin(xyz[2]) ;
		return RaDec(newra,newdec) ;
	}
} ;

/**
 * auxiliary comparison function to support sorted outputs.
 * @param left the sky position to the left of the operator
 * @param right the sky position to the right of the operator
*/
bool operator< (const RaDec & left, const RaDec &right)
{
	if ( left.ra < right.ra)
		return true ;
	else if ( left.ra > right.ra)
		return false ;
	else if ( left.dec < right.dec)
		return true ;
	else
		return false ;
}


/**
 * auxiliary comparison function to support sorted outputs.
 * @param left the sky position to the left of the operator
 * @param right the sky position to the right of the operator
*/
bool operator== (const RaDec & left, const RaDec &right)
{
	if ( left.ra == right.ra && left.dec == right.dec)
		return true ;
	else 
		return false ;
}

/**
* The class represents one line in the WDS star catalogue (one binary)
*/
class wdsstar : public string {
public:
	/**
	* Default constructor.
	*/
	wdsstar() : string("")
	{
	}

	/**
	* Constructor
	* @param wdslinebuf one line of the star catalog
	*/
	wdsstar(char *wdslinebuf) : string(wdslinebuf)
	{
	}

	/** returns the designator.
	* @return the first 10 characters of the 2000 coordinats.
	*/
	string name() const
	{
		return string(*this,0,10) ;
	}

	/** returns the Henry Draper Number, when known, else -1.
	* @return the first 10 characters of the 2000 coordinats.
	*/
	int hdn() const
	{
		/** The first three columns, 99-101 in the 1-based enumeration,
		* are not part of the number but the (signed) region indicator.
		* These are skipped, and the field of up to 5 decimals is actually read.
		*/
		const string h(*this,101,5) ;
		const int resul = atoi(h.c_str()) ;
		return ( resul != 0 ) ? resul : -1 ;
	}

	/**
	* @param preg precompiled regular expression for the star name
	* @return true if the name matches the regular expression
	*/
	bool nameMatches(const regex_t * preg) const
	{
		return ( regexec(preg, name().c_str(),0,0,0) ==0 ) ? true : false ;
	}

	/**
	* Converter to a RaDec coordinate pair.
	* This is in essence done by grabbing the J2000 coordinates at the end of the line.
	*/
	operator RaDec() const
	{
		/* grab the 2000 coordinates and use an implicit ctor with the string */

		return coo() ;
	}

	/*
	* Compute the proper motions in the order of Ra(PS), Dec(PS), Ra(SS), Dec(SS)
	* This still maintains the units of 15*[msec/yr]*cos(delta) for RA and
	* [mas/yr] for DEC.
	* @param Pmot on return the four proper motions for both stars, both coordinates.
	*/
	void pm(double Pmot[4]) const
	{
		/* in case they are all blank: use zeros */
		for(int i=0; i < 4; i++)
			Pmot[i]=0 ;
		istringstream nums0( pm(true,true).c_str() ) ; 
		nums0 >> Pmot[0] ;
		istringstream nums1( pm(false,true).c_str() ) ;
		nums1 >> Pmot[1] ;
		istringstream nums2( pm(true,false).c_str() ) ; 
		nums2 >> Pmot[2] ;
		istringstream nums3( pm(false,false).c_str() ) ;
		nums3 >> Pmot[3] ;

		// check whether this is per 100 yrs and needs a factor of 10
		if ( code().find('P') != string::npos)
			for(int i=0; i < 4; i++)
				Pmot[i] *= 10 ;
	}

	// the position angle of the last observation [deg]
	// In the following functions, star indices are C-type and one less than in http://ad.usno.navy.mil/wds/wdsweb_format.txt
	string posLast() const
	{
		return string(*this,42,3) ;
	}

	// the separation angle of the last observation [as]
	string sepLast() const
	{
		return string(*this,52,5) ;
	}

	/** Extract the K magnitude, if available.
	* @param primary if true, the first component's magnitude is returned, else the second's
	* @return the K magnitude. NaN if not available.
	* @since 2006-03-06
	*/
	double Kmag(const bool primary) const
	{
		double ret ;
		/** Initialize with NaN in case it is not found */
		memset(&ret,0xff,sizeof(double)) ;
		/** Inspect the flags field for the K specifier. If the flag indicates the K-magnitude is
		* present, use it; otherwise return NaN.
		*/
		const string notes = code() ;
		if ( notes.find("K") != string::npos)
			ret = mag(primary) ;
		return ret ;
	}

	/** Extract the V magnitude, if available.
	* @param primary if true, the first component's magnitude is returned, else the second's
	* @return the V magnitude. NaN if not available.
	* @since 2006-03-06
	*/
	double Vmag(const bool primary) const
	{
		double ret ;
		/** Initialize with NaN in case it is not found */
		memset(&ret,0xff,sizeof(double)) ;
		/** Inspect the flags field for the B, K and R specifier. If the flag indicates
		* this is the magnitude provided, return NaN; otherwise assume we have a V magnitude.
		*/
		const string notes = code() ;
		if ( notes.find_first_of("BKR") != string::npos)
			;
		else
			ret = mag(primary) ;
		return ret ;
	}

	/** Return the spectral type.
	* @param primary if true, the first component's spectral type is returned, else the second's
	* @return a spectral type in the standard K-to-O fashion
	*/ 
	string spectyp(const bool primary) const
	{
		string ret ;
		/** Get columns 71-79 (1-based) of the original WDS catalog line.
		*/
		const string spe = string(*this,70,9) ;

		/* debugging
		*	cerr << __FILE__ << " " << __LINE__ << " " << spe << endl ;
		*/

		/** Try to split this into the first and second component if there is a slash inside */
		string::size_type pos=spe.find("/")  ;
		if ( pos != string::npos)
		{
			/** There was a slash. For the primary, the spectral type is what is before the slash */
			if( primary)
				return string(spe,0,pos) ;
			else
			{
				/** For the secondary, one has to cancel some piece before the slash,
				* depending on what has been actually not replicated after it.
				*/
				switch( spe.at(pos+1) )
				{
				case 'A' :
				case 'B' :
				case 'C' :
				case 'D' :
				case 'F' :
				case 'G' :
				case 'K' :
				case 'M' :
				case 'N' :
				case 'O' :
				case 'P' :
				case 'R' :
				case 'S' :
				case 'W' :
				case 'a' :
				case 'd' :
				case 'e' :
				case 'g' :
				case 'p' :
				case 's' :
					return string(spe,pos+1) ;
					break ;
				case '1':
				case '2':
				case '3':
				case '4':
				case '5':
				case '6':
				case '7':
				case '8':
				case '9':
					return string(spe,0,1)+string(spe,pos+1) ;
					break ;
				default :
					return string(spe,0,2)+string(spe,pos+1) ;
				}
			}
		}
		else
		{
			/** If we have no slash, we may have a blank and another separate spec for
			* the secondary. To check this, we search for another letter of a star class, which
			* does not appear right at the beginning.
			*/
			pos=spe.find_first_of("ABCDFGKMNOPRSWadeghps",1)  ;
			if ( pos != string::npos)
			{
				/** If there was an explicit blank in the spec, we either return the
				* the first part before the blank for the primary, or the second part
				* for the secondary.
				*/
				if ( primary)
					return string(spe,0,pos) ;
				else
					return string(spe,pos) ;
			}
			else
				return spe ;
		}
	}

protected:
private:

	/**
	* Get a string representation of proper motions directly from the WDS line.
	* This is 15*pm(ra)*cos(delta) with pm(ra) measured in [msec/yr] for ra=true,
	* and pm(dec) in units of [mas/yr] if ra=false.
	* @param ra flags whether the ra or the dec value is requested
	* @param primary flags whether the primary or the secondary star is requested
	* @return the 4-letter string representation [mas/yr]
	*/
	string pm(const bool ra, bool primary) const
	{
		if( ra)
			if( primary)
				return string(*this,80,4) ;
			else
				return string(*this,89,4) ;
		else
			if( primary)
				return string(*this,84,4) ;
			else
				return string(*this,93,4) ;
	}

	/**
	* the 2000 coordinates
	* @return the 18 letters (9 for ra, 9 for dec) at the end of the line.
	*/
	string coo() const
	{
		return string(*this,112,18) ;
	}

	/**
	* the additional flags like N,B,C,D,I etc
	* @return the 4 letters that flag coding conventions
	*/
	string code() const
	{
		return string(*this,107,4) ;
	}

	/**
	* The number for the magnitude of the first or second component
	* @return the magnitude. This may be NaN if the field is not provided by the catalog.
	*/
	double mag(const bool primary) const
	{
		/** Extract either the columns 59-63 (1-based) for the first or 65-69 for the second component */
		string magString;
		if ( primary) 
			magString = string(*this,58,5) ;
		else
			magString = string(*this,64,5) ;
		/** Initialize with NaN in case there is no such information */
		double ret ;
		memset(&ret,0xff,sizeof(double)) ;
		sscanf(magString.c_str(),"%lf",&ret) ;

		/* debugging
		* cerr << __FILE__ << " " << __LINE__ << " " << magString << " " << ret << endl ;
		*/
		return ret ;
	}
} ;

struct posit_less {
	bool operator()(const wdsstar s1, const wdsstar s2) const
	{
		//const RaDec one(s1) ;
		//const RaDec two(s2) ;
		RaDec one = static_cast<wdsstar>(s1) ;
		RaDec two = static_cast<wdsstar>(s2) ;
		return ( one < two ) ;
	}
} ;

/**@}*/

#endif /* PRWDS_C */

/*
$Log: prWds.C,v $
Revision 1.2  2006/03/06 17:47:29  mathar
Added scanning of spectral type and magnitudes.

Revision 1.1  2006/03/05 20:07:39  mathar
First version of a coherent suite of CCfits based star catalogue functions.

Revision 1.2  2005/09/19 10:30:19  mathar
Added more support for more than one WDS catalog file, and for a regular
expression as the star designator to match a series of star names.

Revision 1.1  2005/09/18 20:29:32  mathar
First version of the program (which converts WDS entries to an APES catalog).

*/