#ifndef PI
#define PI 	3.141592653589793
#endif


#ifndef M_PI
#define M_PI		3.14159265358979323846
#endif

#define rad(x)  ((x)*PI/180.)
#define deg(x)  ((x)*180./PI)
#define hdeg(x) ((x)*15.)
#define degh(x) ((x)/15)


static double REST24( double val );
static double REST( double val, double t );
static void cCurrentStartime( Int32 utcadjust, double longitude, Int32 hsecs );
static void cAzAlt( double latitude, double ra, double dec );
static void cRaDec( double latitude, double az, double alt );
static double makeAbsoluteDouble( double val );




/*
static double rad( double d );
static double deg( double r );
*/

static DateTimeType utcdatetime;
static UInt32       utcseconds;

/*//Julian Date at 0 H UT*/
static double JULD;

/*//middle Greewich-startime at 0 H UT*/
static double GMST0;

/*//current middle greenwich startime*/
static double GMST;

/*//local middle startime for current time*/
static double LMST;

/*//real greenwich startime
//( Greenwich Apparent Sideral Time )*/
static double GAST;

/*//real local Startime for current time
//( Local Apparent Sideral Time )*/
static double LAST;

static double azimuth;
static double altitude;

static double rectaszension;
static double declination;





static double makeAbsoluteDouble( double val ) {
	double d;


	char *p;
	
	printDouble(val,&tmptext[0]);
	p = StrChr(tmptext,'.');
	if(p != NULL){
		p[0] = 0;
	}
	p = StrChr(tmptext,',');
	if(p != NULL){
		p[0] = 0;
	}

	strToDouble((const char *)tmptext,&d);	


/*
	modf(val,&d);
*/
	return(d);
}



static double REST24( double val ) {
	double v,i;

	v = val / 24.0;
	modf(v,&i);
	v = i * 24.0;
	v = val - v;
	
	return( v );
}


static double REST( double val, double t ) {
	double v;

	v = val / t;
	v = makeAbsoluteDouble(v) * t;
	v = val - v;
	
	return( v );
}



/*
//this calculates and sets utcdatetime, utcseconds, GMST0,GMST,LMST,GAST,and at least LAST.

//utcadjust must be the Difference between LocalTime and GMT in seconds
//longitude, the longitude in degrees ( 48.xxx ) for the local place.
//for eldob, we can use scope.utcadjust, and scope.longitude for this
//if hsecs not 0 ( normally ) the time will be set for precession of seconds,
//otherwise for secods + ( hsecs * 1/100) of seconds, so if not 0 it should
//be a number of 0 to 99, indicating hundrets of a second.
*/
static void cCurrentStarTime( Int32 utcadjust, double longitude, Int32 hsecs ) {

	double T,M,J,tmpa,tmpb,tmpc,tmpd, tmpe;
	double p,j,m,jd,b,JD,d;
	int CAL;
	double TN,HOUR,MIN,SEC,TE,Omega,L;
	double Lst,Depsi,Deps,eps0,eps,Depsicoseps;
	
	/*get current time in seconds since 12:00 AM since 1.1.1904 */
	utcseconds = TimGetSeconds() + utcadjust;
	/*convert it to a DateTimeType*/
	TimSecondsToDateTime( utcseconds, &utcdatetime );

	J = utcdatetime.year;
	M = utcdatetime.month;
	T = utcdatetime.day;
	
	tmpa = J;
	tmpb = M / 100;
	tmpc = T / 10000;

	p = tmpa + tmpb + tmpc;
	/*now we have the format JJJ.MMTT*/

	CAL = 2;
	if( p < 1582.1015 ){
		CAL = 1;
	}
	
	if( M > 2 ){
		j = J;	
		m = M;
	}
	else {
		j = J - 1;
		m = M + 12;
	}

	jd = floor( 365.25 * j ) + floor( 30.6001 * ( m + 1)) + T + 1720994.5;
	
	b = 2 - floor( j / 100 ) + floor( floor( j / 100 ) / 4 );

	if( CAL < 2 ){
		JD = jd;
	}
	else {
		JD = jd + b;
	}

	JULD = JD;
	/*
	/////////////////////////////////////////////
	// Now we have the julian date for time 0h UT
	/////////////////////////////////////////////
	*/	
	
	TN = ( JD - 2451545 ) / 36525;

	GMST0 = 24110.54841 + 8640184.812866 * TN + 0.093104 * TN * TN - 0.0000062 * TN * TN * TN;

	GMST0 = REST24( GMST0 / 3600.0 );
	
	if( GMST0 < 0 ){
		GMST0 = GMST0 + 24;		
	}

	/*
	/////////////////////////////////////////////////////
	// Now we have the middle Greewich-startime at 0 H UT
	/////////////////////////////////////////////////////
	*/

	HOUR = utcdatetime.hour;
	MIN  = utcdatetime.minute;
	SEC  = utcdatetime.second;
	
	tmpa = HOUR;
	tmpb = MIN / 60;
	tmpc = SEC / 3600;

	/*hundrets of seconds*/
	if( hsecs > 0 ){
		tmpe = ( 1.0 / 360000.0) * hsecs;
		d    = tmpa + tmpb + tmpc + tmpe; 
	}
	else {	
		d    = tmpa + tmpb + tmpc; 
	}

	tmpd = 1.00273790935 * d;
	GMST = REST24( GMST0 + tmpd);
	
	/*
	/////////////////////////////////////////////////////
	// Now we have the current middle greenwich startime.
	/////////////////////////////////////////////////////
	*/
	
	tmpa = longitude;
	tmpa = tmpa / 15;
	
	LMST = REST24(GMST + tmpa + 24);

	/*	
	//////////////////////////////////////////////////////////
	// Now, we have the local middle startime for current time
	//////////////////////////////////////////////////////////
	*/
	
	tmpb = d / 24;
	tmpa = JD + tmpb - 2451545.0;
	TE   = tmpa / 36525;
	
	Omega = 125.04452 - 1934.136261 * TE + 0.0020708 * TE * TE + TE * TE * TE * 1 / 450000;
	
	L     = 280.4665 + 36000.7698 * TE;
	Lst   = 218.3165 + 481267.8813 * TE;
	
	Omega = rad( Omega );
	L     = rad( L );
	Lst   = rad( Lst );
	
	Depsi = -17.20 * sin( Omega ) - 1.32 * sin( 2 * L ) - 0.23 * sin( 2 * Lst ) + 0.21 * sin( 2 * Omega );
	Deps  = 9.20 * cos( Omega ) + 0.57 * cos( 2 * L ) + 0.10 * cos( 2 * Lst ) - 0.09 * cos( 2 * Omega );
	eps0  = 23.43929111 + ( -46.815 * TE - 0.00059 * TE * TE + 0.001813 * TE * TE * TE ) / 3600;
	eps   = eps0 + Deps / 3600;

	eps = rad( eps );
	
	Depsicoseps = Depsi * cos(eps);

	tmpc  = 15.0 * 3600;	
	GAST  = GMST + Depsicoseps / tmpc;

	/*	
	////////////////////////////////////////////////
	//Now, we have the real greenwich startime GAST
	//( Greenwich Apparent Sideral Time )
	////////////////////////////////////////////////
	*/
	
	
	tmpc = 15.0 * 3600;
	LAST = LMST + Depsicoseps / tmpc;
	if( LAST < 0 ) {
		LAST = LAST + 24;
	}

	/*
	////////////////////////////////////////////////////////
	// Now, we have the real local Startime for current time
	////////////////////////////////////////////////////////
	*/
}


/*
//before using this function CalcCurrentStarTime must
//be called before.
//calculates Rektazension/Declination-kooordinates to azimuth/altitude-koordinates for the
//local time, and place.
//This sets the double-variables azimuth, and altitude ( as floatnumber of degrees )
//ra must be given as hours
//and dec as degrees.
//for latitude, for program eldob the scope.latitude can be used
*/
static void cAzAlt( double latitude, double ra, double dec ) {
	double a,d,Theta,Q,f,Zlr,Nnr,at;
	
	a = ra;
	d = dec;
	
	/*//Stundenwinkel*/
	Theta = REST24( LAST - a + 24);
	/*Theta = LAST - a + 24;*/
	Q     = Theta * 15;
	
	f = rad (latitude);
	d = rad (d);
	Q = rad ( Q );
	
	altitude = asin( sin(f) * sin(d) + cos(f) * cos(d) * cos(Q));
	altitude = deg( altitude );
	
	Zlr = sin(Q);

	Nnr = (cos(Q) * sin(f) - tan(d) * cos(f));
	at  = deg( atan( Zlr / Nnr ) );
	
	azimuth = 0;

/*
	if( ( Nnr < 0 ) && ( Zlr != 0 ) ) {
		azimuth = 180.0 + at;
	}

	if( ( Nnr > 0 ) && ( Zlr > 0 ) ) {
		azimuth = at;
	}

	if( ( Nnr > 0 ) && ( Zlr < 0 ) ) {
		azimuth = 360.0 + at;
	}

	if( ( Nnr == 0 ) && ( Zlr == 0 ) ) {
		azimuth = 0.0;
	}
*/


		
	if( ( Nnr < 0 ) && ( Zlr < 0 )){
		azimuth = at;
		return;
	}

	if( ( Nnr < 0 ) && ( Zlr > 0 )){
		azimuth = 360 + at;
		return;
	}

	if( ( Nnr > 0 ) && ( Zlr != 0 )){
		azimuth = 180 + at;
		return;
	}


}




/*
//before using this function CalcCurrentStarTime must
//be called before.
//calculates Azimut/Altitude-kooordinates to Rectazsension/declination-koordinates using        
//local time, and place.
//This sets the double-variables rectazsension, and declination ( as floatnumber of hours/degrees )
//az, and alt must be given as floatnumber of degrees.
//for latitude, for program eldob the scope.latitude can be used
*/
static void cRaDec( double latitude, double az, double alt ) {


	double _lat,_az,_alt,_last,_z,_dec,_ra,_stw,_stw1,_stw2;
			
	_lat 	= rad( latitude );
	_az 	= rad( az );
	_alt 	= rad( alt );
	_last   = rad( LAST * 15 );
	_z		= rad( 90.0 - alt );

	_dec	= asin(  sin( _alt ) * sin( _lat ) +  cos( _alt ) * cos( _lat ) * cos( _az )  );

/*
	_stw 	= asin( -sin( _az ) * cos( _alt )   /  cos( _dec )  );
	_stw 	= asin( sin( _az ) * cos( _alt )   /  cos( _dec )  );
	_stw	= asin(  ( -sin( _z ) * sin( _az ) ) / cos( _dec ) );
	_stw	= acos( ( cos( _z ) * cos( _lat ) - sin( _z ) * sin( _lat ) * cos( _az ) ) / cos( _dec ) );
*/	

	_stw1	= asin( -sin( _az ) * cos( _alt )   /  cos( _dec )  );
	
	_stw2	= cos( _z ) * cos( _lat ) - sin( _z ) * sin( _lat ) * cos( _az );

	if( _stw1 >= 0.0 ){
		_stw	= acos( _stw2 / cos( _dec ) );
	}
	else {
		_stw	= -acos( _stw2 / cos( _dec ) );
	}

	_ra		= _last - _stw;

	rectaszension = deg( _ra ) / 15;
	
	if( rectaszension < 0.0 ){
		rectaszension = 24.0 + rectaszension;
	}
	
	if( rectaszension >= 23.999999 ){
		rectaszension = rectaszension - 24.0;
	}

	declination   = deg( _dec );
	
}