#ifndef STAR_C
#define STAR_C
/**
* @since 2006-08-28
* @author Richard J. Mathar
*/
#include <iostream>
#include <fstream>
#include <stdlib.h>
#include <string.h>
#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 "units.h"
#include "BlackB.h"

#include "Star.h"

#include "Geocen.h"

/**
* default ctor
*/
Star::Star() : ra(0.), dec(0.), pmra(0.), pmdec(0.), bbscale(0.)
{
}

/**
* ctor
* @param[in] f the FITS file with the PS and SS keywords in the PHDU
* @param[in] indx 0 for the PS and 1 for the SS
*/
Star::Star(FITS &f, const int indx)
{
	PHDU & p = f.pHDU() ;
	string hkey = "HIERARCH ESO PRI " ;
	switch(indx)
	{
	case 0 :
		hkey += "PS " ;
		break ;
	default :
		hkey += "SS " ;
		break ;
	}
	/**
	* We read HIERARCH ESO PRI {PS|SS} {RA|DEC} which are in [deg] in the primary FITS header,
	* and the proper motions in [arcsec/yr] from the same place. IF this is a MIDI
	* file, the alternative is to set the proper motions to zero and just to read RA and DEC.
	*/
	try {
		p.readKey(hkey+"RA", ra) ;
		p.readKey(hkey+"DEC", dec) ;
		p.readKey(hkey+"PMA", pmra) ;
		p.readKey(hkey+"PMD", pmdec) ;
	}
	catch(...)
	{
		p.readKey("RA", ra) ;
		p.readKey("DEC", dec) ;
		pmra=pmdec=0. ;
	}

	ra = DEG2RAD(ra) ;
	dec = DEG2RAD(dec) ;

	/* debugging
	* cout << __FILE__ << " " << __LINE__ << " "  << ra << " " << dec << endl ;
	*/

	/**
	* We read HIERARCH ESO PRI {PS|SS} TEMP and MAG also from the primay FITS header.
	* If these are not available, defaults of 8000K and 7 mag are assumed.
	*/
	try
	{
		p.readKey(hkey+"TEMP", T) ;
		if ( T == 0.)
			T=8000. ;
	}
	catch(...)
	{
		T=8000. ;
	}

	try
	{
		p.readKey(hkey+"MAG", Kmag) ;
		if ( Kmag == 0.)
			Kmag= 7. ;
	}
	catch(...)
	{
		Kmag=7. ;
	}

	/** We set the zero point of the Jy scale in the K band at 2.2 micron and 620 Jy.
	* \c jansk becomes the Jy at our particular magnitude.
	*/
	const double jansk = 620.*pow(10.,-0.4*Kmag) ;
	const BlackB bb(T) ;

	/** \c bbscale is the scaling factor to convert a spectral density in W/m^2/Hz
	* of the template black body to W/m^2/Hz at the telescope entrance.
	*/
	bbscale = jansk*1.e-26/bb.densitnu( WN2HZ(MICRON2WN(2.2)) ) ;

	/* debugging
	*cerr << __FILE__ << " " << __LINE__ << " "  << Kmag << " mag " << T << " K" <<  endl ;
	*/
}

/**
* ctor
* @param[in] alpha right ascension [rad]
* @param[in] delta declination [rad]
* @param[in] pmAlpha proper motion in alpha [arcsec/yr]
* @param[in] pmDelta proper motion in delta [arcsec/yr]
*/
Star::Star(const double alpha, const double delta, const double pmAlpha, const double pmDelta)
	: ra(alpha), dec(delta), pmra(pmAlpha), pmdec(pmDelta)
{
}

/**
* Compute a pointing vector with simple spherical astrometry
* @parameter[in] pos the geocentric position of the observation
* @parameter[in] ha the hour angle [rad]
* @return a topocentric vector with three components normalized to 1.
*   component[0] is the N, component[1] the West and component[2] the zenith direction
*/
vector<double> Star::point(const Geocen & pos, const double ha) const
{
	vector<double> b ;
	/*
	* This is (15) of astro-ph/0608273 with a sign flip to account for the sign in (13).
	*/
	b.push_back( -sin(pos.geolat)*cos(dec)*cos(ha)+cos(pos.geolat)*sin(dec) ) ;
	/*
	* This is (14) of astro-ph/0608273 with a sign flip to account for the different az convetion here
	*/
	b.push_back( cos(dec)*sin(ha) ) ;
	/*
	* This is (16) of astro-ph/0608273.
	*/
	b.push_back( cos(pos.geolat)*cos(dec)*cos(ha)+sin(pos.geolat)*sin(dec) ) ;

	/* debugging
	*for(int c=0 ; c < 3; c++)
	*	cout << c << " " << b[c] << endl ;
	*/
	return b ;
}

/**
* Compute the time derivative of pointing vector with simple spherical astrometry
* @parameter[in] pos the geocentric position of the observation
* @parameter[in] ha the hour angle [rad]
* @return a vector with three components normalized to 1/sec
*   component[0] is the N, component[1] the West and component[2] the zenith direction
*/
vector<double> Star::pointDt(const Geocen & pos, const double ha) const
{
	vector<double> b ;
	/** This is the derivative w.r.t \c ha, as derived from point(), post-mutliplied with the
	* derivative of \c ha (which has units of [rad]) w.r.t time, which is 2 pi in 24 hrs.
	*/
	b.push_back( sin(pos.geolat)*cos(dec)*sin(ha) ) ;
	b.push_back( cos(dec)*cos(ha) ) ;
	b.push_back( -cos(pos.geolat)*cos(dec)*sin(ha) ) ;

	for(int c=0 ; c < 3; c++)
		b[c] *= 2.*M_PI/SECPERDAY ;
	return b ;
}

#if 0
/**
* correlated flux
* @remarks not implemented since we only need the integrated flux over some
*  region of wavenumbers covered by the simulated spectrum.
*
* @param[in] wn wavenumber [1/cm]
* @return the equivalent flux (above the atmosphere) at that point [W/m^2/Hz].
*/
double Star::flux(const double wn) const
{
	return 1. ;
}
#endif

/**
* correlated flux
*
* @param[in] wnstart lower wavenumber in the spectrum [1/cm]
* @param[in] wnstop upper wavenumber in the spectrum [1/cm]
* @return the equivalent flux (above the atmosphere) at that point, integrated over the spectral region [W/m^2]
*/
double Star::flux(double wnstart, double wnstop) const
{
	/* initialize a black body of the same temperature. This is done once per call, so
	* about as many times per program run as there are spectral channels.
	*/
	const BlackB bb(T) ;

	/* the wavenumber region in units of Hz
	*/
	double nu[2] ;
	nu[0]= WN2HZ(wnstart) ; nu[1]= WN2HZ(wnstop) ;

	/** The black body flux in this spectral region is bb.fluxnu(nu,false). We call the
	* BlackB::fluxnu() to integrate the black body spectrum between the corresponding frequencies.
	*/
	return bbscale * bb.fluxnu(nu,true) ;
}

// __oOO__
#endif // STAR_C