jd    =input julian day number
ctut1 =input offset of CT from UT1 (s)
sidtm0=output mean sidereal time at 0 hr uninversal time UT1
=mean sidereal time at julian date jd-.5 (in radians)
fract =output ratio between mean sidereal time and universal
time multiplied by twopi divided by 86,400
gmstera=difference between Greenwich mean sidereal time and
Earth rotation angle at epoch (in radians)
(only if 2000 or later IAU model selected)

logic added mar 1980 to select old or new iau definition of
sidereal time.  the new definition was recommended by iau
commissions 4, 19, and 31 in aug 1980 so that there will be
no change in the rate of ut1 when the iau(1976) precession
constant if used, or a change in the value of ut1 when the fk5
equinox is used.

use IAU 1980 expression for sidtm
changed by r.king dec 1980 to final iau 1982 values
Corresponding expressions in seconds of time and julian centuries
relative to J2000 are:
excess fract= 8640184.81284857 +0.186209355 t -1.864121e-05 t^2
sidtm0=
24110.5483974 +8640184.81285379 t +0.09310441 t^2 -6.197924e-06 t^3

use old expression for sidtm
corresponding expressions in seconds of time and julian centuries
relative to J2000 are:
excess fract= 8640184.72800493 + 0.185874259 t
sidtm0= 24110.4708965 + 8640184.72779644 t + 0.09290001 t^2

use IAU 2006 model
This procedure adopts the new formula for true sidereal time (GST) as
Earth rotation angle (ERA) minus Equation of origins (EO) instead of
the old mean sidereal time (GMST) plus projected nutation in longitude.
Instead of returning the rate of GMST, it now returns the rate of ERA.
However, to preserve compatibility, the operation still calculates and
returns GMST at 0h UT and also now the offset of the latter from ERA.
This offset is to be saved by the caller, and the EO is added to it
later to get the offset of GST from the approximation to GMST computed
by updating GMST at 0h UT with the ERA rate.  Thus:
ERA-EO = GMST(0) + UT*ERA' - (GMST(0)-ERA(0)) - EO
Note: the formulation here comes from Capitaine, Wallace, & Chapront
2003, but it follows the expression in terms of angles, rather than
time, the two expressions being inconsistent.
In terms of time seconds and Julian centuries, this formulation is:
excess fract= 8640184.79448082 + 0.185544227 tt -8.833E-8 tt^2
-7.98827E-6 tt^3 -1.227E-8 tt^4  (not actually used)
sidtm0 = 24110.5493771 + 8639877.31737601 tu + 307.4771023 tt
+0.092772113 tt^2 -2.933E-8 tt^3 -1.99707E-6 tt^4 -2.453E-9 tt^5