Code covered by the BSD License

# A MATLAB Script for Terrestrial and Celestial Coordinate Conversion

### David Eagle (view profile)

MATLAB functions and script for performing terrestrial to/from celestial coordinate transformations.

sidtim (tjdh, tjdl, k)
```function gst = sidtim (tjdh, tjdl, k)

% this function computes the greenwich sidereal time
% (either mean or apparent) at julian date tjdh + tjdl.

%      tjdh   = ut1 julian date, high-order part (in)
%      tjdl   = ut1 julian date, low-order part (in)
%               the julian date may be split at any point, but
%               for highest precision, set tjdh to be the integral
%               part of the julian date, and set tjdl to be the
%               fractional part
%      k      = time selection code (in)
%               set k=0 for greenwich mean sidereal time
%               set k=1 for greenwich apparent sidereal time
%      gst    = greenwich (mean or apparent) sidereal time
%               in hours (out)

% note:  see also subroutine setdt to set the value of delta-t
% (delta-t = tt - ut1) to be used here.

% ported from NOVAS 3.1

%%%%%%%%%%%%%%%%%%%%%%%

degcon = 180.0d0 / pi;

% t0 = tdb julian date of epoch j2000.0 (tt)

t0 = 2451545.0d0;

unitx(1) = 1.0d0;
unitx(2) = 0.0d0;
unitx(3) = 0.0d0;

deltat = getdt;

% time argument for precession and nutation components of sidereal
% time is tdb

utjd = tjdh + tjdl;

ttjd = utjd + deltat;

tdbjd = ttjd;

[a, secdif] = novas_times (tdbjd);

tdbjd = ttjd + secdif / 86400.0d0;

% get method/accuracy mode

mode = getmod;

if (mode == 2)

% equinox-based mode
% see usno circular 179, section 2.6.2

% get -1 times the mean or true right ascension of the cio

rcio = eqxra (tdbjd, k);

% get earth rotation angle

theta = erot (tjdh, tjdl);

% combine to obtain sidereal time

gst = mod (theta / 15.0d0 - rcio, 24.0d0);

if (gst < 0.0d0)

gst = gst + 24.0d0;

end

else

% cio-based mode
% see usno circular 179, section 6.5.4

% get earth rotation angle

theta = erot (tjdh, tjdl);

% obtain the basis vectors, in the gcrs, of the celestial
% intermediate system

[rcio, kcio] = cioloc (tdbjd);

if (rcio == 99.0d0)

pause

end

[x, y, z] = ciobas (tdbjd, rcio, kcio);

% compute the direction of the true equinox in the gcrs

w1 = nutate (-tdbjd, unitx);

w2 = preces (tdbjd, w1, t0);

eq = frame (w2, -1);

% compute the hour angle of the equinox wrt the tio meridian
% (near greenwich, but passes through the cip and tio)

haeq = theta - atan2 (eq(1) * y(1) + eq(2) * y(2) + eq(3) * y(3), ...
eq(1) * x(1) + eq(2) * x(2) + eq(3) * x(3)) * degcon;

% for mean sidereal time, obtain the equation of the equinoxes
% and subtract it

if (k == 0)

[a, a, ee, a, a] = etilt (tdbjd);

haeq = haeq - ee / 240.0d0;

end

haeq = mod (haeq, 360.0d0) / 15.0d0;

if (haeq < 0.0d0)

haeq = haeq + 24.0d0;

end

gst = haeq;

end

```