#include "astro.h"


double REST( double val, double t );
void cCurrentStarTime( void );



/*
  This calculates and sets GAST,and LAST.
  YEAR,MONTH,DAY,MINUTE,SECOND,MSECOND
  must contain correct time before.
  Also LONGITUTE must contain correct value.
*/
void cCurrentStarTime( void ) {
	double TN,TE,Omega,L;
	double Lst,Depsi,Deps,eps0,eps,Depsicoseps;
	double GMST0,GMST,LMST;
	int CAL;

	double tmpa,tmpb,tmpc,tmpe;//tmpd
	double p,j,m,jd,b,JD,d;
	
	tmpa = (double)YEAR;
	tmpb = (double)MONTH / 100;
	tmpc = (double)DAY / 10000;

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

	CAL = 2;
	if( p < 1582.1015 ){
		CAL = 1;
	}

	if( MONTH > 2 ){
		j = (double)YEAR;	
		m = (double)MONTH;
	}
	else {
		j = (double)YEAR - 1;
		m = (double)MONTH + 12;
	}

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


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


	/*
	/////////////////////////////////////////////
	// 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 = REST( GMST0 / 3600.0, 24.0 );
	
	if( GMST0 < 0 ){
		GMST0 = GMST0 + 24;		
	}


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

	
	tmpa = (double)HOUR;
	tmpb = (double)MINUTE / 60;
	tmpc = (double)SECOND / 3600;

	/*milliseconds*/

	if( MSECOND > 0 ){
		tmpe = ( 1.0 / 3600000.0) * (double)MSECOND;
		d    = tmpa + tmpb + tmpc + tmpe; 
	}
	else {	
		d    = tmpa + tmpb + tmpc; 
	}

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




	/*	
	//////////////////////////////////////////////////////////
	// Now, we have the local middle startime for current time
	//////////////////////////////////////////////////////////
	*/
	

//	tmpb = d / 24;
//	tmpa = JD + tmpb - 2451545.0;
//	TE   = tmpa / 36525;
	TE	 = ( JD + ( d / 24 ) - 2451545.0 ) / 36525;


	
	Omega = rad( 125.04452 - 1934.136261 * TE + 0.0020708 * TE * TE + TE * TE * TE * 1 / 450000 );	
	L     = rad( 280.4665 + 36000.7698 * TE );
	Lst   = rad( 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   = rad( eps0 + Deps / 3600 );
//	eps = rad( eps );

	
	Depsicoseps = Depsi * cos(eps);



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

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

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


double REST( double val, double t ) {
	double x;
	modf( val / t ,&x);
	return( val - ( x * t ));
}




/*
  Calculates Rektazension/Declination-kooordinates to azimuth/altitude-koordinates
  for the local time, and place.
  This sets KOOR_AZ, KOOR_ALT ( as floatnumbers of degrees )
  KOOR_RA must be set as hours
  KOOR_DEC dec as degrees.
  the current local startime ( LAST ), must be set before
  by calling cCurrentStartTime
  Also LATITUDE must contain correct value.
*/
void cAzAlt( void ) {

	double a,d,Theta,Q,f,Zlr,Nnr,at,altitude,tmpa,tmpb,tmpc;
		
	a = KOOR_RA;
	d = KOOR_DEC;

	/*Stundenwinkel*/
	
	tmpa  = LAST - a + 24.0;
	tmpb = tmpa / 24.0;
	modf(tmpb,&tmpc);
	tmpb = tmpc * 24.0;
	Theta = tmpa - tmpb;

	Q     = Theta * 15.0;
	
	f = rad ( LATITUDE );
	d = rad ( d );
	Q = rad ( Q );
	
	altitude = asin( sin(f) * sin(d) + cos(f) * cos(d) * cos(Q));
	altitude = deg( altitude );
	KOOR_ALT = altitude;

	Zlr = sin(Q);

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

	KOOR_AZ = 0;

/*q3*/
	if( ( Nnr < 0.0 ) && ( Zlr < 0.0 )){
		KOOR_AZ = at;
		return;
	}

/*q4*/
	if( ( Nnr < 0.0 ) && ( Zlr > 0.0 )){
		KOOR_AZ = ( 360.0 + at );
		return;
	}

/*q1 oder q2*/
	if( ( Nnr > 0.0 ) && ( Zlr != 0.0 )){
		KOOR_AZ = ( 180.0 + at );
		return;
	}

}







/*
  Calculates Azimut/Altitude-kooordinates to Rectaszension/Declination-koordinates
  for the local time, and place.
  This KOOR_RA, KOOR_DEC ( as floatnumbers of Hours/degrees )
  KOOR_AZ, and KOOR_ALT must be set as floatnumber of degrees.
  The current startime LAST must be set before by calling cCurrentStarTime, also
  the LATITUDE.
 */
void cRaDec( void ) {

	double _lat,_az,_alt,_last,_z,_dec,_ra,_stw,_stw1,_stw2;
			
	_lat 	= rad( LATITUDE );
	_az 	= rad( KOOR_AZ );
	_alt 	= rad( KOOR_ALT );
	_last   = rad( LAST * 15 );
	_z		= rad( 90.0 - KOOR_ALT );

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

	_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;

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

	KOOR_DEC   = deg( _dec );
}





/*
	converts a RA-Koordinate ( floatnumber of hours )
	to a long value, representing the value as 1/10 arcseconds
*/
long dHourslArc10( double val ) {
	double x,y;
	x = val * 15 * 36000;
	modf( x , &y );	
	return((long)y);
}


/*
	converts a AZ/ALT, or DEC-Koordinate ( floatnumber of degrees )
	to a long value, representing the value as 1/10 arcseconds.
*/
long dDegreeslArc10( double val ) {
	double x,y;
	x = val * 36000;
	modf( x , &y );	
	return((long)y);
}



/*
	converts a long, representing a RA-Koordinate as 1/10 of arcseconds
	to its aequevalent floatnumber of hours, - stored as double.
*/
double lArc10dHours( long val ) {
	double x;
	x = (double)val / (double)36000.0 / (double)15.0;
	return(x);
}




/*
	converts a AZ/ALT, or DEC-Koordinate, represented as 1/10 of arcseconds,
	to its aequevalent flotnumber of degrees, - stored as double.
*/
double lArc10dDegrees( long val ) {
	double x;
	x = (double)val / (double)36000.0;
	return(x);
}



/*
	val must representing a value of 1/10 arcseconds.
	returns a number, representing
	the 360 degrees-sector, which is fullfilled by the koordinate.
	0 = 0         -  3.239.999
	1 = 3.240.000 -  6.479.999
	2 = 6.480.000 -  9.719.999
	3 = 9.720.000 - 12.959.999
*/
int arc10Sector( long val ) {
	
	if( ( val >= 0 ) && ( val < 3240000 ) ){
		return(0);
	}

	if( ( val >= 3240000 ) && ( val < 6480000 ) ){
		return(1);
	}

	if( ( val >= 6480000 ) && ( val < 9720000 ) ){
		return(2);
	}

	return(3);
}




/*
	val must representing a floatnumber of degrees.
	returns a number, representing
	the 360 degrees-sector, which is fullfilled by the koordinate.
	0 = 000 - 089.999999
	1 = 090 - 179.999999
	2 = 180 - 269.999999
	3 = 270 - 359.999999
*/
int degreesSector( double val ) {

	if( ( val >= 0.0 ) && ( val < 90.0 ) ){
		return(0);
	}

	if( ( val >= 90.0 ) && ( val < 180.0 ) ){
		return(1);
	}

	if( ( val >= 180.0 ) && ( val < 270.0 ) ){
		return(2);
	}

	return(3);

}









void testAstro(void)
{
	char tmp[42];
	int ct;
		
	LATITUDE	= 48.766667;
	LONGITUDE	=  9.183333;
	
	YEAR		= 2003;
	MONTH		= 4;
	DAY			= 11;
	HOUR		= 20;
	MINUTE		= 14;
	SECOND		= 2;
	MSECOND		= 10;

	KOOR_RA		= 19.132456;
	KOOR_DEC	= 64.324567;

	ct = 0;


	while(ct < 2 ){
	
	cCurrentStarTime();
	cAzAlt();

	ct++;
	}

	cRaDec();

	
	sprintf(tmp,"%016.6f", KOOR_RA);
	lcd_puts(0,0,tmp);

	sprintf(tmp,"%016.6f",KOOR_DEC);
	lcd_puts(1,0,tmp);

	sprintf(tmp,"%016.6f",KOOR_AZ);
	lcd_puts(2,0,tmp);

	sprintf(tmp,"%016.6f",KOOR_ALT);
	lcd_puts(3,0,tmp);


}