// Copyright (c) CNES 2008 // // This software is part of CelestLab, a CNES toolbox for Scilab // // This software is governed by the CeCILL license under French law and // abiding by the rules of distribution of free software. You can use, // modify and/ or redistribute the software under the terms of the CeCILL // license as circulated by CEA, CNRS and INRIA at the following URL // 'http://www.cecill.info'. function [pos2,vel2] = CL__fr_fastConvert(frame1,frame2,jd, pos1,vel1) // Conversion of positions in various frames using fast algorithms. // // Calling Sequence // [pos2,vel2] = CL__fr_fastConvert(frame1,frame2,jd,pos1,vel1) // // Description // //

Converts position and velocity from frame 1 to frame 2 at given dates, using fast algorithms.

//

Angular velocity vector is supposed to be zero because all the frames are inertial or quasi-inertial.

//

The models used are very simple and efficient.

//

Available frames are : "GCRS", "EOD", "MOD", "CIRS" and "Veis"

//

//

Accuracy (arcsec):

//

Period MOD <-> EOD EOD <-> GCRS GCRS <-> CIRS CIRS <-> Veis

//

2000-2100 0.002 0.0005 0.9 5e-005

//

1000-3000 2 0.6 150 0.0005

// // // // Parameters // frame1 : (string) Name of frame 1 (1x1) // frame2 : (string) Name of frame 2 (1x1) // jd : Julian date (TT) (1xN) // pos1 : Position in frame 1 (3xN). // vel1 : (optional) Velocity in frame 1 (3xN). Can be [] and in that case vel2 will be []. Default is [] // pos2 : Position in frame 2 (3xN) // vel2 : Velocity in frame 2 (3xN) // // Authors // CNES - DCT/SB // // Examples // // Convert position only // jd = [21915,21916]; // TT time scale ! // pos1 = rand(3,2); // pos2 = CL__fr_fastConvert("EOD","GCRS",jd,pos1,[]) // // // Convert position and velocity // jd = [21915,21916]; // TT time scale ! // pos1 = rand(3,2); // vel1 = rand(3,2); // [pos2, vel2] = CL__fr_fastConvert("EOD","GCRS",jd,pos1,vel1) // Internal function : Frame transformation of a rotation about one axis of a given angle. // (Fast implementation) // naxe = 1 for rotation about X or 3 for rotation about Z (1x1) // ang = (1xN) // direction = 1 or -1 (1x1) // pos1 = (3xN) can be empty, in that case pos2 will be empty // vel1 = (3xN) can be empty, in that case vel2 will be empty // pos2 = (3xN) or empty // vel2 = (3xN) or empty function [pos2,vel2] = CL__fr_fastRot(naxe, ang, direction, pos1, vel1) pos2 = []; vel2 = []; cosa = cos(direction * ang); sina = sin(direction * ang); // Compute position if (pos1 <> []) if (naxe == 1) pos2 = [ pos1(1,:); .. cosa .* pos1(2,:) + sina .* pos1(3,:) ; .. -sina .* pos1(2,:) + cosa .* pos1(3,:) ]; elseif (naxe == 3) pos2 = [ cosa .* pos1(1,:) + sina .* pos1(2,:) ; .. -sina .* pos1(1,:) + cosa .* pos1(2,:) ; .. pos1(3,:) ]; end end // Compute velocity // vel2 = M * vel1 // (omega ignored) if (vel1 <> []) if (naxe == 1) vel2 = [ vel1(1,:); .. cosa .* vel1(2,:) + sina .* vel1(3,:) ; .. -sina .* vel1(2,:) + cosa .* vel1(3,:) ]; elseif (naxe == 3) vel2 = [ cosa .* vel1(1,:) + sina .* vel1(2,:); .. -sina .* vel1(1,:) + cosa .* vel1(2,:) ; .. vel1(3,:) ]; end end endfunction // MOD --> EOD // Rotation about X axis of IAU2006 obliquity troncated (See CL__iers_obl06) function [pos2,vel2] = CL__fr_fastMOD2EOD(TT,pos1,vel1,direction) epsa = (84381.406 - 46.836769 * TT) * (1/3600)*(%pi/180); [pos2,vel2] = CL__fr_fastRot(1, epsa, direction, pos1, vel1); endfunction // EOD --> GCRS // Rotation about Z axis of IAU2006 psib troncated at second degree (See CL__iers_pfw06) // Rotation about X axis of IAU2006 -phib troncated at second degree (See CL__iers_pfw06) // Rotation about Z axis of IAU2006 -gamb troncated at second degree (See CL__iers_pfw06) function [pos2,vel2] = CL__fr_fastEOD2GCRS(TT,pos1,vel1,direction) TT2 = TT.*TT; psib = (-0.041775 + 5038.481484 * TT + 1.5584175 * TT2) * (1/3600)*(%pi/180); phib = (84381.412819 - 46.811016 * TT + 0.0511268 * TT2) * (1/3600)*(%pi/180); gamb = (-0.052928 + 10.556378 * TT + 0.4932044 * TT2) * (1/3600)*(%pi/180); if (direction == 1) [pos2,vel2] = CL__fr_fastRot(3, psib, direction, pos1, vel1); [pos2,vel2] = CL__fr_fastRot(1, -phib, direction, pos2, vel2); [pos2,vel2] = CL__fr_fastRot(3, -gamb, direction, pos2, vel2); else [pos2,vel2] = CL__fr_fastRot(3, -gamb, direction, pos1, vel1); [pos2,vel2] = CL__fr_fastRot(1, -phib, direction, pos2, vel2); [pos2,vel2] = CL__fr_fastRot(3, psib, direction, pos2, vel2); end endfunction // GCRS --> CIRS // See Model CPNd (Concise CIO based precession-nutation formulations. N. Capitaine and P. T. Wallace) function [pos2,vel2] = CL__fr_fastGCRS2CIRS(TT,pos1,vel1,direction) // Delaunay variables (Table C.1) F = 1.6279050815 + 8433.4661569164 * TT; D = 5.1984665887 + 7771.3771455937 * TT; om = 2.1824391966 - 33.7570459536 * TT; // X,Y (Table G1) (conversion from microarcsec to radians) X = (2004191898 * TT - 6844318 * sin(om) - 523908 * sin(2*F-2*D+2*om)) * (1/3600)*(%pi/180)*1e-6; Y = (-22407275 * TT.*TT + 9205236 * cos(om) + 573033 * cos(2*F-2*D+2*om)) * (1/3600)*(%pi/180)*1e-6; // M (Table 1, Method C) // M*pos1 pos2 = [ pos1(1,:) - direction * X.*pos1(3,:) ; ... pos1(2,:) - direction * Y.*pos1(3,:) ; ... direction * X .* pos1(1,:) + direction * Y .* pos1(2,:) + pos1(3,:)]; vel2 = []; if (vel1 <> []) // M*vel1 vel2 = [ vel1(1,:) - direction * X.*vel1(3,:) ; ... vel1(2,:) - direction * Y.*vel1(3,:) ; ... direction * X .* vel1(1,:) + direction * Y .* vel1(2,:) + vel1(3,:)]; end endfunction // CIRS --> Veis // Rotation about Z axis of angle ERA (IAU2000) minus angle TSV // (See CL__iers_era00 and CL__iers_sidTimeVeis) // (s prime neglected) function [pos2,vel2] = CL__fr_fastCIRS2Veis(TT,pos1,vel1,direction) ang = 0.0111817446159 + 5.58629142944e-08 * TT ; [pos2,vel2] = CL__fr_fastRot(3, ang, direction, pos1, vel1); endfunction // Code: if (argn(2) < 1 | argn(2) > 6) CL__error("Invalid number of input arguments"); end if (~exists("vel1","local")); vel1 = []; end; // If output velocity not requested: dont compute it if (argn(1) == 1); vel1 = []; end; // ICRS <=> GCRS if (frame1 == "ICRS"); frame1 = "GCRS"; end if (frame2 == "ICRS"); frame2 = "GCRS"; end frames = ["MOD", "EOD", "GCRS", "CIRS", "Veis"]; I1 = find(frame1 == frames); I2 = find(frame2 == frames) if (I1 ==[] | I2 == []) CL__error("Unrecognized frame names"); end if (frame1 == frame2) pos2 = pos1; vel2 = vel1; return; // <-- RETURN! end // Number of centuries since epoch 2000 12h (TT time scale) TT = (jd - 2451545) / 36525; // func_ref = list of frame conversion functions // "MOD" --> "EOD" --> "GCRS" --> "CIRS" --> "Veis" func_ref = list(CL__fr_fastMOD2EOD, CL__fr_fastEOD2GCRS, CL__fr_fastGCRS2CIRS, CL__fr_fastCIRS2Veis); // func_tab = array of frame conversion numbers (corresponding to func_ref) // direction = direction of frame conversion (1 if normal direction, -1 otherwise) if (I1 < I2) direction = 1; func_tab = I1:(I2-1); else direction = -1; func_tab = (I1-1):-1:I2; end // Computation [pos2,vel2] = func_ref(func_tab(1))(TT,pos1,vel1,direction); for k = 2 : size(func_tab,2) [pos2,vel2] = func_ref(func_tab(k))(TT,pos2,vel2,direction); end endfunction // Not used but stored for possible future use: //// GCRS --> CIRS //// See Model CPNc (Concise CIO based precession-nutation formulations. N. Capitaine and P. T. Wallace) //function [pos2,vel2] = CL__fr_fastGCRS2CIRSc(TT,pos1,vel1,direction) // // Delaunay variables (Table C.1) // l = 2.3555557435 + 8328.6914257191 * TT; // lp = 6.2400601269 + 628.3019551714 * TT; // F = 1.6279050815 + 8433.4661569164 * TT; // D = 5.1984665887 + 7771.3771455937 * TT; // om = 2.1824391966 - 33.7570459536 * TT; // // TT2 = TT.*TT; // TT3 = TT2.*TT; // // // X,Y (Table F1) (conversion from microarcsec TTo radians) // X = -17251 +... // 2004191898 * TT ... // -429783 * TT2 ... // -198618 * TT3 ... // -6844318 * sin(om) ... // -3310 * TT .* sin(om) ... // +205833 * TT .* cos(om) ... // +82169 * sin(2*om) ... // +2521 * sin(2*D) ... // +5096 * sin(2*F-2*D+om) ... // -523908 * sin(2*F-2*D+2*om) ... // +12814 * TT .* cos(2*F-2*D+2*om) ... // -15407 * sin(2*F+om) ... // -90552 * sin(2*F+2*om) ... // -8585 * sin(lp-2*F+2*D-2*om) ... // +58707 * sin(lp) ... // -20558 * sin(lp+2*F-2*D+2*om) ... // -4911 * sin(l-2*F-2*om) ... // -6245 * sin(l-2*D) ... // +28288 * sin(l) ... // +2512 * sin(l+om) ... // -11992 * sin(l+2*F+2*om); // // // Y = -5530 ... // -25896 * TT ... // -22407275 * TT2 ... // +9205236 * cos(om) ... // +153042 * TT .* sin(om) ... // -89618 * cos(2*om) ... // -6918 * cos(2*F-2*D+om) ... // +573033 * cos(2*F-2*D+2*om) ... // +11714 * TT .* sin(2*F-2*D+2*om) ... // +20070 * cos(2*F+om) ... // +97847 * cos(2*F+2*om) ... // -9593 * cos(lp-2*F+2*D-2*om) ... // +7387 * cos(lp) ... // +22438 * cos(lp+2*F-2*D+2*om) ... // +2555 * cos(l-2*F-2*D-2*om) ... // -5331 * cos(l-2*F-2*om) ... // +3144 * cos(l-om) ... // -3324 * cos(l+om) ... // +2636 * cos(l+2*F+om) ... // +12903 * cos(l+2*F+2*om); // // X = X * (1/3600)*(%pi/180)*1e-6; // Y = Y * (1/3600)*(%pi/180)*1e-6; // s = (3809 * TT - 72574 * TT3 - 2641 * sin(om)) * (1/3600)*(%pi/180)*1e-6 - X.*Y / 2; // // // M (Table 1, Method B) // vel2 = []; // if (direction == 1) // // M * pos1 // pos2 = [ (1-X.^2/2) .* pos1(1,:) + (-s-X.*Y/2) .* pos1(2,:) + (-X) .* pos1(3,:) ; ... // (s-X.*Y/2) .* pos1(1,:) + (1-Y.^2/2) .* pos1(2,:) + (-Y-s.*X) .* pos1(3,:) ; ... // (X) .* pos1(1,:) + (Y) .* pos1(2,:) + (1-(X.^2+Y.^2)/2) .* pos1(3,:)]; // if (vel1 <> []) // // M * vel1 // vel2 = [(1-X.^2/2) .* vel1(1,:) + (-s-X.*Y/2) .* vel1(2,:) + (-X) .* vel1(3,:) ; ... // (s-X.*Y/2) .* vel1(1,:) + (1-Y.^2/2) .* vel1(2,:) + (-Y-s.*X) .* vel1(3,:) ; ... // (X) .* vel1(1,:) + (Y) .* vel1(2,:) + (1-(X.^2+Y.^2)/2) .* vel1(3,:)]; // end // // else // // M * pos1 // pos2 = [ (1-X.^2/2) .* pos1(1,:) + (s-X.*Y/2) .* pos1(2,:) + (X) .* pos1(3,:) ; ... // (-s-X.*Y/2) .* pos1(1,:) + (1-Y.^2/2) .* pos1(2,:) + (Y) .* pos1(3,:) ; ... // (-X) .* pos1(1,:) + (-Y-s.*X) .* pos1(2,:) + (1-(X.^2+Y.^2)/2) .* pos1(3,:)]; // if (vel1 <> []) // // M * vel1 // vel2 = [ (1-X.^2/2) .* vel1(1,:) + (s-X.*Y/2) .* vel1(2,:) + (X) .* vel1(3,:) ; ... // (-s-X.*Y/2) .* vel1(1,:) + (1-Y.^2/2) .* vel1(2,:) + (Y) .* vel1(3,:) ; ... // (-X) .* vel1(1,:) + (-Y-s.*X) .* vel1(2,:) + (1-(X.^2+Y.^2)/2) .* vel1(3,:)]; // end // end //endfunction