#include <math.h>
#include "zeitgleichung.h"

double In2Pi( double x ) {

	int n = (int) ( x / ( 2 * M_PI ) );
	x -= n * 2 * M_PI;
	if ( x < 0 ) x += ( 2 * M_PI );
 
	return x;
}

double BerechneZeitgleichung( double *Deklination, double T ) {

	double M, L, Ekliptik, RA, RA_Mittel, dRA;
	
	// M - Winkelweg seit letztem Periheldurchgang, wenn die Erde sich auf einer Kreisbahn bewegen würde
	M = In2Pi( 2 * M_PI * ( 0.993133 + 99.997361 * T ) );
	
	L = In2Pi( 2 * M_PI * ( 0.7859453 + M / ( 2 * M_PI ) + ( 6893.0 * sin( M ) + 72.0 * sin( 2.0 * M ) + 6191.2 * T ) / 1296.0e3 ) );
	
	// Ekliptik - Neigung der Erdachse 
	// J. Laskar, "Secular terms of classical planetary theories using the results of general theory," Astronomy and Astrophysics 157, 59070 (1986).
	Ekliptik = ( M_PI / 180 ) * ( 23.43929111 + ( -46.8150 * T - 0.00059 * T * T + 0.001813 * T * T * T ) / 3600.0 );
	
	//Projektion auf die Äquatorebene
	RA = atan( tan( L ) * cos( Ekliptik ) );
	
	if ( RA < 0.0 ) RA += M_PI;
	if ( L > M_PI ) RA += M_PI;
	
	RA = 12.0 * RA / M_PI;
	*Deklination = asin( sin( Ekliptik ) * sin( L ) );
	
	RA_Mittel = 18.71506921 + 2400.0513369 * T + ( 2.5862e-5 - 1.72e-9 * T ) * T * T;
	
	// Damit 0 <= RA_Mittel < 24
	RA_Mittel = 12.0 * In2Pi( M_PI * RA_Mittel / 12.0 ) / M_PI;
	
	dRA = RA_Mittel - RA;
	if ( dRA < -12.0 ) dRA += 24.0;
	if ( dRA > 12.0 ) dRA -= 24.0;
	
	dRA *= 1.0027379;
	
	return dRA ;
}