Performance_Work_Statemen_060420024t__updated_ICs_CLIN0001.docx

The National Geodetic Survey

Program to Update Initial Conditions and Compute Force Model Parameters for GNSS Orbits

I. SCOPE

The National Geodetic Survey (NGS) is a program office under the National Ocean Service (NOS). NGS provides the framework for all positioning activities in the Nation. The foundational elements - latitude, longitude, elevation, shoreline information and their changes over time - contribute to informed decision making and impact a wide range of important activities including mapping and charting, navigation, flood risk determination, transportation, land use and ecosystem management. NGS' authoritative spatial data, models, and tools are vital for the protection and management of natural and manmade resources and support the economic prosperity and environmental health of the Nation.

In providing this service for the Nation, NGS routinely processes vast amounts of data from Global Navigation Satellite Systems (GNSS) to compute GNSS satellite orbits. Historically, NGS has had capabilities to use only the U.S. Global Positioning System (GPS) data. Adding capability for multi-GNSS processing is an important part of NSRS modernization, allowing NGS to provide orbits for all available constellations. To add this crucial capability, NGS has developed its own multi-GNSS processing software suite, including libraries that can support precise orbit determination.

Currently, NOAA’s production orbit determination software uses a series of old Fortran codes to compute satellite orbits, and only for the GPS constellation. This contract will modernize one portion of the orbit production process, switching over to code written in C++ that leverages a larger suite of software libraries used in NOAA. The specific transition of the code to C++ is CLIN0001.

II. GLOSSARY OF TERMS AND REFERENCES

G-File - A model for satellite initial conditions (IC) and non-gravitational force parameters

· ASCII file

· Format description in Appendix A of this PWS

· Additional historic documentation at https://geodesy.noaa.gov/GRD/GPS/DOC/arc/g_file.html

T-File - Tabular ephemeris

· Binary file

· Format description in Appendix B of this PWS

· Read/write capabilities already provided by BinaryEphemeris class

NRM - normal equations

· Binary file

· Format description https://geodesy.noaa.gov/GRD/GPS/DOC/gpscom/neqfile.html

· Write capabilities already provided by NrmFile::writeBinaryFile()

SP3 - Special Product version 3 for precise ephemeris

· ASCII file

· Format description https://files.igs.org/pub/data/format/sp3c.txt

· Read and interface capabilities already provided by Ephemeris class

III. DESCRIPTION OF WORK

The expected outcome of this contract will result in a new C++ program to form and solve a linearized system of equations (Appendix C). Inputs to the program include a precise Standard Product 3 “sp3” orbit file, to be used as the observations (O), and a “T-file” containing a set of a priori estimates for satellite initial conditions, a model of computed satellite orbits (C), and partial derivatives (Aij). File format descriptions are provided in the appendices. The program will compute pre-fit residuals (observation minus computed, or O-Cs) and construct a system of equations (design matrix A) using the partial derivatives (A = [Aij]; see Appendix C) and a set of a priori parameter constraints. The pre-fit residuals and design matrix will be used to form a normal equation system, which will be written out in a binary format to a file on disk. The system of equations will then be solved using a constrained least-squares algorithm. The computing requirement for solving the system of equations will be a fast, numerically stable method. The resulting solution vector shall contain adjustments to the initial conditions and estimates for non-gravitational force model parameters, so that the updated initial conditions and newly estimated force model parameters are consistent with the precise orbit file. The program shall output a G-file with the updated initial conditions and force model parameter estimates.

To support the work, NGS will provide C++ object-oriented libraries. The contractor will use existing libraries, and also develop new C++ libraries and objects as needed to accomplish the goals of the program. New C++ libraries shall use the object-oriented approach demonstrated in the other NGS code. At minimum, a new C++ library (proposed name “SatIC”) shall be written, to read and write initial conditions in the G-file format.

The new code shall be capable of being compiled and run on an HPC environment, such as Mississippi State High Performance Computing systems. The M-PAGES code (C++ libraries) is already compiled and operational in this Linux environment. NOAA/NGS will assist the awardee with securing system account(s), so that their development and testing can be conducted in this Linux environment.

The design of the new C++ program shall be:

Inputs

· sp3 file: an ascii file containing the precise satellite data (O)

· T-file-input: a binary format file containing the computed satellite orbit models (C), and partial derivatives (Aij).

Existing NOAA NGS Libraries/Classes to support code development

Key steps in the program Format descriptions are either in the appendices or linked in the Glossary/References section

1. Read input raw sp3 file using the Ephemeris class from the NGS library suite

2. Read input raw T-file using the BinaryEphemeris from the NGS library suite

3. Retrieve satellite XYZ from input raw sp3 file at the epochs defined within the T-file, using the Ephemeris class methods from the NGS library suite

4. Form observed minus computed (O-C) prefit residuals: O-C = interpolated_sp3_XYZ - Tfile_XYZ

5. Form design matrix for observation equations using partials from BinaryEphemeris from the NGS library suite and a set of parameter constraints

6. Write out normal equations as NRM files, one per satellite, using the NrmFile class

7. Solve the constrained least squares problem using a fast, numerically stable solver; note that the provided NGS C++ code does include a least squares solver

8. Solution to system of equations = adjustments to satelliteIC and a set of non-gravitational force model parameters.

9. Load the updated satellite IC and non-gravitational force model parameters into the BinaryEphemeris object

10. Write the BinaryEphemeris to a G-file.

Outputs

· NRM files

· Output G-File

The following additional materials will be provided to the awardee:

· All C++ libraries mentioned above, as well as programs which use these libraries for similar purposes and makefiles which successfully compile those programs

· HTML documentation for C++ libraries

· NGS coding standards document

· Example input files (sp3, T)

· Example output files (NRMs, G) that correspond to the input files

· Pseudocode for the beginning for the new C++ program, to demonstrate some of the methods available through objects in the provided C++ libraries

IV. TECHNICAL COMPETENCIES

A technically acceptable contractor will demonstrate the following areas of expertise:

· Understanding of all the file formats used as input and output

· Ability to write new code in C++ using provided libraries

· Proficiency with compiling and running code in a Linux environment, including an HPC environment

· Knowledge of coding linear algebraic algorithms in C++, including the efficient construction and manipulation of large matrices, and matrix multiplication and factorization methods.

· Similarly, knowledge and previous experience for construction and solving of sets of normal equations in C++ code.

V. TASKS AND EVALUATION

Under the guidance of the NGS Project Manager, perform the following tasks:

1. Modify and extend existing BinaryEphemeris class to meet the needs of this contract

a. Add getters and setters to provide interface to the private attributes (do not simply change the attributes to public)

b. Class functionality shall be verified through the unit tests (see /Lib/Ephemeris/UnitTest/unit_test.cpp)

2. Modify and extend existing SatCoords class to meet the needs of this contract

3. Write a new C++ class called SatIC with the following functionality. This new class should be similar in function and behavior to BinaryEphemeris, and use BinaryEphemeris objects to hold appropriate data

a. Read G-File

b. Write G-File

c. Data are contained in private attributes

d. Getters and setters to access and modify protected attributes

4. Write new C++ program, using specified libraries, that will adjust the initial conditions of satellite orbits using an input precise orbit file as the data

a. Compute the difference between an a priori tabular ephemeris (T-File), and a sp3 ephemeris (treated as data).

b. Form matrices: O-C prefit residuals and a design matrix (A)

c. Use A and O-C to form and solve normal equations with parameter constraints

d. Apply the solved parameters (output adjustments to satellite initial conditions and non-gravitational force model parameters)

e. Create G-File.

5. Provide a Makefile for the new C++ program, in the style of other makefiles used in the example programs, so that the program can be compiled from the command line in the Linux shell using “make”, without the use of an IDE.

Successful completion of the CLIN0001 task shall be achieved when the following conditions are met:

1) Using the NGS example input files, the contractor will use the new code to produce output G-files that match the NGS example output G-files to the following levels:

a) within 1 millimeter for position.

b) within 1e-3 millimeters per second for velocity

c) and within 1e-4 (unitless) for force model parameter estimates.

2) The processing time for the full GPS constellation (32 satellites) using the contractor’s new program shall be less than 5 seconds for any set of sp3 and T-files.

VI. IT SECURITY REQUIREMENTS AND GOVERNMENT FURNISHED MATERIALS

NOAA will provide access to the contractor to work a NOAA HPC environment, with compiled M-PAGES libraries and source code on disk.

VII. DELIVERABLES

VIII. TECHNICAL POINT OF CONTACT

The Technical Point of Contact is Andria Bilich, who will be responsible for technical guidance to the contractor and should a dispute arise shall inform the Contracting Officer. This is not a personal services contract and the Technical Point of Contact or no other Federal employee will personally “supervise” the contractor.

APPENDIX A. G-FILE FORMAT DESCRIPTION

The G-file consists of two formatted header lines, one or more comment lines, and initial conditions and force parameters for one or more satellites. Shown below is a partial listing of a typical G file:

24 113 12 00 00 GPST J2000 IAU76 BERNE

15 X Y Z XDOT YDOT ZDOT DRAD YRAD BRAD DCOS DSIN YCOS YSIN BCOS BSIN

reformat_llsolv(01/04/26) GENERATED ON 24/05/01 20:04

END

PRN 02

17880.73992726 -11643.34016308 -15155.52556841 2.880018516541 1.389444240844 2.286403370505 1.02671483 0.00135917 -0.00287208 0.00715075 0.02669477 0.00000000 -0.00000000 -0.00317177 -0.00744675

PRN 32

8073.92271087 14514.69891274 -20559.13048555 -2.467441002448 2.838424248527 1.004838276901 0.99703423 0.00083512 -0.01449323 0.00000000 -0.00000000 0.00000000 -0.00000000 -0.00164743 -0.01180184

END

The required first line gives the epoch of the ICs in GPST or UTC (year, day-of-year, hours, minutes, seconds), followed by the time type (GPST or UTC), the inertial frame for the ICs (J2000 or B1950), the precession constant used (IAU76 or IAU68), then the model for direct solar radiation accelerations (here, BERNE stands for the Bernese model described in Springer et al. [1998]).

The required second line gives the number of ICs, followed by labels for the ICs and force model parameters contained within the G-file. Note that the same order of parameters will be used for each satellite.

The header may include one or more optional comment lines.

The header is terminated by END.

In the rest of the file, satellite-specific blocks begin with ‘Gnn’ or ‘PRNnn’, where nn is the PRN and the spacecraft body name. After that follow the data for each satellite, including the position (X, Y, Z), velocity (XDOT, YDOT, ZDOT), and non-gravitational force model parameters (DRAD, YRAD, BRAD, DCOS, DSIN, YCOS, YSIN, BCOS, BSIN).

A description of the non-gravitational force model parameters is provided in the following table.

cf. Beutler et al., Man. Geod. 19, 367, 1994

APPENDIX B. T-FILE FORMAT DESCRIPTION

Read/write capabilities for this format are provided in the BinaryEphemeris class.

Record 1 header, mdyic, hmsic, mdyb, hmsb, mdye, hmse, sdelte, nepche, ut1utc, xwob, ywob

Record 2 comt, naveph, nintre, (sname(j), (avic(i,j),i=1,nics),j=1,naveph)

Record 3 (repeated for all epochs) yye( i, j, k ), j= 1, nintre), k= 1, nsveph ) yye(par,svs) values at this epoch 3060 x 8-byte double where par = 3 + 3 x ics.

The yye(par) elements are provided rowwise according to the following table, so that element1=X, element2=Y, element3=Z, element4=dX/dXic, … element48 = dZ/dBsin.

This set of 48 elements is repeated for each satellite in the file.

See the table in Appendix A for definitions of the model parameters (X, Y, X, XDOT, YDOT, ZDOT, DRAD, YRAD, BRAD, DCOS, DSIN, YCOS, YSIN, BCOS, BSIN).

Appendix C: Observation equations and Partial Derivatives

Notation and terminology Below we use X, Y, Z to signify numeric Cartesian coordinates of a satellite with respect to an inertial coordinate system at a single epoch of time. The true coordinates of any satellite cannot be known exactly. However, we may infer the coordinates of a satellite based on observational data or theoretical calculations. Both observations and theoretical calculations are subject to errors and uncertainty.

· We will designate an observation of the X coordinate of a satellite as

· Similarly, we will refer to a theoretically computed value of the X coordinate of a satellite as .

· We use O and C (without subscripts) represent the total collection of values and at multiple epochs.

· A represents a matrix.

· “Partials” is shorthand for “partial derivatives”.

· Delta terms ( such as are the terms we need to determine (corrections to the initial condition state vector, and estimates of force model parameters).

Equations Below, we provide detailed terms for the X coordinate. The same form of equations apply to the Y and Z coordinate, naturally using the appropriate and partials instead.

At a single epoch, for a single satellite, the observation equations are written as:

15 terms of the form 15 terms of the form where are the “observed” satellite positions obtained from the sp3 file. are the computed satellite positions obtained from the T-file, and the 15 delta terms ( such as are the terms we need to determine (corrections to the initial condition state vector, and estimates of force model parameters). Note that although the observations and computed values are available at multiple epochs, but there is only one set of 15 that does not depend on time. We want to estimate these delta terms using data from multiple epochs within a day, using constrained least squares or other solution type.

In an unconstrained least squares problem, the partials form the design matrix A in the adjustment equation

The vector e represents the unknown errors that we seek to minimize. The dimension of the A matrix is For a constrained least-squares problem, we also incorporate parameter constraints of the form where 𝛏 are a priori values for the adjustments, typically 𝛏=0 and 𝛎 represents errors in the a priori values. Typically, tight constraints will only be applied to either two or four force model parameters depending on the satellite type, but the code should provide the option to constrain any (or all) of the parameters. These constraints may be treated as pseudo-observations and concatenated to the system of equations A, forming a new augmented system Aaug.

The normal equation system is constructed as where represents the matrix transpose of A and W is a weight matrix (that must be specified) that takes into account the errors in the parameter constraints as well as the errors in the observations, and and are the augmented observations and model vectors, i.e., The matrix in the constrained case is similarly an augmented version of the design matrix, which includes a dimensional identity matrix that associates the a priori values to the corresponding elements in the vector.

The solution is obtained by solving the system for the best estimate of .

Once these are estimated, computation of the corrected, post-adjustment X coordinates (X’) are computed from the a priori value of X plus the terms of the observation equation:

X’ = X + + + + Y’ = Y + 15 dY terms Z’ = Z + 15 dZ terms where the estimated are indicated by a prime (‘) with each term.

The success of the estimation will be assessed based on the fit of the model X’ to the observations O. That is, the success will be measured by the root-mean-square of the post-fit residuals O - X’.

Although the above equations detail only the X-coordinate, the Y- and Z-coordinates must also be included in the adjustment.

The final product of this estimation is to create a new G-file with updated satellite initial conditions and estimated force model parameters. The updated SV IC and force model params which are written to the new G-file are computed as:

where (no a priori value) where the naught (0) terms are values from the input T-file.

image1.png