satposition.cc
上传用户:rrhhcc
上传日期:2015-12-11
资源大小:54129k
文件大小:11k
源码类别:

通讯编程

开发平台:

Visual C++

  1. /* -*-  Mode:C++; c-basic-offset:8; tab-width:8; indent-tabs-mode:t -*- */
  2. /*
  3.  * Copyright (c) 1999 Regents of the University of California.
  4.  * All rights reserved.
  5.  *
  6.  * Redistribution and use in source and binary forms, with or without
  7.  * modification, are permitted provided that the following conditions
  8.  * are met:
  9.  * 1. Redistributions of source code must retain the above copyright
  10.  *    notice, this list of conditions and the following disclaimer.
  11.  * 2. Redistributions in binary form must reproduce the above copyright
  12.  *    notice, this list of conditions and the following disclaimer in the
  13.  *    documentation and/or other materials provided with the distribution.
  14.  * 3. All advertising materials mentioning features or use of this software
  15.  *    must display the following acknowledgement:
  16.  *      This product includes software developed by the MASH Research
  17.  *      Group at the University of California Berkeley.
  18.  * 4. Neither the name of the University nor of the Research Group may be
  19.  *    used to endorse or promote products derived from this software without
  20.  *    specific prior written permission.
  21.  *
  22.  * THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
  23.  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  24.  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  25.  * ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
  26.  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  27.  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
  28.  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  29.  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  30.  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  31.  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  32.  * SUCH DAMAGE.
  33.  *
  34.  * Contributed by Tom Henderson, UCB Daedalus Research Group, June 1999
  35.  * Modified to use period_ and offer isascending(), Lloyd Wood, March 2000. */
  36. #ifndef lint
  37. static const char rcsid[] =
  38.     "@(#) $Header: /cvsroot/nsnam/ns-2/satellite/satposition.cc,v 1.9 2005/08/22 05:08:34 tomh Exp $";
  39. #endif
  40. #include "satposition.h"
  41. #include "satgeometry.h"
  42. #include <stdio.h>
  43. #include <stdlib.h>
  44. #include <math.h>
  45. static class TermSatPositionClass : public TclClass {
  46. public:
  47. TermSatPositionClass() : TclClass("Position/Sat/Term") {}
  48. TclObject* create(int argc, const char*const* argv) {
  49. if (argc == 5) {
  50. float a, b;
  51. sscanf(argv[4], "%f %f", &a, &b);
  52. return (new TermSatPosition(a, b));
  53. } else
  54. return (new TermSatPosition);
  55. }
  56. } class_term_sat_position;
  57. static class PolarSatPositionClass : public TclClass {
  58. public:
  59. PolarSatPositionClass() : TclClass("Position/Sat/Polar") {}
  60. TclObject* create(int argc, const char*const* argv) {
  61. if (argc == 5) {
  62. float a = 0, b = 0, c = 0, d = 0, e = 0;
  63. sscanf(argv[4], "%f %f %f %f %f", &a, &b, &c, &d, &e);
  64. return (new PolarSatPosition(a, b, c, d, e));
  65. }
  66. else {
  67. return (new PolarSatPosition);
  68. }
  69. }
  70. } class_polarsatposition;
  71. static class GeoSatPositionClass : public TclClass {
  72. public:
  73. GeoSatPositionClass() : TclClass("Position/Sat/Geo") {}
  74. TclObject* create(int argc, const char*const* argv) {
  75. if (argc == 5) 
  76. return (new GeoSatPosition(double(atof(argv[4]))));
  77. else 
  78. return (new GeoSatPosition);
  79. }
  80. } class_geosatposition;
  81. static class SatPositionClass : public TclClass {
  82. public:
  83. SatPositionClass() : TclClass("Position/Sat") {}
  84. TclObject* create(int, const char*const*) {
  85. printf("Error:  do not instantiate Position/Satn");
  86. return (0);
  87. }
  88. } class_satposition;
  89. double SatPosition::time_advance_ = 0;
  90. SatPosition::SatPosition() : node_(0)  
  91. {
  92.         bind("time_advance_", &time_advance_);
  93. }
  94. int SatPosition::command(int argc, const char*const* argv) {     
  95. //Tcl& tcl = Tcl::instance();
  96. if (argc == 2) {
  97. }
  98. if (argc == 3) {
  99. if(strcmp(argv[1], "setnode") == 0) {
  100. node_ = (Node*) TclObject::lookup(argv[2]);
  101. if (node_ == 0)
  102. return TCL_ERROR;
  103. return TCL_OK;
  104. }
  105. }
  106. return (TclObject::command(argc, argv));
  107. }
  108. /////////////////////////////////////////////////////////////////////
  109. // class TermSatPosition
  110. /////////////////////////////////////////////////////////////////////
  111. // Specify initial coordinates.  Default coordinates place the terminal
  112. // on the Earth's surface at 0 deg lat, 0 deg long.
  113. TermSatPosition::TermSatPosition(double Theta, double Phi)  {
  114. initial_.r = EARTH_RADIUS;
  115. period_ = EARTH_PERIOD; // seconds
  116. set(Theta, Phi);
  117. type_ = POSITION_SAT_TERM;
  118. }
  119. //
  120. // Convert user specified latitude and longitude to our spherical coordinates
  121. // Latitude is in the range (-90, 90) with neg. values -> south
  122. // Initial_.theta is stored from 0 to PI (spherical)
  123. // Longitude is in the range (-180, 180) with neg. values -> west
  124. // Initial_.phi is stored from 0 to 2*PI (spherical)
  125. //
  126. void TermSatPosition::set(double latitude, double longitude)
  127. {
  128. if (latitude < -90 || latitude > 90)
  129. fprintf(stderr, "TermSatPosition:  latitude out of bounds %fn",
  130.    latitude);
  131. if (longitude < -180 || longitude > 180)
  132. fprintf(stderr, "TermSatPosition: longitude out of bounds %fn",
  133.     longitude);
  134. initial_.theta = DEG_TO_RAD(90 - latitude);
  135. if (longitude < 0)
  136. initial_.phi = DEG_TO_RAD(360 + longitude);
  137. else
  138. initial_.phi = DEG_TO_RAD(longitude);
  139. }
  140. coordinate TermSatPosition::coord()
  141. {
  142. coordinate current;
  143. current.r = initial_.r;
  144. current.theta = initial_.theta;
  145. current.phi = fmod((initial_.phi + 
  146.     (fmod(NOW + time_advance_,period_)/period_) * 2*PI), 2*PI);
  147. #ifdef POINT_TEST
  148. current = initial_; // debug option to stop earth's rotation
  149. #endif
  150. return current;
  151. }
  152. /////////////////////////////////////////////////////////////////////
  153. // class PolarSatPosition
  154. /////////////////////////////////////////////////////////////////////
  155. PolarSatPosition::PolarSatPosition(double altitude, double Inc, double Lon, 
  156.     double Alpha, double Plane) : next_(0), plane_(0) {
  157. set(altitude, Lon, Alpha, Inc);
  158.         bind("plane_", &plane_);
  159.         if (Plane) 
  160. plane_ = int(Plane);
  161. type_ = POSITION_SAT_POLAR;
  162. }
  163. //
  164. // Since it is easier to compute instantaneous orbit position based on a
  165. // coordinate system centered on the orbit itself, we keep initial coordinates
  166. // specified in terms of the satellite orbit, and convert to true spherical 
  167. // coordinates in coord().
  168. // Initial satellite position is specified as follows:
  169. // initial_.theta specifies initial angle with respect to "ascending node"
  170. // i.e., zero is at equator (ascending)-- this is the $alpha parameter in Otcl
  171. // initial_.phi specifies East longitude of "ascending node"  
  172. // -- this is the $lon parameter in OTcl
  173. // Note that with this system, only theta changes with time
  174. //
  175. void PolarSatPosition::set(double Altitude, double Lon, double Alpha, double Incl)
  176. {
  177. if (Altitude < 0) {
  178. fprintf(stderr, "PolarSatPosition:  altitude out of 
  179.    bounds: %fn", Altitude);
  180. exit(1);
  181. }
  182. initial_.r = Altitude + EARTH_RADIUS; // Altitude in km above the earth
  183. if (Alpha < 0 || Alpha >= 360) {
  184. fprintf(stderr, "PolarSatPosition:  alpha out of bounds: %fn", 
  185.     Alpha);
  186. exit(1);
  187. }
  188. initial_.theta = DEG_TO_RAD(Alpha);
  189. if (Lon < -180 || Lon > 180) {
  190. fprintf(stderr, "PolarSatPosition:  lon out of bounds: %fn", 
  191.     Lon);
  192. exit(1);
  193. }
  194. if (Lon < 0)
  195. initial_.phi = DEG_TO_RAD(360 + Lon);
  196. else
  197. initial_.phi = DEG_TO_RAD(Lon);
  198. if (Incl < 0 || Incl > 180) {
  199. // retrograde orbits = (90 < Inclination < 180)
  200. fprintf(stderr, "PolarSatPosition:  inclination out of 
  201.     bounds: %fn", Incl); 
  202. exit(1);
  203. }
  204. inclination_ = DEG_TO_RAD(Incl);
  205. // XXX: can't use "num = pow(initial_.r,3)" here because of linux lib
  206. double num = initial_.r * initial_.r * initial_.r;
  207. period_ = 2 * PI * sqrt(num/MU); // seconds
  208. }
  209. //
  210. // The initial coordinate has the following properties:
  211. // theta: 0 < theta < 2 * PI (essentially, this specifies initial position)  
  212. // phi:  0 < phi < 2 * PI (longitude of ascending node)
  213. // Return a coordinate with the following properties (i.e. convert to a true
  214. // spherical coordinate):
  215. // theta:  0 < theta < PI
  216. // phi:  0 < phi < 2 * PI
  217. //
  218. coordinate PolarSatPosition::coord()
  219. {
  220. coordinate current;
  221. double partial;  // fraction of orbit period completed
  222. partial = 
  223.     (fmod(NOW + time_advance_, period_)/period_) * 2*PI; //rad
  224. double theta_cur, phi_cur, theta_new, phi_new;
  225. // Compute current orbit-centric coordinates:
  226. // theta_cur adds effects of time (orbital rotation) to initial_.theta
  227. theta_cur = fmod(initial_.theta + partial, 2*PI);
  228. phi_cur = initial_.phi;
  229. // Reminder:  theta_cur and phi_cur are temporal translations of 
  230. // initial parameters and are NOT true spherical coordinates.
  231. //
  232. // Now generate actual spherical coordinates,
  233. // with 0 < theta_new < PI and 0 < phi_new < 360
  234. assert (inclination_ < PI);
  235. // asin returns value between -PI/2 and PI/2, so 
  236. // theta_new guaranteed to be between 0 and PI
  237. theta_new = PI/2 - asin(sin(inclination_) * sin(theta_cur));
  238. // if theta_new is between PI/2 and 3*PI/2, must correct
  239. // for return value of atan()
  240. if (theta_cur > PI/2 && theta_cur < 3*PI/2)
  241. phi_new = atan(cos(inclination_) * tan(theta_cur)) + 
  242. phi_cur + PI;
  243. else
  244. phi_new = atan(cos(inclination_) * tan(theta_cur)) + 
  245. phi_cur;
  246. phi_new = fmod(phi_new + 2*PI, 2*PI);
  247. current.r = initial_.r;
  248. current.theta = theta_new;
  249. current.phi = phi_new;
  250. return current;
  251. }
  252. //
  253. // added by Lloyd Wood, 27 March 2000.
  254. // allows us to distinguish between satellites that are ascending (heading north)
  255. // and descending (heading south).
  256. //
  257. bool PolarSatPosition::isascending()
  258. {
  259. double partial = (fmod(NOW + time_advance_, period_)/period_) * 2*PI; //rad
  260. double theta_cur = fmod(initial_.theta + partial, 2*PI);
  261. if ((theta_cur > PI/2)&&(theta_cur < 3*PI/2)) {
  262. return 0;
  263. } else {
  264. return 1;
  265. }
  266. }
  267. int PolarSatPosition::command(int argc, const char*const* argv) {     
  268. Tcl& tcl = Tcl::instance();
  269.         if (argc == 2) {
  270. }
  271. if (argc == 3) {
  272. if (strcmp(argv[1], "set_next") == 0) {
  273. next_ = (PolarSatPosition *) TclObject::lookup(argv[2]);
  274. if (next_ == 0) {
  275. tcl.resultf("no such object %s", argv[2]);
  276. return (TCL_ERROR);
  277. }
  278. return (TCL_OK);
  279. }
  280. }
  281. return (SatPosition::command(argc, argv));
  282. }
  283. /////////////////////////////////////////////////////////////////////
  284. // class GeoSatPosition
  285. /////////////////////////////////////////////////////////////////////
  286. GeoSatPosition::GeoSatPosition(double longitude) 
  287. {
  288. initial_.r = EARTH_RADIUS + GEO_ALTITUDE;
  289. initial_.theta = PI/2;
  290. set(longitude);
  291. type_ = POSITION_SAT_GEO;
  292. period_ = EARTH_PERIOD;
  293. }
  294. coordinate GeoSatPosition::coord()
  295. {
  296. coordinate current;
  297. current.r = initial_.r;
  298. current.theta = initial_.theta;
  299. double fractional = 
  300.     (fmod(NOW + time_advance_, period_)/period_) *2*PI; // rad
  301. current.phi = fmod(initial_.phi + fractional, 2*PI);
  302. return current;
  303. }
  304. //
  305. // Longitude is in the range (0, 180) with negative values -> west
  306. //
  307. void GeoSatPosition::set(double longitude)
  308. {
  309. if (longitude < -180 || longitude > 180)
  310. fprintf(stderr, "GeoSatPosition:  longitude out of bounds %fn",
  311.     longitude);
  312. if (longitude < 0)
  313. initial_.phi = DEG_TO_RAD(360 + longitude);
  314. else
  315. initial_.phi = DEG_TO_RAD(longitude);
  316. }