physicslib.cpp File Reference

#include "physics/physicslib.h"
#include "support/timelib.h"
#include "support/datalib.h"

Include dependency graph for physicslib.cpp:

Macros
#define	MAXDEGREE   360
#define	ASTEP   1

Functions
rvector	gravity_accel (posstruc pos, int model, uint32_t degree)
	Spherical harmonic gravitational vector. More...
rvector	gravity_vector (svector pos, int model, uint32_t degree)
double	gravity_potential (double lambda, double phi, double r, int model, uint32_t degree)
rvector	gravity_accel2 (posstruc pos, int model, uint32_t degree)
	Calculates geocentric acceleration vector from chosen model. More...
int32_t	gravity_params (int model)
	Gravitational model parameters. More...
double	nplgndr (uint32_t l, uint32_t m, double x)
	Legendre polynomial. More...
svector	groundstation (locstruc &satellite, locstruc &groundstation)
	Update Ground Station data. More...
void	hardware_init_eci (cosmosstruc *cinfo, locstruc &loc)
	Initialize Hardware. More...
void	simulate_hardware (cosmosstruc *cinfo, vector< locstruc > &locvec)
	Simulate Hardware data - multiple. More...
void	simulate_hardware (cosmosstruc *cinfo, locstruc &loc)
	Simulate Hardware data - single. More...
void	initialize_imu (uint16_t index, devspecstruc &devspec, locstruc &loc)
	Initialize IMU. More...
void	simulate_imu (int index, cosmosstruc *cinfo, locstruc &loc)
	Simulate IMU. More...
void	att_accel (physicsstruc &physics, locstruc &loc)
	Attitude acceleration. More...
int32_t	pos_accel (physicsstruc &physics, locstruc &loc)
	Position acceleration. More...
double	msis00_density (posstruc pos, float f107avg, float f107, float magidx)
	Calculate atmospheric density. More...
void	orbit_init_tle (int32_t mode, double dt, double utc, cosmosstruc *cinfo)
void	orbit_init_eci (int32_t mode, double dt, double utc, cartpos ipos, cosmosstruc *cinfo)
void	orbit_init_shape (int32_t mode, double dt, double utc, double altitude, double angle, double hour, cosmosstruc *cinfo)
void	propagate (cosmosstruc *cinfo, double utc)
void	gauss_jackson_setup (gj_handle &gjh, uint32_t order, double utc, double &dt)
	Prepare for Gauss-Jackson integration. More...
void	gauss_jackson_init_tle (gj_handle &gjh, uint32_t order, int32_t mode, double dt, double utc, cosmosstruc *cinfo, locstruc &loc)
	Initialize Gauss-Jackson orbit using Two Line Elements. More...
void	gauss_jackson_init_eci (gj_handle &gjh, uint32_t order, int32_t mode, double dt, double utc, cartpos ipos, qatt iatt, physicsstruc &physics, locstruc &loc)
	Initialize Gauss-Jackson orbit using ECI state vector. More...
void	gauss_jackson_init_stk (gj_handle &gjh, uint32_t order, int32_t mode, double dt, double utc, stkstruc &stk, physicsstruc &physics, locstruc &loc)
void	gauss_jackson_init (gj_handle &gjh, uint32_t order, int32_t mode, double dt, double utc, double altitude, double angle, double hour, locstruc &iloc, physicsstruc &physics, locstruc &loc)
locstruc	gauss_jackson_converge_orbit (gj_handle &gjh, physicsstruc &physics)
void	gauss_jackson_converge_hardware (gj_handle &gjh, physicsstruc &physics)
vector< locstruc >	gauss_jackson_propagate (gj_handle &gjh, physicsstruc &physics, locstruc &loc, double tomjd)
int	orbit_init (int32_t mode, double dt, double utc, string ofile, cosmosstruc *cinfo)
	Initialize orbit from orbital data. More...
int	orbit_propagate (cosmosstruc *cinfo, double utc)
	Load TLE's from file. More...
int	update_eci (cosmosstruc *cinfo, double utc, cartpos pos)
int32_t	triangularize (cosmosstruc *cinfo)

Variables
static double	ftl [2 *360+1]
static int	cmodel = -1
static double	coef [360+1][360+1][2]
static double	spmm [360+1]
static double	lastx = 10.
static uint16_t	lastm = 65535
static double	initialutc
static string	orbitfile
static stkstruc	stkhandle
static locstruc	sloc [MAXGJORDER+2]
static double	vc [360+1][360+1]
	Data structures for spherical harmonic expansion. More...
static double	wc [360+1][360+1]

Macro Definition Documentation

#define MAXDEGREE   360

#define ASTEP   1

Function Documentation

void gauss_jackson_init_tle	(	gj_handle &	gjh,
		uint32_t	order,
		int32_t	mode,
		double	dt,
		double	utc,
		cosmosstruc *	cinfo,
		locstruc &	loc
	)

Initialize Gauss-Jackson orbit using Two Line Elements.

Initializes Gauss-Jackson structures using starting time and position from a Two Line Element set.

Parameters

gjh	Reference to Gauss-Jackson variables.
order	the order at which the integration will be performed (must be even)
mode	Type of physical modelling. Zero is direct.
dt	Step size in seconds
utc	Initial step time as UTC in Modified Julian Days
cinfo	Reference to cosmosstruc to use.
loc	Initial location.

{
    uint32_t i;

    cinfo->node.phys.dt = dt;
    cinfo->node.phys.mode = mode;
    gauss_jackson_setup(gjh, order, utc, cinfo->node.phys.dt);
    cinfo->node.phys.dtj = cinfo->node.phys.dt/86400.;

    utc -= (order/2.)*dt/86400.;
    pos_clear(loc);
    lines2eci(utc,cinfo->tle,loc.pos.eci);
    loc.pos.eci.pass++;
    pos_eci(&loc);

    // Initial attitude
    cinfo->node.phys.ftorque = rv_zero();
    switch (cinfo->node.phys.mode)
    {
    case 0:
    //  loc.att.icrf.utc = loc.utc;
    //  loc.att.icrf.s = q_eye();
    loc.att.icrf.a = rv_zero();
    //  att_icrf(&loc);
    //  break;
    case 1:
        loc.att.lvlh.utc = loc.utc;
        loc.att.lvlh.s = q_eye();
        loc.att.lvlh.v = rv_zero();
        att_lvlh2icrf(&loc);
        break;
    case 2:
        loc.att.lvlh.utc = loc.utc;
        loc.att.lvlh.s = q_change_around_y(-DPI2);
        loc.att.lvlh.v = rv_zero();
        att_lvlh2icrf(&loc);
        break;
    case 3:
    case 4:
    case 5:
    case 6:
    case 7:
    case 8:
    case 9:
    case 10:
    case 11:
    case 12:
        loc.att.icrf.s = q_drotate_between_rv(cinfo->faces[abs(cinfo->pieces[cinfo->node.phys.mode-2].face_idx[0])].normal.to_rv(),rv_smult(-1.,loc.pos.icrf.s));
        loc.att.icrf.v = rv_zero();
        att_icrf2lvlh(&loc);
        break;
    }


    pos_accel(cinfo->node.phys, loc);
    // Initialize hardware
    hardware_init_eci(cinfo, loc);
    att_accel(cinfo->node.phys, loc);
    //  groundstations(cinfo,&loc);

    gjh.step[gjh.order2].sloc = loc;

    // Position at t0-dt
    for (i=gjh.order2-1; i<gjh.order2; --i)
    {
        gjh.step[i].sloc = gjh.step[i+1].sloc;
        gjh.step[i].sloc.utc -= dt / 86400.;
        lines2eci(gjh.step[i].sloc.utc,cinfo->tle,gjh.step[i].sloc.pos.eci);
        gjh.step[i].sloc.pos.eci.pass++;
        pos_eci(&gjh.step[i].sloc);

        gjh.step[i].sloc.att.lvlh = gjh.step[i+1].sloc.att.lvlh;
        att_lvlh2icrf(&gjh.step[i].sloc);

        pos_accel(cinfo->node.phys, gjh.step[i].sloc);
        // Initialize hardware
        hardware_init_eci(cinfo,gjh.step[i].sloc);
        att_accel(cinfo->node.phys, gjh.step[i].sloc);
    }

    for (i=gjh.order2+1; i<=gjh.order; i++)
    {
        gjh.step[i].sloc = gjh.step[i-1].sloc;
        gjh.step[i].sloc.utc += dt / 86400.;
        lines2eci(gjh.step[i].sloc.utc,cinfo->tle,gjh.step[i].sloc.pos.eci);
        gjh.step[i].sloc.pos.eci.pass++;
        pos_eci(&gjh.step[i].sloc);

        gjh.step[i].sloc.att.lvlh = gjh.step[i-1].sloc.att.lvlh;
        att_lvlh2icrf(&gjh.step[i].sloc);

        pos_accel(cinfo->node.phys, gjh.step[i].sloc);
        // Initialize hardware
        hardware_init_eci(cinfo,gjh.step[i].sloc);
        att_accel(cinfo->node.phys, gjh.step[i].sloc);
    }

    loc = gauss_jackson_converge_orbit(gjh, cinfo->node.phys);
    gauss_jackson_converge_hardware(gjh, cinfo->node.phys);

    cinfo->node.phys.utc = loc.utc;

    cinfo->timestamp = currentmjd();
}

int32_t triangularize

(

cosmosstruc *

cinfo

)

{
    cinfo->node.phys.triangles.clear();

    // Initialize with existing faces broken into triangles
    for (uint16_t i=0; i<cinfo->pieces.size(); ++i)
    {
        for (uint16_t j=0; j<cinfo->pieces[i].face_cnt; ++j)
        {
            facestruc tface = cinfo->faces[cinfo->pieces[i].face_idx[j]];
            cinfo->node.phys.vertices.push_back(tface.com);
            cinfo->node.phys.vertices.push_back(cinfo->vertexs[tface.vertex_idx[tface.vertex_cnt-1]]);
            for (uint16_t k=0; k<tface.vertex_cnt; ++k)
            {
                cinfo->node.phys.vertices.push_back(cinfo->vertexs[tface.vertex_idx[k]]);
                trianglestruc ttriangle;
                ttriangle.pidx = j;
                ttriangle.normal = tface.normal;

                ttriangle.tidx[0] = cinfo->node.phys.vertices.size() - (k+2);
                ttriangle.tidx[1] = cinfo->node.phys.vertices.size() - 2;
                ttriangle.tidx[2] = cinfo->node.phys.vertices.size() - 1;

                ttriangle.com = cinfo->node.phys.vertices[ttriangle.tidx[0]].cross(cinfo->node.phys.vertices[ttriangle.tidx[1]]);
                ttriangle.area = cinfo->node.phys.vertices[ttriangle.tidx[0]].area(cinfo->node.phys.vertices[ttriangle.tidx[1]]);
                ttriangle.heat = cinfo->pieces[i].heat * ttriangle.area / cinfo->pieces[i].area;

                cinfo->node.phys.triangles.push_back(ttriangle);
            }
        }
    }

    // Loop, breaking triangles into ever smaller triangles, until we are below required resolution (5mm squared)
    bool modified = true;
    while (modified)
    {
        modified = false;
        for (size_t i=0; i<cinfo->node.phys.triangles.size(); ++i)
        {
            trianglestruc tface = cinfo->node.phys.triangles[i];
            if (tface.area > 2.5e-5)
            {
                cinfo->node.phys.vertices.push_back(tface.com);
                trianglestruc ttriangle = tface;
                ttriangle.tidx[0] = cinfo->node.phys.vertices.size() - 1;

                ttriangle.tidx[1] = tface.tidx[2];
                ttriangle.tidx[2] = tface.tidx[0];
                ttriangle.com = cinfo->node.phys.vertices[ttriangle.tidx[0]].cross(cinfo->node.phys.vertices[ttriangle.tidx[1]]);
                ttriangle.area = cinfo->node.phys.vertices[ttriangle.tidx[0]].area(cinfo->node.phys.vertices[ttriangle.tidx[1]]);
                ttriangle.heat = tface.heat * ttriangle.area / tface.area;
                cinfo->node.phys.triangles.push_back(ttriangle);

                ttriangle.tidx[1] = tface.tidx[0];
                ttriangle.tidx[2] = tface.tidx[1];
                ttriangle.com = cinfo->node.phys.vertices[ttriangle.tidx[0]].cross(cinfo->node.phys.vertices[ttriangle.tidx[1]]);
                ttriangle.area = cinfo->node.phys.vertices[ttriangle.tidx[0]].area(cinfo->node.phys.vertices[ttriangle.tidx[1]]);
                ttriangle.heat = tface.heat * ttriangle.area / tface.area;
                cinfo->node.phys.triangles.push_back(ttriangle);

                ttriangle.tidx[1] = tface.tidx[1];
                ttriangle.tidx[2] = tface.tidx[2];
                ttriangle.com = cinfo->node.phys.vertices[ttriangle.tidx[0]].cross(cinfo->node.phys.vertices[ttriangle.tidx[1]]);
                ttriangle.area = cinfo->node.phys.vertices[ttriangle.tidx[0]].area(cinfo->node.phys.vertices[ttriangle.tidx[1]]);
                ttriangle.heat = tface.heat * ttriangle.area / tface.area;
                cinfo->node.phys.triangles.push_back(ttriangle);

                modified = true;
            }
        }
    }

    // For each triangle, identify what is in its field of view in 10 degree resolution grid
    for (size_t i=0; i<cinfo->node.phys.triangles.size(); ++i)
    {
        trianglestruc tface = cinfo->node.phys.triangles[i];
        if (cinfo->node.phys.triangles[i].triangleindex.empty())
        {
            cinfo->node.phys.triangles[i].triangleindex.resize(9);
            for (size_t j=0; j<9; ++j)
            {
                float cel = cos(RADOF(90.-(10.*j+5.)));
                cinfo->node.phys.triangles[i].triangleindex[j].resize((int)(cel*18+.5));
            }
        }

        for (size_t j=0; j<cinfo->node.phys.triangles.size(); ++j)
        {
            if (j != i)
            {
                double sep = cinfo->node.phys.triangles[i].com.separation(cinfo->node.phys.triangles[j].com);
                if (sep > DPI2)
                {
                    float az;
                    float el;
                    Vector topo;
                    body2topo(cinfo->node.phys.triangles[i].com, cinfo->node.phys.triangles[j].com, topo);
                    topo2azel(topo, az, el);
                    float alt = (int16_t)((DPI2 - el) / 9.) + RADOF(.5);
                    uint16_t rowi = 9 * alt / DPI2;
                    uint16_t coli = 9 * cos(alt) * az / DPI2;
                    if (cinfo->node.phys.triangles[i].triangleindex[rowi][coli])
                    {
                        cinfo->node.phys.triangles[i].triangleindex[rowi][coli] = j;
                    }
                    else
                    {

                    }
                }
            }
        }
    }

    return cinfo->node.phys.triangles.size();
}

Variable Documentation

double ftl[2 *360+1]

static

int cmodel = -1

static

double coef[360+1][360+1][2]

static

double spmm[360+1]

static

double lastx = 10.

static

uint16_t lastm = 65535

static

double initialutc

static

string orbitfile

static

stkstruc stkhandle

static

locstruc sloc[MAXGJORDER+2]

static

double vc[360+1][360+1]

static

Data structures for spherical harmonic expansion.

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

double wc[360+1][360+1]

static