Code covered by the BSD License  

Highlights from
A MATLAB Script for Terrestrial and Celestial Coordinate Conversion

A MATLAB Script for Terrestrial and Celestial Coordinate Conversion

by

 

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

demo_tercel.m
% 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');

Contact us