satposition.cc

Go to the documentation of this file.

/* -*-  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.8 2000/06/21 17:44:10 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

        if (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);
        } else 
                printf("ERROR:  inclination_ > PI\n");
        
        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);
}

---