/* -*- Mode:C++; c-basic-offset:8; tab-width:8; indent-tabs-mode:t -*- */
/*
* Copyright (c) 1999 Regents of the University of California.
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
* 1. Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
* 3. All advertising materials mentioning features or use of this software
* must display the following acknowledgement:
* This product includes software developed by the MASH Research
* Group at the University of California Berkeley.
* 4. Neither the name of the University nor of the Research Group may be
* used to endorse or promote products derived from this software without
* specific prior written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
* ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
* OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
* HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
* SUCH DAMAGE.
*
* Contributed by Tom Henderson, UCB Daedalus Research Group, June 1999
* Modified to use period_ and offer isascending(), Lloyd Wood, March 2000. */
#ifndef lint
static const char rcsid[] =
"@(#) $Header: /nfs/jade/vint/CVSROOT/ns-2/satellite/satposition.cc,v 1.9 2005/08/22 05:08:34 tomh Exp $";
#endif
#include "satposition.h"
#include "satgeometry.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
static class TermSatPositionClass : public TclClass {
public:
TermSatPositionClass() : TclClass("Position/Sat/Term") {}
TclObject* create(int argc, const char*const* argv) {
if (argc == 5) {
float a, b;
sscanf(argv[4], "%f %f", &a, &b);
return (new TermSatPosition(a, b));
} else
return (new TermSatPosition);
}
} class_term_sat_position;
static class PolarSatPositionClass : public TclClass {
public:
PolarSatPositionClass() : TclClass("Position/Sat/Polar") {}
TclObject* create(int argc, const char*const* argv) {
if (argc == 5) {
float a = 0, b = 0, c = 0, d = 0, e = 0;
sscanf(argv[4], "%f %f %f %f %f", &a, &b, &c, &d, &e);
return (new PolarSatPosition(a, b, c, d, e));
}
else {
return (new PolarSatPosition);
}
}
} class_polarsatposition;
static class GeoSatPositionClass : public TclClass {
public:
GeoSatPositionClass() : TclClass("Position/Sat/Geo") {}
TclObject* create(int argc, const char*const* argv) {
if (argc == 5)
return (new GeoSatPosition(double(atof(argv[4]))));
else
return (new GeoSatPosition);
}
} class_geosatposition;
static class SatPositionClass : public TclClass {
public:
SatPositionClass() : TclClass("Position/Sat") {}
TclObject* create(int, const char*const*) {
printf("Error: do not instantiate Position/Sat\n");
return (0);
}
} class_satposition;
double SatPosition::time_advance_ = 0;
SatPosition::SatPosition() : node_(0)
{
bind("time_advance_", &time_advance_);
}
int SatPosition::command(int argc, const char*const* argv) {
//Tcl& tcl = Tcl::instance();
if (argc == 2) {
}
if (argc == 3) {
if(strcmp(argv[1], "setnode") == 0) {
node_ = (Node*) TclObject::lookup(argv[2]);
if (node_ == 0)
return TCL_ERROR;
return TCL_OK;
}
}
return (TclObject::command(argc, argv));
}
/////////////////////////////////////////////////////////////////////
// class TermSatPosition
/////////////////////////////////////////////////////////////////////
// Specify initial coordinates. Default coordinates place the terminal
// on the Earth's surface at 0 deg lat, 0 deg long.
TermSatPosition::TermSatPosition(double Theta, double Phi) {
initial_.r = EARTH_RADIUS;
period_ = EARTH_PERIOD; // seconds
set(Theta, Phi);
type_ = POSITION_SAT_TERM;
}
//
// Convert user specified latitude and longitude to our spherical coordinates
// Latitude is in the range (-90, 90) with neg. values -> south
// Initial_.theta is stored from 0 to PI (spherical)
// Longitude is in the range (-180, 180) with neg. values -> west
// Initial_.phi is stored from 0 to 2*PI (spherical)
//
void TermSatPosition::set(double latitude, double longitude)
{
if (latitude < -90 || latitude > 90)
fprintf(stderr, "TermSatPosition: latitude out of bounds %f\n",
latitude);
if (longitude < -180 || longitude > 180)
fprintf(stderr, "TermSatPosition: longitude out of bounds %f\n",
longitude);
initial_.theta = DEG_TO_RAD(90 - latitude);
if (longitude < 0)
initial_.phi = DEG_TO_RAD(360 + longitude);
else
initial_.phi = DEG_TO_RAD(longitude);
}
coordinate TermSatPosition::coord()
{
coordinate current;
current.r = initial_.r;
current.theta = initial_.theta;
current.phi = fmod((initial_.phi +
(fmod(NOW + time_advance_,period_)/period_) * 2*PI), 2*PI);
#ifdef POINT_TEST
current = initial_; // debug option to stop earth's rotation
#endif
return current;
}
/////////////////////////////////////////////////////////////////////
// class PolarSatPosition
/////////////////////////////////////////////////////////////////////
PolarSatPosition::PolarSatPosition(double altitude, double Inc, double Lon,
double Alpha, double Plane) : next_(0), plane_(0) {
set(altitude, Lon, Alpha, Inc);
bind("plane_", &plane_);
if (Plane)
plane_ = int(Plane);
type_ = POSITION_SAT_POLAR;
}
//
// Since it is easier to compute instantaneous orbit position based on a
// coordinate system centered on the orbit itself, we keep initial coordinates
// specified in terms of the satellite orbit, and convert to true spherical
// coordinates in coord().
// Initial satellite position is specified as follows:
// initial_.theta specifies initial angle with respect to "ascending node"
// i.e., zero is at equator (ascending)-- this is the $alpha parameter in Otcl
// initial_.phi specifies East longitude of "ascending node"
// -- this is the $lon parameter in OTcl
// Note that with this system, only theta changes with time
//
void PolarSatPosition::set(double Altitude, double Lon, double Alpha, double Incl)
{
if (Altitude < 0) {
fprintf(stderr, "PolarSatPosition: altitude out of \
bounds: %f\n", Altitude);
exit(1);
}
initial_.r = Altitude + EARTH_RADIUS; // Altitude in km above the earth
if (Alpha < 0 || Alpha >= 360) {
fprintf(stderr, "PolarSatPosition: alpha out of bounds: %f\n",
Alpha);
exit(1);
}
initial_.theta = DEG_TO_RAD(Alpha);
if (Lon < -180 || Lon > 180) {
fprintf(stderr, "PolarSatPosition: lon out of bounds: %f\n",
Lon);
exit(1);
}
if (Lon < 0)
initial_.phi = DEG_TO_RAD(360 + Lon);
else
initial_.phi = DEG_TO_RAD(Lon);
if (Incl < 0 || Incl > 180) {
// retrograde orbits = (90 < Inclination < 180)
fprintf(stderr, "PolarSatPosition: inclination out of \
bounds: %f\n", Incl);
exit(1);
}
inclination_ = DEG_TO_RAD(Incl);
// XXX: can't use "num = pow(initial_.r,3)" here because of linux lib
double num = initial_.r * initial_.r * initial_.r;
period_ = 2 * PI * sqrt(num/MU); // seconds
}
//
// The initial coordinate has the following properties:
// theta: 0 < theta < 2 * PI (essentially, this specifies initial position)
// phi: 0 < phi < 2 * PI (longitude of ascending node)
// Return a coordinate with the following properties (i.e. convert to a true
// spherical coordinate):
// theta: 0 < theta < PI
// phi: 0 < phi < 2 * PI
//
coordinate PolarSatPosition::coord()
{
coordinate current;
double partial; // fraction of orbit period completed
partial =
(fmod(NOW + time_advance_, period_)/period_) * 2*PI; //rad
double theta_cur, phi_cur, theta_new, phi_new;
// Compute current orbit-centric coordinates:
// theta_cur adds effects of time (orbital rotation) to initial_.theta
theta_cur = fmod(initial_.theta + partial, 2*PI);
phi_cur = initial_.phi;
// Reminder: theta_cur and phi_cur are temporal translations of
// initial parameters and are NOT true spherical coordinates.
//
// Now generate actual spherical coordinates,
// with 0 < theta_new < PI and 0 < phi_new < 360
assert (inclination_ < PI);
// asin returns value between -PI/2 and PI/2, so
// theta_new guaranteed to be between 0 and PI
theta_new = PI/2 - asin(sin(inclination_) * sin(theta_cur));
// if theta_new is between PI/2 and 3*PI/2, must correct
// for return value of atan()
if (theta_cur > PI/2 && theta_cur < 3*PI/2)
phi_new = atan(cos(inclination_) * tan(theta_cur)) +
phi_cur + PI;
else
phi_new = atan(cos(inclination_) * tan(theta_cur)) +
phi_cur;
phi_new = fmod(phi_new + 2*PI, 2*PI);
current.r = initial_.r;
current.theta = theta_new;
current.phi = phi_new;
return current;
}
//
// added by Lloyd Wood, 27 March 2000.
// allows us to distinguish between satellites that are ascending (heading north)
// and descending (heading south).
//
bool PolarSatPosition::isascending()
{
double partial = (fmod(NOW + time_advance_, period_)/period_) * 2*PI; //rad
double theta_cur = fmod(initial_.theta + partial, 2*PI);
if ((theta_cur > PI/2)&&(theta_cur < 3*PI/2)) {
return 0;
} else {
return 1;
}
}
int PolarSatPosition::command(int argc, const char*const* argv) {
Tcl& tcl = Tcl::instance();
if (argc == 2) {
}
if (argc == 3) {
if (strcmp(argv[1], "set_next") == 0) {
next_ = (PolarSatPosition *) TclObject::lookup(argv[2]);
if (next_ == 0) {
tcl.resultf("no such object %s", argv[2]);
return (TCL_ERROR);
}
return (TCL_OK);
}
}
return (SatPosition::command(argc, argv));
}
/////////////////////////////////////////////////////////////////////
// class GeoSatPosition
/////////////////////////////////////////////////////////////////////
GeoSatPosition::GeoSatPosition(double longitude)
{
initial_.r = EARTH_RADIUS + GEO_ALTITUDE;
initial_.theta = PI/2;
set(longitude);
type_ = POSITION_SAT_GEO;
period_ = EARTH_PERIOD;
}
coordinate GeoSatPosition::coord()
{
coordinate current;
current.r = initial_.r;
current.theta = initial_.theta;
double fractional =
(fmod(NOW + time_advance_, period_)/period_) *2*PI; // rad
current.phi = fmod(initial_.phi + fractional, 2*PI);
return current;
}
//
// Longitude is in the range (0, 180) with negative values -> west
//
void GeoSatPosition::set(double longitude)
{
if (longitude < -180 || longitude > 180)
fprintf(stderr, "GeoSatPosition: longitude out of bounds %f\n",
longitude);
if (longitude < 0)
initial_.phi = DEG_TO_RAD(360 + longitude);
else
initial_.phi = DEG_TO_RAD(longitude);
}
syntax highlighted by Code2HTML, v. 0.9.1