fix no more working teletext / ac3
[enigma2.git] / lib / dvb / rotor_calc.cpp
1 #include <cmath>
2
3 /*----------------------------------------------------------------------------*/
4 double factorial_div( double value, int x)
5 {
6         if(!x)
7                 return 1;
8         else
9         {
10                 while( x > 1)
11                 {
12                         value = value / x--;
13                 }
14         }
15         return value;
16 }
17
18 /*----------------------------------------------------------------------------*/
19 double powerd( double x, int y)
20 {
21         int i=0;
22         double ans=1.0;
23
24         if(!y)
25                 return 1.000;
26         else
27         {
28                 while( i < y)
29                 {
30                         i++;
31                         ans = ans * x;
32                 }
33         }
34         return ans;
35 }
36
37 /*----------------------------------------------------------------------------*/
38 double SIN( double x)
39 {
40         int i=0;
41         int j=1;
42         int sign=1;
43         double y1 = 0.0;
44         double diff = 1000.0;
45
46         if (x < 0.0)
47         {
48                 x = -1 * x;
49                 sign = -1;
50         }
51
52         while ( x > 360.0*M_PI/180)
53         {
54                 x = x - 360*M_PI/180;
55         }
56
57         if( x > (270.0 * M_PI / 180) )
58         {
59                 sign = sign * -1;
60                 x = 360.0*M_PI/180 - x;
61         }
62         else if ( x > (180.0 * M_PI / 180) )
63         {
64                 sign = sign * -1;
65                 x = x - 180.0 *M_PI / 180;
66         }
67         else if ( x > (90.0 * M_PI / 180) )
68         {
69                 x = 180.0 *M_PI / 180 - x;
70         }
71
72         while( powerd( diff, 2) > 1.0E-16 )
73         {
74                 i++;
75                 diff = j * factorial_div( powerd( x, (2*i -1)) ,(2*i -1));
76                 y1 = y1 + diff;
77                 j = -1 * j;
78         }
79         return ( sign * y1 );
80 }
81
82 /*----------------------------------------------------------------------------*/
83 double COS(double x)
84 {
85         return SIN(90 * M_PI / 180 - x);
86 }
87
88 /*----------------------------------------------------------------------------*/
89 double ATAN( double x)
90 {
91         int i=0; /* counter for terms in binomial series */
92         int j=1; /* sign of nth term in series */
93         int k=0;
94         int sign = 1; /* sign of the input x */
95         double y = 0.0; /* the output */
96         double deltay = 1.0; /* the value of the next term in the series */
97         double addangle = 0.0; /* used if arctan > 22.5 degrees */
98
99         if (x < 0.0)
100         {
101                 x = -1 * x;
102                 sign = -1;
103         }
104
105         while( x > 0.3249196962 )
106         {
107                 k++;
108                 x = (x - 0.3249196962) / (1 + x * 0.3249196962);
109         }
110
111         addangle = k * 18.0 *M_PI/180;
112
113         while( powerd( deltay, 2) > 1.0E-16 )
114         {
115                 i++;
116                 deltay = j * powerd( x, (2*i -1)) / (2*i -1);
117                 y = y + deltay;
118                 j = -1 * j;
119         }
120         return (sign * (y + addangle) );
121 }
122
123 double ASIN(double x)
124 {
125         return 2 * ATAN( x / (1 + std::sqrt(1.0 - x*x)));
126 }
127
128 double Radians( double number )
129 {
130         return number*M_PI/180;
131 }
132
133 double Deg( double number )
134 {
135         return number*180/M_PI;
136 }
137
138 double Rev( double number )
139 {
140         return number - std::floor( number / 360.0 ) * 360;
141 }
142
143 double calcElevation( double SatLon, double SiteLat, double SiteLon, int Height_over_ocean = 0 )
144 {
145         double  a0=0.58804392,
146                         a1=-0.17941557,
147                         a2=0.29906946E-1,
148                         a3=-0.25187400E-2,
149                         a4=0.82622101E-4,
150
151                         f = 1.00 / 298.257, // Earth flattning factor
152
153                         r_sat=42164.57, // Distance from earth centre to satellite
154
155                         r_eq=6378.14,  // Earth radius
156
157                         sinRadSiteLat=SIN(Radians(SiteLat)),
158                         cosRadSiteLat=COS(Radians(SiteLat)),
159
160                         Rstation = r_eq / ( std::sqrt( 1.00 - f*(2.00-f)*sinRadSiteLat*sinRadSiteLat ) ),
161
162                         Ra = (Rstation+Height_over_ocean)*cosRadSiteLat,
163                         Rz= Rstation*(1.00-f)*(1.00-f)*sinRadSiteLat,
164
165                         alfa_rx=r_sat*COS(Radians(SatLon-SiteLon)) - Ra,
166                         alfa_ry=r_sat*SIN(Radians(SatLon-SiteLon)),
167                         alfa_rz=-Rz,
168
169                         alfa_r_north=-alfa_rx*sinRadSiteLat + alfa_rz*cosRadSiteLat,
170                         alfa_r_zenith=alfa_rx*cosRadSiteLat + alfa_rz*sinRadSiteLat,
171
172                         El_geometric=Deg(ATAN( alfa_r_zenith/std::sqrt(alfa_r_north*alfa_r_north+alfa_ry*alfa_ry))),
173
174                         x = std::fabs(El_geometric+0.589),
175                         refraction=std::fabs(a0+a1*x+a2*x*x+a3*x*x*x+a4*x*x*x*x),
176
177                         El_observed = 0.00;
178
179         if (El_geometric > 10.2)
180                 El_observed = El_geometric+0.01617*(COS(Radians(std::fabs(El_geometric)))/SIN(Radians(std::fabs(El_geometric))) );
181         else
182                 El_observed = El_geometric+refraction;
183
184         if (alfa_r_zenith < -3000)
185                 El_observed=-99;
186
187         return El_observed;
188 }
189
190 double calcAzimuth(double SatLon, double SiteLat, double SiteLon, int Height_over_ocean=0)
191 {
192         double  f = 1.00 / 298.257, // Earth flattning factor
193
194                         r_sat=42164.57, // Distance from earth centre to satellite
195
196                         r_eq=6378.14,  // Earth radius
197
198                         sinRadSiteLat=SIN(Radians(SiteLat)),
199                         cosRadSiteLat=COS(Radians(SiteLat)),
200
201                         Rstation = r_eq / ( std::sqrt( 1 - f*(2-f)*sinRadSiteLat*sinRadSiteLat ) ),
202                         Ra = (Rstation+Height_over_ocean)*cosRadSiteLat,
203                         Rz = Rstation*(1-f)*(1-f)*sinRadSiteLat,
204
205                         alfa_rx = r_sat*COS(Radians(SatLon-SiteLon)) - Ra,
206                         alfa_ry = r_sat*SIN(Radians(SatLon-SiteLon)),
207                         alfa_rz = -Rz,
208
209                         alfa_r_north = -alfa_rx*sinRadSiteLat + alfa_rz*cosRadSiteLat,
210
211                         Azimuth = 0.00;
212
213         if (alfa_r_north < 0)
214                 Azimuth = 180+Deg(ATAN(alfa_ry/alfa_r_north));
215         else
216                 Azimuth = Rev(360+Deg(ATAN(alfa_ry/alfa_r_north)));
217
218         return Azimuth;
219 }
220
221 double calcDeclination( double SiteLat, double Azimuth, double Elevation)
222 {
223         return Deg( ASIN(SIN(Radians(Elevation)) *
224                                                                                                 SIN(Radians(SiteLat)) +
225                                                                                                 COS(Radians(Elevation)) *
226                                                                                                 COS(Radians(SiteLat)) +
227                                                                                                 COS(Radians(Azimuth))
228                                                                                                 )
229                                                 );
230 }
231
232 double calcSatHourangle( double SatLon, double SiteLat, double SiteLon )
233 {
234         double  Azimuth=calcAzimuth(SatLon, SiteLat, SiteLon ),
235                         Elevation=calcElevation( SatLon, SiteLat, SiteLon ),
236
237                         a = - COS(Radians(Elevation)) * SIN(Radians(Azimuth)),
238
239                         b = SIN(Radians(Elevation)) * COS(Radians(SiteLat))
240                                 -
241                                 COS(Radians(Elevation)) * SIN(Radians(SiteLat)) * COS(Radians(Azimuth)),
242
243 // Works for all azimuths (northern & southern hemisphere)
244                         returnvalue = 180 + Deg(ATAN(a/b));
245
246         if ( Azimuth > 270 )
247         {
248                 returnvalue = ( (returnvalue-180) + 360 );
249                 if (returnvalue>360)
250                         returnvalue = 360 - (returnvalue-360);
251         }
252
253         if ( Azimuth < 90 )
254                 returnvalue = ( 180 - returnvalue );
255
256         return returnvalue;
257 }