+#include <cmath>
+
+/*----------------------------------------------------------------------------*/
+double factorial_div( double value, int x)
+{
+ if(!x)
+ return 1;
+ else
+ {
+ while( x > 1)
+ {
+ value = value / x--;
+ }
+ }
+ return value;
+}
+
+/*----------------------------------------------------------------------------*/
+double powerd( double x, int y)
+{
+ int i=0;
+ double ans=1.0;
+
+ if(!y)
+ return 1.000;
+ else
+ {
+ while( i < y)
+ {
+ i++;
+ ans = ans * x;
+ }
+ }
+ return ans;
+}
+
+/*----------------------------------------------------------------------------*/
+double SIN( double x)
+{
+ int i=0;
+ int j=1;
+ int sign=1;
+ double y1 = 0.0;
+ double diff = 1000.0;
+
+ if (x < 0.0)
+ {
+ x = -1 * x;
+ sign = -1;
+ }
+
+ while ( x > 360.0*M_PI/180)
+ {
+ x = x - 360*M_PI/180;
+ }
+
+ if( x > (270.0 * M_PI / 180) )
+ {
+ sign = sign * -1;
+ x = 360.0*M_PI/180 - x;
+ }
+ else if ( x > (180.0 * M_PI / 180) )
+ {
+ sign = sign * -1;
+ x = x - 180.0 *M_PI / 180;
+ }
+ else if ( x > (90.0 * M_PI / 180) )
+ {
+ x = 180.0 *M_PI / 180 - x;
+ }
+
+ while( powerd( diff, 2) > 1.0E-16 )
+ {
+ i++;
+ diff = j * factorial_div( powerd( x, (2*i -1)) ,(2*i -1));
+ y1 = y1 + diff;
+ j = -1 * j;
+ }
+ return ( sign * y1 );
+}
+
+/*----------------------------------------------------------------------------*/
+double COS(double x)
+{
+ return SIN(90 * M_PI / 180 - x);
+}
+
+/*----------------------------------------------------------------------------*/
+double ATAN( double x)
+{
+ int i=0; /* counter for terms in binomial series */
+ int j=1; /* sign of nth term in series */
+ int k=0;
+ int sign = 1; /* sign of the input x */
+ double y = 0.0; /* the output */
+ double deltay = 1.0; /* the value of the next term in the series */
+ double addangle = 0.0; /* used if arctan > 22.5 degrees */
+
+ if (x < 0.0)
+ {
+ x = -1 * x;
+ sign = -1;
+ }
+
+ while( x > 0.3249196962 )
+ {
+ k++;
+ x = (x - 0.3249196962) / (1 + x * 0.3249196962);
+ }
+
+ addangle = k * 18.0 *M_PI/180;
+
+ while( powerd( deltay, 2) > 1.0E-16 )
+ {
+ i++;
+ deltay = j * powerd( x, (2*i -1)) / (2*i -1);
+ y = y + deltay;
+ j = -1 * j;
+ }
+ return (sign * (y + addangle) );
+}
+
+double ASIN(double x)
+{
+ return 2 * ATAN( x / (1 + std::sqrt(1.0 - x*x)));
+}
+
+double Radians( double number )
+{
+ return number*M_PI/180;
+}
+
+double Deg( double number )
+{
+ return number*180/M_PI;
+}
+
+double Rev( double number )
+{
+ return number - std::floor( number / 360.0 ) * 360;
+}
+
+double calcElevation( double SatLon, double SiteLat, double SiteLon, int Height_over_ocean = 0 )
+{
+ double a0=0.58804392,
+ a1=-0.17941557,
+ a2=0.29906946E-1,
+ a3=-0.25187400E-2,
+ a4=0.82622101E-4,
+
+ f = 1.00 / 298.257, // Earth flattning factor
+
+ r_sat=42164.57, // Distance from earth centre to satellite
+
+ r_eq=6378.14, // Earth radius
+
+ sinRadSiteLat=SIN(Radians(SiteLat)),
+ cosRadSiteLat=COS(Radians(SiteLat)),
+
+ Rstation = r_eq / ( std::sqrt( 1.00 - f*(2.00-f)*sinRadSiteLat*sinRadSiteLat ) ),
+
+ Ra = (Rstation+Height_over_ocean)*cosRadSiteLat,
+ Rz= Rstation*(1.00-f)*(1.00-f)*sinRadSiteLat,
+
+ alfa_rx=r_sat*COS(Radians(SatLon-SiteLon)) - Ra,
+ alfa_ry=r_sat*SIN(Radians(SatLon-SiteLon)),
+ alfa_rz=-Rz,
+
+ alfa_r_north=-alfa_rx*sinRadSiteLat + alfa_rz*cosRadSiteLat,
+ alfa_r_zenith=alfa_rx*cosRadSiteLat + alfa_rz*sinRadSiteLat,
+
+ El_geometric=Deg(ATAN( alfa_r_zenith/std::sqrt(alfa_r_north*alfa_r_north+alfa_ry*alfa_ry))),
+
+ x = std::fabs(El_geometric+0.589),
+ refraction=std::fabs(a0+a1*x+a2*x*x+a3*x*x*x+a4*x*x*x*x),
+
+ El_observed = 0.00;
+
+ if (El_geometric > 10.2)
+ El_observed = El_geometric+0.01617*(COS(Radians(std::fabs(El_geometric)))/SIN(Radians(std::fabs(El_geometric))) );
+ else
+ El_observed = El_geometric+refraction;
+
+ if (alfa_r_zenith < -3000)
+ El_observed=-99;
+
+ return El_observed;
+}
+
+double calcAzimuth(double SatLon, double SiteLat, double SiteLon, int Height_over_ocean=0)
+{
+ double f = 1.00 / 298.257, // Earth flattning factor
+
+ r_sat=42164.57, // Distance from earth centre to satellite
+
+ r_eq=6378.14, // Earth radius
+
+ sinRadSiteLat=SIN(Radians(SiteLat)),
+ cosRadSiteLat=COS(Radians(SiteLat)),
+
+ Rstation = r_eq / ( std::sqrt( 1 - f*(2-f)*sinRadSiteLat*sinRadSiteLat ) ),
+ Ra = (Rstation+Height_over_ocean)*cosRadSiteLat,
+ Rz = Rstation*(1-f)*(1-f)*sinRadSiteLat,
+
+ alfa_rx = r_sat*COS(Radians(SatLon-SiteLon)) - Ra,
+ alfa_ry = r_sat*SIN(Radians(SatLon-SiteLon)),
+ alfa_rz = -Rz,
+
+ alfa_r_north = -alfa_rx*sinRadSiteLat + alfa_rz*cosRadSiteLat,
+
+ Azimuth = 0.00;
+
+ if (alfa_r_north < 0)
+ Azimuth = 180+Deg(ATAN(alfa_ry/alfa_r_north));
+ else
+ Azimuth = Rev(360+Deg(ATAN(alfa_ry/alfa_r_north)));
+
+ return Azimuth;
+}
+
+double calcDeclination( double SiteLat, double Azimuth, double Elevation)
+{
+ return Deg( ASIN(SIN(Radians(Elevation)) *
+ SIN(Radians(SiteLat)) +
+ COS(Radians(Elevation)) *
+ COS(Radians(SiteLat)) +
+ COS(Radians(Azimuth))
+ )
+ );
+}
+
+double calcSatHourangle( double SatLon, double SiteLat, double SiteLon )
+{
+ double Azimuth=calcAzimuth(SatLon, SiteLat, SiteLon ),
+ Elevation=calcElevation( SatLon, SiteLat, SiteLon ),
+
+ a = - COS(Radians(Elevation)) * SIN(Radians(Azimuth)),
+
+ b = SIN(Radians(Elevation)) * COS(Radians(SiteLat))
+ -
+ COS(Radians(Elevation)) * SIN(Radians(SiteLat)) * COS(Radians(Azimuth)),
+
+// Works for all azimuths (northern & southern hemisphere)
+ returnvalue = 180 + Deg(ATAN(a/b));
+
+ if ( Azimuth > 270 )
+ {
+ returnvalue = ( (returnvalue-180) + 360 );
+ if (returnvalue>360)
+ returnvalue = 360 - (returnvalue-360);
+ }
+
+ if ( Azimuth < 90 )
+ returnvalue = ( 180 - returnvalue );
+
+ return returnvalue;
+}
projects / vuplus_dvbapp / commitdiff