% demo_tercel.m January 9, 2013
% this script demonstrates how to interact with the tercel.m
% and celter.m functions which transform vectors between the
% itrs (terrestrial) to gcrs (celestial) coordinate systems
% ported from Fortran version of NOVAS 3.1
% input
% tjdh = tdb or tt julian date, higher part
% tjdl = tdb or tt julian date, lower part
% xp = conventionally-defined x coordinate of cip
% with respect to itrs pole, arc seconds
% yp = conventionally-defined y coordinate of cip
% with respect to itrs pole, arc seconds
% vec1 = position vector, geocentric equatorial rectangular
% coordinates, referred to itrs axes
% output
% vec2 = position vector, geocentric equatorial rectangular
% coordinates, referred to gcrs axes (double precision)
% NOTE: this script computes the accumulated precession in right
% ascension using the analytic equation in capitaine et al. (2003),
% astronomy and astrophysics 412, 567-586, eq. (42).
% (see functions cioloc.m and eqxra.m)
% Orbital Mechanics with Matlab
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;
global imode lmode
global inutate psicor epscor
% initialize mode flags
imode = 2;
lmode = 2;
% initialize nutation algorithm
inutate = 1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% example from USNO/AA Technical Note 2010-04
% "Testing Coordinate Frame Transformations"
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% higher part of UT1 julian date
tjdh = 2400000.5d0;
% lower part of UT1 julian date
tjdl = 54195.49999916581146d0;
% polar coordinates (arc seconds)
xp = +0.0349282d0;
yp = +0.4833163d0;
% position vector
vec1(1) = 1.0d0;
vec1(2) = 0.0d0;
vec1(3) = 0.0d0;
% define and set the value of delta-t (TT-UT1, seconds)
deltat = 65.25607389d0;
setdt(deltat);
% celestial pole offsets
dx = 0.1750d0;
dy = -0.2259d0;
% "complete" julian date
tjd = tjdh + tjdl;
% compute celestial pole offsets
[psicor, epscor] = celpol(tjd, 2, dx, dy);
%%%%%%%%%%%%%%%%%%%%%%
% set calculation mode
%%%%%%%%%%%%%%%%%%%%%%
% set mode = 0 for cio-based method, full accuracy
% set mode = 1 for cio-based method, reduced accuracy
% set mode = 2 for equinox-based method, full accuracy
% set mode = 3 for equinox-based method, reduced accuracy
setmod(2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute position vector, geocentric equatorial rectangular
% coordinates, referred to gcrs axes (celestial system)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
vec2 = tercel (tjdh, tjdl, xp, yp, vec1);
% print results
clc; home;
fprintf('\ndemo_tercel\n');
fprintf('\nrotates a vector from the terrestrial to the celestial system');
fprintf('\n-------------------------------------------------------------\n');
mode = getmod;
if (mode == 0)
fprintf('\ncio-based method, full accuracy\n');
end
if (mode == 1)
fprintf('\ncio-based method, reduced accuracy\n');
end
if (mode == 2)
fprintf('\nequinox-based method, full accuracy\n');
end
if (mode == 3)
fprintf('\nequinox-based method, reduced accuracy\n');
end
[cdstr, utstr] = jd2str(tjd);
fprintf('\nUT1 calendar date ');
disp(cdstr);
fprintf('\nUT1 time ');
disp(utstr);
fprintf('\nUT1 Julian date %14.8f \n', tjd);
fprintf('\nposition vector, geocentric equatorial rectangular');
fprintf('\ncoordinates, referred to ITRS axes (terrestrial system)\n');
fprintf ('\n %+16.14e %+16.14e %+16.14e \n', vec1(1), vec1(2), vec1(3));
fprintf('\nposition vector, geocentric equatorial rectangular');
fprintf('\ncoordinates, referred to GCRS axes (celestial system)\n');
fprintf ('\n %+16.14e %+16.14e %+16.14e \n', vec2(1), vec2(2), vec2(3));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute transformation matrix
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
v1 = [1 0 0];
v2 = [0 1 0];
v3 = [0 0 1];
row1 = tercel (tjdh, tjdl, xp, yp, v1);
row2 = tercel (tjdh, tjdl, xp, yp, v2);
row3 = tercel (tjdh, tjdl, xp, yp, v3);
fprintf('\nterrestrial-to-celestial transformation matrix\n');
fprintf ('\n %+16.14e %+16.14e %+16.14e', row1(1), row1(2), row1(3));
fprintf ('\n %+16.14e %+16.14e %+16.14e', row2(1), row2(2), row2(3));
fprintf ('\n %+16.14e %+16.14e %+16.14e', row3(1), row3(2), row3(3));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% perform celestial-to-terrestrial transformation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
vec3 = celter (tjdh, tjdl, xp, yp, vec2);
fprintf('\n\nposition vector, geocentric equatorial rectangular');
fprintf('\ncoordinates, referred to ITRS axes (terrestrial system)\n');
fprintf ('\n %+16.14e %+16.14e %+16.14e \n', vec3(1), vec3(2), vec3(3));
fprintf('\n\n');