/** @addtogroup tools
*/

/**@{*/

#ifndef PRSTARCAT_C
#define PRSTARCAT_C

/*
* UL - PRIMA project
*
* "@(#) $Id: prStarCat.C,v 1.8 2006/09/06 13:14:43 mathar Exp $"
*
*/

/** @file prStarCat.C
* The file \c prStarCat.C contains the classes, \c Simbad, \c Star
* and \c starCat. The class \c Star is a collection of astrometric
* properties of a single star with functions to convert these into
* FITS, APES, OIFITS and XML representations. The class \c starCat is a
* vector of these
* (a star catalog), and the class \c Simbad a helper class which scans
* an output of a query to Simbad to generate instances of \c Star.
*/

#include <math.h>
#include <time.h>
#include <iostream>
#include <iomanip>
#include <sstream>
#include <libgen.h>
#include <unistd.h>
#include <algorithm>

// STL headers
#include <string>
#include <map>
#include <vector>
#include <fstream>

#include <CCfits>

#include "slalib.h"
#include "units.h"
#include "prWds.C"

#ifdef HAVE_NOVAS
#include "novas.h"
#endif

using namespace std ;
using namespace CCfits ;

// maximum line length of the Simbad ASCII input (email) file
#define LL 256

/** Name of the binary FITS table with the star catalogue.
*/
#define PRSTARCAT_BINT "STARS"

/** If this macro is defined, RA and DEC columns are generated in
* the FITS files. The very disadvantage of doing so is that editing
* these FITS files manually (say, with \c fv) could easily create
* discrepancies to the COOR columns.
* #define PRSTARCAT_FITS_DEG
*/

/**
* If this macro is defined, we store internally pmra*cos(delta), else the bare 
* proper motion (time derivative of the coordinate, of alpha). The sober
* convention is to include the cosine factor because this is the more natural
* scaling. The conversion to XML, IPHASE and other slalib based formats must
* explicitly remove the cosine factor, whereas the factor is already absorbed
* into the native definitions of APES, WDS and NOVAS, for example.
*/
#define PRSTARCAT_PM_WITHD


/** One file generated by an automated response of the Simbad server.
* An instance contains the information from sending an e-mail with a request
* for a single star's information to the automated server.
* The methods are a set of ASCII scanners.
* @since 2005-08-29
* @author Richard J. Mathar
*/
class Simbad {

public:

	/** \c emaill is the series of lines in the standardized e-mail
	* response, typically startin with the \c From: line.
	*/
	vector<string> emaill ;

	/** This constructor creates an instance from an e-mail file.
	* @param[in] simbfname the UNIX file name of the e-mail.
	* These files ought be a simple ASCII copy of the e-mail response
	* send to  smbmail@simbad.u-strasbg.fr  following the template e-mail
	* body (of course with different e-mail and ID list):
	* <code>
	* SIMBAD interrogation
	* mailaddr: mathar@mpia.de
	* commands:
	* bibyears 1995/2003
	* maxdata
	* idlist:
	* Eta Carinae
	* </code>
	* The file name typically stars with a line
	* <code>
	* From smbmail@simbad.u-strasbg.fr  Mon Aug 29 16:37:32 2005
	* </code>
	* It hase a line like
	* <code>
	*        Your identifier .. is translated ..
	* </code>
	* and would finish as
	* <code>
	* Exit of SIMBADP
	* </code>
	* @author Richard J. Mathar
	* @since 2003-02-07 as mioSimbad2Fits.C in the CMM mio module of MIDI
	*/
	Simbad(char *simbfname)
	{
			// read e-mail file up to the first line after the separator with the dashes
		// cerr << " open " << simbfname << endl ;
		ifstream emfile(simbfname) ;
		if( !emfile)
		{
			cerr << " Cannot open " << simbfname << endl ;
			exit(EXIT_FAILURE) ;
		}
		bool found_eq=false ;	// already passed by the line that starts with the four '=' ? Scanning starts thereafter
		int found_delim =0 ;	// already passed by delimiter lines that are about 70 '-'s ?
		for(;;)
		{
			char buf[LL+1] ;	// e-mail lines are not longer than 80 chars by experience
			emfile.getline(buf,LL) ;
			if ( strncmp(buf,"====",4) == 0)	// wait until we find the line with (at least) 4 starting equals
				found_eq = true ;
			else if ( strncmp(buf,delim,strlen(delim)) == 0)	// delimiters: at least 30 dashes
			{
				emaill.push_back(delim) ;
				found_delim++ ;
			}
			else if ( found_eq )
			{
					// there might be broken lines that we re-wrap if they appear in the first section
				if( found_delim == 0 && strchr(buf,'=') == NULL) 	// missing equal sign indicates such a case
				{
					emaill[emaill.size()-1] += buf ;
				}
				else
					emaill.push_back(buf) ;
			}
			if( strncmp(buf,"Exit of SIMBADP",15) == 0 )
				break ;
			if( emfile.eof() )
				break ;
		}
		emfile.close() ;
	}

	/** Return the name of the object, having replaced blanks by underscores. The name will be drawn from the
	* first 16 characters of the first string line after the first delimiter, deleting all trailing blanks,
	* and replacing all remaining blanks by underscores.
	* @return a name like HD_10361 or SV*_ZI_942
	*/
	string object() const
	{
		// skip first header part
		int found_delim =0 ;	// already passed by delimiter lines that are about 70 '-'s ?
		vector<string>::const_iterator i=emaill.begin() ;
		for( ; i != emaill.end() ; i++)
		{
			if ( strcmp((*i).c_str(),delim) == 0 )	// delimiters: at least 30 dashes
			{
				found_delim++ ;
			}
			else if ( found_delim >= 1 )	// look at first line after the delimiter and take the first name from the left
			{
				string ret(*i,0,16) ; 		// only maximum of 16 characters available : VLT-SPE-ESO-15000-2764
					// delete trailing blanks
				// cerr << __FILE__ << " " << __LINE__ << " " << ret << "|" << endl ;
				string::size_type pos=ret.find_last_not_of(" ")  ;
				if( pos != string::npos)
					ret.erase(pos+1) ;
					// replace remaining 'in-name' blanks by underscores
#if 0
				// cerr << __FILE__ << " " << __LINE__ << " " << ret << endl ;
				while( (pos=ret.find_first_of(" ")) != string::npos)
					ret[pos] = '_' ;
#endif
				// cerr << __FILE__ << " " << __LINE__ << " " << ret << endl ;
				return ret;
			}
		}
		return string("") ;	// didn't find delimiter: nothing returned
	}

	/** Return the bibliographic reference.
	* The bibliographic reference is found in the first line, last string.
	* @return the reference code
	*/
	string biblio() const
	{
		string::size_type pos = emaill[0].find_last_of(" ") ;
		return string(emaill[0],pos) ;	// didn't find delimiter: nothing returned
	}

	/** Return the position by looking into the first line; return value is the equinox, RA and DEC with their errors
	* @param[out] ra right ascension in the format HH MM SS.sss
	* @param[out] dec declination in the format +-DD HH MM
	* @param[out] frame the reference system, ICRS or FK5 for example
	* @param[out] ra_err error in RA [mas]
	* @param[out] dec_err error in DEC [mas]
	* @return equinox [a]
	*/
	string radec(string & ra, string &dec, string & frame, double &ra_err, double &dec_err) const
	{
		// string & icrs = emaill[0] ;
		istringstream icrsl(emaill[0].c_str()) ;
		string skip ,ret ;
		icrsl >> frame >> ret ;
			// skip the equal sign
		icrsl >> skip ;
		ra.clear() ;
		icrsl >> skip ;	 //hours
		ra += skip ;
		ra += ":" ;
		icrsl >> skip ;	 // minutes
		ra += skip ;
		ra += ":" ;
		icrsl >> skip ;	 // seconds
		ra += skip ;

		dec.clear() ;
		icrsl >> skip ;	// hours
		dec += skip ;
		dec += ":" ;
		icrsl >> skip ;	 // minutes
		dec += skip ;
		dec += ":" ;
		icrsl >> skip ;	 // seconds
		dec += skip ;

		memset(&ra_err,0xff,sizeof(double)) ;
		memset(&dec_err,0xff,sizeof(double)) ;
		// find the error ellipse that is embedded in [...]
		string::size_type pmpos = emaill[0].find("[") ;
		ret = string(emaill[0],pmpos+1) ;
		istringstream err(ret.c_str()) ;
		//  evaluate the error ellipse along ra and dec
		double erra, errb , errang;
		memset(&erra,0xff,sizeof(double)) ;
		err >> erra >> errb >> errang;
		errang = DEG2RAD(errang) ;
		ra_err = hypot(sin(errang)*erra,cos(errang)*errb ) ;
		dec_err = hypot(cos(errang)*erra,sin(errang)*errb ) ;

		return ret ;
	}

	/*
	* @param[out] ra proper motion multiplied by cos(delta) [mas/a]
	* @param[out] dec proper motion [mas/a]
	* @param[out] ra_err error in proper motion multiplied by cos(delta) [mas/a]
	* @param[out] dec_err error in proper motion along declination [mas/a]
	* @since 2006-03-04 includes the errors on return
	*/
	void prop(double & ra, double &dec, double & ra_err, double & dec_err) const
	{
		// insert NaN's where information not available
		memset(&ra,0xff,sizeof(double)) ;
		memset(&dec,0xff,sizeof(double)) ;
		memset(&ra_err,0xff,sizeof(double)) ;
		memset(&dec_err,0xff,sizeof(double)) ;
		vector<string>::const_iterator i=emaill.begin() ;
		for( ; i != emaill.end() ; i++)
		{
			string::size_type pmpos ;
			if ( strcmp((*i).c_str(),delim) == 0 )	// delimiters: at least 30 dashes
				break ;
			else if ( (pmpos=(*i).find("pm =")) != string::npos )	// search for "pm =" in a line
			{
				string ret(*i,pmpos+4) ; // explicitly skip the "pm ="

				istringstream icrsl(ret.c_str()) ;
					// get proper motions in mas/a
				icrsl >> ra >> dec ;
					// find the error ellipse that is embedded in [...]
				pmpos=(*i).find("[") ;
				ret = string(*i,pmpos+1) ;
				istringstream err(ret.c_str()) ;
					//  evaluate the error ellipse along ra and dec
				double erra, errb , errang;
				memset(&erra,0xff,sizeof(double)) ;
				err >> erra >> errb >> errang;
				errang = DEG2RAD(errang) ;
				ra_err = hypot(sin(errang)*erra,cos(errang)*errb ) ;
				dec_err = hypot(cos(errang)*erra,sin(errang)*errb ) ;
				return ;
			}
		}
	}

	/**
	* read the K magnitude
	* @return the value of the magnitude in the K-band, or nan if not found.
	*/
	double Kmag() const
	{
		double k ;
		memset(&k,0xff,sizeof(double)) ;
		vector<string>::const_iterator i=emaill.begin() ;
		for( ; i != emaill.end() ; i++)
		{
			string::size_type pmpos ;
			if ( (pmpos=(*i).find("JP11  ")) == 0 )	// search for "JP11" in a line but skip the ones in the alias list
				break ;
		}
		if( i != emaill.end() )
			i += 5 ;	 // skip I: M: line with dashes 3 lines with JP11 and one more line with JP11
		if( i != emaill.end() )
		{
			string::size_type pmpos ;
			string ret(*i,39) ; // explicitly skip two vertical bars
			istringstream icrsl(ret.c_str()) ;
			icrsl >> k ;
		}
		if( isnan(k) )
		{
			if( i != emaill.end() )
				i += 3 ;	 // skip I: M: line with dashes 3 lines with JP11 and one more line with JP11
			if( i != emaill.end() )
			{
				string::size_type pmpos ;
				string ret(*i,39) ; // explicitly skip two vertical bars
				istringstream icrsl(ret.c_str()) ;
				icrsl >> k ;
			}
		}
		return k ;
	}

	/**
	* read the V magnitude
	* @return the value of the magnitude in the visible, or nan if not found.
	*/
	double Vmag() const
	{
		double v ;
		memset(&v,0xff,sizeof(double)) ;
		vector<string>::const_iterator i=emaill.begin() ;
		for( ; i != emaill.end() ; i++)
		{
			string::size_type pmpos ;
			if ( (pmpos=(*i).find("UBV  ")) == 0 )	// search for "UBV" in a line but skip the ones in the star alias list
				break ;
		}
		if( i != emaill.end() )
			i += 2 ;	 // skip I: M: line with dashes 3 lines with JP11 and one more line with JP11
		if( i != emaill.end() )
		{
			string::size_type pmpos ;
			string ret(*i,9) ; // explicitly skip one vertcial bar and mayby the folloing (E) or (J)
			istringstream icrsl(ret.c_str()) ;
			icrsl >> v ;
		}
		return v ;
	}

	/**
	* @param[out] palx the parallax in arcsec
	* @param[out] palx_err the parallax error in arcsec
	*/
	void plx(double & palx, double &palx_err) const
	{
		// insert NaN's where information not available
		memset(&palx,0xff,sizeof(double)) ;
		memset(&palx_err,0xff,sizeof(double)) ;
		vector<string>::const_iterator i=emaill.begin() ;
		for( ; i != emaill.end() ; i++)
		{
			string::size_type pmpos ;
			if ( strcmp((*i).c_str(),delim) == 0 )	// delimiters: at least 30 dashes
				break ;
			else if ( (pmpos=(*i).find("plx =")) != string::npos )	// search for "plx =" in a line
			{
				string ret(*i,pmpos+5) ; // explicitly skip the "plx ="

				istringstream icrsl(ret.c_str()) ;
					// get parallax in mas
				icrsl >> palx ;

					// find the error that is embedded in [...]
				pmpos=(*i).find("[") ;
				ret = string(*i,pmpos+1) ;
				istringstream err(ret.c_str()) ;
					//  evaluate the error ellipse along ra and dec
				err >> palx_err ;

					// convert them to arcsec
				palx /= 1000. ;
				palx_err /= 1000. ;
				return ;
			}
		}
	}

	/** Retrieve the spectral type
	* @return a string representing the spectral type of the object
	*/
	string spectr() const
	{
		vector<string>::const_iterator i=emaill.begin() ;
		for( ; i != emaill.end() ; i++)
		{
			string::size_type pmpos ;
			if ( strcmp((*i).c_str(),delim) == 0 )	// delimiters: at least 30 dashes
				break ;
			else if ( (pmpos=(*i).find("sp type =")) != string::npos )	// search for "sp type =" in a line
			{
				string ret(*i,pmpos+10) ; // explicitly skip the "sp type =" and the next blank
				return ret;
			}
		}
		return string("") ;
	}

	/** Retrieve various star catalogue numbers
	* @param[out] hip Hipparcos
	* @param[out] hd Henry Draper
	* @param[out] sao Smithsonian 
	* @param[out] wds Washington Double star
	* @return ADS
	* @see <a href="http://ad.usno.navy.mil/star/star_cats_rec.shtml">catalogues</a>
	*/
	string catalog(int & hip, int & hd, int & sao, string & wds) const
	{
		hip=hd=-1 ;
		vector<string>::const_iterator i=emaill.begin() ;
		string ret("") ;

		// scroll forward to the names
		for( ; i != emaill.end() ; i++)
		{
			string::size_type pmpos ;
			if ( strcmp((*i).c_str(),delim) == 0 )	// delimiters: at least 30 dashes
				break ;
		}
		if( i != emaill.end() )
			i ++ ;	// skip the delimiter line
		for( ; i != emaill.end() ; i++)
		{
			string::size_type pmpos ;
			if ( strcmp((*i).c_str(),delim) == 0 )	// delimiters: at least 30 dashes
				break ;
			if ( (pmpos=(*i).find("HD ")) != string::npos )	// search for Henry Draper number
			{
				string id(*i,pmpos+3) ; // explicitly skip the "HD" and the next blank
				hd = atoi( id.c_str()) ;
			}
			if ( (pmpos=(*i).find("HIP ")) != string::npos )	// search for Hipparcos number
			{
				string id(*i,pmpos+4) ; // explicitly skip the "HIP" and the next blank
				hip = atoi( id.c_str()) ;
			}
			if ( (pmpos=(*i).find("SAO ")) != string::npos )	// search for Hipparcos number
			{
				string id(*i,pmpos+4) ; // explicitly skip the "SAO" and the next blank
				sao = atoi( id.c_str()) ;
			}
			if ( (pmpos=(*i).find("WDS ")) != string::npos )	// Washington Double Star des
			{
				string id(*i,pmpos+4) ; // explicitly skip the "WDS" and the next blank
				string::size_type idend = id.find(" ") ;	// up to the end of string
				if ( idend != string::npos)
					id.resize(idend) ;
				wds = id ;
			}
			if ( (pmpos=(*i).find("ADS ")) != string::npos )	// Aitken Double S designator
			{
				string id(*i,pmpos+4) ; // explicitly skip the "ADS" and the next blank
				string::size_type idend = id.find("  ") ;	// try to keep the ABC etc component id if present
				if ( idend != string::npos)
					id.resize(idend) ;
				ret = id ;
			}
		}
		return ret ;
	}

	/** Estimate the effective temperature.
	* http://ads.harvard.edu/cgi-bin/bbrowse?book=hsaa&chap=2&page=0068
	* @param[out] V_K the magnitude excess V minus K
	* @return the effective black body temperature [K]
	*/
	double Teff(double & V_K) const
	{
		string styp = spectr() ;
		double T=6000. ;
		switch( *styp.c_str() )
		{
		case 'O':
			//V_J=-0.73 ;
			V_K=-0.94 ;
			switch( *(styp.c_str()+1) )
			{
			case '5':
			case '6':
			case '7':
				T=38000. ;
				break ;
			case '8':
			case '9':
				T=35000. ;
				break ;
			}
			break;
		case 'B':
			switch( *(styp.c_str()+1) )
			{
			case '0':
				T=30000. ;
				//V_J= -0.7;
				V_K= -0.93;
				break ;
			case '1':
				T=24200. ;
				//V_J= -0.61 ;
				V_K= -0.88 ;
				break ;
			case '2':
				T=22100. ;
				//V_J= -0.55 ;
				V_K= -0.74 ;
				break ;
			case '3':
				T=18800. ;
				//V_J= -0.45 ;
				V_K= -0.61 ;
				break ;
			case '4':
				T=17408. ;	// log-log scale interpolation
	/* Maple V
	xa := 3 ;
	xb := 5 ;
	Ta := 18800 ;
	Tb := 16400 ;
	xc := 4 ;
	# T=a*x**b; log T = log a+b*log x
	logTa := log(Ta) ;
	logTb := log(Tb) ;
	logxa := log(xa) ;
	logxb := log(xb) ;
	b := (logTb-logTa)/(logxb-logxa) ;
	loga := logTa-b*logxa ;
	a := exp(loga) ;
	evalf(a*xc^b) ;
	quit
	*/
				break ;
			case '5':
				T=16400. ;
				//V_J= -0.35 ;
				V_K= -0.47 ;
				break ;
			case '6':
				T=15400. ;
				//V_J= -0.30 ;
				V_K= -0.41 ;
				break ;
			case '7':
				T=14500. ;
				//V_J= -0.25 ;
				V_K= -0.35 ;
				break ;
			case '8':
				T=13400. ;
				//V_J= -0.17 ;
				V_K= -0.24 ;
				break ;
			case '9':
				T=12400. ;
				//V_J= -0.09 ;
				V_K= -0.14 ;
				break ;
			}
			break;
		case 'A':
			switch( *(styp.c_str()+1) )
			{
			case '0':
				T=10800. ;
				//V_J= -0.01 ;
				V_K= -0.03 ;
				break ;
			case '2':
				T=9730. ;
				//V_J= 0.11 ;
				V_K= +0.13 ;
				break ;
			case '3':
				T=9222. ;
				break ;
			case '5':
				T=8620. ;
				//V_J= 0.27 ;
				V_K= +0.36 ;
				break ;
			case '7':
				T=8190. ;
				//V_J= 0.35 ;
				V_K= +0.46 ;
				break ;
			}
			break;
		case 'F':
			switch( *(styp.c_str()+1) )
			{
			case '0':
				T=7240. ;
				//V_J= 0.58 ;
				V_K= +0.79 ;
				break ;
			case '2':
				T=6930. ;
				//V_J= 0.68 ;
				V_K= +0.93 ;
				break ;
			case '5':
				T=6540. ;
				//V_J= 0.79 ;
				V_K= +1.07 ;
				break ;
			case '7':
				T=6295. ;
				break ;
			case '8':
				T=6200. ;
				//V_J= 0.96 ;
				V_K= +1.27 ;
				break ;
			}
			break;
		case 'G':
			switch( *(styp.c_str()+1) )
			{
			case '0':
				T=5920. ;
				//V_J= 1.03 ;
				V_K= +1.35 ;
				break ;
			case '2':
				T=5780. ;
				//V_J= 1.10 ;
				V_K= +1.44 ;
				break ;
			case '5':
				T=5610. ;
				//V_J= 1.14 ;
				V_K= +1.49 ;
				break ;
			case '8':
				T=5490. ;
				//V_J= 1.24 ;
				V_K= +1.63 ;
				break ;
			case '9':
				T=5460. ;	// interpolated
				break ;
			}
			break;
		case 'K':
			switch( *(styp.c_str()+1) )
			{
			case '0':
				T=5240. ;
				//V_J= 1.38 ;
				V_K= +1.83 ;
				break ;
			case '2':
				T=4780. ;
				//V_J= 1.57 ;
				V_K= +2.15 ;
				break ;
			case '3':
				T=4656. ;
				break ;
			case '4':
				T=4571. ;
				break ;
			case '5':
				T=4506. ;
				//V_J= 2.04 ;
				V_K= +2.75 ;
				break ;
			case '6':
				T=4453. ;
				break ;
			case '7':
				T=4410. ;
				//V_J= 2.36 ;
				V_K= +3.21 ;
				break ;
			case '8':
				T=4372. ;
				break ;
			}
			break;
		case 'M':
			switch( *(styp.c_str()+1) )
			{
			case '0':
				T=3920. ;
				//V_J= 2.71 ;
				V_K= +3.60 ;
				break ;
			case '1':
				T=3680. ;
				//V_J= 3.06 ;
				V_K= +3.95 ;
				break ;
			case '2':
				T=3500. ;
				//V_J= 3.37 ;
				V_K= +4.27 ;
				break ;
			case '3':
				T=3360. ;
				//V_J= 3.66 ;
				V_K= +4.57 ;
				break ;
			case '4':
				T=3230. ;
				//V_J= 3.97 ;
				V_K= +4.87 ;
				break ;
			case '5':
				T=3120. ;
				//V_J= 4.28 ;
				V_K= +5.17 ;
				break ;
			case '6':
				//V_J= 4.63 ;
				V_K= +5.58 ;
				break ;
			case '7':
				//V_J= 5.20 ;
				V_K= +6.18 ;
				break ;
			case '8':
				T=2660. ;
				//V_J= 5.80 ;
				V_K= +6.75 ;
				break ;
			}
			break;
		}
		return T ;
	}

	/** Retrieve the radial velocity.
	* @param[out] vel the radial velocity [km/s]
	* @return the string BARYCENT for the reference center.
	*/
	string veloc(double &vel) const
	{
		memset(&vel,0xff,sizeof(double)) ;
		vector<string>::const_iterator i=emaill.begin() ;
		for( ; i != emaill.end() ; i++)
		{
			string::size_type pmpos ;
			if ( strcmp((*i).c_str(),delim) == 0 )	// delimiters: at least 30 dashes
				break ;
			else if ( (pmpos=(*i).find("rv =")) != string::npos )	// search for "rv =" in a line
			{
				string ret(*i,pmpos+5) ; // explicitly skip the "rv =" and the next blank
				istringstream icrsl(ret.c_str()) ;
				char v_or_z ;
				icrsl >> v_or_z >> vel ;
					// convert to velocity if redshift give http://iram.fr/ARN/may95/node4.html
				if( v_or_z == 'z' )
					vel *= 299792.458 ;
				return "BARYCENT" ;
			}
		}
		return string("") ;
	}
protected:
private:
	const static char delim[] ;	// delimiter in the e-mail fields

} ; /* class Simbad */
const char Simbad::delim[]="------------------------------" ;
#undef LL

#define CHARM2_LL 471
/** The CHARM2 catalog of ftp://cdsarc.u-strasbg.fr/pub/cats/J/A+A/431/773/
* @see <a href="http://www.eso.org/~arichich/">home page A Richichi</a>
* @since 2007-02-15
* @author Richard J. Mathar
*/
class Charm2 {

public:
	/* internal charm2 index */
	vector<int> N ;

	vector<string> ids ;

	/* the right ascensions in degrees */
	vector<double> RA ;
	/* the declinations in degrees */
	vector<double> DE ;

	/* proper motion in RA multiplied by cos(delta) [as/yr] */
	vector<double> pmRA ;

	/* proper motion in DEC [as/yr] */
	vector<double> pmDE ;

	/* magnitude in V */
	vector<double> Vmag ;

	/* magnitude in K */
	vector<double> Kmag ;

	/* spectral class */
	vector<string> SpType ;

	/* parallax [mas] */
	vector<double> Plx ;

	/* radial velocity [km/s] */
	vector<double> RV ;

	/*
	* @author Richard J. Mathar
	* @since 2007-02-15
	*/
	Charm2(char *fname)
	{
		ifstream emfile(fname) ;
		if( !emfile)
		{
			cerr << " Cannot open " << fname << endl ;
			exit(EXIT_FAILURE) ;
		}
		for(;;)
		{
			char buf[CHARM2_LL+1] ;	// lines of the length documented in the ReadMe file
			emfile.getline(buf,CHARM2_LL) ;
			int seq ;
			char decsign,Spt[14],id1to4[100] ;
			double ra[3],dec[3],pm[2],mag[2],plx,rv ;

			Spt[13]='\0' ;
			id1to4[99]='\0' ;
			/* get N, skip bytes 5-30, get ids, skip bytes 130-332, get RAh, RAm, RAs, DE-, DEd, DEm, DeS, pmRA, pmDE
			* skip bytes 377-383, get Vmag, Kmag, skip bytes 395-407, get SpType, skip 421-425, get Plx, RV
			*/
			sscanf(buf,"%4d",&seq) ;
			sscanf(buf+30,"%99c",id1to4) ;
			sscanf(buf+332,"%2le",&ra[0]) ;
			sscanf(buf+335,"%2le",&ra[1]) ;
			sscanf(buf+338,"%6le",&ra[2]) ;
			sscanf(buf+345,"%1c",&decsign) ;
			sscanf(buf+346,"%2le",&dec[0]) ;
			sscanf(buf+349,"%2le",&dec[1]) ;
			sscanf(buf+352,"%5le",&dec[2]) ;
			sscanf(buf+358,"%8le",&pm[0]) ;
			sscanf(buf+367,"%9le",&pm[1]) ;
			sscanf(buf+383,"%5le",&mag[0]) ;
			sscanf(buf+389,"%5le",&mag[1]) ;
			sscanf(buf+407,"%13c",Spt) ;
			sscanf(buf+425,"%6le",&plx) ;
			sscanf(buf+432,"%6le",&rv) ;

			/** debugging
			cerr << __FILE__ << " " << __LINE__ << buf << endl ;
			cerr << __FILE__ << " " << __LINE__ << " seq " << seq << endl ;
			cerr << __FILE__ << " " << __LINE__ << " id " << id1to4 << endl ;
			cerr << __FILE__ << " " << __LINE__ << " rah " << ra[0] << endl ;
			cerr << __FILE__ << " " << __LINE__ << " ram " << ra[1] << endl ;
			cerr << __FILE__ << " " << __LINE__ << " ras " << ra[2] << endl ;
			cerr << __FILE__ << " " << __LINE__ << " decsign " << decsign << endl ;
			cerr << __FILE__ << " " << __LINE__ << " dech " << dec[0] << endl ;
			cerr << __FILE__ << " " << __LINE__ << " decm " << dec[1] << endl ;
			cerr << __FILE__ << " " << __LINE__ << " decs " << dec[2] << endl ;
			cerr << __FILE__ << " " << __LINE__ << " pmra " << pm[0] << endl ;
			cerr << __FILE__ << " " << __LINE__ << " pmdec " << pm[1] << endl ;
			cerr << __FILE__ << " " << __LINE__ << " mags " << mag[0] << " " << mag[1] << endl ;
			cerr << __FILE__ << " " << __LINE__ << " " << seq << " Spt|" << Spt << "|" << endl ;
			cerr << __FILE__ << " " << __LINE__ << " plx " << plx << endl ;
			cerr << __FILE__ << " " << __LINE__ << " rv " << rv << endl ;
			*/

			N.push_back(seq) ;
			RA.push_back(15.*(ra[0]+ra[1]/60.+ra[2]/3600.)) ;
			if ( decsign == '-')
				DE.push_back(-dec[0]-dec[1]/60.-dec[2]/3600.) ;
			else
				DE.push_back(dec[0]+dec[1]/60.+dec[2]/3600.) ;
			pmRA.push_back(pm[0]) ;
			pmDE.push_back(pm[1]) ;
			Vmag.push_back(mag[0]) ;
			Kmag.push_back(mag[1]) ;
			SpType.push_back(Spt) ;
			Plx.push_back(plx) ;
			RV.push_back(rv) ;
			if( emfile.eof() )
				break ;
		}
		emfile.close() ;
	}

protected:
private:

} ; /* class Charm2 */
#undef CHARM2_LL

/** The properties of a single star.
* @since 2006-03-01
* @author Richard J. Mathar
* @todo write some kind of XERCES based XML scanner (as ctor)
*   This would be the inverse functionality of the xml() member function.
* @todo write some kind of APES input scanner (as ctor)
*   This would be the inverse functionality of the apes() member function.
*/
class Star {
public:
	double equinox ;

	/** Frame, ICRS or FK5 */
	string frame ;

	/** RA in degrees. For computational reasons, one would prefer [rad], but since
	* this here is only a sort of catalogue, we stay with the simpler [deg] to simplify
	* the I/O.
	*/
	double ra, ra_err ;

	/** DEC in degrees */
	double dec , dec_err ;

	/** proper motion in alpha, [mas/a]
#ifdef PRSTARCAT_PM_WITHD
	* This includes a factor cos(delta) where delta is the declination.
#endif
	*/
	double pmra , pmra_err ;

	/** proper motion in delta, [mas/a]
	*/
	double pmdec , pmdec_err ;

	/** parallax [mas] */
	double px, px_err ;

	/** radial velocity [km/s] 
	*/
	double rv ;

	/** radial velocity type. Strings like "LSR", "GEOCENTR".
	*/
	string veltyp ;

	/** radial velocity definition. Strings like "OPTICAL", "RADIO".
	*/
	string veldef ;

	/** effective temperature [K] */
	double Teff ;

	/** magnitude in V */
	double VMag ;

	/** magnitude in K */
	double KMag ;

	/** spectral type */
	string styp ;

	/** orbital period [d] */
	double P ;

	/* orbital semi-major axis [as] */
	double a ;

	/** orbital inclination [deg] */
	double inclin ;

	/* Longitude of the ascending Node [deg] */
	double Omega ;

	/* argument of the perihelion [deg] */
	double omega ;

	/* time of the perihelion passage MJD [d] */
	double T0 ;

	/** eccentricity [1] */
	double eccen ;

	/** usage flag (T or R) */
	string flag ;

	/** Henry Draper catalog number */
	int hdn ;
	/** Hipparcos catalog number */
	int hip ;
	/** FK5 or FK6 catalog number */
	int fk6 ;
	/** SAO catalogue number */
	int sao ;
	/** local PRIMA catalog number */
	int indx ;

	/** WDS identifier */
	string wds ;

	/** Aitken Double Star catalogue */
	string ads ;

	/** default ctor. Nothing known. All parameters set to zero.
	*/
	Star() : hdn(-1), hip(-1), fk6(-1), sao(-1), indx(-1)
	{
	}

	/** Constructor from a Simbad e-mail.
	* @param simb the collection of Simbad -email lines for a single star.
	* @param duplic true if a pseudo automatic classification of target and reference stars is attempted.
	* 	An alternating classification as target and reference is derived from
	*	using classifying each second star as target.
	*/
	Star(const Simbad & simb, bool duplic=false) : hdn(-1), hip(-1), fk6(-1), sao(-1), indx(-1), px_err(0.)
	{

		/** Construct coordinates (RA, DEC, reference frame) . Convert RA and DEC errors
		* from mas to degrees.
		*/
		string raepp, decepp ;
		string nox = simb.radec(raepp,decepp,frame,ra_err,dec_err) ;
		ra_err= MAS2DEG(ra_err) ;
		dec_err= MAS2DEG(dec_err) ;

		int raHr, raMin, decDeg, decMin ;
		double raSec, decSec ;
		sscanf(raepp.c_str(),"%d:%d:%lf",&raHr, &raMin, &raSec) ;
		sscanf(decepp.c_str(),"%d:%d:%lf",&decDeg, &decMin, &decSec) ;

		ra = 15.*(raHr + raMin/60. + raSec/3600.) ;
		dec = abs(decDeg) + decMin/60. + decSec/3600. ;
		if ( decDeg < 0 )
			dec *= -1. ;

		/** construct proper motions and remove the factor cos(delta) from the ra component
		* if the internal definition is different.
		*/
		simb.prop(pmra,pmdec,pmra_err,pmdec_err) ;
#ifndef PRSTARCAT_PM_WITHD
		pmra /= cos(DEG2RAD(dec)) ;
		pmra_err /= cos(DEG2RAD(dec)) ;
#endif


		/** Parallax. Convert from [as] to [mas].
		*/
		simb.plx(px,px_err) ;
		px *= 1.e3 ;
		px_err *= 1.e3 ;

		/** radial velocity.
		*/
		veltyp=simb.veloc(rv) ;

		/** spectral type */
		styp = simb.spectr() ;

		/** K and visible magnitude. If the K-band magnitude is not given
		* by Simbad directly, we try to accesss the V-K excess and extrapolate
		* the K-band magintude starting from the V
		*/
		KMag = simb.Kmag() ;
		VMag = simb.Vmag() ;
		double V_K ;	// V-K excess
		Teff = simb.Teff(V_K) ;
		if ( isnan(KMag) )
		{
			if( !isnan(VMag) )
				KMag = VMag -V_K ;
		}

		equinox = 2000. ;

		/** catalogue name and index assignments */
		ads = simb.catalog(hip,hdn,sao,wds) ;

		if ( duplic)
			flag = ( no % 2 == 0) ? "T" : "R" ;

		target = simb.object() ;
		if ( duplic)
		{
			target += " " ;
			target += char('A'+ no %2) ;
		}

		no++ ;
	}

	/** Constructor from star of the CHARM2 catalog
	* @param charm2 the collection of CHARM2 lines
	* @param c2idx the index of the catalog. This is not the index assigend
	* 	within the catalog but the sequential number, starting at 0, in the Charm2 vector.
	* @since 2007-02-15
	* @author Richard J. Mathar
	* @todo convert the Charm::ids into catalog cross-references (hdn, sao, ads,...)
	*/
	Star(const Charm2 & charm2, const int c2idx) : equinox(2000.), frame("FK5"), veltyp(""), veldef(""), flag("")
	{
		ra = charm2.RA[c2idx] ;
		memset(&ra_err,0xff,sizeof(double)) ;

		dec = charm2.DE[c2idx] ;
		memset(&dec_err,0xff,sizeof(double)) ;

#ifdef PRSTARCAT_PM_WITHD
		pmra = 1.e3*charm2.pmRA[c2idx] ;
#else
		pmra = 1.e3*charm2.pmRA[c2idx]/cos(DEG2RAD(dec)) ;
#endif
		memset(&pmra_err,0xff,sizeof(double)) ;

		pmdec = 1.e3*charm2.pmDE[c2idx] ;
		memset(&pmdec_err,0xff,sizeof(double)) ;

		px = charm2.Plx[c2idx] ;
		memset(&px_err,0xff,sizeof(double)) ;

		rv = charm2.RV[c2idx] ;

		memset(&Teff,0xff,sizeof(double)) ;

		VMag = charm2.Vmag[c2idx] ;
		KMag = charm2.Kmag[c2idx] ;

		styp = charm2.SpType[c2idx] ;

		memset(&P,0xff,sizeof(double)) ;
		memset(&a,0xff,sizeof(double)) ;
		memset(&inclin,0xff,sizeof(double)) ;
		memset(&Omega,0xff,sizeof(double)) ;
		memset(&omega,0xff,sizeof(double)) ;
		memset(&T0,0xff,sizeof(double)) ;
		memset(&eccen,0xff,sizeof(double)) ;

		indx= charm2.N[c2idx] ;
#if 0
		/** Henry Draper catalog number */
		int hdn ;

		/** SAO catalogue number */
		int sao ;

		/** Aitken Double Star catalogue */
		string ads ;
#endif
	}

	/** Constructor from row in the standard FITS representation.
	* @param e the HDU with the star collection
	* @param row the row in this HDU with the particular star, 1-based
	*/
	Star(const ExtHDU & ext, long row)
	{
		ext.makeThisCurrent() ;
		valarray<double> dinp(1) ;
		valarray<float> finp(1) ;
		valarray<long> iinp(1) ;
		//valarray<vec> sinp(1) ;
		vector<string> sinp(1) ;

		/** switch to 1-based indexing */
		row++ ;

		ext.column("COOR").read(sinp,row,row) ;
		coordsInv(sinp[0]) ;

		ext.column("RA_ERR").read(finp,row,row) ;
		ra_err = finp[0] ;
		ext.column("DEC_ERR").read(finp,row,row) ;
		dec_err = finp[0] ;

		ext.column("EQUINOX").read(finp,row,row) ;
		equinox = finp[0] ;

		ext.column("RADECSYS").read(sinp,row,row) ;
		frame = sinp[0] ;

		ext.column("PMRA").read(finp,row,row) ;
		pmra = finp[0] ;

		ext.column("PMRA_ERR").read(finp,row,row) ;
		pmra_err = finp[0] ;

		ext.column("PMDEC").read(finp,row,row) ;
		pmdec = finp[0] ;

		ext.column("PMDEC_ERR").read(finp,row,row) ;
		pmdec_err = finp[0] ;

		ext.column("PX").read(finp,row,row) ;
		px = finp[0] ;
		ext.column("PX_ERR").read(finp,row,row) ;
		px_err = finp[0] ;

		ext.column("RV").read(finp,row,row) ;
		rv = finp[0] ;
		ext.column("RVTYPE").read(sinp,row,row) ;
		veltyp = sinp[0] ;
		ext.column("RVDEF").read(sinp,row,row) ;
		veldef = sinp[0] ;

		ext.column("TEFF").read(finp,row,row) ;
		Teff = finp[0] ;

		ext.column("VMAG").read(finp,row,row) ;
		VMag = finp[0] ;

		ext.column("KMAG").read(finp,row,row) ;
		KMag = finp[0] ;

		ext.column("STYP").read(sinp,row,row) ;
		styp = sinp[0] ;

		ext.column("P").read(finp,row,row) ;
		P = finp[0] ;

		ext.column("A").read(finp,row,row) ;
		a = finp[0] ;

		ext.column("I").read(finp,row,row) ;
		inclin = finp[0] ;

		ext.column("OMEGG").read(finp,row,row) ;
		Omega = finp[0] ;

		ext.column("T0").read(dinp,row,row) ;
		T0 = dinp[0] ;

		ext.column("OMEG").read(finp,row,row) ;
		omega = finp[0] ;

		ext.column("E").read(finp,row,row) ;
		eccen = finp[0] ;

		ext.column("BFLAG").read(sinp,row,row) ;
		flag = sinp[0] ;

		ext.column("HDN").read(iinp,row,row) ;
		hdn = iinp[0] ;
		ext.column("HIC").read(iinp,row,row) ;
		hip = iinp[0] ;
		ext.column("FK6").read(iinp,row,row) ;
		fk6 = iinp[0] ;

		ext.column("ADS").read(sinp,row,row) ;
		ads = sinp[0] ;

#ifdef PRSTARCAT_FITS_DEG
#error not yet implemented
#endif

		ext.column("SAO").read(iinp,row,row) ;
		sao = iinp[0] ;

		ext.column("WDS").read(sinp,row,row) ;
		wds = sinp[0] ;

		ext.column("INDX").read(iinp,row,row) ;
		indx = iinp[0] ;

		/** For the name, we use a WDS type of name augmented by the flag */
		target = wdsName()+flag ;
	}

	/** Constructor from a WDS line.
	* @param[in] w the WDS star to be used as a source of information
	* @param[in] useprim indicates whether we generate a star from the implicit first or implicit second entry
	* @todo complete the transformation of the other information contained in the WDS lines
	* @since 2006-03-05
	* @author Richard J. Mathar
	*/
	Star (const wdsstar w, bool useprim) : equinox(2000.), frame("FK5"),
		px(0.), rv(0.), veltyp(""), veldef(""), Teff(0.), VMag(0.), KMag(0.), styp(""), P(0.),
		a(0.), inclin(0.), Omega(0.), omega(0.), T0(0.), eccen(0.), hip(-1), fk6(-1), sao(-1), indx(0), ads(""), target("")
	{
		/* debugging
		* cerr << __FILE__ << " " << __LINE__ << " " << w << endl ;
		*/
		wds = w.name() ;
		memset(&ra_err,0xff,sizeof(double)) ;
		memset(&dec_err,0xff,sizeof(double)) ;
		memset(&pmra_err,0xff,sizeof(double)) ;
		memset(&pmdec_err,0xff,sizeof(double)) ;
		memset(&px_err,0xff,sizeof(double)) ;
		/** Assume that the first WDS component is PRIMA's reference, the second the target */
		flag = ( useprim) ? "R" : "T" ;

		/** Access the coordinates in RA and DEC */
		RaDec wradec=w ;

		/* debugging
		*cerr << __FILE__ << " " << __LINE__ << " " << wradec.ra.wdsra << endl ;
		*cerr << __FILE__ << " " << __LINE__ << " " << wradec.dec.wdsdec << endl ;
		*/

		/** If this is the second component, move a little bit further
		* according to the two components of the position vector mentioned.
		*/
		if ( !useprim)
			wradec = wradec.move(w.posLast(),w.sepLast()) ;

		/** convert the WDS coordinates to degrees */
		ra=RAD2DEG(wradec.ra.rad) ;
		dec=RAD2DEG(wradec.dec.rad) ;

		/** Get the proper motions individually for the two components */
		double wPmot[4] ;
		w.pm(wPmot) ;
		if ( useprim)
		{
			pmra = wPmot[0] ;
			pmdec = wPmot[1] ;
		}
		else
		{
			pmra = wPmot[2] ;
			pmdec = wPmot[3] ;
		}
#ifndef PRSTARCAT_PM_WITHD
		/** remove the cos(delta) factor from the proper motion in alpha if the internal definition
		* does not use it.
		*/
		pmra /= cos(wradec.dec.rad) ;
#endif

		hdn=w.hdn() ;

		/** retrieve the V or K magnitude (only one of which will be shown in the WDS entry)
		*/
		VMag = w.Vmag(useprim) ;
		KMag = w.Kmag(useprim) ;

		/** check out the spectral types, if available */
		styp = w.spectyp(useprim) ;
	}

#ifdef HAVE_NOVAS
	/** Convert the coordinate system from ICRS to J2000 FK5 coordinates.
	* The current version calls the Fortran77 NOVAS2.9b subroutine CATRAN().
	* @since 03/12/2006
	* @author Richard J. Mathar
	*/
	void tofk5()
	{
		if ( frame == "ICRS" )
		{
			int it = 5 ;
			/* Dummy arguments according to the documentation, but their addresess
			* cannot be set to 0, as it seems
			*/
			double date1, date2 ;  
			/** Convert RA argument from degrees to hours.
			*/
			double ra1 = ra/15. ;
			double dec1 = dec ;
#ifndef PRSTARCAT_PM_WITHD
			double pmra1 = pmra*cos(DEG2RAD(dec)) ;
#else
			double pmra1 = pmra ;
#endif
			double pmdec1 = pmdec ;
			double parx1 = px ;
			double rv1 = rv ;
			if ( isnan(ra1) || isnan(dec1) || isnan(pmra1) || isnan(pmdec1) || isnan(parx1) || isnan(rv1) )
				cerr << "Skipping " << target << ": undef. coordinates, proper motions, parallax or rad vel\n" ;
			else
			{
				double ra2 ;
				double dec2 ;
				double pmra2 ;
				double pmdec2 ;
				double parx2 ;
				double rv2 ;
				catran_(&it,&date1,&ra1,&dec1,&pmra1,&pmdec1,&parx1,&rv1,
					&date2,&ra2,&dec2,&pmra2,&pmdec2,&parx2,&rv2) ;
				ra = ra2*15. ;
				dec = dec2 ;
#ifndef PRSTARCAT_PM_WITHD
				pmra = pmra2/cos(DEG2RAD(dec2)) ;
#else
				pmra = pmra2 ;
#endif
				pmdec = pmdec2 ;
				px =  parx2 ;
				rv =  rv2 ;
				/* debugging
				* cerr << ra1 << " " << dec1 << " " << pmra1 << " " << pmdec1 << endl ;
				* cerr << ra2 << " " << dec2 << " " << pmra2 << " " << pmdec2 << endl ;
				*/
				frame = "FK5" ;
				equinox = 2000. ;
			}
		}
	}
#endif

	/**  Convert into a HH MM SS.sss and/or +-DD MM SS.sss string format.
	* @param use_ra specify whether the output should contain the RA value.
	* @param use_dec specify whether the output should contain the DEC value.
	*    The use of the two flags \c use_ra and \use_dec is independent. If
	*    both are specified, the two strings are concatenated, separated by a blank.
	* @param colon specify whether the hours/degrees and minutes should be followed
	*    by a colon each.
	* @param secprec indicate how many decimals follow after the dot of the
	*    seconds string.
	* @return strings of the hexdecimal format for star coordinates
	*/
	string coords(bool use_ra, bool use_dec, bool colon, int secprec, bool apes=false) const
	{
		char coor[2][33] ;
		char fmt[30] ;
		if ( use_ra )
		{
			if ( apes )
				sprintf(fmt,"%%02dh %%02dm %%0%d.%df",3+secprec,secprec) ;
			else if ( colon )
				sprintf(fmt,"%%02d:%%02d:%%0%d.%df",3+secprec,secprec) ;
			else
				sprintf(fmt,"%%02d %%02d %%0%d.%df",3+secprec,secprec) ;
			double raSec =ra/15. ;
			int raHr= int(raSec) ;
			raSec -= raHr ; raSec *= 60 ;
			int raMin = int(raSec) ;
			raSec -= raMin ; raSec *= 60 ;
			sprintf(coor[0],fmt,raHr,raMin,raSec) ;
		}
		if ( use_dec )
		{
			if ( apes )
				sprintf(fmt,"%%+03dd %%02dm %%0%d.%df",3+secprec,secprec) ;
			else if ( colon)
				sprintf(fmt,"%%+03d:%%02d:%%0%d.%df",3+secprec,secprec) ;
			else
				sprintf(fmt,"%%+03d %%02d %%0%d.%df",3+secprec,secprec) ;
			double decSec = fabs(dec) ;
			int decDeg= int(decSec) ;
			decSec -= decDeg ; decSec *= 60 ;
			int decMin = int(decSec) ;
			decSec -= decMin ; decSec *= 60 ;
			if ( dec >=0.)
				sprintf(coor[1],fmt,decDeg,decMin,decSec) ;
			else
				sprintf(coor[1],fmt,-decDeg,decMin,decSec) ;
		}
		if ( use_ra && use_dec)
		{
			string resul(coor[0]) ;
			resul += " " + string(coor[1]) ;
			return resul ;
		}
		else if ( use_ra )
			return string(coor[0]) ;
		else if ( use_dec )
			return string(coor[1]) ;
		else
			return string() ;
	}

	/**  Convert into a HHMM.m+-DDMM string format in the WDS style.
	* @return strings of the hexdecimal format for star coordinates
	*/
	string wdsName() const
	{
		char coor[2][7] ;

		double raMin =ra/15. ;
		const int raHr= int(raMin) ;
		raMin -= raHr ; raMin *= 60 ;
		sprintf(coor[0],"%02d%04.1f",raHr,raMin) ;

		double decMin = fabs(dec) ;
		const int decDeg= int(decMin) ;
		decMin -= decDeg ; decMin *= 60 ;
		if ( dec >=0.)
			sprintf(coor[1],"%+03d%02.0f",decDeg,decMin) ;
		else
			sprintf(coor[1],"%+03d%02.0f",-decDeg,decMin) ;

		string resul(coor[0]) ;
		resul += string(coor[1]) ;
		return resul ;
	}

	/**  Convert into a JHHMMSS.ss+-DDMMSS.s string format in the IAU style.
	* The idea is to stay within the 16 character format of the OIFITS tabular format.
	* @return strings of the hexdecimal format for star coordinates
	* @since 2006-05-17
	* @mathar 2006-05-17
	*/
	string iauName() const
	{
		char coor[2][11] ;

		double raSec =ra/15. ;
		const int raHr= int(raSec) ;
		raSec -= raHr ; raSec *= 60 ;
		const int raMin= int(raSec) ;
		raSec -= raMin ; raSec *= 60 ;
		/* According to the rules: truncation, not rounding
		* sprintf(coor[0],"%02d%02d%05.2f",raHr,raMin,raSec) ;
		*/
		const int rafullSec= int(raSec) ;
		raSec -= rafullSec ; raSec *= 10 ;
		sprintf(coor[0],"%02d%02d%02d.%1d",raHr,raMin,rafullSec,(int)raSec) ;

		double decSec = fabs(dec) ;
		const int decDeg= int(decSec) ;
		decSec -= decDeg ; decSec *= 60 ;
		const int decMin= int(decSec) ;
		decSec -= decMin ; decSec *= 60 ;
		/* According to the rules: truncation, not rounding */
		if ( dec >=0.)
			// sprintf(coor[1],"%+03d%02d%04.1f",decDeg,decMin,decSec) ;
			sprintf(coor[1],"%+03d%02d%02d",decDeg,decMin,(int)decSec) ;
		else
			//sprintf(coor[1],"%+03d%02d%04.1f",-decDeg,decMin,decSec) ;
			sprintf(coor[1],"%+03d%02d%02d",-decDeg,decMin,(int)decSec) ;

		string resul(coor[0]) ;
		resul += string(coor[1]) ;
		return resul ;
	}

	/** Decompose the RA and/or DEC string into the internal two values.
	* @param loc a string of the format HH MM SS[.sss] +-DD MM SS[.sss]
	*	The fractional seconds are optional, and the minutes specs may
	*	be surrounded by colons. This is some inverse functionality 
	*       compared to coords().
	* @since 2006-03-03
	* @mathar 2006-03-03
	*/
	void coordsInv(string loc)
	{
		/** unify the representation by replacing the colons by blanks */
		string::size_type pos ;
		while ( (pos = loc.find(":")) != string::npos)
			loc.replace(pos,1," ") ;

		int raHr, raMin, decDeg, decMin ;
		double raSec, decSec ;

		/** @todo handle the case of one coordinate only
		* figure out (by looking for a + or - sign) whether this contains a DEC....
		* pos = loc.find(":")) != string::npos)
		*/
		sscanf(loc.c_str(),"%d %d %lf %d %d %lf",&raHr, &raMin, &raSec, &decDeg, &decMin, &decSec) ;
		ra = 15.*(raHr + raMin/60. + raSec/3600.) ;
		dec = fabs(decDeg) + decMin/60. + decSec/3600. ;
		if( decDeg < 0 )
			dec *= -1. ;
	}

	/** place a representation equivalent to the PRIMA simulator's input
	* catalogue into the output stream.
	* @param os output stream (cerr, cout,....)
	*/
	void xml(ostream & os) const
	{
		// cerr << __FILE__ << " " << __LINE__ << " " << target << endl ;
		os << "\t<source name=\"" << target << "\"" ;

		/** convert ra from degrees to HH:MM:SS.ssss and dec from degrees to +-DD:MM.SS.ssss
		*/
		os << " ra=\"" << coords(true,false,true,4) << "\"" ;
		os << " dec=\"" << coords(false,true,true,4) << "\"" << endl  ;

		os << "\tKmag=\"" << KMag << "\"" ;

		/** convert proper motion to as/100y. Also remove the cos(delta) factor which is not
		* used within the simulator.
		*/
#ifdef PRSTARCAT_PM_WITHD
		os << " PmotRA=\"" << pmra/10./cos(DEG2RAD(dec)) << "\"" ;
#else
		os << " PmotRA=\"" << pmra/10. << "\"" ;
#endif

		os << " equinox=\"" << "J" << setprecision(5) << equinox << "\"" ;

		/** convert proper motion to as/100y
		*/
		os << " PmotDec=\"" << pmdec/10. << "\"" ;

		/** convert parallax to arcsec
		*/
		os << " parallax=\"" << px/1.e3 << "\"" ;

		os << " vrad=\"" << rv << "\">" << endl ;
		os << "\t<sed lambda0=\"1.52\" dlambda=\"0.01\" nlambda=\"125\">" << endl ;
		os << "\t\t<blackbody temperature=\"" << Teff << "\"/>" << endl ;
		os << "\t</sed>" << endl ;
		os << "\t</source>" << endl ;
	}

	/** Transform the catalogue into the APES format.
	* Place a representation equivalent to the
	* \htmlonly
	* <a href="http://obswww.unige.ch/~segransa/apes/prototype.html">APES</a>
	* \endhtmlonly
	* \latexonly
	* APES \cite{vlt157510003}
	* \endlatexonly
	* input catalogue into the output stream. The supported keys are
	* type code alpha delta equinox mu_alpha mu_delta parallax T0 Period ecc a inc omega OMEGA
	*	SP_type Teff magV magK coord_syst
	* @param os output stream (cerr, cout,....)
	*/
	void apes(ostream & os) const
	{
		os << flag << "\t" ;
		os << target << "\t" ;

		/** convert ra from degrees to HHh MMm SS.sss and dec from degrees to +-DD:MM.SS.ssss
		*/
		os << coords(true,false,false,3,true) << "\t" ;
		os << coords(false,true,false,3,true) << "\t" ;
		os << setprecision(4) << setiosflags(ios_base::fixed) << equinox << "\t" ;

		/** The proper motions in APES are mas/a, not arcsec/a as in the current version of the
		* APES detailed description.
		* The APES catalog also includes a cos(delta) (D Segransan, priv comm 2006-03-06).
		*/
		os << setiosflags(ios_base::scientific) ;
#ifdef PRSTARCAT_PM_WITHD
		os << setprecision(3) << pmra << "\t" ;
#else
		os << setprecision(3) << pmra*cos(DEG2RAD(dec)) << "\t" ;
#endif
		os << setprecision(3) << pmdec << "\t" ;

		os << px << "\t" ;
		os << T0 << "\t" ;
		os << P << "\t" ;
		os << eccen << "\t" ;

		/* convert semi-major axis to muas */
		os << int(1.e6*a) << "\t" ;
		os << inclin << "\t" ;
		os << omega << "\t" ;
		os << Omega << "\t" ;
		os << styp << "\t" ;
		os << Teff << "\t" ;
		os << VMag << "\t" ;
		os << KMag << "\t" ;
		os << frame  ;

		os << endl ;
	}

	/** Transform the catalogue into the \c iphase format.
	* @param os output stream (cerr, cout,....)
	* @todo iphase requires a FK4 or FK5 frame. This transformation is just missing at the moment.
	*/
	void iphase(ostream & os) const
	{
		os << setw(6) << setfill('0') << indx << " " ;

		/** convert ra from degrees to HH MM SS.sss and dec from degrees to +-DD MM SS.sss
		*/
		os << coords(true,false,false,3,false) << " " ;
		os << coords(false,true,false,3,false) << " " ;

		// os << setiosflags(ios_base::scientific) ;
		os << setiosflags(ios_base::fixed | ios_base::right) ;
		/** convert proper motions in RA to seconds/yr. Iphase proper motions
		* are not used with the factor cos(delta), consistent with the rest of the SLALIB.
		*/
		os << setw(9) << setfill(' ') ;
		if ( finite(pmra) )
#ifdef PRSTARCAT_PM_WITHD
			os << setprecision(4) << pmra/1.e3/15./cos(DEG2RAD(dec)) << " " ;
#else
			os << setprecision(4) << pmra/1.e3/15. << " " ;
#endif
		else
			os << setprecision(4) << 0. << " " ;

		/** convert proper motions in DEC to arcseconds/yr */
		os << setw(7);
		if ( finite(pmdec) )
			os << setprecision(3) << pmdec/1.e3 << " " ;
		else
			os << setprecision(3) << 0. << " " ;

		os << setw(8) << setprecision(1) << setiosflags(ios_base::fixed) << equinox << " " ;

		/** convert parallax to arcseconds */
		os << setw(8) << setprecision(3)  ;
		if ( finite(px) )
			os << px/1.e3 << " " ;
		else
			os << 0. << " " ;

		os << setw(7) << setprecision(1)  ;
		if ( finite(rv) )
			os << rv << " " ;
		else
			os << 0. << " " ;

		os << endl ;
	}

	/** Transform the star coordinates into gnuplot format.
	* @param[in] lst Paranal sidereal time [rad]
	* @param[in] threeD if true generate splot commands, else plot commands
	* @param[in,out] os output stream (cerr, cout,....)
	*/
	void gp(ostream & os, const double lst, bool threeD, const float minElev) const
	{
		/* Paranal GPS geographical latitude [rad] */
		const double phi= -0.429843732 ;

		double ha =lst - DEG2RAD(ra) ;
		os << "# ha " << RAD2DEG(ha) << " deg " << 24.*RAD2DEG(ha)/360. << " hrs " << endl ;
		/* azimuth is N=0, E=pi/2 */
		double az, el ;
		slaDe2h(ha,DEG2RAD(dec),phi,& az,&el) ;
		os << "# az " << RAD2DEG(az) << " deg, el " << RAD2DEG(el) << " deg " << endl ;
		if ( threeD)
		{
			/* Cartesian vector with x towards E, y torwards N, z to zenith.
			* Switch to the slaDcs2c convention with the spherical azimuth
			* zero along x, pi/2 along y */
			az = slaDranrm(M_PI_2-az) ;
			double xyz[3] ;
			slaDcs2c(az,el,xyz) ;
			/* select a point style depending on the quadrant of the azimuth */
			const int pt = (int)(az/M_PI_2) ;
			os << "set label \"" << target  << "\" at " << xyz[0] << "," << xyz[1] << "," << xyz[2] 
				<< " point pt " << pt+1 << " tc palette z " << endl ;

			/* and place another one at the ground */
			if ( el > minElev )
				os << "set label \"" << ""  << "\" at " << xyz[0] << "," << xyz[1] 
					<< ",0 point pt " << pt+1 << endl;
		}
		else
		{
			os << "set label \"" << target  << "\" at " << RAD2DEG(az) << "," << RAD2DEG(el) << " point pt 2" << endl ;
		}
	}

protected:
private:
	string target ;

	static int no ;
} ;
int Star::no =0 ;

/** Ordering operator
* Stars are ordered with respect to right ascension (major order) and declination.
* @return true if the \c left star is ordered first in an ordered listing.
* @since 2006-03-03
* @author Richard J. Mathar
*/
bool operator< (const Star &left, const Star &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 ;
}

/** Ordering operator with respect to declination.
* Stars are alternatively ordered numerically along increasing declination.
* @return true if the \c left star is ordered first in an ordered listing.
* @since 2006-03-09
* @author Richard J. Mathar
*/
struct StarDecCmp : public binary_function<const Star & , const Star & , bool> {
	bool operator() (const Star &left, const Star &right)
	{
		if ( left.dec < right.dec )
			return true ;
		else if ( left.dec > right.dec )
			return false ;
		else if ( left.ra < right.ra )
			return true ;
		else
			return false ;
	}
} ;

/** A star catalogue.
* The star catalogue is a list of entries of the Star class.
*/
class starCat : public vector<Star> {
public:
	/** default ctor. Empty list generated. */
	starCat()
	{
	}

	/** Creator from a standard FITS file 
	* @param fits a FITS instance with the standard PRIMA layout of the STARS table
	*/
	starCat(const FITS & fits)
	{
		const ExtHDU & e = fits.extension(PRSTARCAT_BINT) ;
		const long nrows = e.rows() ;
		for(int r=0 ; r < nrows ; r++)
		{
			Star star(e,r) ;
			add(star) ;
		}
	}

	/** Creator from a WDS listing
	* @param wdsvec a vector of WDS lines
	*/
	starCat(const vector<wdsstar> * wdsvec) : vector<Star>(0)
	{
		/* debugging
		* cerr << __FILE__ << " " << __LINE__ << " " << wdsvec->size() << endl ;
		*/

		/** Loop over all entries in the WDS vector. For each entry, generate
		* the two components and add them individually to the catalog.
		*/
		for(int r=0 ; r < wdsvec->size() ; r++)
		{
			/* debugging
			* cerr << __FILE__ << " " << __LINE__ << " " << r << endl ;
			* cerr << __FILE__ << " " << __LINE__ << " " << (*wdsvec)[r] << endl ;
			*/
			Star starT( (*wdsvec)[r],true) ;
			add(starT) ;
			Star starR( (*wdsvec)[r],false) ;
			add(starR) ;
		}
	}

	/** Add a single star to the list
	* @param star the instance to be appended.
	*/
	void add(Star &star)
	{
		push_back(star) ;
	}

	/** Convert into a FITS binary table.
	* This conversion function defines the columns in the \c STARS binary
	* table of the FITS files.
	* @return a pointer to the FITS representation of the generated file.
	* Note that the binary table STARS is generated here, but not yet filled
	* with the actual data, which is done in starCat::outfits(). 
	* @author Richard J. Mathar
	*/
	operator FITS*()
	{
		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("COOR") ; colUns.push_back("") ; colFmt.push_back("32A") ;
		colNams.push_back("EQUINOX") ; colUns.push_back("") ; colFmt.push_back("E") ;
		colNams.push_back("RADECSYS") ; colUns.push_back("") ; colFmt.push_back("4A") ;
#ifdef PRSTARCAT_FITS_DEG
		colNams.push_back("RA") ; colUns.push_back("deg") ; colFmt.push_back("D") ;
		colNams.push_back("DEC") ; colUns.push_back("deg") ; colFmt.push_back("D") ;
#endif
		colNams.push_back("RA_ERR") ; colUns.push_back("deg") ; colFmt.push_back("E") ;
		colNams.push_back("DEC_ERR") ; colUns.push_back("deg") ; colFmt.push_back("E") ;

		colNams.push_back("PMRA") ; colUns.push_back("mas/a") ; colFmt.push_back("E") ;
		colNams.push_back("PMRA_ERR") ; colUns.push_back("mas/a") ; colFmt.push_back("E") ;
		colNams.push_back("PMDEC") ; colUns.push_back("mas/a") ; colFmt.push_back("E") ;
		colNams.push_back("PMDEC_ERR") ; colUns.push_back("mas/a") ; colFmt.push_back("E") ;
		colNams.push_back("PX") ; colUns.push_back("mas") ; colFmt.push_back("E") ;
		colNams.push_back("PX_ERR") ; colUns.push_back("mas") ; colFmt.push_back("E") ;
		colNams.push_back("RV") ; colUns.push_back("km/s") ; colFmt.push_back("E") ;
		colNams.push_back("RVTYPE") ; colUns.push_back("") ; colFmt.push_back("8A") ;
		colNams.push_back("RVDEF") ; colUns.push_back("") ; colFmt.push_back("8A") ;
		colNams.push_back("TEFF") ; colUns.push_back("K") ; colFmt.push_back("E") ;
		colNams.push_back("VMAG") ; colUns.push_back("1") ; colFmt.push_back("E") ;
		colNams.push_back("KMAG") ; colUns.push_back("1") ; colFmt.push_back("E") ;
		colNams.push_back("STYP") ; colUns.push_back("") ; colFmt.push_back("6A") ;
		colNams.push_back("P") ; colUns.push_back("d") ; colFmt.push_back("E") ;
		colNams.push_back("A") ; colUns.push_back("arcsec") ; colFmt.push_back("E") ;
		colNams.push_back("I") ; colUns.push_back("deg") ; colFmt.push_back("E") ;
		colNams.push_back("OMEGG") ; colUns.push_back("deg") ; colFmt.push_back("E") ;
		colNams.push_back("T0") ; colUns.push_back("d") ; colFmt.push_back("D") ;
		colNams.push_back("OMEG") ; colUns.push_back("deg") ; colFmt.push_back("E") ;
		colNams.push_back("E") ; colUns.push_back("1") ; colFmt.push_back("E") ;
		colNams.push_back("BFLAG") ; colUns.push_back("") ; colFmt.push_back("1A") ;
		colNams.push_back("HDN") ; colUns.push_back("") ; colFmt.push_back("J") ;
		colNams.push_back("HIC") ; colUns.push_back("") ; colFmt.push_back("J") ;
		colNams.push_back("FK6") ; colUns.push_back("") ; colFmt.push_back("I") ;
		colNams.push_back("ADS") ; colUns.push_back("") ; colFmt.push_back("8A") ;
		colNams.push_back("SAO") ; colUns.push_back("") ; colFmt.push_back("J") ;
		colNams.push_back("WDS") ; colUns.push_back("") ; colFmt.push_back("10A") ;
		colNams.push_back("INDX") ; colUns.push_back("1") ; colFmt.push_back("I") ;
		fits->addTable(PRSTARCAT_BINT,0,colNams,colFmt,colUns) ;
			ExtHDU & e = fits->extension(PRSTARCAT_BINT) ;
			e.column("HDN").addNullValue(-1L) ;
			e.column("HIC").addNullValue(-1L) ;
			e.column("FK6").addNullValue(short(-1)) ;
			e.column("SAO").addNullValue(-1L) ;
			e.column("INDX").addNullValue(short(-1)) ;
		return fits ;
	}


	/** Append (add) another star catalog.
	* @param oth the additional star catalog to be appended.
	* @since 2006-03-05
	* @author Richard J. Mathar
	*/
	starCat & operator += (const starCat & oth)
	{
		vector<Star>::iterator last = end() ;
		insert( last,  oth.begin(), oth.end()) ;
		return *this ;
	}

	/** Generate the PRIMA standard FITS format.
	* @param fits the output file which had been created by operator FITS*().
	* @since 2006-03-03
	* @author Richard J. Mathar
	*/
	void outFits(FITS *fits)
	{

		for(int s =0 ; s < int(size()) ; )
		{
			ExtHDU & e = fits->extension(PRSTARCAT_BINT) ;
			e.makeThisCurrent() ;

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

			e.column("EQUINOX").write(& at(s).equinox,1L,nrows+1) ;

			asc[0] = at(s).frame ;
			e.column("RADECSYS").write(asc,nrows+1) ;

			float tmp = at(s).KMag ;
			e.column("KMAG").write(&tmp,1L,nrows+1) ;
			tmp = at(s).VMag ;
			e.column("VMAG").write(&tmp,1L,nrows+1) ;

			tmp = at(s).Teff ;
			e.column("TEFF").write(&tmp,1L,nrows+1) ;

			asc[0] = at(s).styp ;
			e.column("STYP").write(asc,nrows+1) ;

			asc[0] = at(s).ads ;
			e.column("ADS").write(asc,nrows+1) ;

#ifdef PRSTARCAT_FITS_DEG
			e.column("RA").write(& at(s).ra,1L,nrows+1) ;
			e.column("DEC").write(& at(s).dec,1L,nrows+1) ;
#endif
			tmp = at(s).ra_err ;
			e.column("RA_ERR").write(& tmp,1L,nrows+1) ;
			tmp = at(s).dec_err ;
			e.column("DEC_ERR").write(& tmp,1L,nrows+1) ;

			asc[0] = at(s).coords(true,true,false,6) ;
			e.column("COOR").write(asc,nrows+1) ;

			tmp = at(s).pmra ;
			e.column("PMRA").write(&tmp,1L,nrows+1) ;
			tmp = at(s).pmra_err ;
			e.column("PMRA_ERR").write(&tmp,1L,nrows+1) ;
			tmp = at(s).pmdec ;
			e.column("PMDEC").write(&tmp,1L,nrows+1) ;
			tmp = at(s).pmdec_err ;
			e.column("PMDEC_ERR").write(&tmp,1L,nrows+1) ;

			tmp = at(s).rv ;
			e.column("RV").write(&tmp,1L,nrows+1) ;

			asc[0] = at(s).veltyp ;
			e.column("RVTYPE").write(asc,nrows+1) ;
			asc[0] = at(s).veldef ;
			e.column("RVDEF").write(asc,nrows+1) ;

			tmp = at(s).px ;
			e.column("PX").write(&tmp,1L,nrows+1) ;
			tmp = at(s).px_err ;
			e.column("PX_ERR").write(&tmp,1L,nrows+1) ;

			e.column("HIC").write(& at(s).hip,1L,nrows+1) ;
			e.column("SAO").write(& at(s).sao,1L,nrows+1) ;
			e.column("HDN").write(& at(s).hdn,1L,nrows+1) ;
			e.column("FK6").write(& at(s).fk6,1L,nrows+1) ;

			asc[0] = at(s).wds ;
			e.column("WDS").write(asc,nrows+1) ;

			asc[0] = at(s).flag ;
			e.column("BFLAG").write(asc,nrows+1) ;

			s++ ;
			e.column("INDX").write(&s,1L,nrows+1) ;
		}
	}

	/** Generate the OIFITS format on standard output.
	* This contains a \c OI_TARGET table as described in
	* \latexonly
	* \cite{PaulsSPIE5491} and \url{http://www.mrao.cam.ac.uk/research/OAS/oi_data/oifits.html}.
	* \endlatexonly
	* \htmlonly
	* <a href="http://www.mrao.cam.ac.uk/~jsy1001/exchange/">OIFITS</a>.
	* \endhtmlonly
	* @since 2006-03-03
	* @author Richard J. Mathar
	*/
	void oifits()
	{
		FITS *fits = new FITS("-",Write) ;
		fits->pHDU().makeThisCurrent() ;
		fits->pHDU().addKey("CREATOR",__FILE__,"Program to create this file (Richard J. Mathar)") ;
		fits->pHDU().writeComment("OIFITS http://www.mrao.cam.ac.uk/research/OAS/oi_data/oifits.html") ;
		fits->pHDU().writeDate() ;

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

		colNams.push_back("TARGET_ID") ; colUns.push_back("1") ; colFmt.push_back("I") ;
		colNams.push_back("TARGET") ; colUns.push_back("") ; colFmt.push_back("16A") ;
		colNams.push_back("RAEPO") ; colUns.push_back("deg") ; colFmt.push_back("D") ;
		colNams.push_back("DECEPO") ; colUns.push_back("deg") ; colFmt.push_back("D") ;
		colNams.push_back("EQUINOX") ; colUns.push_back("") ; colFmt.push_back("E") ;
		colNams.push_back("RA_ERR") ; colUns.push_back("deg") ; colFmt.push_back("D") ;
		colNams.push_back("DEC_ERR") ; colUns.push_back("deg") ; colFmt.push_back("D") ;
		colNams.push_back("SYSVEL") ; colUns.push_back("m/s") ; colFmt.push_back("D") ;
		colNams.push_back("VELTYP") ; colUns.push_back("") ; colFmt.push_back("8A") ;
		colNams.push_back("VELDEF") ; colUns.push_back("") ; colFmt.push_back("8A") ;
		colNams.push_back("PMRA") ; colUns.push_back("deg/a") ; colFmt.push_back("D") ;
		colNams.push_back("PMDEC") ; colUns.push_back("deg/a") ; colFmt.push_back("D") ;
		colNams.push_back("PMRA_ERR") ; colUns.push_back("deg/a") ; colFmt.push_back("D") ;
		colNams.push_back("PMDEC_ERR") ; colUns.push_back("deg/a") ; colFmt.push_back("D") ;
		colNams.push_back("PARALLAX") ; colUns.push_back("deg") ; colFmt.push_back("E") ;
		colNams.push_back("PARA_ERR") ; colUns.push_back("deg") ; colFmt.push_back("E") ;
		colNams.push_back("SPECTYP") ; colUns.push_back("") ; colFmt.push_back("16A") ;

		fits->addTable("OI_TARGET",size(),colNams,colFmt,colUns) ;
		ExtHDU & ext = fits->extension("OI_TARGET") ;
		ext.column("TARGET_ID").addNullValue(short(-1)) ;
		ext.addKey("OI_REVN",1,"Revision number of the table definition") ;

		for(int s =0 ; s < int(size()) ; s++)
		{
			ExtHDU & ext = fits->extension("OI_TARGET") ;
			ext.makeThisCurrent() ;

			long nrows= s ;
			vector<string> asc(1) ;

			short tmpI = at(s).indx ;
			ext.column("TARGET_ID").write(& tmpI,1L,nrows+1) ;

			// asc[0] = at(s).target ;
			/** As the target's name we use a IAU type of identifier */
			asc[0] = at(s).iauName()+at(s).flag ;
			ext.column("TARGET").write(asc,nrows+1) ;

			ext.column("RAEPO").write(& at(s).ra,1L,nrows+1) ;
			ext.column("DECEPO").write(& at(s).dec,1L,nrows+1) ;

			float tmpE = at(s).equinox ;
			ext.column("EQUINOX").write(& tmpE,1L,nrows+1) ;

			ext.column("RA_ERR").write(& at(s).ra_err,1L,nrows+1) ;
			ext.column("DEC_ERR").write(& at(s).dec_err,1L,nrows+1) ;

			/** For the radial velocities we switch from km/s to m/s */
			double tmpD = 1.e3* at(s).rv ;
			ext.column("SYSVEL").write(&tmpD,1L,nrows+1) ;
			asc[0] = at(s).veltyp ;
			ext.column("VELTYP").write(asc,nrows+1) ;
			asc[0] = at(s).veldef ;
			ext.column("VELDEF").write(asc,nrows+1) ;

			/** Convert the proper motions from mas to degrees on output */
			ext = fits->extension("OI_TARGET") ;
			tmpD = MAS2DEG(at(s).pmra) ;
			ext.column("PMRA").write(& tmpD,1L,nrows+1) ;
			tmpD = MAS2DEG(at(s).pmra_err) ;
			ext.column("PMRA_ERR").write(& tmpD,1L,nrows+1) ;
			tmpD = MAS2DEG(at(s).pmdec) ;
			ext.column("PMDEC").write(&tmpD,1L,nrows+1) ;
			tmpD = MAS2DEG(at(s).pmdec_err) ;
			ext.column("PMDEC_ERR").write(&tmpD,1L,nrows+1) ;

			/** Convert the parallax from mas to degrees on output */
			tmpE = MAS2DEG(at(s).px) ;
			ext.column("PARALLAX").write(&tmpE,1L,nrows+1) ;
			tmpE = MAS2DEG(at(s).px_err) ;
			ext.column("PARA_ERR").write(&tmpE,1L,nrows+1) ;

			asc[0] = at(s).styp ;
			ext.column("SPECTYP").write(asc,nrows+1) ;
		}
		delete fits ;
	}

	/** Generate the APES format on the output stream.
	* @param[in] the outputstream to write to.
	* @since 2006-03-04
	* @author Richard J. Mathar
	*/
	void apes( ostream & os) const
	{
		/** generate the header lines */
		os << "type\tcode\talpha\tdelta\tequinox\tmu_alpha\tmu_delta\tparallax\tT0\tPeriod\tecc\ta\tinc\tomega\tOMEGA\t" ;
		os << "SP_type\tTeff\tmagV\tmagK\tcoord_syst\t" << endl ;
		os << "----\t----\t-----\t-----\t-------\t--------\t--------\t--------\t--\t------\t---\t-\t---\t-----\t-----\t" ;
		os << "-------\t----\t----\t----\t----------\t" << endl ;
		for(int s=0 ; s < size() ; s++)
			at(s).apes(os) ;
	}

	/** Generate the XML format on the output stream.
	* @param[in] the outputstream to write to.
	* @since 2006-03-03
	* @author Richard J. Mathar
	*/
	void xml(ostream & os) const
	{
		for(int s=0 ; s < size() ; s++)
			at(s).xml(os) ;
	}

	/** Translate the stars in the ICRS frame to FK5.
	* @since 2006-03-12
	* @author Richard J. Mathar
	* The work is done by calling Star::tofk5() for all stars in the list.
	*/
	void tofk5()
	{
		for(int s=0 ; s < size() ; s++)
			at(s).tofk5() ;
	}

	/** Generate the iphase format on the output stream.
	* @parameter os the output stream to write to.
	* @since 2006-03-06
	* @author Richard J. Mathar
	* The work is done by calling Star::iphase() for all stars in the list.
	*/
	void iphase(ostream & os) const
	{
		const time_t now = time(NULL) ;
		os << "! generated by " << __FILE__ << " " << ctime(&now) << endl ;
		for(int s=0 ; s < size() ; s++)
			at(s).iphase(os) ;
		os << "END" << endl ;
	}

	/** Generate a gnuplot command file on the output stream.
	* @paramp[in] os the output stream to write to.
	* @paramp[in] utc the UTC time stamp of the observation
	* @paramp[in] threeD if true, splot command are generated, else plot commands.
	* @since 2006-03-08
	* @author Richard J. Mathar
	*/
	void gp(ostream & os, string utc, bool threeD) const
	{
		const time_t now = time(NULL) ;
		os << "# generated by " << __FILE__ << " " << ctime(&now) << endl ;
		const float minElev = 30. ;

		if ( threeD)
		{
			os << "set xrange [-1:1]" << endl ;
			os << "set label \"E\" at 1.1,0,0" << endl ;
			os << "set label \"W\" at -1.1,0,0" << endl ;
			os << "set yrange [-1:1]" << endl ;
			os << "set label \"N\" at 0,1.1,0" << endl ;
			os << "set label \"S\" at 0,-1.1,0" << endl ;
			os << "set zrange [0:1]" << endl ;
			os << "set zlabel \"z\" " << endl ;
			os << "set ticslevel 0" << endl ;
			os << "set parametric" << endl ;
			os << "set angles degrees" << endl ;
			os << "set urange [0:360]" << endl ;
			os << "set vrange [0:90]" << endl ;
			// os << "set size 1,0.5" << endl ;
			os << "unset xtics" << endl ;
			os << "unset ytics" << endl ;
			os << "unset ztics" << endl ;
		}
		else
		{
			os << "set xrange [-10:370]" << endl ;
			os << "set xlabel \"azimuth (deg), N=0, E=90\" " << endl ;
			os << "set yrange [0:90]" << endl ;
			os << "set ylabel \"elevation (deg)\" " << endl ;
			os << "set label \"N\" at 0,0" << endl ;
			os << "set label \"E\" at 90,0" << endl ;
			os << "set label \"S\" at 180,0" << endl ;
			os << "set label \"W\" at 270,0" << endl ;
			os << "set label \"N\" at 360,0" << endl ;
		}

		/** Decompose the UTC string into year, month, day etc */
		int year, month, day , hrs, min, sec, status=0;
		sscanf(utc.c_str(),"%d-%d-%dT%d:%d:%d",&year,&month,&day,&hrs,&min,&sec) ;

		/** compute modified Julian date */
		double mjd ;
		slaCldj(year,month,day,&mjd,&status) ;
		if ( status )
			cerr << "Invalid year " << year << ", month " << month << " or day " << day << " in " << utc << endl ;

		double fracdays ;
		slaDtf2d(hrs,min,sec,&fracdays,&status) ;
		if ( status )
			cerr << "Invalid hour " << hrs << ", minute " << min << " or second " << sec << " in " << utc << endl ;

		/** convert modified Julian date to Greenwich sidereal time [rad], and then Paranal sidereal time */
		const double lst = slaDranrm(slaGmsta(mjd,fracdays) - 1.228803656)  ;

		os << "# for " << utc ;
		os << ", MJD " << mjd+fracdays << " LST " << lst*12*3600./M_PI << " s = " << lst*12/M_PI << " h" << endl ;

		for(int s=0 ; s < size() ; s++)
			at(s).gp(os,lst,threeD,DEG2RAD(minElev)) ;

		if ( threeD)
		{
			/* Plot two hula-hoops at the horizon and the minimum elevation */
			os << "splot cos(u),sin(u),0 wi lines lt 0 notitle" << endl ;
			os << "replot 0,cos(v),sin(v) wi lines lt 0 notitle" << endl ;
			os << "replot cos(u)*cos(" << minElev << "),sin(u)*cos(" << minElev << "),sin(" << minElev << ")" ;
			os << " wi lines lt 0 title \"" << utc << "\"" << endl ;
		}
		else
		{
			os << "plot " << minElev << " title \"" << utc << "\"" << endl ;
		}
	}
protected:
private:
} ;

#undef PRSTARCAT_BINT

#endif /* PRSTARCAT_C */

/**@}*/

/*
$Log: prStarCat.C,v $
Revision 1.8  2006/09/06 13:14:43  mathar
Corrected units of star catalog output to OIFITS format.

Revision 1.7  2006/03/12 20:38:56  mathar
Implemented functions tofk5() to convert ICRS to FK5 frames.

*/