Skip to content

Commit

Permalink
Support for array arithmetics
Browse files Browse the repository at this point in the history 
Signed-off-by: Adam.Dybbroe <a000680@c20671.ad.smhi.se>
  • Loading branch information
Adam.Dybbroe committed Jul 14, 2016
1 parent d22cb3f commit 7108093
Show file tree
Hide file tree
Showing 3 changed files with 213 additions and 60 deletions.
26 changes: 26 additions & 0 deletions doc/source/moon_calculations.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,21 @@ Computing the phase of the moon
>>> 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
----------------------------------

Expand All @@ -29,6 +44,17 @@ Computing the position of the moon
>>> 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
--------
Expand Down
147 changes: 90 additions & 57 deletions pyorbital/moon_phase.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,7 +33,7 @@
SMALL_FLOAT = 1.0e-12


def seqJulian(dobj):
def Julian(dobj):
"""Returns the number of julian days for the specified date-time.
Input = datetime object.
"""
Expand All @@ -46,52 +46,112 @@ def seqJulian(dobj):
dobj = np.array(dobj, 'datetime64[us]')

day = (
dobj - dobj.astype('datetime64[M]')).astype('f') / (24. * 3600. * 1000) + 1.0
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.less(year, 1582)
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, -10)
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 Julian(dobj):
"""Returns the number of julian days for the specified date-time.
Input = datetime object.
"""
# 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.))
# 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 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
# 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
# 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):
Expand All @@ -102,9 +162,9 @@ def sun_position(jday):
# int i;

if isinstance(jday, collections.Sequence):
jday = np.array(jday)
jday = np.array(jday, 'float64')
elif not isinstance(jday, np.ndarray):
jday = np.array([jday], 'f')
jday = np.array([jday], 'float64')

n__ = 360. / 365.2422 * jday
i__ = np.divide(n__, 360.).astype('i')
Expand Down Expand Up @@ -132,46 +192,19 @@ def sun_position(jday):

return sunpos

# def moon_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 moon_position(jday, lsun):
"""Get the moon position"""

if isinstance(jday, collections.Sequence):
jday = np.array(jday)
jday = np.array(jday, 'float64')
elif not isinstance(jday, np.ndarray):
jday = np.array([jday], 'f')
jday = np.array([jday], 'float64')

if isinstance(lsun, collections.Sequence):
lsun = np.array(lsun)
lsun = np.array(lsun, 'float64')
elif not isinstance(lsun, np.ndarray):
lsun = np.array([lsun], 'f')
lsun = np.array([lsun], 'float64')

ms_ = 0.985647332099 * jday - 3.762863
ms_ = np.where(np.less(ms_, 0), ms_ + 360.0, ms_)
Expand Down
100 changes: 97 additions & 3 deletions pyorbital/tests/test_moon.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,8 @@
import unittest
import numpy as np
from datetime import datetime, timedelta
from pyorbital.moon_phase import moon_phase
from pyorbital.moon_phase import (moon_phase, Julian,
sun_position, moon_position)
from pyorbital import planets

LAT, LON = 58.5875, 16.1875
Expand All @@ -53,6 +54,13 @@
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,
Expand All @@ -61,6 +69,17 @@
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"""
Expand All @@ -87,6 +106,78 @@ 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
Expand All @@ -103,7 +194,10 @@ def test_moon_phase(self):

self.assertNumpyArraysEqual(phase, MOON_PHASE, 6)

def test_moon_position(self):
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)
Expand All @@ -121,7 +215,7 @@ def test_moon_position(self):
self.assertAlmostEqual(rasc, 24.15641340, 5)
self.assertAlmostEqual(decl, 12.9278790, 5)

def test_moon_positions(self):
def test_planets_moon_positions(self):

moon = planets.Moon(self.start_time)
rasc, decl, alt, azi = moon.topocentric_position(LONARR, LATARR)
Expand Down

0 comments on commit 7108093

Please sign in to comment.