#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <stdarg.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <fcntl.h>
#include <ctype.h>
#include <time.h>
#include <sys/time.h>
#include <sys/timeb.h>
#include <unistd.h>
#include <malloc.h>
#include <math.h>
#include <termios.h>
#include <dirent.h>
#include <signal.h>



#define STARDAYSECONDS	86164.0
#define HOURSECONDS		3600.0
#define HOUR24SECONDS   86400.0


#ifndef PI
#define PI 	3.141592653589793
#endif

#ifndef M_PI
/*
# define M_PI		3.14159265358979323846
*/
# define M_PI		3.1415926535897932384626433832795029L
#endif



#define rad(x)  	((x)*PI/180.0)
#define deg(x)  	((x)*180.0/PI )
#define hdeg(x) 	((x)*15.0)
#define degh(x) 	((x)/15.0)
#define hrad(x)     rad(hdeg(x))

#define arctohs(a,s) ((a)/(s))
#define hstoarc(a,s) ((a)*(s))

#define hourstosecs(x)		((x)*3600.0)
#define hourstomillis(x)	((x)*3600.0*1000.0)
#define secstohours(x)		((x)/3600.0)
#define millistohours(x)	((x)/3600.0/1000.0)

#define circle_radius(x,y)	(sqrt(((x)*(x))+((y)*(y))));
#define circle_circumf(r)   ((r)*2.0*PI)
#define circle_angle(b,u)	(((b)/(u))*360.0)
	

#define ARC0 	0.0
#define ARC45   162000.0
#define ARC90	324000.0
#define ARC135	486000.0
#define ARC180	648000.0
#define ARC225  810000.0
#define ARC270	972000.0
#define ARC315	1134000.0
#define ARC360	1296000.0



unsigned short GMTYEAR;
unsigned short GMTMONTH;
unsigned short GMTDAY;
unsigned short GMTHOUR;
unsigned short GMTMINUTE;
unsigned short GMTSECOND;
unsigned short GMTMSSECOND;



double LAST;
double GAST;
double JULD;
double GMST0;
double GMST;
double LMST;
double LATITUDE;
double LONGITUDE;

double KOOR_RA;
double KOOR_DEC;
double KOOR_AZ;
double KOOR_ALT;



#include "util.c"


/*------------------------------------------------------------------------------------------*/
/*------------------------------------------------------------------------------------------*/
/*------------------------------------------------------------------------------------------*/
/*------------------------------------------------------------------------------------------*/
/*------------------------------------------------------------------------------------------*/


double makeAbsoluteDouble( double val ) {
	double d;
	modf(val, &d);
	return(d);
}



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

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



/*
  This calculates and sets GMST0,GMST,LMST,GAST,and at least LAST.
*/
void cCurrentStarTime( void ) {

	double tmpa,tmpb,tmpc,tmpd,tmpe;
	double p,j,m,jd,b,JD,d;
	int CAL;
	double TN,TE,Omega,L;
	double Lst,Depsi,Deps,eps0,eps,Depsicoseps;
	
	tmpa = (double)GMTYEAR;
	tmpb = (double)GMTMONTH / 100;
	tmpc = (double)GMTDAY / 10000;

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

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

	jd = floor( 365.25 * j ) + floor( 30.6001 * ( m + 1)) + (double)GMTDAY + 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 = REST( GMST0 / 3600.0, 24.0 );
	
	if( GMST0 < 0 ){
		GMST0 = GMST0 + 24;		
	}

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

	
	tmpa = (double)GMTHOUR;
	tmpb = (double)GMTMINUTE / 60;
	tmpc = (double)GMTSECOND / 3600;

	/*milliseconds*/

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

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

	/*	
	//////////////////////////////////////////////////////////
	// 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
	////////////////////////////////////////////////////////
	*/
}








/*
  Calculates Rektazension/Declination-kooordinates to azimuth/altitude-koordinates
  for the local time, and place.
  This sets KOOR_AZ, KOOR_ALT ( as doublenumbers 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
*/
void cAzAlt( void ) {

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

	/*//Stundenwinkel*/
	Theta = REST( LAST - a + 24.0, 24.0 );
	/*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 );
	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 ) && ( Zlr < 0 )){
		KOOR_AZ = at;
		return;
	}

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

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

}




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

	double _lat,_az,_alt,_last,_z,_dec,_ra,_stw,_stw1,_stw2, rectaszension,declination;
			
	_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;

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

	declination   = deg( _dec );

	KOOR_RA  = rectaszension;
	KOOR_DEC = declination;
}





char str1[84];
char str2[84];

int main(int argc, char *argv[])
{
	struct tm * t;
	time_t secs;
				
	LATITUDE	= 48.766667;
	LONGITUDE	=  9.183333;

	secs = time(NULL);
	t = localtime(&secs);
	
	GMTYEAR		= t->tm_year + 1900;
	GMTMONTH	= t->tm_mon + 1;
	GMTDAY		= t->tm_mday;
	GMTHOUR		= t->tm_hour;
	GMTMINUTE	= t->tm_min;
	GMTSECOND	= t->tm_sec;
	GMTMSSECOND = 10;


	cCurrentStarTime();

	KOOR_RA  = 19.132456;
	KOOR_DEC = 64.324567;

	fprintf(stdout,"\nLAST: %f\n",LAST);

/*
	printRa(KOOR_RA,&str1[0]);printDec(KOOR_DEC,&str2[0]);	
	fprintf(stdout,"%020.6f  %020.6f  %s  %s\n", KOOR_RA, KOOR_DEC, str1, str2 );	
	cAzAlt();

	printAz(KOOR_AZ,&str1[0]);printAlt(KOOR_ALT,&str2[0]);	
	fprintf(stdout,"%020.6f  %020.6f  %s  %s\n", KOOR_AZ, KOOR_ALT, str1, str2 );

	cRaDec();

	printRa(KOOR_RA,&str1[0]);printDec(KOOR_DEC,&str2[0]);	
	fprintf(stdout,"%020.6f  %020.6f  %s  %s\n", KOOR_RA, KOOR_DEC, str1, str2 );
*/
	exit(0);
}