#179

pending completion

Build # #179

Build Type

push

coveralls-python

Committed by

web-flow

Commit Message

<a href="https://github.com/USDA-ARS-NWRC/smrf/commit/<a class=hub.com/USDA-ARS-NWRC/smrf/commit/a56cd2425c51792bf551f981dc8226e182044cb7">a56cd2425<a href="https://github.com/USDA-ARS-NWRC/smrf/commit/a56cd2425c51792bf551f981dc8226e182044cb7">">Merge </a><a class="double-link" href="https://github.com/USDA-ARS-NWRC/smrf/commit/<a class="double-link" href="https://github.com/USDA-ARS-NWRC/smrf/commit/bebb0d901ed7e57909a0934b6e4902d642404154">bebb0d901</a>">bebb0d901</a><a href="https://github.com/USDA-ARS-NWRC/smrf/commit/a56cd2425c51792bf551f981dc8226e182044cb7"> into 042fbff88">042fbff88</a>

Run Details

343 of 343 new or added lines in 3 files covered. (100.0%)

2704 of 3820 relevant lines covered (70.79%)

0.71 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

91.41

/smrf/envphys/sunang.py

import logging
from math import asin, atan, cos, fmod, sin, sqrt, tan

import numpy as np
import pytz

JULIAN_CENTURY = 36525                # days in Julian century
DEGS_IN_CIRCLE = 3.6e2        # degrees in circle
TOLERANCE = 1.1920928955e-07
M_PI = 3.14159265358979323846
M_PI_2 = 2 * M_PI
MAX_DECLINATION = 23.6


def sunang(date_time, latitude, longitude, truncate=True):
    """
    Calculate the sun angle (the azimuth and zenith angles of
    the sun's position) for a given geodetic location for a single
    date time and coordinates. The function can take either latitude
    longitude position as a single point or numpy array.

    Args:
        date_time: python datetime object
        latitude: value or np.ndarray (in degrees)
        longitude: value or np.ndarray (in degrees)
        truncate: True will replicate the IPW output precision, not applied
            if position is an array

    Returns
        cosz - cosine of the zenith angle, same shape as input position
        azimuth - solar azimuth, same shape as input position
        rad_vec - Earth-Sun radius vector
    """

    # Calculate the ephemeris parameters
    declination, omega, rad_vec = ephemeris(date_time)

    # calculate the sun angle
    rad_latitude = latitude * M_PI /180
    rad_longitude = longitude * M_PI / 180
    mu, azimuth = sunpath(rad_latitude, rad_longitude, declination, omega)

    azimuth = azimuth * 180 / M_PI

    if truncate and not isinstance(mu, np.ndarray):
        mu = float('{0:.6f}'.format(mu))
        azimuth = float('{0:.5g}'.format(azimuth))
        rad_vec = float('{0:.5g}'.format(rad_vec))

    return mu, azimuth, rad_vec


def sunang_thread(queue, date, lat, lon):
    """
    See sunang for input descriptions

    Args:
        queue: queue with cosz, azimuth
        date: loop through dates to accesss queue, must be same as rest of queues

    """

    if 'cosz' not in queue.keys():
        raise ValueError('queue must have cosz key')
    if 'azimuth' not in queue.keys():
        raise ValueError('queue must have cosz key')

    log = logging.getLogger(__name__)

    for t in date:

        log.debug('%s Calculating sun angle' % t)

        cosz, azimuth, rad_vec = sunang(t.astimezone(pytz.utc), lat, lon)

        queue['cosz'].put([t, cosz])
        queue['azimuth'].put([t, azimuth])


def sunpath(latitude, longitude, declination, omega):
    """
    Sun angle from solar declination and longtitude

    Args:
        latitude: in radians
        longitude: in radians
        declination: solar declination (radians)
        omega: solar longitude (radians)

    Returns
        cosz: cosine of solar zenith
        azimuth: solar azimuth in radians
    """

    if isinstance(latitude, np.ndarray):
        if np.any(np.abs(latitude) > M_PI_2):
            raise ValueError("latitude array not between -90 and +90 degrees")

        if np.any(np.abs(longitude) > M_PI):
            raise ValueError("longitude array not between -180 and +180 degrees")

    else:
        if np.abs(latitude) > M_PI_2:
            raise ValueError("latitude ({} deg) not between -90 and +90 degrees".format(latitude*180/M_PI))

        if np.abs(longitude) > M_PI:
            raise ValueError("longitude ({} deg) not between -180 and +180 degrees".format(longitude*180/M_PI))

    if np.abs(declination) > M_PI*MAX_DECLINATION/180:
        raise ValueError("declination ({} deg) > maximum declination ({} deg)".format(declination*180/M_PI, MAX_DECLINATION))

    if np.abs(omega) > M_PI:
        raise ValueError("longitude of sun ({} deg) not between -180 and +180 degrees".format(omega*180/M_PI))

    cosz, az = rotate(np.sin(declination), omega, np.sin(latitude), longitude)

    return cosz, az

def dsign(a, b):
    """
    modified from /usr/src/lib/libF77/d_sign.c
    """
    
    x = a if a >= 0 else -a
    y = x if b >= 0 else -x
    return y

# original ibm 360/44 fortran ivf - vislab - wilson - 29jul79
# translated, modified and reduced by dozier - ucsb - 11/81

def ephemeris(dt):
    """
    Calculates radius vector, declination, and apparent longitude
    of sun, as function of the given date and time.

    The routine is adapted from:

    W. H. Wilson, Solar ephemeris algorithm, Reference 80-13, 70
        pp., Scripps Institution of Oceanography, University of
        California, San Diego, La Jolla, CA, 1980.

    Args:
        dt: date time python object

    Returns:
        declin: solar declination angle, in radians
        omega: sun longitude, in radians
        r: Earth-Sun radius vector
    """
        
    one = 1.0
    degrd = atan(one) / 45.0

    # Convert time to seconds since midnight
    gmts = 3600.0 * dt.hour + 60.0 * dt.minute + dt.second

    # p51 = gmts/3600/24*360

    # The int() is required for compatability with IPW as the divide by 100 keeps
    # the value as an integer where python converts to a float...
    p51 = gmts / 10.0 / 24.0
    p22 = int(((dt.year - 1900) * JULIAN_CENTURY - 25) / 100) + yearday(dt.year, dt.month, dt.day) - 0.5
    p23 = (p51 / DEGS_IN_CIRCLE + p22) / JULIAN_CENTURY
    p22 = p23 * JULIAN_CENTURY

    # mean longitude - p24
    p11 = 279.69668
    p12 = 0.9856473354
    p13 = 3.03e-4
    p24 = p11 + fmod(p12 * p22, DEGS_IN_CIRCLE) + p13 * p23 * p23
    p24 = fmod(p24, DEGS_IN_CIRCLE)

    # mean anomaly - p25
    p11 = 358.47583
    p12 = 0.985600267
    p13 = -1.5e-4
    p14 = -3.e-6
    p25 = p11 + fmod(p12 * p22, DEGS_IN_CIRCLE) + p23 * p23 * (p13 + p14 * p23)
    p25 = fmod(p25, DEGS_IN_CIRCLE)
 
    # eccentricity - p26
    p11 = 0.01675104
    p12 = -4.18e-5
    p13 = -1.26e-7
    p26 = p11 + p23 * (p12 + p13 * p23)
    p11 = p25 * degrd
    p12 = p11

    while True:
        p13 = p12
        p12 = p11 + p26 * sin(p13)
        if abs((p12 - p13) / p12) < TOLERANCE:
            break

    p13 = p12 / degrd

    # true anomaly - p27`
    p27 = 2.0 * atan(sqrt((1.0 + p26) / (1.0 - p26)) * tan(p13 / 2.0 * degrd)) / degrd
    if dsign(1.0, p27) != dsign(1.0, sin(p13 * degrd)):
        p27 += 1.8e2
    if p27 < 0.0:
        p27 += DEGS_IN_CIRCLE
 
    # radius vector - r
    r = 1.0 - p26 * cos(p13 * degrd)

    # aberration - p29
    p29 = -20.47 / r / 3600

    # mean obliquity - p43
    p11 = 23.452294
    p12 = -0.0130125
    p13 = -1.64e-6
    p14 = 5.03e-7
    p43 = p11 + p23 * (p12 + p23 * (p13 + p14 * p23))
 
    #  mean ascension - p45
    p11 = 279.6909832
    p12 = 0.98564734
    p13 = 3.8707
    p13 = 3.8708e-4
    p45 = p11 + fmod(p12 * p22, DEGS_IN_CIRCLE) + p13 * p23 * p23
    p45 = fmod(p45, DEGS_IN_CIRCLE)

    # nutation and longitude pert
    p11 = 296.104608
    p12 = 1325 * DEGS_IN_CIRCLE
    p13 = 198.8491083
    p14 = 0.00919167
    p15 = 1.4388e-5
    p28 = p11 + fmod(p12 * p23, DEGS_IN_CIRCLE) + p23 * (p13 + p23 * (p14 + p15 * p23))
    p28 = fmod(p28, DEGS_IN_CIRCLE)

    # mean elongation of moon - p30
    p11 = 350.737486
    p12 = 1236 * DEGS_IN_CIRCLE
    p13 = 307.1142167
    p14 = 1.436e-3
    p30 = p11 + fmod(p12 * p23, DEGS_IN_CIRCLE) + p23 * (p13 + p14 * p23)
    p30 = fmod(p30, DEGS_IN_CIRCLE)

    # moon long of ascending node - p31
    p11 = 259.183275
    p12 = -5 * DEGS_IN_CIRCLE
    p13 = -134.142008
    p14 = 2.0778e-3
    p31 = p11 + fmod(p12 * p23, DEGS_IN_CIRCLE) + p23 * (p13 + p14 * p23)
    p31 = fmod(p31, DEGS_IN_CIRCLE)

    # mean long of moon - p31
    p11 = 270.434164
    p12 = 1336 * DEGS_IN_CIRCLE
    p13 = 307.8831417
    p14 = -1.1333e-3
    p32 = p11 + fmod(p12 * p23, DEGS_IN_CIRCLE) + p23 * (p13 + p14 * p23)
    p32 = fmod(p32, DEGS_IN_CIRCLE)

    # moon perturbation of sun long - p33
    p33 = 6.454 * sin(p30 * degrd) + 0.013 * sin(3 * p30 * degrd) + \
        0.177 * sin((p30 + p28) * degrd) - \
        0.424 * sin((p30 - p28) * degrd)
    p33 = p33 / 3600

    # nutation of long - p34
    p34 = -(17.234 - 0.017 * p23) * sin(p31 * degrd) + \
        0.209 * sin(2 * p31 * degrd) - \
        0.204 * sin(2 * p32 * degrd)
    p34 = p34 - 1.257 * sin(2 * p24 * degrd) + 0.127 * sin(p28 * degrd)
    p34 = p34 / 3600

    # nutation in obliquity - p34
    p35 = 9.214 * cos(p31 * degrd) + 0.546 * cos(2 * p24 * degrd) - \
        .09 * cos(2 * p31 * degrd) + 0.088 * cos(2 * p32 * degrd)
    p35 = p35 / 3600

    # inequalities of long period - p36
    p36 = 0.266 * sin((31.8 + 119 * p23) * degrd) + \
        ((1.882 - 0.016 * p23) * degrd) * \
        sin((57.24 + 150.27 * p23) * degrd)
    p36 = p36 + 0.202 * sin((315.6 + 893.3 * p23) * degrd) + \
        1.089 * p23 * p23 + 6.4 * sin((231.19 + 20.2 * p23) * degrd)
    p36 = p36 / 3600

    # apparent longitude - p41
    p41 = p27 - p25 + p24 + p29 + p33 + p36 + p34

    p43 += p35

    # apparent right ascension - p44
    p44 = atan(tan(p41 * degrd) * cos(p43 * degrd)) / degrd
    if dsign(one, p44) != dsign(one, sin(p41 * degrd)):
        p44 += 1.8e2
    if p44 < 0:
        p44 += DEGS_IN_CIRCLE

    # equation of time - p46
    p46 = p45 - p44
    if p46 > 1.8e2:
        p46 -= DEGS_IN_CIRCLE

    # declination - p47
    p47 = asin(sin(p41 * degrd) * sin(p43 * degrd))
    declin = p47

    # hour angle in degrees - p48
    p48 = p51 + p46 - 1.8e2
    if p48 > 180:
        p48 -= DEGS_IN_CIRCLE
    elif p48 < -180:
        p48 += DEGS_IN_CIRCLE
    omega = -p48 * degrd

    return declin, omega, r

def rotate(mu, azm, mu_r, lam_r):
    """
    Calculates new spherical coordinates if system rotated about
    origin.  Coordinates are right-hand system.  All angles are in
    radians.

    Args:
        mu: cosine of angle theta from z-axis in old coordinate 
            system, sin(declination)
        azm: azimuth (+ccw from x-axis) in old coordinate system, 
            hour angle of sun (long. where sun is vertical)
        mu_r: cosine of angle theta_r of rotation of z-axis, sin(latitude)
        lam_r: azimuth (+ccw) of rotation of x-axis, longitude

    Returns:
        muPrime: cosine of the solar zenith
        aPrime: solar azimuth in radians
    """

    # Check input values: mu = cos(theta) mu_r = cos(theta-sub-r)
    if mu < -1.0 or mu > 1.0:
        raise ValueError("rotate: mu = cos(theta) = {} is not between -1 and +1".format(mu))

    if isinstance(mu_r, np.ndarray):
        if np.any(mu_r < -1.0) or np.any(mu_r > 1.0):
            raise ValueError("rotate, mu_r array is not between -1 and +1")
    else:
        if mu_r < -1.0 or mu_r > 1.0:
            raise ValueError("rotate: mu_rotate = cos(theta-sub-r) = {} is not between -1 and +1".format(mu_r))

    if np.abs(azm) > (2 * np.pi):
        raise ValueError("rotate: azimuth ({} deg) is not between -360 and 360".format(180 * (azm * 360) / np.pi))

    if isinstance(lam_r, np.ndarray):
        if np.any(np.abs(lam_r) > (2 * np.pi)):
            raise ValueError("rotate, lam_r array is not between -360 and 360")
    else:
        if np.abs(lam_r) > (2 * np.pi):
            raise ValueError("rotate: lam_r ({} deg) = axis rotation, is not between -360 and 360".format(180 * lam_r / np.pi))

    # difference between azimuth and rotation angle of x-axis
    omega = lam_r - azm

    #sine of angle theta (cosine is an argument)
    #sine of angle theta-sub-r (cosine is an argument)
    sinTheta = np.sqrt((1. - mu) * (1. + mu))
    sinThr = np.sqrt((1. - mu_r) * (1. + mu_r))
        
    #cosine of difference between azimuth and aziumth of rotation
    cosOmega = np.cos(omega)

    #Output:
    #azimuth and cosine of angle in rotated axis system.
    #(bug if trig routines trigger error)
    aPrime = -np.arctan2(sinTheta * np.sin(omega),
                mu_r * sinTheta * cosOmega - mu * sinThr)
    muPrime = sinTheta * sinThr * cosOmega + mu * mu_r
    
    return muPrime, aPrime


def yearday(year, month, day):
    """
    yearday returns the yearday for the given date.  yearday is
    the 'day of the year', sometimes called (incorrectly) 'julian day'.

    Args:
        year
        month
        day

    Returns:
        yday: day of year
        
    """

    assert(year >= 0)
    assert(1 <= month and month <= 12)
    assert(1 <= day and day <= numdays(year, month))

    yday = day

    # Add in number of days for preceeding months.
    for mon in range(month-1, 0, -1):
        yday += numdays(year, mon)

    return yday


def numdays(year, month):
    """
    numdays returns the number of days in the given month of
    the given year.

    Args:
        year
        month
    
    Returns:
        ndays: number of days in month
    """

    NDAYS = list([31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31])

    assert(year >= 0)
    assert(1 <= month and month <= 12)

    ndays = NDAYS[month-1]

    # Check for leap year for February
    if ((month == 2) and leapyear(year)):
        ndays += 1

    return ndays


def leapyear(year):
    """
    leapyear determines if the given year is a leap year or not.
    year must be positive, and must not be abbreviated; i.e.
    89 is 89 A.D. not 1989.

    Args:
        year

    Returns:
        True if a leap year, False if not a leap year
    """

    assert(year >= 0)

    if ((year % 4 == 0 and year % 100 != 0) or year % 400 == 0):
        return True
    else:
        return False

1	import logging	1✔
2	from math import asin, atan, cos, fmod, sin, sqrt, tan	1✔
3
4	import numpy as np	1✔
5	import pytz	1✔
6
7	JULIAN_CENTURY = 36525 # days in Julian century	1✔
8	DEGS_IN_CIRCLE = 3.6e2 # degrees in circle	1✔
9	TOLERANCE = 1.1920928955e-07	1✔
10	M_PI = 3.14159265358979323846	1✔
11	M_PI_2 = 2 * M_PI	1✔
12	MAX_DECLINATION = 23.6	1✔
13
14
15	def sunang(date_time, latitude, longitude, truncate=True):	1✔
16	"""
17	Calculate the sun angle (the azimuth and zenith angles of
18	the sun's position) for a given geodetic location for a single
19	date time and coordinates. The function can take either latitude
20	longitude position as a single point or numpy array.
21
22	Args:
23	date_time: python datetime object
24	latitude: value or np.ndarray (in degrees)
25	longitude: value or np.ndarray (in degrees)
26	truncate: True will replicate the IPW output precision, not applied
27	if position is an array
28
29	Returns
30	cosz - cosine of the zenith angle, same shape as input position
31	azimuth - solar azimuth, same shape as input position
32	rad_vec - Earth-Sun radius vector
33	"""
34
35	# Calculate the ephemeris parameters
36	declination, omega, rad_vec = ephemeris(date_time)	1✔
37
38	# calculate the sun angle
39	rad_latitude = latitude * M_PI /180	1✔
40	rad_longitude = longitude * M_PI / 180	1✔
41	mu, azimuth = sunpath(rad_latitude, rad_longitude, declination, omega)	1✔
42
43	azimuth = azimuth * 180 / M_PI	1✔
44
45	if truncate and not isinstance(mu, np.ndarray):	1✔
46	mu = float('{0:.6f}'.format(mu))	1✔
47	azimuth = float('{0:.5g}'.format(azimuth))	1✔
48	rad_vec = float('{0:.5g}'.format(rad_vec))	1✔
49
50	return mu, azimuth, rad_vec	1✔
51
52
53	def sunang_thread(queue, date, lat, lon):	1✔
54	"""
55	See sunang for input descriptions
56
57	Args:
58	queue: queue with cosz, azimuth
59	date: loop through dates to accesss queue, must be same as rest of queues
60
61	"""
62
63	if 'cosz' not in queue.keys():	1✔
64	raise ValueError('queue must have cosz key')	×
65	if 'azimuth' not in queue.keys():	1✔
66	raise ValueError('queue must have cosz key')	×
67
68	log = logging.getLogger(__name__)	1✔
69
70	for t in date:	1✔
71
72	log.debug('%s Calculating sun angle' % t)	1✔
73
74	cosz, azimuth, rad_vec = sunang(t.astimezone(pytz.utc), lat, lon)	1✔
75
76	queue['cosz'].put([t, cosz])	1✔
77	queue['azimuth'].put([t, azimuth])	1✔
78
79
80	def sunpath(latitude, longitude, declination, omega):	1✔
81	"""
82	Sun angle from solar declination and longtitude
83
84	Args:
85	latitude: in radians
86	longitude: in radians
87	declination: solar declination (radians)
88	omega: solar longitude (radians)
89
90	Returns
91	cosz: cosine of solar zenith
92	azimuth: solar azimuth in radians
93	"""
94
95	if isinstance(latitude, np.ndarray):	1✔
96	if np.any(np.abs(latitude) > M_PI_2):	1✔
97	raise ValueError("latitude array not between -90 and +90 degrees")	×
98
99	if np.any(np.abs(longitude) > M_PI):	1✔
100	raise ValueError("longitude array not between -180 and +180 degrees")	×
101
102	else:
103	if np.abs(latitude) > M_PI_2:	1✔
104	raise ValueError("latitude ({} deg) not between -90 and +90 degrees".format(latitude*180/M_PI))	×
105
106	if np.abs(longitude) > M_PI:	1✔
107	raise ValueError("longitude ({} deg) not between -180 and +180 degrees".format(longitude*180/M_PI))	×
108
109	if np.abs(declination) > M_PI*MAX_DECLINATION/180:	1✔
110	raise ValueError("declination ({} deg) > maximum declination ({} deg)".format(declination*180/M_PI, MAX_DECLINATION))	×
111
112	if np.abs(omega) > M_PI:	1✔
113	raise ValueError("longitude of sun ({} deg) not between -180 and +180 degrees".format(omega*180/M_PI))	×
114
115	cosz, az = rotate(np.sin(declination), omega, np.sin(latitude), longitude)	1✔
116
117	return cosz, az	1✔
118
119	def dsign(a, b):	1✔
120	"""
121	modified from /usr/src/lib/libF77/d_sign.c
122	"""
123
124	x = a if a >= 0 else -a	1✔
125	y = x if b >= 0 else -x	1✔
126	return y	1✔
127
128	# original ibm 360/44 fortran ivf - vislab - wilson - 29jul79
129	# translated, modified and reduced by dozier - ucsb - 11/81
130
131	def ephemeris(dt):	1✔
132	"""
133	Calculates radius vector, declination, and apparent longitude
134	of sun, as function of the given date and time.
135
136	The routine is adapted from:
137
138	W. H. Wilson, Solar ephemeris algorithm, Reference 80-13, 70
139	pp., Scripps Institution of Oceanography, University of
140	California, San Diego, La Jolla, CA, 1980.
141
142	Args:
143	dt: date time python object
144
145	Returns:
146	declin: solar declination angle, in radians
147	omega: sun longitude, in radians
148	r: Earth-Sun radius vector
149	"""
150
151	one = 1.0	1✔
152	degrd = atan(one) / 45.0	1✔
153
154	# Convert time to seconds since midnight
155	gmts = 3600.0 * dt.hour + 60.0 * dt.minute + dt.second	1✔
156
157	# p51 = gmts/3600/24*360
158
159	# The int() is required for compatability with IPW as the divide by 100 keeps
160	# the value as an integer where python converts to a float...
161	p51 = gmts / 10.0 / 24.0	1✔
162	p22 = int(((dt.year - 1900) * JULIAN_CENTURY - 25) / 100) + yearday(dt.year, dt.month, dt.day) - 0.5	1✔
163	p23 = (p51 / DEGS_IN_CIRCLE + p22) / JULIAN_CENTURY	1✔
164	p22 = p23 * JULIAN_CENTURY	1✔
165
166	# mean longitude - p24
167	p11 = 279.69668	1✔
168	p12 = 0.9856473354	1✔
169	p13 = 3.03e-4	1✔
170	p24 = p11 + fmod(p12 * p22, DEGS_IN_CIRCLE) + p13 * p23 * p23	1✔
171	p24 = fmod(p24, DEGS_IN_CIRCLE)	1✔
172
173	# mean anomaly - p25
174	p11 = 358.47583	1✔
175	p12 = 0.985600267	1✔
176	p13 = -1.5e-4	1✔
177	p14 = -3.e-6	1✔
178	p25 = p11 + fmod(p12 * p22, DEGS_IN_CIRCLE) + p23 * p23 * (p13 + p14 * p23)	1✔
179	p25 = fmod(p25, DEGS_IN_CIRCLE)	1✔
180
181	# eccentricity - p26
182	p11 = 0.01675104	1✔
183	p12 = -4.18e-5	1✔
184	p13 = -1.26e-7	1✔
185	p26 = p11 + p23 * (p12 + p13 * p23)	1✔
186	p11 = p25 * degrd	1✔
187	p12 = p11	1✔
188
189	while True:	1✔
190	p13 = p12	1✔
191	p12 = p11 + p26 * sin(p13)	1✔
192	if abs((p12 - p13) / p12) < TOLERANCE:	1✔
193	break	1✔
194
195	p13 = p12 / degrd	1✔
196
197	# true anomaly - p27`
198	p27 = 2.0 * atan(sqrt((1.0 + p26) / (1.0 - p26)) * tan(p13 / 2.0 * degrd)) / degrd	1✔
199	if dsign(1.0, p27) != dsign(1.0, sin(p13 * degrd)):	1✔
200	p27 += 1.8e2	×
201	if p27 < 0.0:	1✔
202	p27 += DEGS_IN_CIRCLE	1✔
203
204	# radius vector - r
205	r = 1.0 - p26 * cos(p13 * degrd)	1✔
206
207	# aberration - p29
208	p29 = -20.47 / r / 3600	1✔
209
210	# mean obliquity - p43
211	p11 = 23.452294	1✔
212	p12 = -0.0130125	1✔
213	p13 = -1.64e-6	1✔
214	p14 = 5.03e-7	1✔
215	p43 = p11 + p23 * (p12 + p23 * (p13 + p14 * p23))	1✔
216
217	# mean ascension - p45
218	p11 = 279.6909832	1✔
219	p12 = 0.98564734	1✔
220	p13 = 3.8707	1✔
221	p13 = 3.8708e-4	1✔
222	p45 = p11 + fmod(p12 * p22, DEGS_IN_CIRCLE) + p13 * p23 * p23	1✔
223	p45 = fmod(p45, DEGS_IN_CIRCLE)	1✔
224
225	# nutation and longitude pert
226	p11 = 296.104608	1✔
227	p12 = 1325 * DEGS_IN_CIRCLE	1✔
228	p13 = 198.8491083	1✔
229	p14 = 0.00919167	1✔
230	p15 = 1.4388e-5	1✔
231	p28 = p11 + fmod(p12 * p23, DEGS_IN_CIRCLE) + p23 * (p13 + p23 * (p14 + p15 * p23))	1✔
232	p28 = fmod(p28, DEGS_IN_CIRCLE)	1✔
233
234	# mean elongation of moon - p30
235	p11 = 350.737486	1✔
236	p12 = 1236 * DEGS_IN_CIRCLE	1✔
237	p13 = 307.1142167	1✔
238	p14 = 1.436e-3	1✔
239	p30 = p11 + fmod(p12 * p23, DEGS_IN_CIRCLE) + p23 * (p13 + p14 * p23)	1✔
240	p30 = fmod(p30, DEGS_IN_CIRCLE)	1✔
241
242	# moon long of ascending node - p31
243	p11 = 259.183275	1✔
244	p12 = -5 * DEGS_IN_CIRCLE	1✔
245	p13 = -134.142008	1✔
246	p14 = 2.0778e-3	1✔
247	p31 = p11 + fmod(p12 * p23, DEGS_IN_CIRCLE) + p23 * (p13 + p14 * p23)	1✔
248	p31 = fmod(p31, DEGS_IN_CIRCLE)	1✔
249
250	# mean long of moon - p31
251	p11 = 270.434164	1✔
252	p12 = 1336 * DEGS_IN_CIRCLE	1✔
253	p13 = 307.8831417	1✔
254	p14 = -1.1333e-3	1✔
255	p32 = p11 + fmod(p12 * p23, DEGS_IN_CIRCLE) + p23 * (p13 + p14 * p23)	1✔
256	p32 = fmod(p32, DEGS_IN_CIRCLE)	1✔
257
258	# moon perturbation of sun long - p33
259	p33 = 6.454 * sin(p30 * degrd) + 0.013 * sin(3 * p30 * degrd) + \	1✔
260	0.177 * sin((p30 + p28) * degrd) - \
261	0.424 * sin((p30 - p28) * degrd)
262	p33 = p33 / 3600	1✔
263
264	# nutation of long - p34
265	p34 = -(17.234 - 0.017 * p23) * sin(p31 * degrd) + \	1✔
266	0.209 * sin(2 * p31 * degrd) - \
267	0.204 * sin(2 * p32 * degrd)
268	p34 = p34 - 1.257 * sin(2 * p24 * degrd) + 0.127 * sin(p28 * degrd)	1✔
269	p34 = p34 / 3600	1✔
270
271	# nutation in obliquity - p34
272	p35 = 9.214 * cos(p31 * degrd) + 0.546 * cos(2 * p24 * degrd) - \	1✔
273	.09 * cos(2 * p31 * degrd) + 0.088 * cos(2 * p32 * degrd)
274	p35 = p35 / 3600	1✔
275
276	# inequalities of long period - p36
277	p36 = 0.266 * sin((31.8 + 119 * p23) * degrd) + \	1✔
278	((1.882 - 0.016 * p23) * degrd) * \
279	sin((57.24 + 150.27 * p23) * degrd)
280	p36 = p36 + 0.202 * sin((315.6 + 893.3 * p23) * degrd) + \	1✔
281	1.089 * p23 * p23 + 6.4 * sin((231.19 + 20.2 * p23) * degrd)
282	p36 = p36 / 3600	1✔
283
284	# apparent longitude - p41
285	p41 = p27 - p25 + p24 + p29 + p33 + p36 + p34	1✔
286
287	p43 += p35	1✔
288
289	# apparent right ascension - p44
290	p44 = atan(tan(p41 * degrd) * cos(p43 * degrd)) / degrd	1✔
291	if dsign(one, p44) != dsign(one, sin(p41 * degrd)):	1✔
292	p44 += 1.8e2	1✔
293	if p44 < 0:	1✔
294	p44 += DEGS_IN_CIRCLE	1✔
295
296	# equation of time - p46
297	p46 = p45 - p44	1✔
298	if p46 > 1.8e2:	1✔
299	p46 -= DEGS_IN_CIRCLE	×
300
301	# declination - p47
302	p47 = asin(sin(p41 * degrd) * sin(p43 * degrd))	1✔
303	declin = p47	1✔
304
305	# hour angle in degrees - p48
306	p48 = p51 + p46 - 1.8e2	1✔
307	if p48 > 180:	1✔
308	p48 -= DEGS_IN_CIRCLE	×
309	elif p48 < -180:	1✔
310	p48 += DEGS_IN_CIRCLE	1✔
311	omega = -p48 * degrd	1✔
312
313	return declin, omega, r	1✔
314
315	def rotate(mu, azm, mu_r, lam_r):	1✔
316	"""
317	Calculates new spherical coordinates if system rotated about
318	origin. Coordinates are right-hand system. All angles are in
319	radians.
320
321	Args:
322	mu: cosine of angle theta from z-axis in old coordinate
323	system, sin(declination)
324	azm: azimuth (+ccw from x-axis) in old coordinate system,
325	hour angle of sun (long. where sun is vertical)
326	mu_r: cosine of angle theta_r of rotation of z-axis, sin(latitude)
327	lam_r: azimuth (+ccw) of rotation of x-axis, longitude
328
329	Returns:
330	muPrime: cosine of the solar zenith
331	aPrime: solar azimuth in radians
332	"""
333
334	# Check input values: mu = cos(theta) mu_r = cos(theta-sub-r)
335	if mu < -1.0 or mu > 1.0:	1✔
336	raise ValueError("rotate: mu = cos(theta) = {} is not between -1 and +1".format(mu))	×
337
338	if isinstance(mu_r, np.ndarray):	1✔
339	if np.any(mu_r < -1.0) or np.any(mu_r > 1.0):	1✔
340	raise ValueError("rotate, mu_r array is not between -1 and +1")	×
341	else:
342	if mu_r < -1.0 or mu_r > 1.0:	1✔
343	raise ValueError("rotate: mu_rotate = cos(theta-sub-r) = {} is not between -1 and +1".format(mu_r))	×
344
345	if np.abs(azm) > (2 * np.pi):	1✔
346	raise ValueError("rotate: azimuth ({} deg) is not between -360 and 360".format(180 * (azm * 360) / np.pi))	×
347
348	if isinstance(lam_r, np.ndarray):	1✔
349	if np.any(np.abs(lam_r) > (2 * np.pi)):	1✔
350	raise ValueError("rotate, lam_r array is not between -360 and 360")	×
351	else:
352	if np.abs(lam_r) > (2 * np.pi):	1✔
353	raise ValueError("rotate: lam_r ({} deg) = axis rotation, is not between -360 and 360".format(180 * lam_r / np.pi))	×
354
355	# difference between azimuth and rotation angle of x-axis
356	omega = lam_r - azm	1✔
357
358	#sine of angle theta (cosine is an argument)
359	#sine of angle theta-sub-r (cosine is an argument)
360	sinTheta = np.sqrt((1. - mu) * (1. + mu))	1✔
361	sinThr = np.sqrt((1. - mu_r) * (1. + mu_r))	1✔
362
363	#cosine of difference between azimuth and aziumth of rotation
364	cosOmega = np.cos(omega)	1✔
365
366	#Output:
367	#azimuth and cosine of angle in rotated axis system.
368	#(bug if trig routines trigger error)
369	aPrime = -np.arctan2(sinTheta * np.sin(omega),	1✔
370	mu_r * sinTheta * cosOmega - mu * sinThr)
371	muPrime = sinTheta * sinThr * cosOmega + mu * mu_r	1✔
372
373	return muPrime, aPrime	1✔
374
375
376	def yearday(year, month, day):	1✔
377	"""
378	yearday returns the yearday for the given date. yearday is
379	the 'day of the year', sometimes called (incorrectly) 'julian day'.
380
381	Args:
382	year
383	month
384	day
385
386	Returns:
387	yday: day of year
388
389	"""
390
391	assert(year >= 0)	1✔
392	assert(1 <= month and month <= 12)	1✔
393	assert(1 <= day and day <= numdays(year, month))	1✔
394
395	yday = day	1✔
396
397	# Add in number of days for preceeding months.
398	for mon in range(month-1, 0, -1):	1✔
399	yday += numdays(year, mon)	1✔
400
401	return yday	1✔
402
403
404	def numdays(year, month):	1✔
405	"""
406	numdays returns the number of days in the given month of
407	the given year.
408
409	Args:
410	year
411	month
412
413	Returns:
414	ndays: number of days in month
415	"""
416
417	NDAYS = list([31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31])	1✔
418
419	assert(year >= 0)	1✔
420	assert(1 <= month and month <= 12)	1✔
421
422	ndays = NDAYS[month-1]	1✔
423
424	# Check for leap year for February
425	if ((month == 2) and leapyear(year)):	1✔
426	ndays += 1	1✔
427
428	return ndays	1✔
429
430
431	def leapyear(year):	1✔
432	"""
433	leapyear determines if the given year is a leap year or not.
434	year must be positive, and must not be abbreviated; i.e.
435	89 is 89 A.D. not 1989.
436
437	Args:
438	year
439
440	Returns:
441	True if a leap year, False if not a leap year
442	"""
443
444	assert(year >= 0)	1✔
445
446	if ((year % 4 == 0 and year % 100 != 0) or year % 400 == 0):	1✔
447	return True	1✔
448	else:
449	return False	1✔

USDA-ARS-NWRC / smrf / #179

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous