work on rotor support, add usals stuff (rotor stuff not completed yet),

only send diseqc when needed

1 files changed, 257 insertions, 0 deletions

diff --git a/lib/dvb/rotor_calc.cpp b/lib/dvb/rotor_calc.cpp
new file mode 100644
index 00000000..41322360
--- /dev/null
+++ b/lib/dvb/rotor_calc.cpp
@@ -0,0 +1,257 @@
+#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;
+}