source: titan/titan/rotorcalc.h @ 15272

Last change on this file since 15272 was 12610, checked in by nit, 10 years ago

[titan] add calculation for usals

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