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 }