---

Orbit.cxx

Go to the documentation of this file.

// $Header: /nfs/slac/g/glast/ground/cvs/flux/src/Orbit.cxx,v 1.4 2001/05/14 16:37:50 burnett Exp $

#include "flux/Orbit.h"

#include "CLHEP/Random/RandFlat.h"
#include "flux/GPS.h"

#include <cmath>
#include <cassert>


Orbit::Orbit(double asclon, double alt, double inc, double phs
           , double pitch , double yaw, double roll) :
// constructor
m_ascendingLon(asclon),
m_altitude(alt),
m_period(1.6612e-4 * pow(alt+6371.2, 1.5)),// Kepler's Law
m_startphase(phs)

{
    inclination( inc );
}

void    Orbit::inclination ( double inc )
{
    m_inclination = inc;
    m_sini = (sin(m_inclination*M_2PI/360.));
    m_cosi = (cos(m_inclination*M_2PI/360.));
}

void    Orbit::ascendingLon ( double asc )
{
    m_ascendingLon = asc;
}

void    Orbit::startphase ( double phs )
{
    m_startphase = phs;
}

double Orbit::latitude(double time) const {
    // latitude as a function of time (in minutes)
    return (360./M_2PI) * asin(m_sini * sin(M_2PI*time/m_period + startphase()));
}

double Orbit::longitude(double time) const {
    // longitude as a function of time, taking into account the starting
    // longitude and the eastward rotation of the earth

    double phase = M_2PI * time / m_period + startphase();
    double lon = (360./M_2PI) * atan2(m_cosi*sin(phase), cos(phase)) -
        360.* time /1440. + m_ascendingLon; // orbital motion - earth rotation
    lon = fmod(lon, 360.);   // Fold into the range [0, 360)
    if (lon < 0.) lon += 360.;
    return lon;
}

double Orbit::pitch (double time) const {
    return 0.;
}

double Orbit::yaw (double time) const {
    return 0.;
}

double Orbit::roll (double time) const {
    return 0.;
}

double Orbit::phase (double time) const {
    double p = period();
    return (p == 0) ? 0. : time/p * 2*M_PI + startphase();
}

std::pair<double,double> Orbit::coords(double time) const {
    return std::make_pair(latitude(time), longitude(time));
}

---