#include "elastro.h" #include "_calc.c" /* the following variables will be set by calling cCurrentStarTime*/ /*//the time in seconds when function is called and the adjusted time*/ time_t cCSTcurseconds; long cCSTmillisecs; time_t cCSTadjcurseconds; /*//Julian Date at 0 H UT*/ double JULD; /*//middle Greewich-startime at 0 H UT*/ double GMST0; /*//current middle greenwich startime*/ double GMST; /*//local middle startime for current time*/ double LMST; /*//real greenwich startime //( Greenwich Apparent Sideral Time )*/ double GAST; /*//real local Startime for current time //( Local Apparent Sideral Time )*/ double LAST; 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 //for summertime at Stuttgart this is f.Expl. -7200, -3600 for wintertime //longitude, the longitude in degrees ( 48.xxx ) for the local place. //also the latitude. //seconds must be the number seconds for the current local time. //if seconds <= 0 current time will be calculated automatically, ( to a precession of //milliseconds ) otherwise this value ( additionalyy hsecs ) is used. //if hsecs not 0 ( normally ) the time will be set for precession of seconds, //otherwise for secods + ( hsecs * 1/1000) of seconds, so if not 0 it should //be a number of 0 to 999, indicating milliseconds. //( for unix you can use ftime for getting seconds, and milliseconds of current time ) //function will be set the timevariables above, and the LAST value will be addionally //returned. */ double cCurrentStarTime( double longitude , double latitude, long seconds, long hsecs, long utcadjust ) { 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; long utcseconds; time_t sec; struct tm *t; long ilhsecs; struct timeb tpa; if( seconds <= 0) { ftime(&tpa); cCSTcurseconds = tpa.time; utcseconds = tpa.time + utcadjust; } else { utcseconds = seconds + utcadjust; cCSTcurseconds = seconds; } cCSTadjcurseconds = utcseconds; sec = (time_t)utcseconds; t = localtime(&sec); SEC = t->tm_sec; MIN = t->tm_min; HOUR = t->tm_hour; T = t->tm_mday; M = t->tm_mon + 1; J = t->tm_year + 1900; 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 = REST( GMST0 / 3600.0, 24.0 ); if( GMST0 < 0 ){ GMST0 = GMST0 + 24; } /* ///////////////////////////////////////////////////// // Now we have the middle Greewich-startime at 0 H UT ///////////////////////////////////////////////////// */ tmpa = HOUR; tmpb = MIN / 60; tmpc = SEC / 3600; /*milliseconds*/ if(seconds <= 0 ){ cCSTmillisecs = tpa.millitm; if( tpa.millitm > 0 ){ tmpe = ( 1.0 / 3600000.0) * (double)tpa.millitm; d = tmpa + tmpb + tmpc + tmpe; } else { d = tmpa + tmpb + tmpc; } } else { if( hsecs > 0 ){ cCSTmillisecs = hsecs; tmpe = ( 1.0 / 3600000.0) * (double)hsecs; d = tmpa + tmpb + tmpc + tmpe; } else { cCSTmillisecs = 0; 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 //////////////////////////////////////////////////////// */ return( LAST ); } /* //calculates Rektazension/Declination-kooordinates to azimuth/altitude-koordinates for the //local time, and place. //This sets the struct AZALT ( as floatnumbers of degrees ) //ra must be given as hours //and dec as degrees. //stt must be the current startime ( LAST ) */ void cAzAlt( double latitude, double ra, double dec, double stt, struct AZALT *azalt ) { double a,d,Theta,Q,f,Zlr,Nnr,at,altitude; a = ra; d = dec; /*//Stundenwinkel*/ Theta = REST( stt - a + 24.0, 24.0 ); /*Theta = stt - 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 ); azalt->alt = altitude; Zlr = sin(Q); Nnr = (cos(Q) * sin(f) - tan(d) * cos(f)); at = deg( atan( Zlr / Nnr ) ); azalt->az = 0; /* if( ( Nnr < 0 ) && ( Zlr != 0 ) ) { azalt->az = 180.0 + at; } Q1 if( ( Nnr > 0 ) && ( Zlr > 0 ) ) { azalt->az = at; } Q2 if( ( Nnr > 0 ) && ( Zlr < 0 ) ) { azalt->az = 360.0 + at; } if( ( Nnr == 0 ) && ( Zlr == 0 ) ) { azalt->az = 0.0; } */ /*Q3*/ if( ( Nnr < 0 ) && ( Zlr < 0 )){ azalt->az = at; return; } /*q4*/ if( ( Nnr < 0 ) && ( Zlr > 0 )){ azalt->az = 360 + at; return; } /*q1 oder q2*/ if( ( Nnr > 0 ) && ( Zlr != 0 )){ azalt->az = 180 + at; return; } } /* //calculates Azimut/Altitude-kooordinates to Rectaszension/Declination-koordinates for the //local time, and place. //This sets the struct RADEC ( as floatnumbers of Hours/degrees ) //az, and alt must be given as floatnumber of degrees. //stt must be the current startime ( LAST ) */ void cRaDec( double latitude, double az, double alt, double stt, struct RADEC * radec ) { double _lat,_az,_alt,_last,_z,_dec,_ra,_stw,_stw1,_stw2, rectaszension,declination; _lat = rad( latitude ); _az = rad( az ); _alt = rad( alt ); _last = rad( stt * 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 ); radec->ra = rectaszension; radec->dec = declination; } /* While refraction, -the seeing altitude is a greater value then the real altitude. For the zenit it is 0. For the horizon, and 1013 hectopascal airpressure, and +10 degrees ( Celsius ) temperature, it is 35'24" ( 2124 arcseconds ). The value of refraction increases with increasing airpressure, and decreasing temperature. It decreases with decreasing airpressure, and increasing temperature. Example: average values for +10 degees /1013 hPa virtual altitude real altitude difference(arcsecs) 90 90 0 (0) 60 59.99 0.01(36) 30 29.97 0.03(108) 20 19.96 0.04(144) 10 9.91 0.09(324) 5 4.84 0.16(576) 4 3.8 0.2 (720) 3 2.76 0.24(864) 2 1.7 0.3 (1080) 1 0.59 0.41(1476) 0.5 0.01 0.49(1764) 0 -0.59 0.59(2124) For our refraction calculations we will not need any greater, time-eating, searching-algorithm's. So we have to produce a list at programstart, which values will be get during running. */ void * cRfValues = NULL; /*This should called whenever a program is end. It makes sure, that memory is freed at time. */ void cRfFreeValues( void ){ if(cRfValues != NULL){ free(cRfValues); cRfValues = NULL; } } /* Produces list of refractionvalues, which will be used by cAltRfToScope, and cAltRfFromScope. temperature should be the current temperature in degrees C airpressue the current airpressue in millibars. if init set to 0, no list will be produces, and no calculation of refraction will be done. */ int cRfInitValues( double temperature, double airpressure, int init ) { double v,vr; short _v; long dummy1,dummy2; short *p; long _ct; double ct; cRfFreeValues(); if( ! init ){ return(1); } cRfValues = malloc( 324002 * sizeof( short ) ); if(cRfValues == NULL){ return(0); } p = (short *)cRfValues; _ct = 0; while( _ct <= 324000 ){ p[_ct] = 0; _ct++; } dummy1 = 0L; dummy2 = 0L; _ct = 0; ct = 0.0; while(_ct < 324000 ){ refract ( airpressure, temperature, rad( ct / 3600.0 ), &v); vr = deg(v) * 3600.0; _v = (short)makeAbsoluteDouble( vr - ct ); p[_ct] = _v; ct = ct + 1.0; _ct++; } return(1); } /* adjust alt-value (arcseconds) with calculated refraction returnes new alt-value of arcseconds, useable for the scope. */ double cAltRfToScope( double arcsec ) { double d; long l; unsigned short *p; unsigned short val; if( cRfValues == NULL ){ return(arcsec); } if( (arcsec >= 324000.0 ) || ( arcsec < 0.0 ) ){ return(arcsec); } p = (unsigned short *)cRfValues; l = (long)makeAbsoluteDouble( arcsec ); val = p[l]; d = arcsec + val; return( d ); } /* recalculates refraction from a alt-value (arcseconds ) which contains calculated refraction returnes recalculated real sky-value of arcseonds */ double cAltRfFromScope( double arcsec ) { double d; long l; unsigned short *p; unsigned short val; if( cRfValues == NULL ){ return(arcsec); } if( (arcsec >= 324000.0 ) || ( arcsec < 0.0 ) ){ return(arcsec); } p = (unsigned short *)cRfValues; l = (long)makeAbsoluteDouble( arcsec ); val = p[l]; d = arcsec - val; return( d ); }