Cosmos::Physics Namespace Reference

Classes
class	ElectricalPropagator
class	GaussJacksonPositionPropagator
class	InertialAttitudePropagator
class	InertialPositionPropagator
class	IterativeAttitudePropagator
class	IterativePositionPropagator
class	LVLHAttitudePropagator
class	Propagator
class	State
class	Structure
class	ThermalPropagator

Functions
int32_t	PhysCalc (locstruc *loc, physicsstruc *phys)
int32_t	PhysSetup (physicsstruc *phys)
int32_t	AttAccel (locstruc *loc, physicsstruc *phys)
	Attitude acceleration. More...
int32_t	PosAccel (locstruc *loc, physicsstruc *phys)
	Position acceleration. More...
double	Msis00Density (posstruc pos, float f107avg, float f107, float magidx)
	Calculate atmospheric density. More...
Vector	GravityAccel (posstruc pos, uint16_t model, uint32_t degree)
	Spherical harmonic gravitational vector. More...
int32_t	GravityParams (int16_t model)
double	nplgndr (uint32_t l, uint32_t m, double x)
	Legendre polynomial. More...
double	Nplgndr (uint32_t l, uint32_t m, double x)

Variables
static const uint8_t	GravityPGM2000A = 1
static const uint8_t	GravityEGM2008 = 2
static const uint8_t	GravityPGM2000A_NORM = 3
static const uint8_t	GravityEGM2008_NORM = 4
static const uint16_t	maxdegree = 360
	Data structures for spherical harmonic expansion. More...
static double	vc [maxdegree+1][maxdegree+1]
static double	wc [maxdegree+1][maxdegree+1]
static double	coef [maxdegree+1][maxdegree+1][2]
static double	ftl [2 *maxdegree+1]
static double	spmm [maxdegree+1]

Function Documentation

int32_t Cosmos::Physics::PhysCalc	(	locstruc *	loc,
		physicsstruc *	phys
	)

Calculate dynamic physical attributes Calculate various derived physical quantities that change, like heat, power generation, torque and drag

Parameters

loc	Pointer to locstruc
phys	Pointer to ::physstruc

Returns

Zero, or negative error.

        {
            Vector unitv = Quaternion(loc->att.geoc.s).irotate(Vector(loc->pos.geoc.v).normalize());
            Vector units = Quaternion(loc->att.icrf.s).irotate(Vector(loc->pos.icrf.s).normalize());
            Vector geov(loc->pos.geoc.v);
            double speed = geov.norm();
            double density;
            if (loc->pos.geod.s.h < 10000. || std::isnan(loc->pos.geod.s.h))
                density = 1.225;
            else
                density = 1000. * Msis00Density(loc->pos, 150., 150., 3.);
            double adrag = density * 1.1 * speed * speed;

            // External panel effects
            phys->powgen = 0.F;
            phys->adrag.clear();
            phys->atorque.clear();
            phys->rdrag.clear();
            phys->rtorque.clear();
            for (trianglestruc& triangle : phys->triangles)
            {
                if (triangle.external)
                {
                    // Atmospheric effects
                    double vdot = unitv.dot(triangle.normal);
                    double ddrag = 0.;
                    if (vdot > 0. && phys->mass > 0.F)
                    {
                        ddrag = adrag * vdot / phys->mass;
                    }
                    Vector dtorque = ddrag * triangle.twist;
                    phys->atorque += dtorque;
                    Vector da = ddrag * triangle.shove;
                    phys->adrag += da;

                    // Solar effects
                    double sdot = units.dot(triangle.normal);
                    if (loc->pos.sunradiance && sdot > 0)
                    {
                        ddrag = loc->pos.sunradiance * sdot / (3e8*phys->mass);
                        dtorque = ddrag * triangle.twist;
                        phys->rtorque += dtorque;
                        da = ddrag * triangle.shove;
                        phys->rdrag += da;
                    }
                }
            }
        }

int32_t Cosmos::Physics::PhysSetup

(

physicsstruc *

phys

)

Calculate static physical attributes Calculate various derived physical quantities, like Center of Mass and Moments of Inertia

Parameters

loc	Pointer to locstruc
phys	Pointer to ::physstruc

Returns

Zero, or negative error.

        {
            // Calculate Center of mass and internal area
            phys->area = 0.;
            phys->com.clear();
            phys->mass = 0.;
            for (trianglestruc triangle : phys->triangles)
            {
                if (triangle.external <= 1)
                {
                    phys->area += triangle.area;
                }
                if (!triangle.external)
                {
                    phys->area += triangle.area;
                }
                phys->com += triangle.com * triangle.mass;
                phys->mass += triangle.mass;
            }
            if (phys->mass == 0.F)
            {
                return GENERAL_ERROR_TOO_LOW;
            }
            phys->com /= phys->mass;

            // Calculate Principal Moments of Inertia WRT COM
            phys->moi.clear();
            for (trianglestruc triangle : phys->triangles)
            {
                phys->moi.x += (phys->com.x - triangle.com.x) * (phys->com.x - triangle.com.x) * triangle.mass;
                phys->moi.y += (phys->com.y - triangle.com.y) * (phys->com.y - triangle.com.y) * triangle.mass;
                phys->moi.z += (phys->com.z - triangle.com.z) * (phys->com.z - triangle.com.z) * triangle.mass;
            }

            return 0;
        }

int32_t Cosmos::Physics::AttAccel	(	locstruc *	loc,
		physicsstruc *	phys
	)

Attitude acceleration.

Calculate the torque forces on the specified satellite at the specified location/

Parameters

physics	Pointer to structure specifying satellite.
loc	Structure specifying location.

        {
            //    rvector ue, ta, tv;
            //    rvector ttorque;
            Vector ue, ta, tv;
            Vector ttorque;
            rmatrix mom;

            att_extra(loc);

            ttorque = phys->ctorque;

            // Now calculate Gravity Gradient Torque
            // Unit vector towards earth, rotated into body frame
            //    ue = irotate((loc->att.icrf.s),rv_smult(-1.,loc->pos.eci.s));
            //    normalize_rv(ue);
            ue = Quaternion(loc->att.icrf.s).irotate(-1. * Vector(loc->pos.eci.s)).normalize();

            //    phys->gtorque = rv_smult((3.*GM/pow(loc->pos.geos.s.r,3.)),rv_cross(ue,rv_mult(phys->moi,ue)));
            phys->gtorque = (3. * GM / pow(loc->pos.geos.s.r,3.)) * ue.cross(phys->moi * ue);

            //    ttorque = rv_add(ttorque,phys->gtorque);
            ttorque += phys->gtorque;

            // Atmospheric and solar torque
            //  ttorque = rv_add(ttorque,phys->atorque);
            //  ttorque = rv_add(ttorque,phys->rtorque);

            // Torque from rotational effects

            // Moment of Inertia in Body frame
            mom = rm_diag(phys->moi.to_rv());
            // Attitude rate in Body frame
            tv = irotate(loc->att.icrf.s,loc->att.icrf.v);

            // Torque from cross product of angular velocity and angular momentum
            //    phys->htorque = rv_smult(-1., rv_cross(tv,rv_add(rv_mmult(mom,tv),phys->hmomentum)));
            //    ttorque = rv_add(ttorque,phys->htorque);
            phys->htorque = -1. * tv.cross(Vector(rv_mmult(mom, tv.to_rv())) + phys->hmomentum);
            ttorque += phys->htorque;

            // Convert torque into accelerations, doing math in Body frame

            // I x alpha = tau, so alpha = I inverse x tau
            //    ta = rv_mmult(rm_inverse(mom),ttorque);
            ta = Vector(rv_mmult(rm_inverse(mom),ttorque.to_rv()));

            // Convert body frame acceleration back to other frames.
            loc->att.icrf.a = irotate(q_conjugate(loc->att.icrf.s), ta.to_rv());
            loc->att.topo.a = irotate(q_conjugate(loc->att.topo.s), ta.to_rv());
            loc->att.lvlh.a = irotate(q_conjugate(loc->att.lvlh.s), ta.to_rv());
            loc->att.geoc.a = irotate(q_conjugate(loc->att.geoc.s), ta.to_rv());
            loc->att.selc.a = irotate(q_conjugate(loc->att.selc.s), ta.to_rv());
            return 0;
        }

int32_t Cosmos::Physics::PosAccel	(	locstruc *	loc,
		physicsstruc *	phys
	)

Position acceleration.

Calculate the linear forces on the specified sattelite at the specified location/

Parameters

phys	Pointer to structure specifying satellite.
loc	Structure specifying location.

        {
            int32_t iretn;
            double radius;
            Vector ctpos, da, tda;
            cartpos bodypos;

            radius = length_rv(loc->pos.eci.s);

            loc->pos.eci.a = rv_zero();

            // Earth gravity
            // Calculate Geocentric acceleration vector

            if (radius > REARTHM)
            {
                // Start with gravity vector in ITRS

                da = GravityAccel(loc->pos,GravityEGM2008_NORM,12);

                // Correct for earth rotation, polar motion, precession, nutation

                da = Matrix(loc->pos.extra.e2j) * da;
            }
            else
            {
                // Simple 2 body
                da = -GM/(radius*radius*radius) * Vector(loc->pos.eci.s);
            }
            loc->pos.eci.a = rv_add(loc->pos.eci.a, da.to_rv());

            // Sun gravity
            // Calculate Satellite to Sun vector
            ctpos = rv_sub(rv_smult(-1., loc->pos.extra.sun2earth.s), loc->pos.eci.s);
            radius = ctpos.norm();
            da = GSUN/(radius*radius*radius) * ctpos;
            loc->pos.eci.a = rv_add(loc->pos.eci.a, da.to_rv());

            // Adjust for acceleration of frame
            radius = length_rv(loc->pos.extra.sun2earth.s);
            da = rv_smult(GSUN/(radius*radius*radius), loc->pos.extra.sun2earth.s);
            tda = da;
            loc->pos.eci.a = rv_add(loc->pos.eci.a, da.to_rv());

            // Moon gravity
            // Calculate Satellite to Moon vector
            bodypos.s = rv_sub( loc->pos.extra.sun2earth.s, loc->pos.extra.sun2moon.s);
            ctpos = rv_sub(bodypos.s, loc->pos.eci.s);
            radius = ctpos.norm();
            da = GMOON/(radius*radius*radius) * ctpos;
            loc->pos.eci.a = rv_add(loc->pos.eci.a, da.to_rv());

            // Adjust for acceleration of frame due to moon
            radius = length_rv(bodypos.s);
            da = rv_smult(GMOON/(radius*radius*radius),bodypos.s);
            tda -= da;
            loc->pos.eci.a = rv_sub(loc->pos.eci.a, da.to_rv());

            /*
        // Jupiter gravity
        // Calculate Satellite to Jupiter vector
        jplpos(JPL_EARTH,JPL_JUPITER, loc->pos.extra.tt,(cartpos *)&bodypos);
        ctpos = rv_sub(bodypos.s, loc->pos.eci.s);
        radius = length_rv(ctpos);

        // Calculate acceleration
        da = rv_smult(GJUPITER/(radius*radius*radius),ctpos);
        // loc->pos.eci.a = rv_add( loc->pos.eci.a,da);
        */


            Quaternion iratt = Quaternion(loc->att.icrf.s).conjugate();
            if (phys->adrag.norm() > 0.)
            {
                da = iratt.irotate(phys->adrag);
                loc->pos.eci.a = rv_add(loc->pos.eci.a, da.to_rv());
            }
            // Solar drag
            if (phys->rdrag.norm() > 0.)
            {
                da = iratt.irotate(phys->rdrag);
                loc->pos.eci.a = rv_add(loc->pos.eci.a, da.to_rv());
            }
            // Fictitious drag
            if (phys->fdrag.norm() > 0.)
            {
                da = iratt.irotate(phys->fdrag);
                loc->pos.eci.a = rv_add(loc->pos.eci.a, da.to_rv());
            }

            loc->pos.eci.pass++;
            iretn = pos_eci(loc);
            if (iretn < 0)
            {
                return iretn;
            }
            if (std::isnan( loc->pos.eci.a.col[0]))
            {
                loc->pos.eci.a.col[0] = 0.;
            }
            return 0;
        }

double Cosmos::Physics::Msis00Density	(	posstruc	pos,
		float	f107avg,
		float	f107,
		float	magidx
	)

Calculate atmospheric density.

Calculate atmospheric density at indicated Latitute/Longitude/Altitude using the NRLMSISE-00 atmospheric model.

Parameters

pos	Structure indicating position
f107avg	Average 10.7 cm solar flux
f107	Current 10.7 cm solar flux
magidx	Ap daily geomagnetic index

Returns

Density in kg/m3

        {
            struct nrlmsise_output output;
            struct nrlmsise_input input;
            struct nrlmsise_flags flags = {
            {0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1},
            {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.},
            {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.}};
            int year, month;
            double day, doy;
            static double lastmjd = 0.;
            static double lastdensity = 0.;
            static double lastperiod = 0.;

            if (lastperiod != 0.)
            {
                if (fabs(lastperiod) > (pos.extra.utc-lastmjd))
                {
                    return (lastdensity*(1.+(.001*(pos.extra.utc-lastmjd)/lastperiod)));
                }
            }

            mjd2ymd(pos.extra.utc,year,month,day,doy);
            input.doy = static_cast <int32_t>(doy);
            input.g_lat = pos.geod.s.lat*180./DPI;
            input.g_long = pos.geod.s.lon*180./DPI;
            input.alt = pos.geod.s.h / 1000.;
            input.sec = (doy - input.doy)*86400.;;
            input.lst = input.sec / 3600. + input.g_long / 15.;
            input.f107A = f107avg;
            input.f107 = f107;
            input.ap = magidx;
            gtd7d(&input,&flags,&output);

            if (lastdensity != 0. && lastdensity != output.d[5])
                lastperiod = (pos.extra.utc-lastmjd)*.001*output.d[5]/(output.d[5]-lastdensity);
            lastmjd = pos.extra.utc;
            lastdensity = output.d[5];
            return((double)output.d[5]);
        }

Vector Cosmos::Physics::GravityAccel	(	posstruc	pos,
		uint16_t	model,
		uint32_t	degree
	)

Spherical harmonic gravitational vector.

Calculates a spherical harmonic expansion of the chosen model of indicated order and degree for the requested position. The result is returned as a geocentric vector calculated at the epoch.

Parameters

pos	a posstruc providing the position at the epoch
model	Model to use for coefficients
degree	Order and degree to calculate

Returns

A ::Vector pointing toward the earth

See also

pgm2000a_coef.txt

        {
            uint32_t il, im;
            double tmult;
            double ratio, rratio, xratio, yratio, zratio;
            Vector accel;
            double fr;

            // Zero out vc and wc
            memset(vc,0,sizeof(vc));
            memset(wc,0,sizeof(wc));

            // Load Params
            GravityParams(model);

            // Calculate cartesian Legendre terms
            vc[0][0] = REARTHM/pos.geos.s.r;
            wc[0][0] = 0.;
            ratio = vc[0][0] / pos.geos.s.r;
            rratio = REARTHM * ratio;
            xratio = pos.geoc.s.col[0] * ratio;
            yratio = pos.geoc.s.col[1] * ratio;
            zratio = pos.geoc.s.col[2] * ratio;
            vc[1][0] = zratio * vc[0][0];
            wc[1][0] = 0.;
            for (il=2; il<=degree+1; il++)
            {
                vc[il][0] = (2*il-1)*zratio * vc[il-1][0] / il - (il-1) * rratio * vc[il-2][0] / il;
                wc[il][0] = 0.;
            }
            for (im=1; im<=degree+1; im++)
            {
                vc[im][im] = (2*im-1) * (xratio * vc[im-1][im-1] - yratio * wc[im-1][im-1]);
                wc[im][im] = (2*im-1) * (xratio * wc[im-1][im-1] + yratio * vc[im-1][im-1]);
                if (im <= degree)
                {
                    vc[im+1][im] = (2*im+1) * zratio * vc[im][im];
                    wc[im+1][im] = (2*im+1) * zratio * wc[im][im];
                }
                for (il=im+2; il<=degree+1; il++)
                {
                    vc[il][im] = (2*il-1) * zratio * vc[il-1][im] / (il-im) - (il+im-1) * rratio * vc[il-2][im] / (il-im);
                    wc[il][im] = (2*il-1) * zratio * wc[il-1][im] / (il-im) - (il+im-1) * rratio * wc[il-2][im] / (il-im);
                }
            }

            //  dr = dlon = dlat = 0.;

            accel.clear();
            for (im=0; im<=degree; im++)
            {
                for (il=im; il<=degree; il++)
                {
                    if (im == 0)
                    {
                        accel[0] -= coef[il][0][0] * vc[il+1][1];
                        accel[1] -= coef[il][0][0] * wc[il+1][1];
                        accel[2] -= (il+1) * (coef[il][0][0] * vc[il+1][0]);
                    }
                    else
                    {
                        fr = ftl[il-im+2] / ftl[il-im];
                        accel[0] -= .5 * (coef[il][im][0] * vc[il+1][im+1] + coef[il][im][1] * wc[il+1][im+1] - fr * (coef[il][im][0] * vc[il+1][im-1] + coef[il][im][1] * wc[il+1][im-1]));
                        accel[1] -= .5 * (coef[il][im][0] * wc[il+1][im+1] - coef[il][im][1] * vc[il+1][im+1] + fr * (coef[il][im][0] * wc[il+1][im-1] - coef[il][im][1] * vc[il+1][im-1]));
                        accel[2] -= (il-im+1) * (coef[il][im][0] * vc[il+1][im] + coef[il][im][1] * wc[il+1][im]);
                    }
                }
            }
            tmult = GM / (REARTHM*REARTHM);
            accel[0] *= tmult;
            accel[2] *= tmult;
            accel[1] *= tmult;

            return (accel);
        }

int32_t Cosmos::Physics::GravityParams

(

int16_t

model

)

        {
            static int16_t cmodel = -1;

            int32_t iretn = 0;
            uint32_t il, im;
            double norm;
            uint32_t dil, dim;
            double dummy1, dummy2;

            // Calculate factorial
            if (ftl[0] == 0.)
            {
                ftl[0] = 1.;
                for (il=1; il<2*maxdegree+1; il++)
                {
                    ftl[il] = il * ftl[il-1];
                }
            }

            // Load Coefficients
            if (cmodel != model)
            {
                coef[0][0][0] = 1.;
                coef[0][0][1] = 0.;
                coef[1][0][0] = coef[1][0][1] = 0.;
                coef[1][1][0] = coef[1][1][1] = 0.;
                string fname;
                FILE *fi;
                switch (model)
                {
                case GravityEGM2008:
                case GravityEGM2008_NORM:
                    fname = get_cosmosresources();
                    if (fname.empty())
                    {
                        return GENERAL_ERROR_EMPTY;
                    }
                    fname += "/general/egm2008_coef.txt";
                    fi = fopen(fname.c_str(),"r");

                    if (fi==nullptr)
                    {
                        cout << "could not load file " << fname << endl;
                        return iretn;
                    }

                    for (il=2; il<101; il++)
                    {
                        for (im=0; im<= il; im++)
                        {
                            iretn = fscanf(fi,"%u %u %lf %lf\n",&dil,&dim,&coef[il][im][0],&coef[il][im][1]);
                            if (iretn && model == GravityEGM2008_NORM)
                            {
                                norm = sqrt(ftl[il+im]/((2-(im==0?1:0))*(2*il+1)*ftl[il-im]));
                                coef[il][im][0] /= norm;
                                coef[il][im][1] /= norm;
                            }
                        }
                    }
                    fclose(fi);
                    cmodel = model;
                    break;
                case GravityPGM2000A:
                case GravityPGM2000A_NORM:
                default:
                    iretn = get_cosmosresources(fname);
                    if (iretn < 0)
                    {
                        return iretn;
                    }
                    fname += "/general/pgm2000a_coef.txt";
                    fi = fopen(fname.c_str(),"r");
                    for (il=2; il<361; il++)
                    {
                        for (im=0; im<= il; im++)
                        {
                            iretn = fscanf(fi,"%u %u %lf %lf %lf %lf\n",&dil,&dim,&coef[il][im][0],&coef[il][im][1],&dummy1,&dummy2);
                            if (iretn && model == GravityPGM2000A_NORM)
                            {
                                norm = sqrt(ftl[il+im]/((2-(il==im?1:0))*(2*il+1)*ftl[il-im]));
                                coef[il][im][0] /= norm;
                                coef[il][im][1] /= norm;
                            }
                        }
                    }
                    fclose(fi);
                    cmodel = model;
                    break;
                }
            }
            return 0;
        }

double Cosmos::Physics::Nplgndr	(	uint32_t	l,
		uint32_t	m,
		double	x
	)

Variable Documentation

const uint8_t Cosmos::Physics::GravityPGM2000A = 1

static

const uint8_t Cosmos::Physics::GravityEGM2008 = 2

static

const uint8_t Cosmos::Physics::GravityPGM2000A_NORM = 3

static

const uint8_t Cosmos::Physics::GravityEGM2008_NORM = 4

static

const uint16_t Cosmos::Physics::maxdegree = 360

static

Data structures for spherical harmonic expansion.

Coefficients for real and imaginary components of expansion. Of order and rank maxdegree

double Cosmos::Physics::vc[maxdegree+1][maxdegree+1]

static

double Cosmos::Physics::wc[maxdegree+1][maxdegree+1]

static

double Cosmos::Physics::coef[maxdegree+1][maxdegree+1][2]

static

double Cosmos::Physics::ftl[2 *maxdegree+1]

static

double Cosmos::Physics::spmm[maxdegree+1]

static