// 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 [varargout] = CL_mod_moonSunG50(cjd,body) // Moon and Sun positions in Gamma50 (Veis) reference frame - DEPRECATED // // Calling Sequence // [u_moon,r_moon,u_sun,r_sun] = CL_mod_moonSunG50(cjd [,body="b"]) // [u_moon,r_moon] = CL_mod_moonSunG50(cjd, body="m") // [u_sun,r_sun] = CL_mod_moonSunG50(cjd, body="s") // // Description // //

This function is deprecated.

//

Replacement functions: CL_eph_moon and // CL_eph_sun

//

// //

Computes the Moon and/or Sun positions relative to Gamma50 (Veis) reference frame using // a simplified version of BROWN theories for the moon and NEWCOMB for the Sun.

//

The Moon or Sun position is defined by a unit vector u_body // and a distance r_body from the Earth centre to the centre of the body.

//

Either the Moon position, the Sun position or both can be computed.

//

//

Note: The time scale is TT and not UTC as wrongly indicated in previous versions.

// // // // Parameters // cjd : Modified (1950.0) julian day (Time scale: TT) (1xN) // body : (string, optional) Type can be 'm' or 'moon' for the moon,'s' or 'sun' for the Sun (or 'b' or 'both' for both). // u_body : Unit vector from Earth to body [rx;ry;rz] (3xN) // r_body : Distance from Earth to body [m] (1xN) // // Bibliography // 1 CNES - MSLIB FORTRAN 90, Volume S (ms_pos_soleil_lune) // // Authors // CNES - DCT/SB // // Examples // // Sun position in J2000 frame on the 25th of october 2008 and 2009 // cjd = CL_dat_cal2cjd([2008 2009],[10 10],[25 25]); // [r_sun,rs] = CL_mod_moonSunG50(cjd,'s'); // pos_sun_G50 = CL_dMult(r_sun,rs); // M = CL_fr_G502J2000Mat(cjd); // pos_sun_J2000 = M*pos_sun_G50; // // // Sun right ascension and declination at the same dates : // pos_sph = CL_co_car2sph(pos_sun_J2000); // alpha_sun = pos_sph(1,:); // delta_sun = pos_sph(2,:); // //outputs decision // Declarations: // Code: CL__warnDeprecated(); // deprecated function [lhs,rhs] = argn(); if ~(rhs == 1 | rhs == 2) CL__error("Invalid number of input arguments"); end if (~exists('body','local')); body = "both"; end compute_moon = %f; compute_sun = %f; if (rhs == 1) compute_moon = %t; compute_sun = %t; else if (body == 'm' | body == 'moon') compute_moon = %t; elseif (body == 's' | body == 'sun') compute_sun = %t; elseif (body == 'b' | body == 'both') compute_moon = %t; compute_sun = %t; else CL__error("Invalid ''body'' value"); end end if( body ~= 'both' & body ~= 'b' & lhs > 2) CL__error("Number of outputs should be 2 for ''moon'' or ''sun'' "); end // -------------------------- // ephemeris computation // -------------------------- t = cjd - 10000; f = 225.768 + 13.2293505.*t; f = pmodulo(CL_deg2rad(f),2*%pi); xl = 185.454 + 13.064992.*t; xl = pmodulo(CL_deg2rad(xl),2*%pi); d = 11.786+ 12.190749.*t; d = pmodulo(CL_deg2rad(d),2*%pi); xlp = 134.003 + 0.9856.*t; xlp = pmodulo(CL_deg2rad(xlp),2*%pi); g = 282.551+ 0.000047.*t; g = pmodulo(CL_deg2rad(g),2.*%pi); e = 23.44223-0.00000035626.*t; e = pmodulo(CL_deg2rad(e),2.*%pi); ce = cos(e); se = sin(e); rot = 0.0000006119022.*cjd; cr = cos(rot); sr = sin(rot); if compute_moon //Position de la lune dl = 10976.*sin(xl) - 2224.*sin(xl-d-d) + 1149.*sin(d+d); dl = dl + 373 .*sin(xl+xl) - 324 .*sin(xlp) - 200 .*sin(f+f); dl = dl - 103 .*sin(xl+xl-d-d) - 100 .*sin(xl+xlp-d-d); dl = dl + 93.*sin(xl+d+d); dl = dl - 80 .*sin(xlp-d-d) + 72 .*sin(xl-xlp) - 61 .*sin(d); dl = dl - 53 .*sin(xl+xlp); dl = dl + 14 .*sin(xl-xlp-d-d) + 19 .*sin(xl-f-f); dl = dl - 19 .*sin(xl-4 .*d); dl = dl + 17 .*sin(3 .*xl) - 27 .*sin(f+f-d-d); dl = dl - 12 .*sin(xlp+d+d); dl = dl - 22 .*sin(xl+f+f) - 15 .*sin(xl+xl-4 .*d); dl = dl + 7 .*sin(xl+xl+d+d) + 9 .*sin(xl-d); dl = dl - 6 .*sin(3 .*xl-d-d); dl = dl + 7 .*sin(4 .*d) + 9 .*sin(xlp+d) + 7 .*sin(xl-xlp+d+d); dl = dl + 5 .*sin(xl+xl-xlp); dl = dl.*0.00001; b = 8950 .*sin(f) + 490 .*sin(xl+f) + 485 .*sin(xl-f); b = b - 302 .*sin(f-d-d); b = b - 97 .*sin(xl-f-d-d) - 81 .*sin(xl+f-d-d); b = b + 57 .*sin(f+d+d); b = b - 14 .*sin(xlp+f-d-d) + 16 .*sin(xl-f+d+d); b = b + 15.*sin(xl+xl-f) + 30 .*sin(xl+xl+f); b = b - 6 .*sin(xlp-f+d+d) - 7 .*sin(xl+xl+f-d-d); b = b + 7 .*sin(xl+f+d+d); b = b.*0.00001; u = 68.341 + 13.176397 .*t; u = pmodulo(CL_deg2rad(u),2.*%pi) + dl; cu = cos(u); su = sin(u); cb = cos(b); sb = sin(b); rx = cu.*cb; ry = su.*cb.*ce-sb.*se; rz = sb.*ce + su.*cb.*se; q = rx.*cr + ry.*sr; ry = ry.*cr - rx.*sr; rx = q; dasr = 5450 .*cos(xl) + 1002 .*cos(xl-d-d) + 825 .*cos(d+d); dasr = dasr + 297 .*cos(xl+xl) + 90 .*cos(xl+d+d); dasr = dasr + 56 .*cos(xlp-d-d); dasr = dasr + 42 .*cos(xl+xlp-d-d) + 34 .*cos(xl-xlp); dasr = dasr - 12 .*cos(xlp) - 29 .*cos(d) - 21 .*cos(xl-f-f); dasr = dasr + 18 .*cos(xl-4 .*d) - 28 .*cos(xl+xlp); dasr = dasr + 11 .*cos(xl+xl-4 .*d) + 18 .*cos(3 .*xl); dasr = dasr - 9 .*cos(xlp+d+d) - 7 .*cos(xl-xlp-d-d); dasr = dasr + 7 .*cos(xl-xlp+d+d); dasr = dasr - 9 .*cos(xl+xl-d-d) + 8 .*cos(xl+xl+d+d); dasr = dasr + 8 .*cos(4 .*d); asrl = 1 + 0.00001.*dasr; rl = 384389.3 ./ asrl; r_moon=[rx;ry;rz]; rl=rl*1000; // conversion en m end if compute_sun // Position du soleil cl = 99972 .*cos(xlp+g) + 1671 .*cos(xlp+xlp+g) - 1678 .*cos(g); cl = cl + 32 .*cos(3 .*xlp+g) + cos(4 .*xlp+g); cl = cl - 4 .*cos(g-xlp) - 2 .*cos(xlp-d+g) + 4 .*cos(f-d); cl = cl - 4 .*cos(xlp+xlp-f+d+g+g) + 2 .*cos(xlp+d+g); cl = cl.*0.00001; sl = 99972 .*sin(xlp+g) + 1671 .*sin(xlp+xlp+g) - 1678 .*sin(g); sl = sl + 32 .*sin(3 .*xlp+g) + sin(4 .*xlp+g); sl = sl - 4 .*sin(g-xlp) - 2 .*sin(xlp-d+g) + 4 .*sin(f-d); sl = sl - 4 .*sin(xlp+xlp-f+d+g+g) + 2 .*sin(xlp+d+g); sl = sl .* 0.00001; q = sqrt(cl.*cl+sl.*sl); cl = cl./q; sl = sl./q; sx = cl; sy = sl.*ce; sz = sl.*se; q = sx.*cr + sy.*sr; sy = sy.*cr - sx.*sr; sx = q; dasr = 1672.2.*cos(xlp)+28 .*cos(xlp+xlp)-0.35.*cos(d); asrs = 1 + 0.00001 .*dasr; rs = 149597870 ./asrs; r_sun=[sx;sy;sz]; rs=rs*1000; // conversion en m end if (body == 'm' | body == 'moon') varargout = list(r_moon,rl); elseif (body == 's' | body == 'sun') varargout = list(r_sun,rs); elseif (body == 'b' | body == 'both') varargout = list(r_moon,rl,r_sun,rs); end endfunction