Compute shadow length using PyEphem

Question

I am using PyEphem and want to calculate the length of a shadow (assume a stick of unit length is planted in the ground). The length will be given by cot(phi), where phi is the solar elevation angle (please correct me if I'm wrong). I'm not sure what field to use on the Sun? In the example below, I'm using the angle alt:

import ephem, math
o = ephem.Observer()
o.lat, o.long = '37.0625', '-95.677068'
sun = ephem.Sun()
sunrise = o.previous_rising(sun, start=ephem.now())
noon = o.next_transit(sun, start=sunrise)
shadow = 1 / math.tan(sun.alt)

Please check my interpretation below:

1. If the tangent is infinite, it indicates the sun is directly overhead and there is no shadow.
2. If the tangent is zero, it indicates that the sun is at the horizon and the shadow is infinitely long.
3. I don't know how to interpret negative results from cot(phi). Can someone help me?

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().

I would appreciate help with this.

Answer 1

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).

3. 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