what field to use on the Sun?
The sun.alt is correct. alt is an altitude above horizon; together with an azimuth east of north they define an apparent position relative to horizon.
Your calculations are almost correct. You've forgot to supply an observer: sun = ephem.Sun(o).
- I don't know how to interpret negative results from cot(phi). Can
someone help me?
The Sun is below horizon in this case.
Finally, I'm confused about how to use
PyEphem to work backwards from a
shadow length to the next time when
the sun will cast a shadow of that
length, given an ephem.Observer().
Here's a script that given a function: g(date) -> altitude computes the next time when the sun will cast a shadow with the same length as right now (an azimuth -- direction of the shadow is not considered):
#!/usr/bin/env python
import math
import ephem
import matplotlib.pyplot as plt
import numpy as np
import scipy.optimize as opt
def main():
# find a shadow length for a unit-length stick
o = ephem.Observer()
o.lat, o.long = '37.0625', '-95.677068'
now = o.date
sun = ephem.Sun(o) #NOTE: use observer; it provides coordinates and time
A = sun.alt
shadow_len = 1 / math.tan(A)
# find the next time when the sun will cast a shadow of the same length
t = ephem.Date(find_next_time(shadow_len, o, sun))
print "current time:", now, "next time:", t # UTC time
####print ephem.localtime(t) # print "next time" in a local timezone
def update(time, sun, observer):
"""Update Sun and observer using given `time`."""
observer.date = time
sun.compute(observer) # computes `sun.alt` implicitly.
# return nothing to remember that it modifies objects inplace
def find_next_time(shadow_len, observer, sun, dt=1e-3):
"""Solve `sun_altitude(time) = known_altitude` equation w.r.t. time."""
def f(t):
"""Convert the equation to `f(t) = 0` form for the Brent's method.
where f(t) = sun_altitude(t) - known_altitude
"""
A = math.atan(1./shadow_len) # len -> altitude
update(t, sun, observer)
return sun.alt - A
# find a, b such as f(a), f(b) have opposite signs
now = observer.date # time in days
x = np.arange(now, now + 1, dt) # consider 1 day
plt.plot(x, map(f, x))
plt.grid(True)
####plt.show()
# use a, b from the plot (uncomment previous line to see it)
a, b = now+0.2, now+0.8
return opt.brentq(f, a, b) # solve f(t) = 0 equation using Brent's method
if __name__=="__main__":
main()
Output
current time: 2011/4/19 23:22:52 next time: 2011/4/20 13:20:01
