% moonpos calculates moon position from Astronomical Algorithms, Jean Meeus
%******************************************************************************
%
% MATLAB (R) is a trademark of The Mathworks (R) Corporation
%
% Function: moonpos
% Filename: moonpos.c
% Programmer: James Tursa
% Version: 1.1
% Date: April 7, 2009
% Copyright: (c) 2009 by James Tursa, All Rights Reserved
%
% This code uses the BSD License:
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:
%
% * Redistributions of source code must retain the above copyright
% notice, this list of conditions and the following disclaimer.
% * Redistributions in binary form must reproduce the above copyright
% notice, this list of conditions and the following disclaimer in
% the documentation and/or other materials provided with the distribution
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
% POSSIBILITY OF SUCH DAMAGE.
%
% moonpos calculates the moon position using the method in the book
% Astronomical Algorithms by Jean Meeus, Chapter 45, pp 307 - 313. The
% algorithm used is accurate to about 10 arcsec in longitude and 4 arcsec
% in latitude according to the book, although the time range is not
% mentioned. The periodic terms used in this algorithm are based on the
% Chapront ELP-2000/82 lunar theory, including the later improvements by
% Chapront. The results of this theory are referred to the mean equinox
% of date, but moonpos converts these to the apparent geocentric position
% of the moon referred to the true equinox of date.
%
% Building:
%
% moonpos requires that a mex routine be built (one time only). This
% process is typically self-building the first time you call the function
% as long as you have the files moonpos.m and moonpos.c in the same
% directory somewhere on the MATLAB path. If you need to manually build
% the mex function, here are the commands:
%
% >> mex -setup
% (then follow instructions to select a C / C++ compiler of your choice)
% >> mex moonpos.c
%
% If you have an older version of MATLAB, you may need to use this command:
%
% >> mex -DDEFINEMWSIZE moonpos.c
%
% Syntax
%
% B = moonpos(A)
% B = moonpos(A,S1)
% B = moonpos(A,S1,S2,...)
%
% Description
%
% A = a Nx1 or 1xN double vector containing the date values. By default,
% the data is expected to be in Serial Date Number format, but this
% can be changed to Julian Day, Modified Julian Day, Reduced Julian
% Day, or Truncated Julian Day with optional strings (see below).
%
% = a Nx6 double date matrix, where each row is a full date vector
% specifying year, month, day, hours, minutes, seconds in that order.
% Dates must be compatible with the MATLAB datenum function.
%
% = a Nx3 double date matrix, where each row is a partial date vector
% specifying year, month, and day in that order. Dates must be
% compatible with the MATLAB datenum function.
%
% = a NxM character array containing dates that are compatible with the
% datenum function.
%
% For the last three cases the datenum function is called automatically
% by moonpos to convert the dates to SDN format.
%
% B = a NxM double vector containing the moon position data. By default a
% Nx3 vector is returned where the three columns are right ascension (rad),
% declination (rad), and distance (km), but this can be changed with
% optional strings (see below). Output is true-of-date from earth
% center.
%
% Sn = Optional string argument(s) that alter input and output settings. You
% can specify more than one optional string. They are processed left to
% right and are case insensitive. If you specify mutually exclusive
% options (such as 'jd' and 'rjd'), the rightmost option will prevail.
%
% 'julian' or 'julian day' or 'jd'
% = Input vector A is in Julian Day format. i.e., the interval of time
% in days and fractions of a day since January 1, 4713 BC Greenwich noon
%
% 'reduced julian day' or 'reduced' or 'rjd'
% = Input vector A is in Reduced Julian Day format, RJD = JD - 2400000.
%
% 'modified julian day' or 'modified' or 'mjd'
% = Input vector A is in Modified Julian Day format, MJD = JD - 2400000.5
%
% 'truncated julian day' or 'truncated' or 'tjd'
% = Input vector A is in Truncated Julian Day format, TJD = JD - 2440000.5
%
% 'radec'
% = Output is Nx3 matrix in right ascension, declination, distance format
% in units of radians, radians, kilometers (this is the default).
% But note that you can change the units of the angles or distance
% (see below), and you can add the x, y, z cartesian output at the
% same time with the 'xyz' option and get a Nx6 matrix result.
%
% 'xyz'
% = Output is Nx3 matrix with x, y, z cartesian coordinates. Note
% that this can be specified in addition to 'radec' if you want
% a Nx6 matrix output with both data sets present. The order in
% which the 'radec' and 'xyz' options appear in the argument list
% determines the order in which the output appears in the result.
%
% 'feet' or 'ft'
% = Output distances will be in feet instead of km.
%
% 'meter' or 'meters' or 'm'
% = Output distances will be in meters instead of km.
%
% 'degree' or 'degrees' or 'deg'
% = Output angles will be in degrees instead of radians.
%
% 'hms'
% = Output right ascension will be in hh.mmsss... format.
%
% 'dms'
% = Output declination will be in dd.mmsss... format.
%
% 'nutationlow'
% = Use the slightly lower precision nutation formulae at the bottom
% of page 132 instead of the full TABLE 21.A terms.
%
% Examples
%
% >> format long
% >> x = moonpos(2448724.5,'jd') % The example on pages 312-313 of the book
%
% You must build the mex routine before you can use moonpos.
% Attempting to do so now ...
%
% Found file moonpos.c in C:\Program Files\MATLAB\R2006a\work\moonpos.c
%
% Now attempting to compile ...
% (If prompted, please press the Enter key and then select any C/C++
% compiler that is available, such as lcc.)
%
% mex('C:\Program Files\MATLAB\R2006a\work\moonpos.c')
%
% moonpos.c
% C:\Program Files\MATLAB\R2006a\work\moonpos.c(202) : error C2065: 'mwSize' : undeclared identifier
% C:\Program Files\MATLAB\R2006a\work\moonpos.c(202) : error C2146: syntax error : missing ';' before identifier 'i'
% :
% (a list of over 200 compile errors)
% :
% C:\Program Files\MATLAB\R2006a\work\moonpos.c(405) : warning C4244: '=' : conversion from 'double' to 'int', possible loss of data
% C:\Program Files\MATLAB\R2006a\work\moonpos.c(408) : error C2100: illegal indirection
%
% C:\PROGRAM FILES\MATLAB\R2006A\BIN\MEX.PL: Error: Compile of 'C:\Program Files\MATLAB\R2006a\work\moonpos.c' failed.
%
% Well, *that* didn't work ... this is often caused on older versions
% of MATLAB that do not have the typedef mwSize defined. So now we will
% attempt to compile it by manually defining mwSize ...
%
% mex('-DDEFINEMWSIZE','C:\Program Files\MATLAB\R2006a\work\moonpos.c')
%
% Success!
% mex moonpos.c build completed ... attempting to call this moonpos.
%
% x =
% 1.0e+005 *
% 0.00002350757329 0.00000240303299 3.68409684525585
%
% There are several things to note about the previous example. First, of
% course, is the long list of errors the first time it was run. This is
% because the MATLAB API for R2006a does not have the typedef for mwSize.
% A try-catch logic accounts for this, however, and then manually defines
% the mwSize for another compile try, and the second time it is a success.
% This process is all automated, so don't worry if you get a bunch of
% errors the very first time you try this function ... they can safely be
% ignored as long as you get the Success message.
%
% The second thing to note is that the return results are in radians, but
% a glance at the book shows the results in degrees and hms & dms. Before
% you start converting this result by hand to do the compare, realize that
% moonpos will do this for you. Just select the 'hms' and 'dms' options:
%
% >> x = moonpos(2448724.5,'jd','hms','dms')
% x =
% 1.0e+005 *
% 0.00008584523366 0.00013460611349 3.68409684525585
%
% The results given in the book are:
% RA = 8h 58m 45.2s
% DEC = 13deg 46' 06"
% DELTA = 368409.7 (corrected from errata, early editions are incorrect)
%
% So the results of moonpos match the book. Also note that the 2nd time
% moonpos was run there were no compile error messages because the mex
% file has already been built.
%
% Several of the values in early editions of the book are incorrect. They
% have been corrected with errata of later editions. The corrections are:
%
% http://www.willbell.com/MATH/AA-First%20Ed.%20Errata.pdf
% http://www.willbell.com/MATH/AA-Second%20Ed.%20Errata.pdf
%
% Sigma r = -16590875
% Delta = 385000.56 - 16590.875 = 368409.7
% pi = arcsine(6378.14 / 368409.7) = 0.991990 = 0deg 59' 31".2
%
% Some more examples:
%
% % Request x,y,z output instead of ra,dec,delta (results will be in km)
% >> x = moonpos(2448724.5,'jd','xyz')
% x =
% 1.0e+005 *
% -2.51640151355659 2.54391554623140 0.87680481376215
%
% % Request x,y,z output in feet instead km
% >> x = moonpos(2448724.5,'jd','xyz','feet')
% x =
% 1.0e+008 *
% -8.25591047754786 8.34617961362006 2.87665621313040
%
% % Request both ra,dec,delta and x,y,z output, with angles in degrees and
% % distances in meters
% >> x = moonpos(2448724.5,'jd','radec','deg','xyz','m')
% x =
% 1.0e+008 *
% Columns 1 through 4
% 0.00000134688474 0.00000013768365 3.68409684525585 -2.51640151355659
% Columns 5 through 6
% 2.54391554623140 0.87680481376215
%
% % Request both ra,dec,delta and x,y,z output, with angles in degrees and
% % distances in meters, but this time list the x,y,z first
% >> x = moonpos(2448724.5,'jd','xyz','radec','deg','m')
% x =
% 1.0e+008 *
% Columns 1 through 4
% -2.51640151355659 2.54391554623140 0.87680481376215 0.00000134688474
% Columns 5 through 6
% 0.00000013768365 3.68409684525585
%
% % Request position for a string date, with right ascension output in
% % hours, minutes, seconds, and declination output in degrees, minutes,
% % seconds. Note that we don't have to specify a time format with this
% % example since the result of the datenum function (called in the
% % background) is a Serial Date Number, which is the default for moonpos.
% >> x = moonpos('01-Jan-2000','hms','dms')
% x =
% 1.0e+005 *
% 0.00014263869309 -0.00008592808456 4.00927160087148
%
% % Plot the radial distance of the moon for the year 2009
% >> t = datenum('01-Jan-2009'):datenum('31-Dec-2009');
% >> x = moonpos(t);
% >> plot(t,x(:,3));
% >> title('2009 Moon Radial Distance from Earth');
% >> ylabel('km');
% >> dateaxis('x',12);
%
%**************************************************************************
function varargout = moonpos(varargin)
mname = mfilename;
disp(' ');
disp(['You must build the mex routine before you can use ' mname '.']);
disp('Attempting to do so now ...');
disp(' ');
mfull = mfilename('fullpath');
cname = [mfull '.c'];
if( isempty(dir(cname)) )
disp(['Cannot find the file ' mname '.c in the same directory as the']);
disp(['file ' mname '.m. Please ensure that they are in the same']);
disp('directory and try again. The following file was not found:');
disp(' ');
disp(cname);
disp(' ');
error(['Unable to compile ' mname '.c']);
else
disp(['Found file ' mname '.c in ' cname]);
disp(' ');
disp('Now attempting to compile ...');
disp('(If prompted, please press the Enter key and then select any C/C++');
disp('compiler that is available, such as lcc.)');
disp(' ');
disp(['mex(''' cname ''')']);
disp(' ');
try
mex(cname,'-output',mfull);
disp('Success!');
disp(['mex ' mname '.c build completed ... attempting to call this ' mname '.']);
disp(' ');
catch
disp(' ');
disp('Well, *that* didn''t work ... this is often caused on older versions');
disp('of MATLAB that do not have the typedef mwSize defined. So now we will');
disp('attempt to compile it by manually defining mwSize ...');
disp(' ');
try
disp(' ');
disp(['mex(''-DDEFINEMWSIZE'',''' cname ''')']);
disp(' ');
mex('-DDEFINEMWSIZE',cname,'-output',mfull);
disp('Success!');
disp(['mex ' mname '.c build completed ... attempting to call this ' mname '.']);
disp(' ');
catch
disp('Hmmm ... That didn''t work either.');
disp(' ');
disp('The mex command failed. This may be because you have already run');
disp('mex -setup and selected a non-C compiler, such as Fortran. If this');
disp('is the case, then rerun mex -setup and select a C/C++ compiler.');
disp(' ');
error(['Unable to compile ' mname '.c']);
end
end
[varargout{1:nargout}] = moonpos(varargin{:});
end
end
