source: titan/titan/rotorcalc.h @ 24101

Last change on this file since 24101 was 16511, checked in by nit, 12 years ago

[titan] delete status.screencalc

File size: 5.0 KB
Line 
1#ifndef ROTORCALC_H
2#define ROTORCALC_H
3
4double 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
16double 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
34double 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
74double COS(double x)
75{
76        return SIN(90 * M_PI / 180 - x);
77}
78
79double 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
113double ASIN(double x)
114{
115        return 2 * ATAN( x / (1 + sqrt(1.0 - x * x)));
116}
117
118double RAD(double nr)
119{
120        return nr * M_PI / 180;
121}
122
123double DEG(double nr)
124{
125        return nr * 180 / M_PI;
126}
127
128double REV(double nr)
129{
130        return nr - floor(nr / 360.0) * 360;
131}
132
133double 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
169double 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
193double 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
202double 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
Note: See TracBrowser for help on using the repository browser.