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