prmred - nicolob/pep-full GitHub Wiki
Type: subroutine
Defined in: prmred.f
Definition: subroutine PRMRED(in0,sepeat,moprnt,iseq1,nstop,jdpad)
ICT(23) binary coded control for atmospheric delay model selection
1 bit=1: use new mapping function code (Mendes & Pavlis, 2002)
=0: use old mapping function code
2 bit=1: use new zenith delay code (Mendes & Pavlis, 2004)
=0: use old zenith delay code
IN =input data and control constants
IOUT =print output
JOUT =additional print output for analiz link to get extra printout
of correlation matrix and parameter adjustments (if zero, no
such printout)
KOUT =additional print output for extended print output in subroutine
SBOUT for SBFN integration of probe or planet (if zero, no
such printout)
LOUT =enabling switch for extra input link print of adjustable
parameter summary and a priori information (if zero, no such)
MOUT =extra "print" data set for summary of program flow, essentially
a digest of important data on iout, designed for speedy
on-line examination via time-sharing terminals. (if zero, no
such output)
NOUT =extra "print" data set for list of observations (with errors
and residuals) processed in compar and prdict links. designed
for on-line processing and/or plotting via "smart" video
terminals. (if zero, no such output)
IPUNCH=punch output
JPUNCH=extra punched output -- new PEP input deck after each
iteration other than the last iteration
KPUNCH
IGRAPH=graphical output
INTERN=internal buffer
IBUF =general purpose buffer
IMAT =saved inverse of coefficient matrix of normal equations and
solution vector calculated in SNORML (output)
(IMAT=0 don't save, IMAT.GT.0 do save)
used for series-by-series adjust
JMAT =same as IMAT but input to SOLFRM if ICT(5).gt.1
IPERT =pertubing planets and Moon positions
JPERT =additional perturbing body positions (asteroids or satellites)
KPERT =positions of asteroid system c-of-m
IPERT2=moon tape used as source of moon for nbody integration,
superseding IPERT
IMAT0 = input saved normal equations and right hand side
= JMAT0+ 100 * JMATA with JMAT0 used as data set number
and with JMATA a flag for type of restoration of
normal equations - 0 normal restore, 1 skip rhs, 2 skip lhs
data sets IMAT0(j),j=1,NUMMT0 are used if NUMMT0 .gt. 0
.and. NUMMT0 .le. maxmt0
IMAT1 = output saved normal equations and right hand side
(series by series)
used for series-by-series adjust (temp)
IMAT2 = output total saved normal equations and right hand side
(must be assigned if restoring normal equations with imbedded
a priori constraint matrices, unless JCT(60)>0)
IMAT3 = output total saved partially pre-reduced normal eqns
also used as input for de-correlation of partials. see ICT(11)
IMAT4 unassigned
LIBHOC= input data set of lunar libration corrections
ICTAT = input integrated CT-AT data set
IENG =data set for gas leak and other man-made-force data
JENG =data set for radiation pressure data
KENG
IEM =Earth-Moon barycenter coordinates,partials
IMN =Moon coordinates,partials
INUT =Earth nutation
ILIB =Moon physical libration, partials
IPLNT(i)=planet NPLNT(i) coordinates,partials (i=1,numpln)
IPLCON=Earth-Moon barcenter,Moon,Earth rotation,Moon rotation,planet
input data not in labeled commons
IVCNT =velocity of center of mass of solar system relative to sun
(INUT,IVCNT not used. Nutation & libration come from Moon or
n-body tape normally.)
IOBS =observation cards (in sysin stream following input data and
control constants if IOBS=5, which is default value. if
ICT(1)=-1 or IOBS input 0, no observation cards read by
program. if IOBS input as 8 or some other number not 0 or 5,
then observation cards on tape)
1 radar observations of time delay and/or Doppler
shift
2 radio tracking observations of time delay and/or
Doppler shift
3 differential radar observation
(feature relative to subradar point if spot name
not blank, bandwidth if spot name is blank)
(If observing sites on Earth, observables 1,2,3 processed
in radar link. If observing sites on non-Earth body,
observables 1,2,3 processed in STRADR link.)
For observations made from the Earth
(processed in OPTIC link)
4 meridian circle observations of right ascension and
declination referred to the true equinox and equator
of date for observed body not Earth,Moon satellite
(geocentric)
4 azimuth and elevation observations of Earth
or Moon satellites (topocentric)
5 photographic observations of right ascension and
declination referred to the mean equinox and equator
of reference epoch (elliptic aberration not removed:
astrometric) (topocentric)
6 same as 5 but with all aberration removed
(astrographic)
For observations made from a non-Earth body
(processed in STOPTC link)
4 look angles in spacecraft coordinate system
5,6 angles relative to star backgound (astometric,
astrographic -- but abberation not quite right?)
13 time delay measurement for signal sent from ground
to a spacecraft to another spacecraft and then
back to the ground
14 time delay measurement for signal sent from ground
to a spacecraft to another spacecraft back to the
first spacecraft and then back to the ground
cols 3-5 Integer (i3 format) giving planet number of observed body
0 Sun 6 Saturn
1 Mercury 7 Uranus
2 Venus 8 Neptune
4 Mars 9 Pluto
5 Jupiter 10 Moon
11,...,30 natural asteroids and satellites
31,... artificial space probes
negative stars
cols 17-20 First four characters (a4 format) of eight character
sending site name (ignored for optical observations)
for satellite mutual events, distinguishes between eclipse
and occultation.
cols 22-25 Four characters (a4 format) giving name spot observed on
given body (if blank, subradar point for rader observations
or center of body for optical observations). However, if
planet number of observed body is negative, then this is
four character star name.
For radar observations of a planet, spot name = '&&&&'
denotes observations away from subradar point with spot
longitude and latitude on observation card for each obs.
For NCODF=13,14 spot name is first four characters of
second spacecraft name.
Also, for ncodf=7-9 spot name is first four characters
of second body name.
cols 26-31 single precision number (e6.0 format) by which error of
first measurement (time delay or right ascension) of every
observation in series is to be multiplied for use in
forming normal equations.
cols 32-37 Single-precision number (e6.0 format) by which error of
second measurement (Doppler shift or declination) of every
observation in series is to be multiplied for use in
forming normal equations.
cols 38-43 Single-precision number (e6.0 format) giving accuracy
criterion for iterations in processing observations in
series (iterations to determine reflection and send times
given receive time for radar observations, iteration to get
body over meridian for meridian cicle observations, etc.).
The units and meaning depend on the observable type. For
most types, the value represents a change of distance in
light-seconds, but it gives a change in the time of transit
or occultation (in seconds) or a change in Earth rotation
phase (in radians) for meridian-circle observations. If
the value is omitted, 0.001 is used as a default.
cols 44-45 Integer (i2 format) giving type of time epoch for
observations in series.
Radar observations: 1 WWV (UTC) epoch
0 UT2 epoch
-1 (A.1-epoch) and (UT1-epoch) are
on observation card
2 same as -1
3 epoch is send time
5 epoch is CT, normal point
pseudo-observations.
Optical observations: 1 epoch in 20th century for which
atomic time exists
0 epoch in 20th century for which
atomic time does not exist
-1 epoch in 1800's, etc.
(needed because only last two
digits of year are on card)
All observations: 10*i+j First two digits of the year are
19+i (for i.ge.0), j interpreted
as above
From this point on, 2 formats are possible, one using full 80 columns,
the other stopping at column 72. For IOBS cards, the format selection
is based on JCT(69); for dummy obs and error overrides, the format
selection can be locally changed by insertion of '#72' or '#80'
cards before any series.
cols 46-52 Single-precision number (f7.1) format) giving offset*1e10
from A.1 atomic time of unit of time delay measurement in
radar observation series.
In the 72 column format this field appears on an optional
2nd header card (present if column 6 contains ':').
cols 71-72 Integer (i2 format) indicating if ephemeris tapes (usually
or 64-65 backward in time) are to be rewound before processing
observations in series.
-1,0 tapes not to be rewound
1,2 tapes to be rewound
if 0,1 next card after header card(s) of series is
observation data
if -1,2 next card after header card(s) of series is
site card(s) then come observation data
cols 73-75 Integer (i3 format) giving tape sequence number for output
or 66-68 observed minus theory and partial derivative tape for
observation series. (Must be positive and non-decreasing
from one observation series to the next.)
cols 76-80 Integer (i5 format) giving series sequence number (must be
or 69-72 positive and increasing from one observation series to the
next for given value of tape sequence number).
In the 72 column format this field is i4.
col 1 Flag for optional further header cards.
cols 2- 8 Same as cols 46-52 of 80 column format above.
cols 9-16 Eight-character name (a8 format) of star catalog used in
reducing the observations in this series.
col 1 Integer (i1 format) giving type of observation
1 radar observation of time delay
2 radar observation of time delay and Doppler shift
3 radar observation of Doppler shift
3 bandwidth if ncodf=3, spot name blank
cols 2-4 Integer (i3 format) giving hour of receive time epoch
cols 5-7 Integer (i3 format) giving minute of receive time epoch
cols 5-15 Single-precision number (f8.4 format) giving second of
receive time epoch.
cols 16-28 Double-precision number (f13.7 format) giving observed time
delay in seconds.
cols 29-35 Single-precision number (e7.0 format) giving error of time
delay measurement in seconds.
cols 36-50 Double-precision number (d15.8 format) giving observed
Doppler shift in cycles per second.
cols 51-57 Single-precision number (e7.0 format) giving error of
Doppler shift measurement in cycles per second.
cols 58-65 Single-precision number (f8.4 format) giving A.1 time minus
epoch (if needed).
cols 66-72 Single-precision number (f7.4 format) giving UT1 time minus
epoch (if needed).
cols 73-74 Integer (i2 format) giving month observation (from 1 to 12)
cols 75-76 Integer (i2 format) giving day of observation (from 1 to
31)
cols 77-78 Integer (i2 format) giving last two digits of year of
observation (see also cols 44-45 on the series header).
col 79 Character (a1 format) indicating observing site (not read,
for labeling purposes only).
col 80 Character (a1 format) indicating observed body (not read,
for labeling purposes only).
col 1 integer (i1 format) giving type of observation
1 or 4 observation of right ascension
2 or 5 observation of right ascension and declination
3 or 6 observation of declination
cols 2-15 give the approximate universal time of transit of the
observed body over the meridian of the observatory (within
twelve hours of the actual time of transit). the center is
observed unless first card of series indicates the limb is
observed.
cols 2-4 integer (i3 format) giving approximate hour of transit
cols 5-7 integer (i3 format) giving approximate minute of transit
cols 8-15 Single-precision number (f8.4 format) giving approximate
second of transit
col 16 blank (not used)
cols 17-30 give the right ascension of center from observation unless
first card of series indicates the limb is observed
cols 17-19 integer (i3 format) giving hour of right ascension
cols 20-22 integer (i3 format) giving minute of right ascension
cols 23-30 Single-precision number (f8.4 format) giving second of
right ascension
cols 31-37 Single-precision number (f7.4 format) giving error of
right ascension in seconds of time
cols 38-50 give the geocentric declination of the center from
observation unless first card of series indicates the limb
is observed
cols 38-40 integer (i3 format) giving degree of declination
cols 41-43 integer (i3 format) giving minute of declination
cols 44-50 Single-precision number (f7.3 format) giving second of
declination
cols 51-57 Single-precision number (f7.3 format) giving error of
declination in seconds of arc
cols 58-66 blank (not used)
col 67 character (a1 format) indicating clamp
cols 68-69 characters (2a1 format) indicating limb of right ascension
and declination (used in computing o-c if published
observation not reduced to center)
cols 70-71 characters (2a1 format) indicating observer
col 72 blank (not used)
cols 73-74 integer (i2 format) giving month of observation
(from 1 to 12)
cols 75-76 integer (i2 format) giving day of observation
(from 1 to 31)
cols 77-78 integer (i2 format) giving last two digits of year of
observation (the century of observation is given on first
card of observation series)
col 79 character (a1 format) indicating observing site (not read,
for labeling purposes only)
col 80 character (a1 format) indicating observed body (not read,
for labeling purposes only)
IOBCON=data for generation of dummy observations and deletion and
error weight alteration of real observations (written in
subroutine dltred of input link)
numobt=number of observation library tapes (must be between 1 and 10
inclusive). then for i=1,numobt, we have
IOBS0(i) = input observation library tapes for first least squares
iteration
IOBS1(i) = output observation library tapes on odd least squares
iterations, input on even iterations
IOBS2(1-10) =output observation library tapes on even least squares
iterations, input on odd iterations greater than first
Observations to be processed are on up to 10 observation library files
plus input cards. New libraries are written with input card data
inserted at appropriate points. Besides information on cards, library
files have observed minus theoretical values of observations computed
in run which produced tapes, data concerned with observations which
stay essentially constant from run to run (e.g., precession-nutation
matrix), and partial derivatives of observations. Partial derivatives
on newly written observation library tapes include all that were on
input library tapes plus others needed in forming normal equations for
given run. Partial derivatives on observation library tapes whose
values do not change signifigantly between iterations are not
recalculated, whereas those whose values might change signifigantly
(e.g., those with respect to planetary shape parameters) are
recalculated before incrementing normal equations. There is option
(ICT(80)=1 with iterat=1) of using observed minus theoretical values
on library tapes in forming normal equations rather than computing the
values. Data set IOBCON contains data conrolling (1) alterations in
error weightings for use in the normal equations of whole observation
series or of individual observations within an observation series, (2)
deletion of specific observations or whole observation series from the
newly written library tapes, (3) overriding constants concerned with
observation series in addition to error weighting constants, and (4)
calculation of theoretical value of dummy observations and calculation
of covariance matrix for dummy observations with or without additional
real observations.
IOBS = observation card data set. first card contains information
about observation series, perhaps followed by additional such
cards. then come cards giving observation data for the
observation series, one card per observation. end of observa-
tion series is signaled by a blank card. then comes next obser-
vation series, etc. end of all the observation series (and end
of data set IOBS) is signaled by two blank cards, one termina-
ting last observation series and an extra one terminating all
the observation series.
IOBS = 0, no observation card data set. in addition IOBS is not read
if ICT(80).gt.0. if ICT(80).eq.0 and IOBS.gt.0, then IOBS is
read on first least squares iteration (iterat=1), but IOBS
will be set equal to zero when
end is reached on iteration iterat=1 if ICT(80).eq.0 so that
IOBS will not be read on subsequent iterations, since the data
on IOBS will have been written on the observation library tapes
to be read on subsequent iterations. ICT(80).lt.0 is signal
that there are no input observation library tapes, program is
to read observation data only from IOBS, rewinding IOBS between
iterations, in which case partial derivatives have to be
recalculated each iteration.
in subroutine compar
iabs1= input observation library data set. if 0, no such input.
iabs2= output observation library data set. if 0, no such output.
there is a loop on jtape from 1 to 10 to define values of iabs1,iabs2.
if iterat=1, then iabs1=IOBS0(jtape)
iabs2=IOBS1(jtape)
if iterat.gt.1 is even,iabs1=IOBS1(jtape)
iabs2=IOBS2(jtape)
if iterat.gt.1 is odd, iabs1=IOBS2(jtape)
iabs2=IOBS1(jtape)
jtape is incremented either (1) if end has been reached on iabs1 with
iabs1.gt.0 and either IOBS=0 or ntapa(2).gt.ntapa(3), or (2) if
iabs1=0 and ntapa(2) makes a jump to a larger number. if ICT(80).ge.0,
for the first least squares analysis iteration, input library tapes
IOBS0 are read(along with cards IOBS if ICT(80).eq.0) and IOBS1 is
written. for the second
iteration, IOBS1 (written in first iteration) is read and IOBS2 is
written. in the third iteration IOBS2 is read and IOBS1 is written.
and so it continues. the original input library tapes IOBS0 are
unchanged, but the input-output library tapes IOBS1,IOBS2 are
alternately read and written, the program eating its tail from
iteration to iteration.
if ICT(80).lt.0, then only observation cards IOBS are read and
iabs2=IOBS1(1) is written from iteration to iteration.
if ICT(80).ge.0 saved partial derivatives are used from input
observation library tapes where possible.
if ICT(80).lt.0 partial derivatives are recalculated each iteration.
Global deletions are performed when Tdlt0, Tdlton, and Tdltof are
specified. Beginning at epoch Tdlt0 (julian ephemeris date) input
data are alternately processed normally for Tdlton days and skipped
for Tdltof days regardless of the other data editing controls in PEP.
This periodic windowing extends both forward and backward in time
from Tdlt0. For input observation cards, the skipped data are given
an overriding weight of zero, but such data from observation libraries
are skipped altogether.
epsa no longer used
ibuf1 = data set for parameter names. if > 0, then names will appear
next to correlations in analiz link printout.
ibuf2 = apriori b matrix
ibuf3 saved right hand sides of normal equations for
series-by-series adjust
ibuf4 partial solutions of series-by-series adjust
ibuf5 multiple parameter set requests, pseudo-iptr's (iptrx's)
ibuf6 parameter(and associated) values (by run) for
multiple parameter set (or multiple epoch) runs
also used for pprne generation (see arrwrt)
ibuf7 temporary data set for partial prereduction of normal eqns
ibuf8 temporary data set for partial prereduction of normal eqns
ibuf9 temporary data set for partial prereduction of normal eqns
ibuf10 unassigned
WGTOBS - Array of weighting factors applied to obslib tapes (IOBS0,
IOBS1, or IOBS2) in forming normal equations (on top of any
error weight information). The sense of weighting is such
that wgtobs>1 increases the weight of the data (as if all
measurement uncertainties were divided by 'wgtobs'). A
value of zero causes the corresponding obslib to be skipped.
WGTMT0 - Array of weighting factors applied to input saved normal
equations (on top of any error weight information).
The sense of weighting is the same as for wgtobs above.
A value of 0 causes the corresponding imat0 to be skipped.
A negative value causes the data on the corresponding imat0
to be deducted from the normal equations and error
statistics, i.e., reduces the count of measurements and
sum-squared residuals and subtracts information from the
normal equations. This is meaningful only if some superset
of the same data is being included with a positive weight.
use of jcal :
jcal is an override vector set in the PEP input stream
jcal(i) = 00 is the default
any digit of jcal that is = 0 is replaced by the corresponding
digit of JCT(2)
jcal (i) is a two digit number
each digit of jcal is independent of the other
the first (ones) digit of jcal(i) refers
to the calculation of corrections
the second (tens) digit of jcal(i) refers to the use of correction
ones digit =
1 : do not calculate the ith correction
2 : use standard logic for correction calculation
3 : calculate the correction
tens digit =
1 : do not use the ith correction
2 : use standard logic for correction use
3 : use the correction regardless of rank
read following variables from above labeled commons
from /INODTA/ and /NRMWGT/
from /PLNDTA/ and /OBSDTA/
from /FCNTRL/, /LCNTRL/ and /APRTBF/
from /PARAM/
from /BDCTRL/
from /DTPARM/
from /MASCON/
maximum number of planets allowed
integer control varible for propagation corrections
initialize remainder of fcntrl labeled common (except for
iterat which is set equal to zero in main program, nparam
which is calculated in subroutine check, and date which was
determined in main program)
HEDING =Input page heading describing program run.
DATE =Calendar date determined from operating system.
ITIMST =Real time at start of program run (binary integer number of
hundreths of seconds, determined from operating system).
NPARAM =Number of parameters to be adjusted in least squares analysis
(calculated in subroutine CHECK and again in ANSET).
ITERAT =Iteration number (initialized and incremented in main prog.)
NPAGE =Page number (initialized in subroutine input, incremented
after each new page heading is written).
Observed minus theory for measurement j deleted in normal equations
if abs(obs-th).ge.EPS(j)*error
j=1 time delay
j=2 Doppler
j=3 right ascension
j=4 declination
j=5 first limb of body crossing for transit or occultation
j=6 second limb of body crossing for transit or occultation
j=7 interferometric differential delay
j=8 interferometric differential delay rate
EPS(9) = Partial convergence criterion for adjustments to parameters.
EPS(10)= Full convergence criterion for adjustments to parameters
least squares iteration terminates if either ITERAT.ge.
ICT(1) or if all adjustments to parameters are less than
EPS(10) times corresponding standard deviations.
EPS(11)= Partial convergence criterion for cleanup of matrix inverse
in subroutine SNORML.
EPS(12)= Full convergence criterion for cleanup of matrix inverse
in subroutine SNORML.
ICT(1)=-99 Just print input parameters and stop.
ICT(1) =-1 No comparison of theory and observation, no least squares
analysis.
ICT(1) = 0 Comparison of theory and observation, no parameters
adjusted.
ICT(1) > 0 Maximum number of least squares iterations involving
integration,comparison of theory and observation,solution
of normal equations,adjustment of parameters.
ICT(2) =-1 Printout supressed for optical observed minus theory.
ICT(2) = 0 No printout supression for observed minus theory.
ICT(2) = 1 Printout supressed for radar and optical observed minus
theory.
ICT(3) =-1 Dummy observations are not written on output observation
library tapes and partials are not calculated.
ICT(3) = 0 Dummy observations are not written on output observation
library tapes but partials are calculated.
ICT(3) = 1 Dummy observations are so written with partials calculated
ICT(3) = 2 Same as 1, but also input observation library tapes with
ntape<0 (dummy data) are treated as real data, i.e.,
ntape -> |ntape|.
ICT(4) Governs the computing/copying of partials on obslib tapes.
Much of the logic is consolidated in routines pcops ('slow'
copying, called from everywhere) and lvt... (merging of L
vectors, called from cmpar1,2,3, shpcnt, trgcnt).
A 'quick' copy, if done, occurs at the beginning of PARTL
and preempts all the rest of the PARTL routines.
Note: some partials (such as satellite/probe parameters
and certain non-motion parameters) are always (re)computed
unless there is a quick copy.
In the following, 'M vector' refers to an L vector read in
from tape.
ICT(4) =-3 Use only saved partial derivatives and do quick copy.
Bug: output tape has the merged L vectors, not the copied
M vectors.
ICT(4) =-2 Use only saved partial derivatives if, for every L vector
set, corresponding M vector is also set. If true, quick
copy made, otherwise behaves like ICT(4) = 0.
ICT(4) =-1 Recompute only those partials that are specified in the
input L vector. Saved partials that are not in L vector
should be copied and not recomputed.
(Not implemented because PARTL has only the merged L vector
and not the original input L vector.)
ICT(4) = 0 Must used saved partial from tape even if the input L
vector for that partial is set. If the L vector is set but
the partial is not on tape, then compute the partial.
ICT(4) = 1 Partial derivatives of observations are computed, no saved
values are used, and the M vectors are ignored.
ICT(5) =-1 No initialization or forming of normal equations from
subset of saved normal equations.
ICT(5) = 0 Initialize normal equations from subset of saved normal
equations, process additional observations and solve
normal equations for ICT(1).gt.0 and iterat=1,2,3,...
ICT(5) = 1 For iterat=1, form normal equations from subset of saved
normal equations and solve them if ICT(1).gt.0 (for
iterat.gt.1, behaves like ICT(5)=-1 if ICT(1).gt.1)
ICT(5) = 2 For iterat=1, restore statistics, scale factors, solution
and inverse of normal equations from saved quantities on
data set jmat if ICT(1).gt.0 (for
iterat.gt.1, behaves like ICT(5)=-1 if ICT(1).gt.1)
ICT(5) = 3 Same as ICT(5)=2, iterat=1 except go to compar link
first to generate an observation library tape
ICT(5) = 4 Restore solution from dataset jmat plus zbar and
fbar-adjoint from partially prereduced normal equations
on imat0. Enables computation of adjusts and sigmas for
"reduced away" parameters, but not covariances; uncertainty
of predict cannot be calculated.
ICT(7) = 0 No overriding error weighting for dummy observation tape
(ntape=negative) in forming normal equations in NRMICT or
NRMFRM or predicting residuals in PRDICT.
ICT(7) = 1 Such overriding is allowed.
ICT( 8)= 0 Series sequencing checked in subroutine compar if
ICT(80).ge.0 (if ICT(80).lt.0, no sequence checking)
ICT( 8)= 1 Series sequencing not checked in subroutine COMPAR
for obs.library tape, but it is for cards IOBS and over-
riding observation series iobcon unless ICT(80).lt.0
ICT( 9)= 0 Perturbing planet data set for all integrations is ipert
ICT( 9)= 1 Perturbing planet data set for individual body integrations
is n-body data set ibody if nbody.gt.0 rather than ipert
(perturbing planet data set for n-body integration is still
ipert)
ICT(10)=-2 No prediction of observed minus theory or de-correlation
of partial derivatives
ICT(10)=-1 Such prediction for radio data only if ICT(11).gt.-2
ICT(10)= 0 Such prediction for radar and optical data if ICT(11).gt.-2
ICT(10)= 1 Such prediction for optical data only if ICT(11).gt.-2
ICT(11)=-3 No prediction, but de-correlation of partial derivatives
using f-bar-adjoint from imat3 as projection matrix (with
output tapes determined by the same rules as prediction)
if ICT(10).gt.-2
ICT(11)=-2 No prediction of observed minus theory or de-correlation,
even if ICT(10).gt.-2
ICT(11)=-1 Prediction of observed minus theory with printout,no tape
if ICT(10).gt.-2
ICT(11)= 0 Prediction of observed minus theory with printout and tape
if ICT(10).gt.-2
ICT(11)= 1 Prediction of observed minus theory with tape,no printout
if ICT(10).gt.-2
ICT(13)= 0 Normal sequencing of observation library data sets in
PRDICT link.
ICT(13).gt.0 & iterat=1 in PRDICT link, input observation library
data set iabs1 is IOBS0(j) for j=1,...,ICT(13) and
IOBS1(j) for j=ICT(13)+1,...,numobt
ICT(14)= 0 No uncertainty of prediction calculated in PRDICT link
ICT(14)= 1 Calculate uncertainty of prediction in PRDICT link
need ICT(10).ge.-1, imat.gt.0 if ICT(5).lt.2 or
jmat.gt.0 if ICT(5)=2
ICT(15)=-1 Nothing done in ANALIZ link even if ICT(1).gt.0.
Program is stopped after call to ANALIZ.
ICT(15)= 0 Everything done in ANALIZ link which other control integers
say to do.
ICT(15)= 1 Normal equations restored in ANALIZ link (if ICT(1).gt.0),
but there is no matrix inversion and solution.
Program is stopped after call to ANALIZ.
ICT(16)= 0 Number of tapes read in PRDICT link is numobt
ICT(16)> 0 Number of tapes read in PRDICT link is ICT(16) which could
be less than or greater than numobt (but not greater than
ten) there must be appropriate data set numbers in
IOBS1(i) or IOBS2(i) as the case may be.
ICT(17)= 0 Form both the left hand side (lhs) and the right hand
side (rhs) of the normal equations in subroutine NRMICT.
ICT(17)= 1 Form the rhs only in subroutine NRMICT.
ICT(17)=-1 Form the lhs only in subroutine NRMICT.
ICT(20) = 0 No printout of partials of observations in COMRIT.
ICT(20).gt.0 Printout of partials in COMRIT for first ICT(20) obs.
ICT(20).lt.0 Printout of partials in COMRIT for first iabs(ICT(20))
observations in each observation series.
ICT(24)=-1 Print out all elevation angles for radio observations if
the effect of the atmosphere or ionosphere is included.
ICT(24)= 0 Do not print out elevation angles for radio observations.
ICT(24)= n Print out all elevation angles .le. abs(n) degrees.
ICT(24)<-1 Print out temp,pressure and humidity at both sites.
ICT(24)=-3 Print out value for utrec in ATMION.
ICT(27)= 0 No center of mass of solar system used in radar and
interferometer calculations.
ICT(27)=-1 Position of Earth relative to center of mass of solar
system used in radar and interferometer calculations in
cislunar space.
ICT(27)=-2 Same as ICT(27)=-1 with velocity also used.
ICT(27)= 1 Same as ICT(27)=-1 plus position of center of mass of solar
system relative to sun used in planetary and mariner-type
space probe radar and interferometer calculations.
ICT(27)= 2 Same as ICT(27)=1 with velocity also used.
ICT(28)=-1 Nutation precesion evaluated midway between receive
and send time for radar and radio tracking data.
ICT(28)= 0 Nutation precesion evaluated midway for radar data,
both ends for radio data.
ICT(28)= 1 Nutation precesion evaluated at both send and receive
times for radar and radio tracking data.
ICT(29)=-1 Derivative of nutation precesion not used for
radar or radio Doppler.
ICT(29)= 0 Derivative of nutation precesion used for radio
Doppler but not radar Doppler.
ICT(29)= 1 Derivative of nutation precesion used for radar
and radio Doppler.
ICT(31)= 0 No checkpoint restart in COMPAR link
ICT(31)= 1,...,10 Checkpoint restart in COMPAR link with jtape=ICT(31)
ICT(31)=-1,...,-10 Checkpoint restart in COMPAR link with
jtape=iabs(ICT(31)).
Restart at start of tape.
(can only have checkpoint restart if ICT(80)=0)
ICT(33)= Value of nseq (series sequence number) at which compar
checkpoint restart occurs if ICT(33).gt.0.
If ICT(33).lt.0, compar checkpoint restart occurs at series
following observation series with nseq=-ICT(33).
ICT(34) Seasonal terms in Earth rotation.
ut2-ut1 before 1956, uses erotat con(14-21)
a1-ut1 after 1956, uses erotat con(14-21)
see JCT(29)
= 0 no terms included
> 0 number of terms included
initial epoch is ercon1(1), also used for polynomial
epoch. see JCT(34)
ICT(35)=-1 stop program if error on n-body,planet,embary,Moon data
set in compar link
ICT(35)= 0 skip optical observation, otherwise stop program if error
on one of above body data sets
ICT(35)= 1 skip all execpt probe observations, which stop program,
if error on body data set
ICT(35)= 2 skip all observations for which error on body data set
ICT(37)=-1 do not skip dummy observation which is below horizon
of observing site
ICT(37)= 0 skip such a dummy observation
ICT(37).gt.0 same as ICT(37)=0 plus increment time of observation by
ICT(37) minutes for first occurance
ICT(38)=-1 do not skip dummy observation which is is occulted by
central body of observed body
ICT(38)= 0 skip such a dummy observation
ICT(38).gt.0 same as ICT(38)=0 plus increment time of observation by
ICT(38) minutes for first occurance
ICT(39)=-2 take nothing from n-body tape, just there for center of
mass of solar system
ICT(39)=-1 take only Moon from n-body tape if there is one
ICT(39)= 0 take everything from n-body tape that is there (compar l.)
ICT(39)= 1 use any individual body tapes which are there
(note: must set iem=0 to override default if no embary tape
ICT(42) = 1 filter input in input stream
ICT(42) = 0 no Kalman-Bucy filter (just maximum liklihood,
least-squares paramter estimation)
ICT(42) > 1 filter input on data set no. ICT(42)
filter input described in subroutine FILTIN
note: old code was removed in 1983 march
ICT(43) = 0 use the old code in the radar and interferometry links
(DELDOP, SBDLDP, MDELDP; INTERF, MNTERF, etc.)
ICT(43) = 1 use the new code in the radar and interferometry links
(RADCTL, RADSB, RADLND, RADMN; FERCTL, FERMN, etc.)
ICT(44) = 1 apriori covariance matrices are in input stream
ICT(44) =-1 apriori covariance matrices are in input stream, with no
right-hand-side contribution, and the right-hand side is
not to be adjusted after each iteration
ICT(44) = 0 no apriori parameter covariance matrix
ICT(44) > 1 apriori covar. matrices are on data set no. ICT(44)
ICT(46)= 0 normal equations solved and coefficient matrix inverted
with double precision arithmetic using fact that it is
symmetric in subroutine SNORML
ICT(46)= 1 same as ICT(46)=0 but extended precision used in
intermediate calculations
ICT(46)= 2 same as ICT(46)=0 but extended precision used for
matrix inversion, and the results are then returned
as double precision
ICT(46)>10 solve subset of the normal equations by eliminating the
ICT(46)-10 smallest eigenvalues and inverting the reduced
matrix that results
ICT(48)=-2 no printout of covariance or correlation matrices
ICT(48)=-1 no printout of inverse of coefficient matrix of normal
equations (covariance matrix)
ICT(48)= 0 printout of inverse of coefficient matrix of normal eqs.
ICT(49)= maximum number of iterations for cleanup of matrix inverse
in subroutine SNORML (if 0, no cleanup)
ICT(49)=-1 no iterative cleanup nor calculation of error matrix
(product of calcuated inverse and orignal matrix minus the
identity matrix)
ICT(66) binary coded control for mars rotation model and code
1 bit=0: turns off new code. use old code for mars rot.
=1: use new code
2 bit=1: print conversion to proper time
4 bit=1: print out i, psi, phi as calculated in spotpl
8 bit=1: print out spot position vector in mars fixed and
1978.0 frame (mars centered)
16 bit=1: print out spot logic vector after call to rotlog.
controls drop through sequence that computes partials
of rv w.r.t. various parameters
32 bit=1: print out spot position vector in 1978.0 and reference
frame (mars centered)
64 bit=1: print out partials of rv w.r.t. various parameters
128 bit=1: print out final partial of position vector in frame
of reference epoch w.r.t. parameters
256 bit=1: print out full rotation matrix
ICT(67): control for "fixing" range observable based on the
assumption that a bit error occurred in the high order
bits that make up the range observable
blim=2**(-abs(ICT(67)))
blim is the maximum allowed range residual (seconds)
blim is also the minimum correction to range allowed
if(ICT(67) is less than zero, then the predicted residual
is used in the computations
this is all done in subroutine rngmod
ICT(76) control for disparate parameter sne combination
= 0 no DPSNEC operations
= 2 flag disparities for adjustable nominals
= 4 flag disparities for all nominals
= 6 copy nominals from imat0(1)
=10+n same as n, but also correct rhs for differences
ICT(77)= 0 no punch output of input data at end of program run
ICT(77)= 1 punch output of input data at end of program run so that
job can restarted at termination point with last adjusted
values of parameters if ICT(1).gt.0
ICT(77).gt.1 write punch output to data set no. ICT(77)
after every iteration plus do same as ICT(77)=1
ICT(78)= 0 no tapes are saved at end of least squares iteration
ICT(78)=positive at end of least squares iterations 1 to ICT(78),
inclusive, there is a pause with typewriter message
'remove and mount tapes as specified' for operator to
follow special instructions (except for last iteration)
ICT(79)=-2 dummy observations not calculated each iteration nor
after convergence
ICT(79)=-1 dummy observations not calculated each iteration, but are
after convergence
ICT(79)= 0 dummy observations calculated each iteration, but not
after convergence
ICT(79)= 1 dummy observations calculated each iteration and after
convergence
ICT(80)=-1 input observation library tapes not used
ICT(80)= 0 input observation library tapes are used, program eats its
tail from iteration to iteration
ICT(80)= 1 input observation library tapes only are used to form the
normal equations for first iteration with no computation
of observed minus theory.
subsequent iterations same as ICT(80)=0
itrwnd(i) = 0 or 1 according to whether data set i is in a rewound
state or not, i=1,ntptop
actually, any positive number indicates non-rewound.
in some cases, values other than 0 or 1 may have
significance as to the contents of the data set.
note that if itrwnd in /fcntrl/ labeled common is given a length
longer than 99, then subroutines xload (referring to extprc) and
oprmsg (referring to typout) must have dimension of itrwnd changed
besides all fortran routines that use the JCT vector.
ntptop should be increased to the new length of itrwnd
JCT(2) = jcal meta default-
any digit of jcal that is = 0 is replaced by the corresponding
digit of JCT(2)
JCT(2) = 11 is the default
if JCT(2).lt.0 propco is shut off
JCT(3): binary coded printout control for propco
1 bit=1: print sumcor
2 bit=1: print uscor vector
4 bit=1: print calcor vector
8 bit=1: not used (formerly, print elevation angle and/or rate
see ICT(24) )
16 bit=1: print each cal as it is added to sumcor
32 bit=1: print cal formation for phase delay dop.
JCT(6): binary coded flags for compar link debug
1 bit: print delays in routine tsav
2 bit: print vectors & magnitudes in uvectr
4 bit: print positions in interpolators
8 bit: print velocities in interpolators
64 bit: print initial distance guess in obsred
128 bit: print nutation, precession and wobble angles
256 bit: print vector partials in cpartl
512 bit: print nutation,precession matrices in prcnut,preces
1024 bit: print time argument passed to integration reader
2048 bit: print relativistic time delay correction
4096 bit: print eccentric telescope mount correction
8192 bit: print CTAT and A1-UT1 corrections
JCT(10): binary coded control for solid Earth tides
1 bit: calculate effect of lunar tides on radio and laser
ranging observables (coded only in radmn, fermn,
and ferstr)
2 bit: calculate the effect of solar tides on these
observables
4 bit: print out love numbers and lag angle used
8 bit: print out radial,longitudinal, and latitudinal
components of tides at each site
64 bit: use IERS tide model
JCT(11): binary coded control for solid body tides on the Moon
1 bit: calculate the effect of Earth-induced tides on radio
and laser ranging observables in radmn and fermn
2 bit: calculate the effect of sun-induced tides on the Moon
4 bit: print out love numbers and lag angle used
8 bit: print out radial, longitudinal and latitudinal
components of tides at each site
JCT(12): controls meaning of dt parameters
-1: dt(i) is constant value of a1-ut1 or wobble
over interval jddt(i) to jddt(i+1)
0: dt(i) is value of a1-ut1 or wobble at jddt(i),
tabular point in piecewise linear function
JCT(14): controls use of cross partials on embary and observed body
tapes (i.e. 406 on embary tape or 306 on mars tape)
does not affect target body partials of non-observed bodies
<=0: don't use (default)
>0: use
JCT(15): controls underflow interrupt during inversion of d-matrix
in partial prereduction
0: underflows generate interrupts which are intercepted by
the fortran error handler. result is set to zero and
traceback is printed for first 15 occurrences.
1: underflow interrupt is disabled. result is set to zero
by hardware. there is no indication of any error, nor
is there any overhead to process the error.
JCT(16): controls time tag used for planetary rotation (bodies other
than Earth)
0: (default) use proper time (omly implemented for mars new
rotation model)
1: use ephemeris time
JCT(17): controls Lorentz contraction calculation for observing and
observed bodies
0: (default) apply contraction
1: omit contraction
JCT(20): binary coded control determines whether PEP stops or
continues after certain errors. 0 --> continue,
1 --> abort.
1 bit: reserved for arithmetic exception in input link
2 bit: errors in multiple parameter solution control
JCT(21): controls precession model (also mean sidereal time)
-1 use pre-1976 IAU model (default if nothing specified)
0 use IAU 1976 (default if non-old precession constant)
1 use IAU 2000 (not implemented)
2 use IAU 2006
JCT(24)=-1 print out all planetocentric elevation angles for radio
observations of a probe within the atmosphere of another
planet
JCT(24)= 0 do not print out planetocentric elevation angles
JCT(24)= n print out elevation angles .le. n degrees
JCT(25) =-2 do not calculate transmitter freq. drift rate corrections
in alsep differenced n-count observable
JCT(25) =-1 do not include transmitter freq. drift corrections in
alsep differenced n-count observable
JCT(25) = 0 include such correction
JCT(26) = 0 if ilib.eq.0 get Moon physical libration from Moon or
nbody tape
JCT(26) = 1 if ilib.eq.0 get Moon physical libration by evaluating
analytical theory
JCT(27) = 0 usual ordering required in input stream
= 1 free ordering of input stream under control of * commands
with the following constraints:
a. title card
b. &nmlst1
c. any *object packets
d. any other * commands (note: *apriori must follow any
input giving parameters to be adjusted.)
note: with * commands the usual &end or blank card
is not needed to end a given portion of the input stream.
JCT(28): binary coded control for corrections to Earth rotation.
1 bit=1, include corrections to nutation to account for
the Earth's elasticity and fluid core.
(ercond(13-16) in prcnut)
2 bit=1, include corrections to a1-ut1 to restore short-
period terms removed by b.i.h. smoothing.
4 bit=1, include corrections to diurnal polar motion.
8 bit=1, include free nutation corrections.
(ercond(11-12) in prcnut)
JCT(30) = 0 no e,f,g (+edot,fdot,gdot) tape output
Earth centered rotating coordinates. e in true equator of
date towards greenwich meridian, g to north, f completes
right hand system
JCT(30).ne.0 such tape output on data set iabs(JCT(30))
unless JCT(40) says there is angle output instead
.lt.0 no frintout of e,f,g
.gt.0 printout of e,f,g at end of dummy observation series
JCT(31) = space defense center satellite number for e,f,g output
JCT(32) = 0 no binary e,f,g output in subroutine efgout for fitting
polynomials in subroutine efgfit
JCT(32).ne.0 such binary output on data set iabs(JCT(32)) if
JCT(30).ne.0
.gt.0 JCT(30) output with JCT(32) output
.lt.0 JCT(30) output is supressed
JCT(33) = 0 use BIH UT1 and wobble if available (before 1969 Jan 3
must use USNO values; before 1956 Jan 17, must use IPMS
wobble and cannot calculate A.1-UT1)
=-1 use USNO values of A.1-UT1 and wobble (available from
1956 Jan 3 to 1975 Feb 26 for A.1-UT1, to 1975 May 7 for
wobble)
.gt.0 read A.1-UT1 and wobble tables from data sets
JCT(33) and JCT(33)+1, respectively.
JCT(34) = 0 no polynomial correction to CT-UT.
= 1 such a polynomial correction, before 1956. (coded for
mercury transit observations only)
coefficients are ercond(17-19),[con(11-13)].
Epoch is hardwired as 1956 Jan 17.0 (see CTUTF).
Note: can logically be used to control the polynomial
correction to A.1-UT1 after 1956(also using ERCOND(17-19)).
However, the current code does not test JCT(34), but
rather CON1(1).gt.0 (see A1UT1F). See ICT(34).
JCT(35) = hour of start of dummy observations each day
JCT(36) = minute of start
JCT(37) = hour of end of dummy observations each day
JCT(38) = minute of end
If JCT(35 to 38) are zero, every time within day is allowed.
Else, dummy observations done just between JCT(35&36) and JCT(37&38).
If JCT(37&38).lt.JCT(35&36), then period of obs includes 0 hr UTC.
These time limits could be the start and end of an optical observatory
observing evening. This feature is not needed for radar predictions.
If JCT(35 to 38) nonzero, best to have ICT(37) & ICT(38) .le.0.
JCT(39)=-1 right ascension and declination rates instead of angles
themselves are calculated in subroutines of OPTIC
(for use in calculating covariance of rate predictions
because the partial derivatives of the rates are coded).
JCT(39)= 0 photographic topcentric right ascension-declination
observations are referred to the mean equinox and equator
of the reference epoch
JCT(39)= 1 photographic topocentric right ascension-declination
observations are referred to the true equinox and equator
of date
JCT(39)= 2 same as JCT(39)=1 with hour angle replacing right asc.
(also applies to transit-circle observations)
JCT(40) = 0 no angle output for optical observatories
JCT(40).ne.0 such tape output on data set iabs(JCT(30))
in format compatable with data general nova computer
JCT(40).gt.0 printout of angles, data set JCT(40) used as intermediate
ebcdic buffer
JCT(40).le.0 no printout of angles
JCT(41) = 0 angle output referred to true equinox and equator of date
for JCT(40) output
JCT(41).gt.0 angle output referred to mean equinox and equator of
1900+JCT(41). for eample, JCT(41)=76 denotes output
referred to mean equinox and equator of 1976.0
JCT(42) = 0 no refraction correction made in computing JCT(40) angles
JCT(42) = 1 refraction correction made in computing JCT(40) angles
(optical frequency correction)
JCT(42) = 2 refraction correction made in computing JCT(40) angles
(radio frequency correction, millstone sub.dell)
JCT(45) = 0 no prediction print and tape output for Arecibo and other
radar sites observing asteroids and Jupiter satellites
and also planets. See also JCT(70)
JCT(45) < 0 such print output (az,el,ra,decl,delay,Doppler in ANGOT)
JCT(45) > 0 such print and arecibo bcd site tape from ANGOT. Output
file is JCT(45) and, if >50, also work file on JCT(45)+10.
JCT(46) = 0 JCT(40) print angle output has range in meters and
range rate in cm/min
JCT(46) = 1 JCT(40) print angle output has one way time delay in
seconds and one way Doppler shift in Hz
JCT(47) = 0 usual JCT(40) print and tape output
JCT(47) > 0 JCT(40) tape output has satellite radius instead
of Doppler rate, print output has special les-8/9
format with 6 frequency Doppler
JCT(47) = 8 les-8 output
JCT(47) = 9 les-9 output
JCT(48) = 0 right ascension and hour angle rates are seconds
of time per second in JCT(40) print output
JCT(48) = 1 right ascension and hour angle rates are seconds
of arc per second in JCT(40) print output
JCT(49) = control for fluid displacement affecting sites + unit number
of external dataset of displacements
control has packed bits stored as control*100, independently
selecting effects to include in the calculated displacement
1: ocean loading
2: atmospheric loading
4: groundwater loading
also a print flag stored as print*10000 to control whether
the results are to be printed out as calculated
JCT(51)< 0 restore normal equations after solution and use disparate
parameter facility to compute postfit rms for each set of
normal equations
JCT(51)= 0 no multiple parameter set solutions
JCT(51)> 0 solve for multiple parameter sets
JCT(52)= 0 no summary printout for multiple parameter set (or
multiple epoch ) runs
JCT(52)> 0 summary printout for multiple parameter set (or
multiple epoch ) runs
JCT(53)= 0 no partial prereduction of normal equations
JCT(53)= 1 partially prereduce normal equations using ne generated
during this run
JCT(53)= 2 partially prereduce normal equations using ne generated
during an earlier run only (total sne)
JCT(55) is a packed-bits control for the planet/system center of mass
correction performed in plcms.
JCT(55) = 0 coordinates computed in pltrp are purely from tabulated
values for the planet itself.
1 bit = 1: get satellite coordinates from interpolation if possible
2 bit = 1: get satellite coordinates from elliptic if necessary
4 bit = 1: use only satellites that are on s-body tape
8 bit = 1: for elliptic method, obtain elements from s-body tape
(if possible). warning: this option should be set as a
rule because the initial conditions in /empcnd/ are
likely to have been overwritten from tape anyway and
may not match the epoch in con1(1).
16 bit = 1: allow calculation of barycenter of planet+observed
satellite instead of planet (not done unless partials
are being calculated)
32 bit = 1: use precessing elliptic orbit instead of stationary
JCT(56)= 0 NUMINT calls SBOUT, not SBFOUT.
JCT(56)= 1 NUMINT calls SBFOUT, not SBOUT.
JCT(56)> 1 NUMINT calls SBFOUT, not SBOUT; SBFOUT restarts filter
by reading ICONOF with ITERAT = 1.
JCT(59)= 0 'equator-equinox' terms applied to interferometric
observations are used as a pseudo-clock polynomial
affecting the observable (just as in radio observations)
JCT(59)> 0 'equator-equinox' terms are used as true clock corrections
(affecting ct epoch + phase and frequency) for
interferometric observations only
JCT(60)=-1 like 0, but print diagnostic information
= 0 copy a priori matrix + mean residual sensitivity vector
to sne datasets and "d" variances to pprsne
= 1 do not write a priori or mean residual info (for downward
compatibility). this can allow restoring sne datasets
with imbedded a priori matrices without having assigned
imat2 (q.v.)
JCT(67)= 0 use true velociies in compar link for observations
= 1 use apparent velocities in compar link (logic contained
in subroutine vlrtrd). also apply light time iteration
corrections to delay observable partials in compar link
(forces computation of velocities).
JCT(68) controls the logical time direction for observation series
except dummy observations (which go from start time to
finish time)
JCT(68)=-1 logical time direction is backward
JCT(68)= 1 logical time direction is forward
JCT(68)= 0 logical time direction obeys usual criteria, i.e.,
forward if nplnt0>30 or ncodf<4 or ncodf>9 or nspot>0
or nspot2>0 or non-Earth observing site
otherwise, backward.
JCT(69)= 0 input observation header cards follow old format(see
subroutine dltred for details)
JCT(69)= 1 input observation header cards follow new format(see
subroutine dltred for details)
JCT(70) packed bits controls for ANGOT
1 abbreviated print to JOUT in ANGOT if JCT(45).ne.0 (else verbose)
2 use mean equator/equinox of reference epoch (else true of date)
4 increase reference epoch by 50 yrs (else use default)
JCT(71)= 0 normal choice of source ct-at values
JCT(71)> 0 radar link uses JCT(71) as selector for ct-at values
in ctatf calling sequence (variable nterm)
JCT(71)= 2 radar link uses analytic annual ct-at term
JCT(71)= 4 radar link uses integrated annual ct-at term
JCT(78)= 0 Unnormalized gravitational potential tesseral harmonic
partial derivatives in SBFN (also zonals)
JCT(78)= 1 Normalized gravitational potential tesseral harmonic
partial derivatives in SBFN (zonals still unnormalized)
JCT(79)= 0 no reintegrations before calculating dummy observations
after orbit fit convergence
JCT(79)= 1 reintegration of motion without partials after orbit fit
convergence to be used in calculation of dummy observations
if iabs(ICT(79))=1
lprm(i)=k parameter PRMTER(k) to be adjusted (i=1,100)(k=1,100). if
k=0 no parameter adjusted. we must have
(1) if j.gt.i and lprm(i)=0, then lprm(j)=0
(2) if j.gt.i and lprm(i),lprm(j) not 0, then
lprm(j).gt.lprm(i)
lem(i) =least-squares adjustment control constants for econd,i=1,30
lmn(i) =least-squares adjustment control constants for mcond,i=1,30
ler(i) =least-squares adjustment control constants for ercond,i=1,30
lmr(i) =least-squares adjustment control constants for mrcond,i=1,30
lpl(i,j)=least-squares adj.control const for pcond(.,j),i=1,30,j=1,
numpln
see subroutine bodred for explanation of lem,lmn,ler,lmr,lpl where
these quantities are read in &nmlst2 namelist as l.
MASS(j) =(mass of body j)/(mass of sun)
1 Mercury
2 Venus
3 Earth+Moon
4 Mars
5 Jupiter
6 Saturn
7 Uranus
8 Neptune
9 Pluto
MASS(10)=(mass of Moon)/(mass of Earth+Moon)
MASS(j),j=11,30 masses of additional bodies (divided by mass of sun
if no nplnt(i)=j (i=1,u_mxpl). if nplnt(i)=j for some i, then
mass is divided by the total mass of the planet-satellite
system npcent(i) or mass of sun if npcent(i)=0.)
By definition, all MASS values are less than 1, and an input value
greater than 1 will be replaced by its inverse, except that values
between 1 and 500 are allowed for body 9, which may thus be
"hijacked" from the solar system and treated as a distant star.
PRMTER(31)=RELFCT=general relativity motion factor
PRMTER(32)=GMVARY=time variation factor for gravitational constant
(unit=1/day)
PRMTER(33)=SUNHAR=second harmonic of gravitational potential of sun
(unit=(sun radius)**2
PRMTER(34-38)=densities (g/cc) of five classes of asteroids
PRMTER(39)=SUNPOE=ratio of gravitational to total energy of the sun
(negative, used for beta, gamma contribution to
principle of equivalence violation)
PRMTER(40)= principle of equivalence parameter deltas (see SBSET).
PRMTER(41)=metric parameter beta, excluding principle of equivalence
violation contribution
PRMTER(42)=metric parameter gamma, excluding p.o.e. violation
PRMTER(43)=metric parameter beta, including principle of equivalence
violation contribution
PRMTER(44)=metric parameter gamma, including p.o.e. violation
PRMTER(45)=mass of 2nd extra asteroid/comet belt (solar masses)
PRMTER(46)=mass of 1st extra asteroid/comet belt (solar masses)
PRMTER(47)=right ascension of the ascending node of the asteroid belt
on the mean equator of the reference epoch measured from
the mean equinox (degrees)
PRMTER(48)=inclination of the asteroid belt to the mean equator of
the reference epoch (degrees)
PRMTER(49)=distance of circular asteroid belt from sun (au)
PRMTER(50)=mass of asteroid belt (solar masses)
PRMTER(51)=AULTSC=value of the astronomical unit in light seconds
(inverse of the speed of light at a large distance
from the sun)
PRMTER(52)=LTVARY=time variation factor for speed of light
(unit=1/day)
PRMTER(53)=RELDEL=general relativity time delay factor
PRMTER(54)=RELDOP=general relativity Doppler shift factor (not
programmed) just used as a switch. if greater than
zero, subroutine radar sets ndop for general rel.
effect on Doppler with factor reldel like time delay
PRMTER(60)= constant plasma factor for interplanetary media
PRMTER(61)= time variable plasma factor for interplanetary media
PRMTER(62)= Earth atmosphere effect factor for delay and Doppler
PRMTER(63)= Earth ionosphere effect factor for delay and Doppler
PRMTER(72)=CTVARY=coefficient of scale variation for time varying
gravitatonal constant
PRMTER(73)=radius of 1st extra asteroid belt (au). see PRMTER(47-49)
PRMTER(74)=inclination of 1st extra asteroid belt (deg). ditto
PRMTER(75)=ascending node of 1st extra asteroid belt (deg). ditto
PRMTER(76)=radius of 2nd extra asteroid belt (au). ditto
PRMTER(77)=inclination of 2nd extra asteroid belt (deg). ditto
PRMTER(78)=ascending node of 2nd extra asteroid belt (deg). ditto
PRMTER(81)= scale factor for atomic time to coordinate time conversion
PRMTER(82)= phase for atomic time to coordinate time conversion
ECINC = mean inclination of ecliptic on equator at epoch (degrees)
SEQINC = inclination of suns equator on ecliptic (degrees)
SEQASC = longitude of ascending node of suns equator on ecliptic
measured from the mean equinox at reference epoch (degrees)
SUNRAD = radius of sun in kilometers
PRMTER(95)= radius of sun in kilometers for transit observations
PRMTER(97)=2440000.5=epoch for time varying gravitational constant
MDSTAU = distance unit for lunar ephemeris on perturbing planet tape
in astronomical units (if 0, assumed to be 1)
MDSTSC = distance unit for lunar ephemeris on Moon tape in light-
seconds (if 0, assumed to be aultsc)
LTVEL = velocity of light in kilometers per second
nbody = 0 no n-body integration or n-body data set
nbody = n.gt.0, n-body integration and/or n-body data set with n
bodies plus Moon if not one of n=9 bodies
nplbdy(j),j=1,nbody, are the planet numbers of the bodies in the
n-body integration in the order which they are integrated
the initial conditions for the integration come from
/empcnd/ labeled common
nplbdy(j)= 3 body is Earth-Moon barycenter relative to the sun with
initial osculating elliptic orbital elements econd(1-6)
nplbdy(j)=10 body is Moon relative to Earth with initial osculating
elliptic orbital elements mcond(1-6)
nplbdy(j)=nplnt(k) for some k=1,...u_mxpl, body is planet nplnt(k)
relative to central body ncentr(k) (the sun if ncentr(k)
is 0) with initial osculating elliptic orbital
elements pcond(1-6,k) (nplnt and ncentr are in /namtim/
labeled common)
jdbdy0=initial epoch of n-body integration is midnight beginning of
day of day with julian day number jdbdy0 (if negative there is
checkpoint restart at the epoch -jdbdy0 unless
jdbdy0.le.-3000000. if jdbdy0=-1 there is checkpoint restart
at record just before end of file of n-body ephemeris tape.)
(if jdbdy0=0 or jdbdy0=-3000000, there is no
numerical integration, comparison of theory and observation
assumes motion is already on tape, but if initial conditions
are adjusted, jdbdy0 is set to positive value on second record
of tape so there will be reintegration on next least squares
iteration, whereas if no initial conditions are adjusted,
jdbdy0 is set to zero in the compar link so there will be no
reintegration on subsequent least squares iterations)
(jdbdy0.le.-3000000 is signal to compar link on first
iteration to override end time on tape by input jdbdy2)
jdbdy1=julian day number of first data record on output tape
jdbdy2 julian day number of last data record on output tape
jdpad if positive is used to extend the range of jdpby1-jdbdy2 at both
ends (provided those are also both positive). in other words,
if jdbdy2>jdbdy1, then jdbdy2 is replaced by jdbdy2+jdpad, and
jdbdy1 is replaced by jdbdy1-jdpad. if jdbdy2<jdbyd1, then
jdbdy2 is replaced by jdbdy2-jdpad, and jdbdy1 by jdbdy1+jdpad.
this same logic applies to all jd1/jd2 pairs in the same run.
must have jdbdy0 lying between jdbdy1 and jdbdy2, unless jdbdy0 is
negative indicating checkpoint restart, in which case -jdbdy0
must be between jdbdy1 and jdbdy2, or unless jdbdy0 is zero
indicating no integration or less than or equal to -3000000
indicating no integration and input jdbdy2 to everride that on
second record of input n-body data set in compar link (needed
because such data set might have been produced by check point
restart going beyound point indicated on second record)
(if jdbdy0=-1, consistency check for -jdbdy0 between
jdbdy1,jdbdy2 is not made)
time direction of output tape is from jdbdy1 to jdbdy2, forward in
time if jdbdy1.lt.jdbdy2, backward in time if jdbdy1.gt.jdbdy2
if jdbdy0.lt.0 & .gt.-3000000, output tape from a previous integration
is read to the time -jdbdy0 (or to the record just before end
of file if jdbdy0=-1) and the coordinates on the tape at
that epoch are used to restart that integration, proceeding
from -jdbdy0 to jdbdy2. integration cannot go in both
directions from epoch -jdbdy0 in restart mode and a
previous integration cannot be restarted at a time in the
first direction of the previous integration if it was in the
two direction mode, since output tape is not written until it
gets to second direction. in the checkpoint restart mode with
jdbdy0 negative, the only data which can be input for the
integration are jdbdy1,jdbdy0,jdbdy2,KBDY(39) because the rest
of the integration data is taken from the second record of the
tape of the previous integration. if KBDY(39).lt.0, the tape
of the previous integration is rewound , whereas if
KBDY(39).ge.0, the restarted integration continues writing
the tape of the previous integration.
if jdbdy0 is positive between jdbdy1 and jdbdy2, we have
(a) if jdbdy0=jdbdy1, numerical integration starts at jdbdy0
and proceeds to jdbdy2, writing output tape as integration
proceeds if KBDY(39).ge.0
(b) if jdbdy0 is strictly between jdbdy1 and jdbdy2, integration
starts at jdbdy0 and goes in direction of jdbdy1 (forward
in time if jdbdy0.lt.jdbdy1, backward in time if jdbdy0.gt.
jdbdy1). the integration polynomial coefficients determined
at jdbdy0 by the starting proceedure are saved and the
numerical integration proceeds towards jdbdy1, writing onto
disk buffer ibuf if KBDY(39).ge.0. when jdbdy1 is reached,
ibuf is backspaced and read to write output tape ibody in
time direction from jdbdy1 to jdbdy0 if KBDY(39).ge.0. then
using the saved starting polynomial coefficients, the
integration goes from jdbdy0 to jdbdy2, writing output
tape ibody as the integration proceeds if KBDY(39).ge.0.
when the integration is completed, we have an output tape
from jdbdy1 to jdbdy2, even though integration went from
jdbdy0 to jdbdy1 and then from jdbdy0 to jdbdy2. these
contortions are gone through in order to simplify the logic
needed to read the tape in the compar link.
JVLBDY = 0 velocity as well as position are on output tape
JVLBDY = m (m.gt.0) velocity as well as position are on output tape
forward from the initial epoch and backward from the epoch
for m days, but only position data is on tape backward
from initial epoch more than m days
INTBDY =-m, m.ge.0, tabular interval is 2**-m days for n-body integr.
INTBDY = m, m.gt.0, tabular interval is m days for n-body integration
This variable is arbitrarily defined to refer to Mercury.
There are 5 tabular points per record on output tape for
Venus through Pluto, but 10 points for Mercury and 40 for
the Moon. Any value other than 2 will cause problems.
If the Earth-Moon barycenter is integrated in the n-body routine, then
kbdy(j) is examined for all j=1,nbody. If, for any j, we have
KBDY(j)>= 0, then the exact equations of motion are integrated,
with the coordinates of the Moon relative to the Earth
being determined from the integration if the Moon is
one of the integrated bodies and from the perturbing
planet data set if it is not. Otherwise, if, for all j,
KBDY(j)< 0, then the Earth-Moon barycenter is treated as a point particle
unless the Moon is also integrated.
If the Moon is integrated here, then all kbdy(j) are treated as 1.
Further, all planets in the integration are forced to be included
as perturbers of the lunar orbit, superseding the values specified
in K(31-39) for the Moon. (see BODRED)
KBDY(19) =-1 Cowell's method used for the Moon if it is one of the
integrated bodies (default)
KBDY(19) = 0 Encke's method with elliptic reference orbit is used for
the Moon if it is one of the integrated bodies
KBDY(19) = 1 mean orbit method is used for the Moon if it is one of
the integrated bodies
KBDY(21) =-1 no relativity effect in n-body integration
KBDY(21) = 0 effect of general relativity included in n-body integr.
(with same relativity motion factor relfct used for
all bodies)
KBDY(21) = 1 same as KBDY(21)=0 except that each body has its
individual relativity motion factor
KBDY(22) =-1 no effect of time variation of gravitational constant
in n-body integration
KBDY(22) =>0 effect of time variation of gravitational constant
included in n-body integration
KBDY(23) =-1 no effect of second harmonic of gravitational potential
of sun included in n-body integration
KBDY(23) =>0 effect of second harmonic of gravitational potential
of sun included in n-body integration
KBDY(25)=-1 effect of any input limited asteroids not included in
n-body integration
KBDY(25)=>0 effect of all input limited asteroids (if any) included in
n-body integration
KBDY(27) = positive, step size for Adams-Moulton or second sum
numerical integration is KBDY(27) days
KBDY(27).le.zero, interval size for Adams-Moulton or second sum
numerical integration is 2**KBDY(27) days
KBDY(30) = number of equations controlling step size for n-body
integration (must equal the number of equations)
interval size controlled throughout nordsieck method, but
only during nordsieck starting proceedure for adams-
moulton method. This input variable is currently ignored,
and the actual number of equations is used instead.
KBDY(38) = 0 printout every tabular point if KBDY(39).le.0 for n-body
integration
KBDY(38) = n.gt.0, printout every n tabular points if KBDY(39).le.0
for n-body integration
KKBDY(71) = 0 do not print accelerations
KKBDY(71) = 1 print accelerations every tabular point on Kout, if any
KKBDY(71) = 64 print high-precision coordinates on Kout, if any
KKBDY(71) = 192 same as 64, but suppress printing of libration angles
jddt0 =-1 dt table refers to wobble only
jddt0 = 0 dt table refers to A.1-UT1 and wobble (used after 1956)
jddt0 = 1 dt table refers to A.1-UT1 (used after 1956)
jddt0 > 1 dt table refers to ET-UT2 (used before 1956)
numdt = number of tabular points in each segment (A.1-UT1, xwob, and
ywob, or just A.1-UT1 or just ET-UT2) of adjustable dt table
if A.1-UT1, xwob, and ywob are respresented by the table,
the maximum number of each is 200. if just A.1-UT1 or
ET-UT2 are represented, the maximum number is 600.
dt0 = value of A.1-UT2 at Julian date jddt0-0.5
dt(i) = value of A.1-UT1 or ET-UT2 at Julian date jddt(i)
dt(i+200)=value of xwob(in arcsecs) at Julian date jddt(i), i=1,numdt
dt(i+400)=value of ywob(in arcsecs) at Julian date jddt(i), i=1,numdt
ldt(i) = 0 dt(i) not adjusted in least squares analysis, i=1,numdt
ldt(i) = 1 dt(i) adjusted in least squares analysis, i=1,numdt
see also JCT(12)
iseq =0 check tape and series sequencing numbers for increasing
order for overriding error weighting and dummy
observtion cards in subroutine prnobs and stop program
if there is a sequencing error
iseq =1 do not stop program if there is such a sequencing error
jcal controls the calculation and use of propagation media
corrections in the radar and interferometry links
jcal (i) is a two digit number
each digit of jcal is independent of the other
the first (ones) digit of jcal(i) refers
to the calculation of corrections
the second (tens) digit of jcal(i) refers to the use of correction
ones digit =
1 : do not calculate the ith correction
2 : use standard logic for correction calculation
3 : calculate the correction
tens digit =
1 : do not use the ith correction
2 : use standard logic for correction use
3 : use the correction regardless of rank
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
globdefs.inc
inodta.inc
plndta.inc
obsdta.inc
aprtbf.inc
bdctrl.inc
dtparm.inc
fcntrl.inc
funcon.inc
lcntrl.inc
nrmwgt.inc
mascon.inc
maxmt0dt.inc
param.inc
prpgat.inc
loadid.inc
timstf.inc
in0
sepeat
moprnt
iseq1
nstop
jdpad
- ...
- ...