diff --git a/doc/_static/moonazimuth_compare.png b/doc/_static/moonazimuth_compare.png new file mode 100644 index 00000000..b355f9b4 Binary files /dev/null and b/doc/_static/moonazimuth_compare.png differ diff --git a/doc/_static/moonheight_compare.png b/doc/_static/moonheight_compare.png new file mode 100644 index 00000000..b8abaa8f Binary files /dev/null and b/doc/_static/moonheight_compare.png differ diff --git a/doc/_static/moonphase_compare.png b/doc/_static/moonphase_compare.png new file mode 100644 index 00000000..7ac0379a Binary files /dev/null and b/doc/_static/moonphase_compare.png differ diff --git a/doc/source/api.rst b/doc/source/api.rst new file mode 100644 index 00000000..96640446 --- /dev/null +++ b/doc/source/api.rst @@ -0,0 +1,25 @@ + +API +--- + +Orbital computations +~~~~~~~~~~~~~~~~~~~~ + +.. automodule:: pyorbital.orbital + :members: + :undoc-members: + +TLE handling +~~~~~~~~~~~~ + +.. automodule:: pyorbital.tlefile + :members: + :undoc-members: + +Astronomical computations +~~~~~~~~~~~~~~~~~~~~~~~~~ + +.. automodule:: pyorbital.astronomy + :members: + :undoc-members: + diff --git a/doc/source/astronomy.rst b/doc/source/astronomy.rst new file mode 100644 index 00000000..9f62757b --- /dev/null +++ b/doc/source/astronomy.rst @@ -0,0 +1,12 @@ +Computing astronomical parameters +--------------------------------- +The astronomy module enables computation of certain parameters of interest for +satellite remote sensing for instance the Sun-zenith angle: + + >>> from pyorbital import astronomy + >>> from datetime import datetime + >>> utc_time = datetime(2012, 5, 15, 15, 45) + >>> lon, lat = 12, 56 + >>> astronomy.sun_zenith_angle(utc_time, lon, lat) + 62.685986438071602 + diff --git a/doc/source/conf.py b/doc/source/conf.py index b753217f..a82d6040 100644 --- a/doc/source/conf.py +++ b/doc/source/conf.py @@ -22,14 +22,16 @@ sys.path.insert(0, os.path.abspath('../../pyorbital')) from pyorbital.version import __version__ -# -- General configuration ----------------------------------------------------- +# -- General configuration ----------------------------------------------- # If your documentation needs a minimal Sphinx version, state it here. #needs_sphinx = '1.0' # Add any Sphinx extension module names here, as strings. They can be extensions # coming with Sphinx (named 'sphinx.ext.*') or your custom ones. -extensions = ['sphinx.ext.autodoc', 'sphinx.ext.doctest', 'sphinx.ext.coverage'] +#extensions = ['sphinx.ext.autodoc', 'sphinx.ext.doctest', 'sphinx.ext.coverage', "sphinxtogithub"] +extensions = [ + 'sphinx.ext.autodoc', 'sphinx.ext.doctest', 'sphinx.ext.coverage'] # Add any paths that contain templates here, relative to this directory. templates_path = ['.templates'] @@ -45,7 +47,7 @@ # General information about the project. project = u'pyorbital' -copyright = u'2012-2015, The Pytroll crew' +copyright = u'2012-2016, The Pytroll crew' # The version info for the project you're documenting, acts as replacement for # |version| and |release|, also used in various other places throughout the @@ -91,7 +93,7 @@ #modindex_common_prefix = [] -# -- Options for HTML output --------------------------------------------------- +# -- Options for HTML output --------------------------------------------- # The theme to use for HTML and HTML Help pages. See the documentation for # a list of builtin themes. @@ -124,7 +126,7 @@ # Add any paths that contain custom static files (such as style sheets) here, # relative to this directory. They are copied after the builtin static files, # so a file named "default.css" will overwrite the builtin "default.css". -html_static_path = ['.static'] +html_static_path = ['../_static'] # If not '', a 'Last updated on:' timestamp is inserted at every page bottom, # using the given strftime format. @@ -171,7 +173,7 @@ htmlhelp_basename = 'pyorbitaldoc' -# -- Options for LaTeX output -------------------------------------------------- +# -- Options for LaTeX output -------------------------------------------- # The paper size ('letter' or 'a4'). #latex_paper_size = 'letter' @@ -210,7 +212,7 @@ #latex_domain_indices = True -# -- Options for manual page output -------------------------------------------- +# -- Options for manual page output -------------------------------------- # One entry per manual page. List of tuples # (source start file, name, description, authors, manual section). diff --git a/doc/source/index.rst b/doc/source/index.rst index 28bccd8e..7975a28a 100644 --- a/doc/source/index.rst +++ b/doc/source/index.rst @@ -21,16 +21,8 @@ TLE files --------- Pyorbital has a module for parsing NORAD TLE-files - >>> from pyorbital import tlefile - >>> tle = tlefile.read('noaa 18', '/path/to/my/tle_file.txt') - >>> tle.inclination - 99.043499999999995 -If no path is given pyorbital tries to read the earth observation TLE-files from celestrak.com - -Computing satellite postion ---------------------------- -The orbital module enables computation of satellite position and velocity at a specific time: +.. _github: http://github.com/pytroll/pyorbital >>> from pyorbital.orbital import Orbital >>> from datetime import datetime @@ -71,49 +63,21 @@ If we take a TLE from one week earlier we get a slightly different result: (104.1539184988462, 79.328272480878141, 838.81555967963391) +>>>>>>> master -Computing astronomical parameters ---------------------------------- -The astronomy module enables computation of certain parameters of interest for satellite remote sensing for instance the Sun-zenith angle: - - >>> from pyorbital import astronomy - >>> from datetime import datetime - >>> utc_time = datetime(2012, 5, 15, 15, 45) - >>> lon, lat = 12, 56 - >>> astronomy.sun_zenith_angle(utc_time, lon, lat) - 62.685986438071602 - -API ---- - -Orbital computations -~~~~~~~~~~~~~~~~~~~~ - -.. automodule:: pyorbital.orbital - :members: - :undoc-members: - -TLE handling -~~~~~~~~~~~~ - -.. automodule:: pyorbital.tlefile - :members: - :undoc-members: - -Astronomical computations -~~~~~~~~~~~~~~~~~~~~~~~~~ - -.. automodule:: pyorbital.astronomy - :members: - :undoc-members: +.. toctree:: + :maxdepth: 2 + tle + satellite_position + astronomy + moon_calculations + api + -.. Contents: - .. toctree:: - :maxdepth: 2 - Indices and tables - ================== - * :ref:`genindex` - * :ref:`modindex` - * :ref:`search` +Indices and tables +================== +* :ref:`genindex` +* :ref:`modindex` +* :ref:`search` diff --git a/doc/source/moon_calculations.rst b/doc/source/moon_calculations.rst new file mode 100644 index 00000000..9ca1c39b --- /dev/null +++ b/doc/source/moon_calculations.rst @@ -0,0 +1,81 @@ +Moon phase and position with Pyorbital +====================================== + +It is possible to calculate the phase and position of the moon at any given +time in the past or future. This is done by an implementation of the +astronimical methods described at Stjaernhimlen_. The calculations have been +compared to the pyephem library, and deviations are small, see below. + + +Computing the phase of the moon +------------------------------- + + >>> from datetime import datetime + >>> from pyorbital.moon_phase import moon_phase + >>> time_t = datetime(2011, 12, 1, 12) + >>> print moon_phase(time_t) + 0.409752921579 + +It works also with datetime sequences: + + >>> import numpy as np + >>> from pyorbital.moon_phase import moon_phase + >>> time_t = np.arange('2005-02', '2005-03', dtype='datetime64[D]') + >>> print moon_phase(time_t) + [ 6.36537786e-01 5.32201222e-01 4.23413367e-01 3.15141455e-01 + 2.13350898e-01 1.24769455e-01 5.62102473e-02 1.34598103e-02 + 4.68042562e-05 1.64235364e-02 5.99725484e-02 1.25824054e-01 + 2.08064456e-01 3.00824090e-01 3.98939205e-01 4.98161316e-01 + 5.95059508e-01 6.86787971e-01 7.70852355e-01 8.44942289e-01 + 9.06852663e-01 9.54487201e-01 9.85925226e-01 9.99530703e-01 + 9.94086558e-01 9.68941138e-01 9.24152715e-01 8.60612554e-01] + + +Computing the position of the moon +---------------------------------- + + >>> from datetime import datetime + >>> from pyorbital import planets + >>> time_t = datetime(2016, 7, 30, 0, 0) + >>> moon = planets.Moon(time_t) + >>> lon = 20.0 + >>> lat = 65.0 + >>> rasc, decl, alt, azi = moon.topocentric_position(lon, lat) + >>> print alt, azi + 6.33389454706 61.6795817556 + +And for an array of longitudes and latitudes: + + >>> import numpy as np + >>> lons = np.arange(100).reshape(10,10) + >>> lats = np.arange(100).reshape(10,10) * 0.9 + >>> rasc, decl, alt, azi = moon.topocentric_position(lons, lats) + >>> print alt.shape + (10, 10) + >>> print alt[4,8] + 14.9564426346 + + +Accuracy +-------- + +In lack of an absolute truth or reference the moon calclations with Pyorbital +have been compared to pyephem. There are indeed deviations, but for the moon +phase they are in general rather small. See image below, where we compare the +moon phase over 416 days starting from December 1st 2015. As seen from the +figure the deviations are within 0.3 %. + + .. image:: _static/moonphase_compare.png + +For moon height and azimuth differences are larger, as seen from the figures +below. We have calculated the position of the moon relative to the City of +Norrköping, Sweden, over four months from March 7, 2012. The deviations are +within 4 degrees for the height, and within 9 degrees for the azimuth. + + + .. image:: _static/moonheight_compare.png + .. image:: _static/moonazimuth_compare.png + + + +.. _`Stjaernhimlen`: http://www.stjarnhimlen.se/comp/ppcomp.html diff --git a/doc/source/satellite_position.rst b/doc/source/satellite_position.rst new file mode 100644 index 00000000..ec5a3a3c --- /dev/null +++ b/doc/source/satellite_position.rst @@ -0,0 +1,17 @@ + +Computing satellite postion +--------------------------- +The orbital module enables computation of satellite position and velocity at a +specific time: + + >>> from pyorbital.orbital import Orbital + >>> from datetime import datetime + >>> orb = Orbital("noaa 18") + >>> now = datetime.utcnow() + >>> # Get normalized position and velocity of the satellite: + >>> orb.get_position(now) + ([0.57529384846822862, 0.77384005228105424, 0.59301408257897559], + [0.031846489698768146, 0.021287993461926374, -0.05854106186659274]) + >>> # Get longitude, latitude and altitude of the satellite: + >>> orb.get_lonlatalt(now) + (-1.1625895579622014, 0.55402132517640568, 847.89381184656702) diff --git a/doc/source/tle.rst b/doc/source/tle.rst new file mode 100644 index 00000000..d365341f --- /dev/null +++ b/doc/source/tle.rst @@ -0,0 +1,12 @@ +TLE files +--------- +Pyorbital has a module for parsing NORAD TLE-files + + >>> from pyorbital import tlefile + >>> tle = tlefile.read('noaa 18', '/path/to/my/tle_file.txt') + >>> tle.inclination + 99.043499999999995 + +If no path is given pyorbital tries to read the earth observation TLE-files +from celestrak.com + diff --git a/pyorbital/moon_phase.py b/pyorbital/moon_phase.py new file mode 100644 index 00000000..fc6c1814 --- /dev/null +++ b/pyorbital/moon_phase.py @@ -0,0 +1,245 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- + +# Copyright (c) 2013, 2016 Adam.Dybbroe + +# Author(s): + +# Adam.Dybbroe + +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. + +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. + +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + +"""Calculations of the phase of the moon +""" + +from numpy import deg2rad, sin, cos, tan, arctan, fabs +import numpy as np +import datetime +import collections + +import math +PI = math.pi +SMALL_FLOAT = 1.0e-12 + + +def Julian(dobj): + """Returns the number of julian days for the specified date-time. + Input = datetime object. + """ + + if isinstance(dobj, datetime.datetime): + dobj = np.array([dobj], 'datetime64[us]') + elif isinstance(dobj, np.ndarray): + dobj = dobj.astype('datetime64[us]') + elif isinstance(dobj, collections.Sequence): + dobj = np.array(dobj, 'datetime64[us]') + + day = ( + dobj - dobj.astype('datetime64[M]')).astype('f') / (24. * 3600. * 1000000) + 1.0 + month = ( + dobj.astype('datetime64[M]') - dobj.astype('datetime64[Y]')).astype('f') + 1 + year = dobj.astype('datetime64[Y]').astype('f') + 1970 + + year = np.where(np.less(month, 3), year - 1, year) + month = np.where(np.less(month, 3), month + 12, month) + + cond1 = np.greater(year, 1582) + cond2 = np.logical_and(np.equal(year, 1582), np.greater(month, 10)) + cond3 = np.logical_and(np.logical_and(np.equal(year, 1582), np.equal(month, 10)), + np.greater(day, 15)) + + a__ = np.divide(year, 100).astype('i') + b__ = np.where(np.logical_or(np.logical_or(cond1, cond2), cond3), + 2 - a__ + a__ / 4, 0) + + c__ = (365.25 * year).astype('i') + e__ = (30.6001 * (month + 1)).astype('i') + return b__ + c__ + e__ + day + 1720994.5 + + +# def xJulian(dobj): +# """Returns the number of julian days for the specified date-time. +# Input = datetime object. +# """ + +# year = dobj.year +# month = dobj.month +# day = dobj.day + (dobj.hour / 24. + dobj.minute / (24 * 60.) + +# dobj.second / (24 * 3600.) + +# dobj.microsecond / (24 * 3600 * 1000000.)) + +# if month < 3: +# year = year - 1 +# month += 12 + +# if (year > 1582 or +# (year == 1582 and month > 10) or +# (year == 1582 and month == 10 and day > 15)): +# a__ = int(year / 100) +# b__ = 2 - a__ + a__ / 4 +# else: +# b__ = 0 + +# c__ = int(365.25 * year) +# e__ = int(30.6001 * (month + 1)) +# return b__ + c__ + e__ + day + 1720994.5 + + +# def xsun_position(jday): +# """Get sun position""" + +# # double n,x,e,l,dl,v; +# # double m2; +# # int i; + +# n__ = 360. / 365.2422 * jday +# i__ = int(n__ / 360.) +# n__ = n__ - i__ * 360.0 +# x__ = n__ - 3.762863 +# if x__ < 0: +# x__ = x__ + 360. + +# x__ = deg2rad(x__) +# e__ = x__ +# while 1: +# dl_ = e__ - .016718 * sin(e__) - x__ +# e__ = e__ - dl_ / (1 - .016718 * cos(e__)) +# if fabs(dl_) < SMALL_FLOAT: +# break + +# v__ = 360. / PI * arctan(1.01686011182 * tan(e__ / 2)) +# sunpos = v__ + 282.596403 +# i__ = int(sunpos / 360.) +# sunpos = sunpos - i__ * 360.0 + +# return sunpos + + +# def xmoon_position(jday, lsun): +# """Get the moon position""" + +# ms_ = 0.985647332099 * jday - 3.762863 +# if ms_ < 0: +# ms_ = ms_ + 360.0 + +# mpos = 13.176396 * jday + 64.975464 +# i__ = int(mpos / 360.0) +# mpos = mpos - i__ * 360.0 +# if mpos < 0: +# mpos = mpos + 360.0 + +# mm_ = mpos - 0.1114041 * jday - 349.383063 +# i__ = int(mm_ / 360.0) +# mm_ = mm_ - i__ * 360.0 + +# ev_ = 1.2739 * sin(deg2rad(2 * (mpos - lsun) - mm_)) +# sms = sin(deg2rad(ms_)) +# ae_ = 0.1858 * sms +# mm_ = mm_ + ev_ - ae_ - 0.37 * sms +# ec_ = 6.2886 * sin(deg2rad(mm_)) +# mpos = mpos + ev_ + ec_ - ae_ + 0.214 * sin(deg2rad(2 * mm_)) +# mpos = 0.6583 * sin(deg2rad(2 * (mpos - lsun))) + mpos + +# return mpos + + +def sun_position(jday): + """Get sun position""" + + # double n,x,e,l,dl,v; + # double m2; + # int i; + + if isinstance(jday, collections.Sequence): + jday = np.array(jday, 'float64') + elif not isinstance(jday, np.ndarray): + jday = np.array([jday], 'float64') + + n__ = 360. / 365.2422 * jday + i__ = np.divide(n__, 360.).astype('i') + n__ = n__ - i__ * 360.0 + x__ = n__ - 3.762863 + x__ = np.where(np.less(x__, 0), x__ + 360., x__) + + x__ = deg2rad(x__) + e__ = x__ + while 1: + dl_ = e__ - .016718 * sin(e__) - x__ + e__ = e__ - dl_ / (1 - .016718 * cos(e__)) + if np.alltrue(np.less(np.fabs(dl_), SMALL_FLOAT)): + break + + v__ = 360. / PI * arctan(1.01686011182 * tan(e__ / 2)) + sunpos = v__ + 282.596403 + i__ = np.divide(sunpos, 360.).astype('i') + sunpos = sunpos - i__ * 360.0 + + if isinstance(sunpos, np.ndarray) and len(sunpos) == 1: + return sunpos[0] + else: + return sunpos + + return sunpos + + +def moon_position(jday, lsun): + """Get the moon position""" + + if isinstance(jday, collections.Sequence): + jday = np.array(jday, 'float64') + elif not isinstance(jday, np.ndarray): + jday = np.array([jday], 'float64') + + if isinstance(lsun, collections.Sequence): + lsun = np.array(lsun, 'float64') + elif not isinstance(lsun, np.ndarray): + lsun = np.array([lsun], 'float64') + + ms_ = 0.985647332099 * jday - 3.762863 + ms_ = np.where(np.less(ms_, 0), ms_ + 360.0, ms_) + + mpos = 13.176396 * jday + 64.975464 + i__ = np.divide(mpos, 360.0).astype('i') + mpos = mpos - i__ * 360.0 + mpos = np.where(np.less(mpos, 0), mpos + 360.0, mpos) + + mm_ = mpos - 0.1114041 * jday - 349.383063 + i__ = np.divide(mm_, 360.0).astype('i') + mm_ = mm_ - i__ * 360.0 + + ev_ = 1.2739 * sin(deg2rad(2 * (mpos - lsun) - mm_)) + sms = sin(deg2rad(ms_)) + ae_ = 0.1858 * sms + mm_ = mm_ + ev_ - ae_ - 0.37 * sms + ec_ = 6.2886 * sin(deg2rad(mm_)) + mpos = mpos + ev_ + ec_ - ae_ + 0.214 * sin(deg2rad(2 * mm_)) + mpos = 0.6583 * sin(deg2rad(2 * (mpos - lsun))) + mpos + + if isinstance(mpos, np.ndarray) and len(mpos) == 1: + return mpos[0] + else: + return mpos + + +def moon_phase(dobj): + """Calculate the phase of the moon""" + + jday = Julian(dobj) - 2444238.5 + lsun = sun_position(jday) + lmoon = moon_position(jday, lsun) + + phase = (1.0 - cos(deg2rad(lmoon - lsun))) / 2.0 + if isinstance(phase, np.ndarray) and len(phase) == 1: + phase = phase[0] + return phase diff --git a/pyorbital/phase b/pyorbital/phase new file mode 100755 index 00000000..368fe443 Binary files /dev/null and b/pyorbital/phase differ diff --git a/pyorbital/phase.c b/pyorbital/phase.c new file mode 100644 index 00000000..b3d7ba07 --- /dev/null +++ b/pyorbital/phase.c @@ -0,0 +1,249 @@ +/* standalone moon phase calculation */ + +#include +#include + +#define PI 3.1415926535897932384626433832795 +#define RAD (PI/180.0) +#define SMALL_FLOAT (1e-12) + +typedef struct { + int year,month,day; + double hour; +} TimePlace; + +void JulianToDate(TimePlace* now, double jd) +{ + long jdi, b; + long c,d,e,g,g1; + + jd += 0.5; + jdi = jd; + printf("jd jdi: %f %d\n", jd, jdi); + if (jdi > 2299160) { + long a = (jdi - 1867216.25)/36524.25; + b = jdi + 1 + a - a/4; + } + else b = jdi; + + printf("b: %d\n", b); + + c = b + 1524; + d = (c - 122.1)/365.25; + e = 365.25 * d; + g = (c - e)/30.6001; + g1 = 30.6001 * g; + now->day = c - e - g1; + now->hour = (jd - jdi) * 24.0; + if (g <= 13) now->month = g - 1; + else now->month = g - 13; + if (now->month > 2) now->year = d - 4716; + else now->year = d - 4715; +} + +double +Julian(int year,int month,double day) +{ + /* + Returns the number of julian days for the specified day. + */ + + int a,b,c,e; + if (month < 3) { + year--; + month += 12; + } + if (year > 1582 || (year == 1582 && month>10) || + (year == 1582 && month==10 && day > 15)) { + a=year/100; + b=2-a+a/4; + } + + c = 365.25*year; + e = 30.6001*(month+1); + return b+c+e+day+1720994.5; +} + +double sun_position(double j) +{ + double n,x,e,l,dl,v; + double m2; + int i; + + n=360/365.2422*j; + i=n/360; + n=n-i*360.0; + x=n-3.762863; + if (x<0) x += 360; + x *= RAD; + e=x; + do { + dl=e-.016718*sin(e)-x; + e=e-dl/(1-.016718*cos(e)); + } while (fabs(dl)>=SMALL_FLOAT); + v=360/PI*atan(1.01686011182*tan(e/2)); + l=v+282.596403; + i=l/360; + l=l-i*360.0; + return l; +} + +double moon_position(double j, double ls) +{ + + double ms,l,mm,n,ev,sms,z,x,lm,bm,ae,ec; + double d; + double ds, as, dm; + int i; + + /* ls = sun_position(j) */ + ms = 0.985647332099*j - 3.762863; + if (ms < 0) ms += 360.0; + l = 13.176396*j + 64.975464; + i = l/360; + l = l - i*360.0; + if (l < 0) l += 360.0; + mm = l-0.1114041*j-349.383063; + i = mm/360; + mm -= i*360.0; + /*n = 151.950429 - 0.0529539*j; + i = n/360; + n -= i*360.0; + */ + ev = 1.2739*sin((2*(l-ls)-mm)*RAD); + sms = sin(ms*RAD); + ae = 0.1858*sms; + mm += ev-ae- 0.37*sms; + ec = 6.2886*sin(mm*RAD); + l += ev+ec-ae+ 0.214*sin(2*mm*RAD); + l= 0.6583*sin(2*(l-ls)*RAD)+l; + return l; +} + +double moon_phase(int year,int month,int day, double hour, int* ip) +{ + /* + Calculates more accurately than Moon_phase , the phase of the moon at + the given epoch. + returns the moon phase as a real number (0-1) + */ + + double j= Julian(year,month,(double)day+hour/24.0)-2444238.5; + double ls = sun_position(j); + double lm = moon_position(j, ls); + + double t = lm - ls; + if (t < 0) t += 360; + *ip = (int)((t + 22.5)/45) & 0x7; + return (1.0 - cos((lm - ls)*RAD))/2; +} + +static void nextDay(int* y, int* m, int* d, double dd) +{ + TimePlace tp; + double jd = Julian(*y, *m, (double)*d); + + jd += dd; + JulianToDate(&tp, jd); + + *y = tp.year; + *m = tp.month; + *d = tp.day; +} + +main() +{ + int y, m, d; + int m0; + int h; + int i; + double step = 1; + int begun = 0; + + double pmax = 0; + double pmin = 1; + int ymax, mmax, dmax, hmax; + int ymin, mmin, dmin, hmin; + double plast = 0; + + printf("tabulation of the phase of the moon for one month\n\n"); + + printf("year: "); fflush(stdout); + scanf("%d", &y); + + printf("month: "); fflush(stdout); + scanf("%d", &m); + + d = 1; + m0 = m; + + /* + { + TimePlace tp; + int ip; + double p; + double jd = Julian(y, m, (double)d); + printf("Julian day = %f\n", jd); + + JulianToDate(&tp, jd); + printf("Year, Month, Day, Hour: %d %d %d %f\n", + tp.year, tp.month, tp.day, tp.hour); + + p = moon_phase(tp.year, tp.month, tp.day, tp.hour, &ip); + printf("Moon phase = %f\n", p); + + return 0; + } + */ + + printf("\nDate Time Phase Segment\n"); + for (;;) { + double p; + int ip; + + for (h = 0; h < 24; h += step) { + + p = moon_phase(y, m, d, h, &ip); + + if (begun) { + if (p > plast && p > pmax) { + pmax = p; + ymax = y; + mmax = m; + dmax = d; + hmax = h; + } + else if (pmax) { + printf("%04d/%02d/%02d %02d:00 (fullest)\n", + ymax, mmax, dmax, hmax); + pmax = 0; + } + + if (p < plast && p < pmin) { + pmin = p; + ymin = y; + mmin = m; + dmin = d; + hmin = h; + } + else if (pmin < 1) { + printf("%04d/%02d/%02d %02d:00 (newest)\n", + ymin, mmin, dmin, hmin); + pmin = 1.0; + } + + if (h == 16) { + printf("%04d/%02d/%02d %02d:00 %5.1f%% (%d)\n", + y, m, d, h, floor(p*1000+0.5)/10, ip); + } + } + else begun = 1; + + plast = p; + + } + nextDay(&y, &m, &d, 1.0); + if (m != m0) break; + } + return 0; +} diff --git a/pyorbital/planets.py b/pyorbital/planets.py new file mode 100644 index 00000000..9a89f324 --- /dev/null +++ b/pyorbital/planets.py @@ -0,0 +1,514 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- + +# Copyright (c) 2011, 2012 + +# Author(s): + +# Adam Dybbroe + +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. + +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. + +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + +"""Astronomy module. +Parts taken from http://www.stjarnhimlen.se/comp/ppcomp.html +Positions of the sun, the moon, the earth and the other planets +""" + +import datetime +import numpy as np + + +class Moon(object): + def __init__(self, dtobj): + # Orbital elements: + self.orbelem = {"N": 0.0, "i": 0.0, "w": 0.0, + "a": 0.0, "e": 0.0, "M": 0.0, + "E": 0.0, 'v': 0.0} + self.datetime = dtobj + self.day_number = day_number(dtobj) + #self.lonsun = 0.0 + self.hour_angle = 0.0 + self.right_ascension = 0.0 + self.declination = 0.0 + + # Pertubations: + self.dlon = 0.0 + self.dlat = 0.0 + self.drad = 0.0 + self.pertubations() + + self.get_lonsun() + self.rdist = 0.0 + self.ecl = obliquity_of_ecliptic(self.day_number) + + # Get the orbital elements of the Sun + orbelem = orbital_elem(self.day_number, body='moon') + + orbelem = eccentric_anomaly(orbelem) + self.rdist, orbelem = distance_and_true_anomaly(orbelem) + + self.orbelem = orbelem + pos = position(orbelem, self.rdist) + self.geocentric_position = pos + # Apply the pertubations: + self.rdist = self.rdist + self.drad + + xg_ = self.geocentric_position['xh'] + yg_ = self.geocentric_position['yh'] + zg_ = self.geocentric_position['zh'] + + # Cartesian to polar coordinates: + lon = np.rad2deg(np.arctan2(yg_, xg_)) + self.dlon + lat = np.rad2deg(np.arctan2(zg_, np.sqrt( + xg_ * xg_ + yg_ * yg_))) + self.dlat + + # Polar to cartesian coordinates: + xg_ = self.rdist * np.cos(np.deg2rad(lon)) + yg_ = self.rdist * np.sin(np.deg2rad(lon)) + zg_ = self.rdist * np.sin(np.deg2rad(lat)) + self.geocentric_position['xh'] = xg_ + self.geocentric_position['yh'] = yg_ + self.geocentric_position['zh'] = zg_ + + # The Moon's parallax, i.e. the apparent size of the (equatorial) + # radius of the Earth, as seen from the Moon: + self.parallax = np.rad2deg(np.arcsin(1. / self.rdist)) + + right_ascension, declination, gdist = ecliptic_to_equatorial(xg_, + yg_, + zg_, + self.ecl) + self.right_ascension = right_ascension + self.declination = declination + + def topocentric_position(self, longitude, latitude): + """The Moon's position, as computed earlier, is geocentric, i.e. as + seen by an imaginary observer at the center of the Earth. Real + observers dwell on the surface of the Earth, though, and they will see + a different position - the topocentric position. This position can + differ by more than one degree from the geocentric position. To compute + the topocentric positions, we must add a correction to the geocentric + position. + """ + + azi, alt_geoc, ha_ = azimuthal_coord(self.datetime, + self.lonsun, + longitude, + latitude, + self.right_ascension, + self.declination) + #print "azi, alt_geoc, ha_: ", azi, alt_geoc, ha_ + + self.hour_angle = ha_ + alt_topoc = alt_geoc - self.parallax * np.cos(np.deg2rad(alt_geoc)) + + # Account for the flattening of the Earth: + gclat = latitude - 0.1924 * np.sin(np.deg2rad(2*latitude)) + rho = 0.99833 + 0.00167 * np.cos(np.deg2rad(2*latitude)) + + g__ = np.arctan(np.tan(np.deg2rad(gclat)) / + np.cos(np.deg2rad(self.hour_angle))) + g__ = np.rad2deg(g__) + + # Now we're ready to convert the geocentric Right Ascention and + # Declination (RA, Decl) to their topocentric values (topRA, topDecl): + + topRA = (self.right_ascension - + self.parallax * rho * np.cos(np.deg2rad(gclat)) * + np.sin(np.deg2rad(ha_)) / + np.cos(np.deg2rad(self.declination)) + ) + topDecl = (self.declination - + self.parallax * rho * np.sin(np.deg2rad(gclat)) * + np.sin(np.deg2rad(g__ - self.declination)) / + np.sin(np.deg2rad(g__)) + ) + + return topRA, topDecl, alt_topoc, azi + + def get_lonsun(self): + """Calculate the lonsun""" + # Ls = Ms + ws Mean Longitude of the Sun (Ns=0) + sun = Sun(self.datetime) + self.lonsun = sun.lonsun + + def pertubations(self): + """Pertubations of the moon, to be added to + its lon,lat and distance""" + + # Ms, Mm Mean Anomaly of the Sun and the Moon + # Nm Longitude of the Moon's node + # ws, wm Argument of perihelion for the Sun and the Moon + # Ls = Ms + ws Mean Longitude of the Sun (Ns=0) + # Lm = Mm + wm + Nm Mean longitude of the Moon + # D = Lm - Ls Mean elongation of the Moon + # F = Lm - Nm Argument of latitude for the Moon + + sun = Sun(self.datetime) + Ms_ = sun.orbelem['M'] + Mm_ = self.orbelem['M'] + Nm_ = self.orbelem['N'] + wm_ = self.orbelem['w'] + Lm_ = Mm_ + wm_ + Nm_ + D__ = Lm_ - sun.lonsun + F__ = Lm_ - Nm_ + + # Add these terms to the Moon's longitude (degrees): + dlon = (-1.274 * np.sin(np.deg2rad(Mm_ - 2*D__)) # (the Evection) + + 0.658 * np.sin(np.deg2rad(2*D__)) # (the Variation) + # (the Yearly Equation) + - 0.186 * np.sin(np.deg2rad(Ms_)) + - 0.059 * np.sin(np.deg2rad(2*Mm_ - 2*D__)) + - 0.057 * np.sin(np.deg2rad(Mm_ - 2*D__ + Ms_)) + + 0.053 * np.sin(np.deg2rad(Mm_ + 2*D__)) + + 0.046 * np.sin(np.deg2rad(2*D__ - Ms_)) + + 0.041 * np.sin(np.deg2rad(Mm_ - Ms_)) + # (the Parallactic Equation) + - 0.035 * np.sin(np.deg2rad(D__)) + - 0.031 * np.sin(np.deg2rad(Mm_ + Ms_)) + - 0.015 * np.sin(np.deg2rad(2*F__ - 2*D__)) + + 0.011 * np.sin(np.deg2rad(Mm_ - 4*D__)) + ) + self.dlon = dlon + + # Add these terms to the Moon's latitude (degrees): + dlat = (-0.173 * np.sin(np.deg2rad(F__ - 2*D__)) + - 0.055 * np.sin(np.deg2rad(Mm_ - F__ - 2*D__)) + - 0.046 * np.sin(np.deg2rad(Mm_ + F__ - 2*D__)) + + 0.033 * np.sin(np.deg2rad(F__ + 2*D__)) + + 0.017 * np.sin(np.deg2rad(2*Mm_ + F__)) + ) + self.dlat = dlat + + # Add these terms to the Moon's distance (Earth radii): + drad = (-0.58 * np.cos(np.deg2rad(Mm_ - 2*D__)) + - 0.46 * np.cos(np.deg2rad(2*D__)) + ) + self.drad = drad + + +class Sun(object): + def __init__(self, dtobj): + # Orbital elements: + self.orbelem = {"N": 0.0, "i": 0.0, "w": 0.0, + "a": 1.0, "e": 0.0, "M": 0.0, + "E": 0.0, 'v': 0.0} + self.datetime = dtobj + self.day_number = day_number(dtobj) + self.lonsun = 0.0 + self.rdist = 0.0 + self.ecl = obliquity_of_ecliptic(self.day_number) + + # Get the orbital elements of the Sun + orbelem = orbital_elem(self.day_number, body='sun') + + orbelem = sun_eccentric_anomaly(orbelem) + self.rdist, orbelem = distance_and_true_anomaly(orbelem) + + self.orbelem = orbelem + + # True longitude: + self.lonsun = self.orbelem['v'] + self.orbelem['w'] + + self.position = {} + + def get_position(self): + """Calculate the position of the Sun. + Equatorial, rectangular, geocentric coordinates""" + + # Convert lonsun,r to ecliptic rectangular geocentric coordinates xs,ys: + xs_ = self.rdist * np.cos(np.deg2rad(self.lonsun)) + ys_ = self.rdist * np.sin(np.deg2rad(self.lonsun)) + # To convert this to equatorial, rectangular, geocentric coordinates, compute: + xe_ = xs_ + ye_ = ys_ * np.cos(np.deg2rad(self.ecl)) + ze_ = ys_ * np.sin(np.deg2rad(self.ecl)) + + # Compute the planet's Right Ascension (RA) and Declination (Dec): + ra_ = np.rad2deg(np.arctan2(ye_, xe_)) + dec = np.rad2deg(np.arctan2(ze_, np.sqrt(xe_ * xe_ + ye_ * ye_))) + + self.position = {'xe': xe_, 'ye': ye_, 'ze': ze_, + 'right_ascension': ra_, 'declination': dec} + + +def day_number(dtobj): + """Get the 'day number' from the date, given by a datetime object + """ + # d = 367*y - 7 * ( y + (m+9)/12 ) / 4 + 275*m/9 + D - 730530 + # dnumber = (367 * dtobj.year - + # 7 * ( dtobj.year + (dtobj.month + 9) / 12 ) / 4 + + # 275 * dtobj.month / 9 + dtobj.day - 730530) + # return dnumber + dtobj.hour/24.0 + + origo = datetime.datetime(1999, 12, 31, 0) + delta_t = dtobj - origo + return delta_t.days + (delta_t.seconds + + delta_t.microseconds/1000000.)/(3600. * 24.) + + +def _gmst(lonsun, utc): + """The Greenwich Mean Sideral Time (GMST) is the LST at Greenwich""" + # GMST0 = 15 * (Ls + 180_degrees) + # GMST = GMST0 + UT + # LST = GMST + local_longitude/15 + gmst0 = ((lonsun + 180.0) % 360) / 15.0 + ut_ = (utc.hour + utc.minute / 60.0 + + (utc.second + utc.microsecond / 1000000.) / 3600.) + return gmst0 + ut_ + + +def _lst(gmst, local_lon): + """The Local Sideral Time (LST). + This is simply the RA of your local meridian""" + return gmst + local_lon / 15.0 + + +def _hour_angle(utc_time, local_lon, lonsun, right_ascension): + """Hour angle at *utc_time* for the given *local_lon* and + *right_ascension*. gmst and lst in hours. the hour angle HA is given in + degrees: + """ + gmst = _gmst(lonsun, utc_time) + return (15 * _lst(gmst, local_lon) - right_ascension) % 360 + + +def azimuthal_coord(utc_time, lonsun, local_lon, local_lat, + right_ascension, declination): + """Calculate the azimuthal coordinates (azimuth and altitude)""" + + ha_ = _hour_angle(utc_time, local_lon, lonsun, right_ascension) + + x__ = np.cos(np.deg2rad(ha_)) * np.cos(np.deg2rad(declination)) + y__ = np.sin(np.deg2rad(ha_)) * np.cos(np.deg2rad(declination)) + z__ = np.sin(np.deg2rad(declination)) + + xhor = (x__ * np.sin(np.deg2rad(local_lat)) - + z__ * np.cos(np.deg2rad(local_lat))) + yhor = y__ + zhor = (x__ * np.cos(np.deg2rad(local_lat)) + + z__ * np.sin(np.deg2rad(local_lat))) + + azi = np.rad2deg(np.arctan2(yhor, xhor)) + 180.0 + # = atan2( zhor, sqrt(xhor*xhor+yhor*yhor) ) + alt = np.rad2deg(np.arcsin(zhor)) + + return azi, alt, ha_ + + +def obliquity_of_ecliptic(dnum): + """compute the obliquity of the ecliptic""" + return 23.4393 - 3.563E-7 * dnum + + +def orbital_elem(dnum, **options): + """Orbital elements of the Sun, the Moon, and the planets""" + + if "body" in options and options["body"] == 'moon': + N__ = 125.1228 - 0.0529538083 * dnum + i__ = 5.1454 + w__ = 318.0634 + 0.1643573223 * dnum + a__ = 60.2666 # (Earth radii) + e__ = 0.054900 + M__ = 115.3654 + 13.0649929509 * dnum + elif "body" in options and options["body"] == 'sun': + N__ = 0.0 + i__ = 0.0 + w__ = 282.9404 + 4.70935E-5 * dnum + a__ = 1.000000 # (AU) + e__ = 0.016709 - 1.151E-9 * dnum + M__ = 356.0470 + 0.9856002585 * dnum + else: + raise IOError("Celestial body not specified or not supported yet!") + + if M__ > 360: + M__ = M__ - 360 * (int(M__)/360) + if N__ > 360: + N__ = N__ - 360 * (int(N__)/360) + + return {"N": N__, "i": i__, "w": w__, + "a": a__, "e": e__, "M": M__} + + +def sun_eccentric_anomaly(orbelem): + """Computation of the eccentric anomaly from the mean anomaly + for ths Sun""" + # Note that the formulae for computing E are not exact; + # however they're accurate enough here. + E__ = (orbelem["M"] + + np.rad2deg(orbelem["e"]) * + np.sin(np.deg2rad(orbelem["M"])) * + (1.0 + orbelem["e"] * np.cos(np.deg2rad(orbelem["M"])))) + + orbelem["E"] = E__ + return orbelem + + +def eccentric_anomaly(orbelem): + """Computation of the eccentric anomaly from the mean anomaly""" + + E__ = (orbelem["M"] + + np.rad2deg(orbelem["e"]) * + np.sin(np.deg2rad(orbelem["M"])) * + (1.0 + orbelem["e"] * np.cos(np.deg2rad(orbelem["M"])))) + + E0_ = E__ + epsilon = 0.001 + delta_e = 2 * epsilon + niter = 0 + while delta_e > epsilon and niter < 10: + #print "E0_ = ",E0_ + E1_ = (E0_ - + (E0_ - np.rad2deg(orbelem["e"]) * np.sin(np.deg2rad(E0_)) + - orbelem["M"]) / + (1 - orbelem["e"] * np.cos(np.deg2rad(E0_)))) + delta_e = abs(E1_ - E0_) + #print "Iteration = %d, dev = %f" % (niter, delta_e) + niter = niter + 1 + E0_ = E1_ + + orbelem["E"] = E1_ + return orbelem + + +def distance_and_true_anomaly(orbelem): + """Computation of the planets distance 'r' and its true anomaly 'v'""" + # xv = r * cos(v) + # yv = r * sin(v) + xv_ = orbelem["a"] * (np.cos(np.deg2rad(orbelem["E"])) - orbelem["e"]) + yv_ = orbelem["a"] * (np.sqrt(1.0 - orbelem["e"] * orbelem["e"]) * + np.sin(np.deg2rad(orbelem["E"]))) + + v__ = np.rad2deg(np.arctan2(yv_, xv_)) + r__ = np.sqrt(xv_ * xv_ + yv_ * yv_) + orbelem["v"] = v__ + + return r__, orbelem + + +def position(orbelem, radius): + """Computation of the planet's position in 3-dimensional space. For the + Moon, this is the geocentric (Earth-centered) position in the ecliptic + coordinate system. For the planets, this is the heliocentric (Sun-centered) + position, also in the ecliptic coordinate system. + """ + + xh_ = radius * (np.cos(np.deg2rad(orbelem["N"])) * + np.cos(np.deg2rad(orbelem["v"] + orbelem["w"])) - + np.sin(np.deg2rad(orbelem["N"])) * + np.sin(np.deg2rad(orbelem["v"] + orbelem["w"])) * + np.cos(np.deg2rad(orbelem["i"]))) + yh_ = radius * (np.sin(np.deg2rad(orbelem["N"])) * + np.cos(np.deg2rad(orbelem["v"] + orbelem["w"])) + + np.cos(np.deg2rad(orbelem["N"])) * + np.sin(np.deg2rad(orbelem["v"] + orbelem["w"])) * + np.cos(np.deg2rad(orbelem["i"]))) + zh_ = radius * (np.sin(np.deg2rad(orbelem["v"] + orbelem["w"])) * + np.sin(np.deg2rad(orbelem["i"]))) + + # Compute the ecliptic longitude and latitude + lonecl = np.rad2deg(np.arctan2(yh_, xh_)) + latecl = np.rad2deg(np.arctan2(zh_, np.sqrt(xh_ * xh_ + yh_ * yh_))) + check = np.sqrt(xh_ * xh_ + yh_ * yh_ + zh_ * zh_) + + #print "sqrt(xh*xh+yh*yh+zh*zh) = r: ", check, radius + return {'xh': xh_, 'yh': yh_, 'zh': zh_, + 'lonecl': lonecl, 'latecl': latecl} + + +def heliocentric_to_geocentric(lonecl, latecl, radius, lonsun, rsun): + """Conversion from heliocentric to geocentric coordinates. (Not needed for + the moon). + Now we have computed the heliocentric (Sun-centered) coordinate + of the planet, and we have included the most important perturbations. We + want to compute the geocentric (Earth-centerd) position. We should convert + the perturbed lonecl, latecl, r to (perturbed) xh, yh, zh + """ + xh_ = radius * np.cos(lonecl) * np.cos(np.deg2rad(latecl)) + yh_ = radius * np.sin(lonecl) * np.cos(np.deg2rad(latecl)) + zh_ = radius * np.sin(np.deg2rad(latecl)) + + # If we are computing the Moon's position, this is already the geocentric + # position, and thus we simply set xg=xh, yg=yh, zg=zh. Otherwise we must + # also compute the Sun's position: convert lonsun, rs (where rs is the r + # computed here) to xs, ys: + + xs_ = rsun * np.cos(np.deg2rad(lonsun)) + ys_ = rsun * np.sin(np.deg2rad(lonsun)) + + xg_ = xh_ + xs_ + yg_ = yh_ + ys_ + zg_ = zh_ + + return {'xg': xg_, 'yg': yg_, 'zg': zg_} + + +def ecliptic_to_equatorial(xg_, yg_, zg_, ecl): + """Convert the rectangular, ecliptic coordinates to rectangular, + equatorial coordinates: Simply rotate the y-z-plane by ecl, the + angle of the obliquity of the ecliptic""" + + xe_ = xg_ + ye_ = (yg_ * np.cos(np.deg2rad(ecl)) - + zg_ * np.sin(np.deg2rad(ecl))) + ze_ = (yg_ * np.sin(np.deg2rad(ecl)) + + zg_ * np.cos(np.deg2rad(ecl))) + + # Compute the planet's Right Ascension (RA) and Declination (Dec): + ra_ = np.rad2deg(np.arctan2(ye_, xe_)) + dec = np.rad2deg(np.arctan2(ze_, np.sqrt(xe_ * xe_ + ye_ * ye_))) + + # Compute the geocentric distance: + # = sqrt(xe*xe+ye*ye+ze*ze) + rg_ = np.sqrt(xg_ * xg_ + yg_ * yg_ + zg_ * zg_) + + return ra_, dec, rg_ + + +if __name__ == "__main__": + now = datetime.datetime.utcnow() + + #dNum = day_number(now) + #orb = orbital_elem(dNum, body='moon') + #orb = eccentric_anomaly(orb) + #rdist, orb = distance_and_true_anomaly(orb) + + #pos = position(orb, rdist) + #xgeo, ygeo, zgeo = pos['xh'], pos['yh'], pos['zh'] + + #ecliptic = obliquity_of_ecliptic(dNum) + + #tup = ecliptic_to_equatorial(xgeo, ygeo, zgeo, ecliptic) + + # Norrköping: Lat N 58° 35′ 15″ Lon E 16° 11′ 15″ + llon = 16 + 11.0/60. + 15./3600. + llat = 58. + 35./60. + 15./3600. + #moon = Moon(now) + #rasc, decl, alt, azi = moon.topocentric_position(17., 58.0) + currtime = datetime.datetime(2012, 1, 7, 12, 0) + endtime = datetime.datetime(2012, 1, 7, 13, 0) + delta_t = datetime.timedelta(seconds=120) + heights = [] + dtimes = [] + while currtime < endtime: + moon = Moon(currtime) + rasc, decl, alt, azi = moon.topocentric_position(llon, llat) + heights.append(alt) + dtimes.append(currtime) + currtime = currtime + delta_t + break + + llon = np.arange(100)/10.0 + 10.0 + llat = np.arange(100)/10.0 + 50.0 + rasc, decl, alt, azi = moon.topocentric_position(llon, llat) diff --git a/pyorbital/pyephem_compare/ephem_compare_moon_phase.py b/pyorbital/pyephem_compare/ephem_compare_moon_phase.py new file mode 100644 index 00000000..c93e608f --- /dev/null +++ b/pyorbital/pyephem_compare/ephem_compare_moon_phase.py @@ -0,0 +1,77 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- + +# Copyright (c) 2015, 2016 Adam.Dybbroe + +# Author(s): + +# Adam.Dybbroe + +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. + +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. + +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + +"""Compare accuracy of the moon_phase function with the ephem library""" + +import matplotlib.pyplot as plt +import ephem +import datetime +from pyorbital.moon_phase import moon_phase +import numpy as np + +start_time = datetime.datetime(2015, 12, 1, 12) +delta_t = datetime.timedelta(hours=1) + +norrk = ephem.Observer() +norrk.lat, norrk.lon = '58.5875', '16.1875' +#norrk.lat, norrk.lon = '0.0', '0.0' +norrk.pressure = 0 +time_t = start_time + +tlist = [] +ephem_pha = [] +pha = [] +for idx in range(10000): + norrk.date = time_t.strftime('%Y/%m/%d %H:%M') + m = ephem.Moon(norrk) + ephem_phase = m.phase + phase = moon_phase(time_t) + + #diff = 100*phase - ephem_phase + tlist.append(time_t) + # print time_t.strftime("%Y-%m-%d %H:%M"), diff + ephem_pha.append(ephem_phase) + pha.append(phase) + time_t = time_t + delta_t + + +# for item in pha: +# print "%f, " % item, +# print + +pha = np.array(pha) * 100 + +ephem_pha = np.array(ephem_pha) + +# Plot the moon phases to compare: + +fig = plt.figure() + +ax = fig.add_subplot(111) +ax.scatter(ephem_pha, ephem_pha - pha) + +ax.set_xlabel('Ephem moon phase (%)', fontsize=20) +ax.set_ylabel('Difference - Ephem-pyorbital (%)', fontsize=20) +ax.set_title('Comparing Moon phase derivation') +ax.grid(True) + +plt.savefig('./moonphase_compare.png') diff --git a/pyorbital/pyephem_compare/ephem_compare_moon_position.py b/pyorbital/pyephem_compare/ephem_compare_moon_position.py new file mode 100644 index 00000000..bb8724d1 --- /dev/null +++ b/pyorbital/pyephem_compare/ephem_compare_moon_position.py @@ -0,0 +1,102 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- + +# Copyright (c) 2015, 2016 Adam.Dybbroe + +# Author(s): + +# Adam.Dybbroe + +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. + +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. + +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + + +"""Compare my derivation of the altitude and azimuth of the moon with that +derived from ephem + +""" + +import ephem +import datetime +import numpy as np +from pyorbital import planets + +if __name__ == "__main__": + + # Norrköping: Lat N 58° 35′ 15″ Lon E 16° 11′ 15″ + llon = 16 + 11.0 / 60. + 15. / 3600. + llat = 58. + 35. / 60. + 15. / 3600. + + currtime = datetime.datetime(2012, 3, 7, 12, 0) + endtime = datetime.datetime(2012, 7, 7, 12, 0) + delta_t = datetime.timedelta(seconds=600) + heights = [] + dtimes = [] + + norrk = ephem.Observer() + norrk.lat, norrk.lon = '58.5875', '16.1875' + norrk.pressure = 0 + #norrk.horizon = '-0:34' + + dtimes = [] + ephem_heights = [] + heights = [] + ephem_azimuths = [] + azimuths = [] + while currtime < endtime: + norrk.date = currtime.strftime('%Y/%m/%d %H:%M') + m = ephem.Moon(norrk) + ephem_heights.append(np.rad2deg(m.alt)) + ephem_azimuths.append(np.rad2deg(m.az)) + + moon = planets.Moon(currtime) + rasc, decl, alt, azi = moon.topocentric_position(llon, llat) + # azi, alt, ha_ = planets.azimuthal_coord(currtime, + # moon.lonsun, + # llon, llat, rasc, decl) + heights.append(alt) + azimuths.append(azi) + dtimes.append(currtime) + currtime = currtime + delta_t + + heights = np.array(heights) + ephem_heights = np.array(ephem_heights) + azimuths = np.array(azimuths) + ephem_azimuths = np.array(ephem_azimuths) + + # Plot the moon positions (altitude above horizon): + import matplotlib.pyplot as plt + + fig = plt.figure() + + ax = fig.add_subplot(111) + ax.scatter(ephem_heights, ephem_heights - heights) + + ax.set_xlabel('Ephem moon height (deg)', fontsize=20) + ax.set_ylabel('Difference - Ephem-pyorbital (deg)', fontsize=20) + ax.set_title('Comparing Moon height derivation') + ax.grid(True) + + plt.savefig('./moonheight_compare.png') + + fig = plt.figure() + ax = fig.add_subplot(111) + x = np.mod((ephem_azimuths - azimuths), 360) + ax.scatter(ephem_azimuths, np.where(np.greater(x, 180), x - 360, x)) + + ax.set_xlabel('Ephem moon azimuth (deg)', fontsize=20) + ax.set_ylabel('Difference - Ephem-pyorbital (deg)', fontsize=20) + ax.set_title('Comparing Moon azimuth derivation') + ax.grid(True) + + plt.savefig('./moonazimuth_compare.png') diff --git a/pyorbital/pyephem_compare/mooncalcorg_comparison.py b/pyorbital/pyephem_compare/mooncalcorg_comparison.py new file mode 100644 index 00000000..56eb4ab2 --- /dev/null +++ b/pyorbital/pyephem_compare/mooncalcorg_comparison.py @@ -0,0 +1,73 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- + +# Copyright (c) 2015 Adam.Dybbroe + +# Author(s): + +# Adam.Dybbroe + +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. + +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. + +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + +"""Compare pyorbital and ephem to data found at mooncalc.org +""" + +import ephem +from datetime import datetime +from pyorbital.moon_phase import moon_phase +from pyorbital import planets + +# lon = 16.183 +# lat = 58.6 + +# From mooncalc.org: +# Moon rising: 11:40 +# Moon peak: 17:07 +# Moon set: 22:47 +# Moon distance: 367697km +# Moon elevation angle: -21.47° +# Moon horizontal angle: 288.28° +# Phase = 31.6% + +lon = 10.39307 +lat = 57.07657 + +# Moon elevation angle: -23.09° +# Moon horizontal angle: 289.98° +# Phase = 31.9% + +#currtime = datetime(2015, 12, 16, 23, 16) +currtime = datetime(2015, 12, 16, 23, 50) + +# Pyorbital: +phase = moon_phase(currtime) +moon = planets.Moon(currtime) +rasc, decl, alt, azi = moon.topocentric_position(lon, lat) +print("pyorbital: phase=%f alt=%f azi=%f" % (phase * 100, alt, azi)) + +# Ephem: +norrk = ephem.Observer() +# norrk.lat, norrk.lon = '58.6', '16.183' +norrk.lat, norrk.lon = '57.07657', '10.39307' +norrk.pressure = 0 +norrk.date = currtime.strftime('%Y/%m/%d %H:%M') +m = ephem.Moon(norrk) +print("ephem: phase=%f alt=%s azi=%s" % + (m.phase, str(m.alt), str(m.az))) + +# pyorbital: phase=31.638283 alt=-21.904673 azi=288.377773 +# ephem: phase=31.796560 alt=-21:32:21.0 azi=288:19:25.5 + +# pyorbital: phase=31.884964 alt=-23.521642 azi=290.057501 +# ephem: phase=32.042305 alt=-23:09:50.4 azi=290:01:13.8 diff --git a/pyorbital/tests/__init__.py b/pyorbital/tests/__init__.py index fb67767b..9ff57b89 100644 --- a/pyorbital/tests/__init__.py +++ b/pyorbital/tests/__init__.py @@ -1,11 +1,13 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- -# Copyright (c) 2014 Martin Raspaud +# Copyright (c) 2014-2016 Martin Raspaud # Author(s): # Martin Raspaud +# Adam Dybbroe + # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by @@ -24,23 +26,36 @@ """ from pyorbital.tests import (test_aiaa, test_tlefile, test_orbital, - test_astronomy, test_geoloc) + test_astronomy, test_geoloc, test_moon) +from pyorbital import planets + import unittest +import doctest + +import os +TRAVIS = os.environ.get("TRAVIS", False) + def suite(): """The global test suite. """ mysuite = unittest.TestSuite() - # Test the documentation strings - #mysuite.addTests(doctest.DocTestSuite(image)) + if not TRAVIS: + # Test sphinx documentation pages: + # mysuite.addTests(doctest.DocFileSuite('../doc/source/index.rst')) + # Test the documentation strings + mysuite.addTests(doctest.DocTestSuite(planets)) + # Use the unittests also mysuite.addTests(test_aiaa.suite()) mysuite.addTests(test_tlefile.suite()) - mysuite.addTests(test_orbital.suite()) mysuite.addTests(test_astronomy.suite()) + mysuite.addTests(test_orbital.suite()) mysuite.addTests(test_geoloc.suite()) - + mysuite.addTests(test_moon.suite()) + return mysuite + if __name__ == '__main__': unittest.TextTestRunner(verbosity=2).run(suite()) diff --git a/pyorbital/tests/test_aiaa.py b/pyorbital/tests/test_aiaa.py index a820e68b..3a286768 100644 --- a/pyorbital/tests/test_aiaa.py +++ b/pyorbital/tests/test_aiaa.py @@ -1,11 +1,11 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- -# Copyright (c) 2011, 2014 SMHI - +# Copyright (c) 2011, 2012, 2014, 2016, 2018 SMHI # Author(s): # Martin Raspaud +# Adam Dybbroe # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by @@ -38,6 +38,7 @@ class LineOrbital(Orbital): + """Read TLE lines instead of file. """ @@ -80,6 +81,7 @@ def get_results(satnumber, delay): class AIAAIntegrationTest(unittest.TestCase): + """Test against the AIAA test cases. """ @@ -131,6 +133,7 @@ def test_aiaa(self): delta_pos = 5e-6 # km = 5 mm delta_vel = 5e-9 # km/s = 5 um/s delta_time = 1e-3 # 1 milisecond + self.assertTrue(abs(res[0] - pos[0]) < delta_pos) self.assertTrue(abs(res[1] - pos[1]) < delta_pos) self.assertTrue(abs(res[2] - pos[2]) < delta_pos) @@ -145,8 +148,7 @@ def test_aiaa(self): def suite(): - """The suite for test_aiaa - """ + """The suite for aiaa integration""" loader = unittest.TestLoader() mysuite = unittest.TestSuite() mysuite.addTest(loader.loadTestsFromTestCase(AIAAIntegrationTest)) diff --git a/pyorbital/tests/test_moon.py b/pyorbital/tests/test_moon.py new file mode 100644 index 00000000..be56d6f0 --- /dev/null +++ b/pyorbital/tests/test_moon.py @@ -0,0 +1,233 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- + +# Copyright (c) 2016 Adam.Dybbroe + +# Author(s): + +# Adam.Dybbroe + +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. + +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. + +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + +"""Unit testing the moon phase and moon position/angles calculations +""" + +import unittest +import numpy as np +from datetime import datetime, timedelta +from pyorbital.moon_phase import (moon_phase, Julian, + sun_position, moon_position) +from pyorbital import planets + +LAT, LON = 58.5875, 16.1875 + +MOON_PHASE = np.array( + [0.409753, 0.458902, 0.507830, 0.556157, 0.603532, + 0.649628, 0.694135, 0.736750, 0.777180, 0.815135, + 0.850327, 0.882472, 0.911289, 0.936505, 0.957855, + 0.975088, 0.987971, 0.996291, 0.999865, 0.998545, + 0.992221, 0.980826, 0.964349, 0.942831, 0.916374, + 0.885145, 0.849373, 0.809354, 0.765447, 0.718071, + 0.667707, 0.614887, 0.560201, 0.504283, 0.447817, + 0.391531, 0.336191, 0.282592, 0.231546, 0.183869, + 0.140348, 0.101725, 0.068656, 0.041694, 0.021252, + 0.007593, 0.000819, 0.000868, 0.007533, 0.020473, + 0.039247, 0.063336, 0.092173, 0.125171, 0.161743, + 0.201318, 0.243355, 0.287348, 0.332827, 0.379362, + 0.426553, 0.474030, 0.521443, 0.568459, 0.614757, + 0.660015, 0.703914, 0.746129, 0.786327, 0.824166, + 0.859293, 0.891348, 0.919965, 0.944774, 0.965412, + 0.981527, 0.992790, 0.998906, 0.999623, 0.994750, + 0.984166, 0.967833, 0.945803, 0.918221, 0.885333, + 0.847479, 0.805085, 0.758660, 0.708778, 0.656070, + 0.601211, 0.544909, 0.487898, 0.430928, 0.374757, + 0.320146, 0.267844, 0.218581, 0.173048, 0.131881]) + +MOON_PHASE_ARR1 = np.array([0.59276087, 0.581796, 0.57076903, 0.55968404, + 0.54854521, 0.5373568, 0.5261231, 0.51484851, + 0.5035375, 0.49219465, 0.48082457, 0.46943203, + 0.45802189, 0.4465991, 0.43516876, 0.42373606, + 0.41230631, 0.40088499, 0.38947769, 0.3780901]) + + +LONARR = np.array([10.000, 11.000, 12.000, 13.000, 14.000, 14.989]) +LATARR = np.array([50.000, 51.000, 52.000, 53.000, 54.000, 54.989]) +RESULT1_ALT = np.array([6.46076487, 6.65664799, 6.82103169, + 6.95359249, 7.05405387, 7.12161456]) +RESULT1_AZI = np.array([112.50462386, 113.42950322, 114.3687665, + 115.32140457, 116.28639318, 117.25189999]) + + +DTOBJ_LIST = [datetime(2016, 3, 1) + timedelta(seconds=10000 * i) + for i in range(20)] +JDAYS_RESULT = np.array([2457448.5, 2457448.61574078, 2457448.73148143, + 2457448.84722221, 2457448.96296299, 2457449.07870364, + 2457449.19444442, 2457449.31018519, 2457449.42592597, + 2457449.54166651, 2457449.65740728, 2457449.77314806, + 2457449.88888884, 2457450.00462961, 2457450.12037039, + 2457450.23611116, 2457450.35185194, 2457450.46759272, + 2457450.58333325, 2457450.69907403]) + + +class TestMoon(unittest.TestCase): + + """Unit testing the moon calculations""" + + def assertNumpyArraysEqual(self, this, that, ndigits=7, msg=''): + """ + modified from http://stackoverflow.com/a/15399475/5459638 + """ + atol = 1. / (10.**ndigits) + + if this.shape != that.shape: + raise AssertionError("Shapes don't match") + if not np.allclose(this, that, atol=atol, rtol=0): + raise AssertionError("Elements don't match!") + + def setUp(self): + """Set up""" + + self.start_time = datetime(2011, 12, 1, hour=12) + self.delta_t = timedelta(hours=12) + return + + def tearDown(self): + """Clean up""" + return + + def test_julian(self): + """Test the conversion from datetime objects to Julian days""" + + dtobj = datetime(2016, 1, 1, 12) + jday = Julian(dtobj) + self.assertEqual(jday, 2457389.0) + + dtobj = datetime(2016, 4, 1, 0) + jday = Julian(dtobj) + self.assertEqual(jday, 2457479.5) + + dtobj = datetime(2016, 6, 15, 18, 30) + jday = Julian(dtobj) + self.assertAlmostEqual(jday, 2457555.270833, 5) + + dtobj = datetime(2011, 9, 30, 22, 30, 30) + jday = Julian(dtobj) + self.assertAlmostEqual(jday, 2455835.4378472, 5) + + dtobj = datetime(1617, 12, 5, 3, 50, 1) + jday = Julian(dtobj) + self.assertAlmostEqual(jday, 2311995.659734, 5) + + dtobj = datetime(300, 10, 17, 1, 15, 15) + jday = Julian(dtobj) + self.assertAlmostEqual(jday, 1830922.552257, 5) + + jdays = Julian(DTOBJ_LIST) + self.assertNumpyArraysEqual(jdays, JDAYS_RESULT, 7) + + def test_sun_position(self): + """Test the derivation of the sun position""" + + jday = 2457389.0 + sunp = sun_position(jday) + self.assertAlmostEqual(sunp, 318.867306265, 7) + + jday = 2457800.0 + sunp = sun_position(jday) + self.assertAlmostEqual(sunp, 4.7408030401432484, 7) + + jday = 2000000.1 + sunp = sun_position(jday) + self.assertAlmostEqual(sunp, 211.79460905279308, 7) + + jday = 2433000.345 + sunp = sun_position(jday) + self.assertAlmostEqual(sunp, 40.830400752277228, 7) + + def test_moon_position(self): + """Test the derivation of the moon position""" + + jday = 2457389.0 + sunp = 318.86730626473212 + lmoon = moon_position(jday, sunp) + self.assertAlmostEqual(lmoon, 110.66501268828912, 7) + + jday = 2457800.0 + sunp = 4.7408030401432484 + lmoon = moon_position(jday, sunp) + self.assertAlmostEqual(lmoon, 123.49626136219904, 7) + + jday = 2000000.1 + sunp = 211.79460905279308 + lmoon = moon_position(jday, sunp) + self.assertAlmostEqual(lmoon, 138.58194894427129, 7) + + jday = 2433000.345 + sunp = 40.830400752277228 + lmoon = moon_position(jday, sunp) + self.assertAlmostEqual(lmoon, 236.15012441309614, 7) + + def test_moon_phase(self): + + time_t = self.start_time + phase = moon_phase(time_t) + self.assertAlmostEqual(phase, 0.409753, 5) + + # Should eventually be possible to take an aray of datetimes! + phase = [] + for i in range(100): + phase.append(moon_phase(time_t)) + time_t = time_t + self.delta_t + + phase = np.array(phase) + + self.assertNumpyArraysEqual(phase, MOON_PHASE, 6) + + phases = moon_phase(DTOBJ_LIST) + self.assertNumpyArraysEqual(phases, MOON_PHASE_ARR1, 6) + + def test_planets_moon_position(self): + + moon = planets.Moon(self.start_time) + rasc, decl, alt, azi = moon.topocentric_position(LON, LAT) + + self.assertAlmostEqual(alt, 6.01043303, 5) + self.assertAlmostEqual(azi, 118.7126239, 5) + self.assertAlmostEqual(rasc, -32.800324, 5) + self.assertAlmostEqual(decl, -9.183704, 5) + + moon = planets.Moon(self.start_time + self.delta_t * 10) + rasc, decl, alt, azi = moon.topocentric_position(LON, LAT) + + self.assertAlmostEqual(alt, -0.8814571, 5) + self.assertAlmostEqual(azi, 62.974637, 5) + self.assertAlmostEqual(rasc, 24.15641340, 5) + self.assertAlmostEqual(decl, 12.9278790, 5) + + def test_planets_moon_positions(self): + + moon = planets.Moon(self.start_time) + rasc, decl, alt, azi = moon.topocentric_position(LONARR, LATARR) + + self.assertNumpyArraysEqual(alt, RESULT1_ALT, 6) + self.assertNumpyArraysEqual(azi, RESULT1_AZI, 6) + + +def suite(): + """The suite for moon phase and angles testing""" + loader = unittest.TestLoader() + mysuite = unittest.TestSuite() + mysuite.addTest(loader.loadTestsFromTestCase(TestMoon)) + + return mysuite diff --git a/pyorbital/tests/test_orbital.py b/pyorbital/tests/test_orbital.py index d5e4be0b..82450091 100644 --- a/pyorbital/tests/test_orbital.py +++ b/pyorbital/tests/test_orbital.py @@ -1,11 +1,13 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- -# Copyright (c) 2012-2014 Martin Raspaud +# Copyright (c) 2012-2014, 2016, 2018 Martin Raspaud # Author(s): # Martin Raspaud +# Adam Dybbroe + # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by @@ -109,7 +111,6 @@ def test_get_next_passes_apogee(self): self.assertTrue(abs(res[0][2] - datetime(2018, 3, 7, 3, 48, 13, 178439)) < timedelta(seconds=0.01)) - def suite(): """The suite for test_orbital """ diff --git a/pyorbital/utils.py b/pyorbital/utils.py new file mode 100644 index 00000000..3c065716 --- /dev/null +++ b/pyorbital/utils.py @@ -0,0 +1,66 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- + +# Copyright (c) 2016 Adam.Dybbroe + +# Author(s): + +# Adam.Dybbroe + +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. + +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. + +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + +""" +""" + + +def JulianToDate(jday): + """Get the current datetime object from the julian day""" + import datetime + + jday = jday + 0.5 + jdi = long(jday) + if (jdi > 2299160): + a__ = long((jdi - 1867216.25) / 36524.25) + b__ = long(jdi + 1 + a__ - a__ / 4) + else: + b__ = jdi + + c__ = long(b__ + 1524) + d__ = long((c__ - 122.1) / 365.25) + e__ = long(365.25 * d__) + g__ = long((c__ - e__) / 30.6001) + g1_ = long(30.6001 * g__) + + day = c__ - e__ - g1_ + fhour = float((jday - jdi) * 24.0) + if (g__ <= 13): + month = g__ - 1 + else: + month = g__ - 13 + + if (month > 2): + year = d__ - 4716 + else: + year = d__ - 4715 + + hour = int(fhour) + minutes_of_hour = 60 * (fhour - int(hour)) + minutes = int(minutes_of_hour) + seconds = 60 * (minutes_of_hour - minutes) + microsec = int((seconds - int(seconds)) * 1000000) + seconds = int(seconds) + + return datetime.datetime(year, month, day, + hour=hour, minute=minutes, + second=seconds, microsecond=microsec) diff --git a/setup.py b/setup.py index 41ab2923..1e92276c 100644 --- a/setup.py +++ b/setup.py @@ -1,11 +1,13 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- -# Copyright (c) 2011-2014 + +# Copyright (c) 2011-2016, 2018 SMHI # Author(s): # Martin Raspaud +# Adam Dybbroe # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by @@ -38,10 +40,10 @@ "Programming Language :: Python", "Topic :: Scientific/Engineering", "Topic :: Scientific/Engineering :: Astronomy"], - url="https://github.com/mraspaud/pyorbital", + url="https://github.com/pytroll/pyorbital", test_suite='pyorbital.tests.suite', - package_dir = {'pyorbital': 'pyorbital'}, - packages = ['pyorbital'], + package_dir={'pyorbital': 'pyorbital'}, + packages=['pyorbital'], install_requires=['numpy>=1.6.0,!=1.14.0'], zip_safe=False, )