1 #include <cmath>
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 }
18 /*----------------------------------------------------------------------------*/
19 double powerd( double x, int y)
20 {
21         int i=0;
22         double ans=1.0;
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 }
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;
46         if (x < 0.0)
47         {
48                 x = -1 * x;
49                 sign = -1;
50         }
52         while ( x > 360.0*M_PI/180)
53         {
54                 x = x - 360*M_PI/180;
55         }
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         }
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 }
82 /*----------------------------------------------------------------------------*/
83 double COS(double x)
84 {
85         return SIN(90 * M_PI / 180 - x);
86 }
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 */
99         if (x < 0.0)
100         {
101                 x = -1 * x;
102                 sign = -1;
103         }
105         while( x > 0.3249196962 )
106         {
107                 k++;
108                 x = (x - 0.3249196962) / (1 + x * 0.3249196962);
109         }
111         addangle = k * 18.0 *M_PI/180;
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 }
123 double ASIN(double x)
124 {
125         return 2 * ATAN( x / (1 + std::sqrt(1.0 - x*x)));
126 }
128 double Radians( double number )
129 {
130         return number*M_PI/180;
131 }
133 double Deg( double number )
134 {
135         return number*180/M_PI;
136 }
138 double Rev( double number )
139 {
140         return number - std::floor( number / 360.0 ) * 360;
141 }
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,
151                         f = 1.00 / 298.257, // Earth flattning factor
153                         r_sat=42164.57, // Distance from earth centre to satellite
160                         Rstation = r_eq / ( std::sqrt( 1.00 - f*(2.00-f)*sinRadSiteLat*sinRadSiteLat ) ),
167                         alfa_rz=-Rz,
172                         El_geometric=Deg(ATAN( alfa_r_zenith/std::sqrt(alfa_r_north*alfa_r_north+alfa_ry*alfa_ry))),
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),
177                         El_observed = 0.00;
179         if (El_geometric > 10.2)
181         else
182                 El_observed = El_geometric+refraction;
184         if (alfa_r_zenith < -3000)
185                 El_observed=-99;
187         return El_observed;
188 }
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
194                         r_sat=42164.57, // Distance from earth centre to satellite
201                         Rstation = r_eq / ( std::sqrt( 1 - f*(2-f)*sinRadSiteLat*sinRadSiteLat ) ),
205                         alfa_rx = r_sat*COS(Radians(SatLon-SiteLon)) - Ra,
207                         alfa_rz = -Rz,
211                         Azimuth = 0.00;
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)));
218         return Azimuth;
219 }
221 double calcDeclination( double SiteLat, double Azimuth, double Elevation)
222 {
228                                                                                                 )
229                                                 );
230 }
232 double calcSatHourangle( double SatLon, double SiteLat, double SiteLon )
233 {
234         double  Azimuth=calcAzimuth(SatLon, SiteLat, SiteLon ),
235                         Elevation=calcElevation( SatLon, SiteLat, SiteLon ),
240                                 -