#include "astro.h"
volatile double LATITUDE;
volatile double LONGITUDE;
volatile double KOOR_RA;
volatile double KOOR_DEC;
volatile double KOOR_AZ;
volatile double KOOR_ALT;
volatile double GAST;
volatile double LAST;
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);
}