4. Apparent AZ & EL:
        Apparent azimuth and elevation of target. Corrected for light-time,
     the gravitational deflection of light, stellar aberration, precession and
     nutation. There is an optional (approximate) correction for atmospheric
     refraction (Earth only). Azimuth measured North(0) -> East(90) ->
     South(180) -> West(270). Elevation is with respect to plane perpendicular
     to local zenith direction. TOPOCENTRIC ONLY. Units: DEGREES

        Labels: Azi_(a-appr)_Elev (airless)
                 Azi_(r-appr)_Elev (refracted)

{{{id=34| %python from sunmodel import * import sunmodel reload(sunmodel) # reload, if the attached .py file has changed spm = SolarPositionModel('horizons', debug=False) def get_angle_data(start, end, step, location=None, utc_offset=None): if location: spm.set_location(location[0], location[1]) if utc_offset: spm.set_utc_offset(utc_offset) return spm.get_angles(start, end, step) /// Connecting to Horizons... Connected }}} {{{id=43| # Return vector of stick given its angular definition def polar_to_cartesian(az, el, length = 1): x = length * float(sin(az)*cos(el)) y = length * float(cos(az)*cos(el)) z = length * float(sin(el)) return [x,y,z] # Takes an angle in degrees and returns a float representing radians def r_from_d(angle): return float(angle/360*2*pi) # Returns vector of the stick given polar definition def stick_vector(az, el, length, base_vector = vector([0,0,0])): return base_vector + vector(polar_to_cartesian(az, el, length)) # Returns vector of sun given azimuth and elevation def sun_vector(az, el): return vector(polar_to_cartesian(az, el)) # Takes four sage Vector objects. The 'normal' vector is normal to the # plane that one would like to cast a shadow on. The 'anchor' vector # is the initial point of the normal vector with respect to the horizon # coordinate system. The 'sun' vector is one that points to the sun from # origin of the horizon plane. The 'stick' vector is the vector pointing # from the horizon origin to the tip of the stick. def get_shadow_vector(normal, anchor, sun, stick): sun_and_stick = sun + stick # Numerator variable num = normal.dot_product(anchor - stick) # Denominator variable den = normal.dot_product(sun_and_stick - stick) u = num / den return stick + u * (sun_and_stick - stick) # Takes normal and anchor vectors to define shadow plane, as well as the # sun vector pointing to the sun to cast projections from both endpoints # of our shadow casting stick defined by the vectors stick_start and # stick_end def get_stick_endpoint_shadows(normal, anchor, sun, stick_start, stick_end): # Check to see if the Sun is below the horizon and return None if # it is because there would be no shadows if sun[2] <= 0: return None shadow_start = get_shadow_vector(normal, anchor, sun, stick_start) shadow_end = get_shadow_vector(normal, anchor, sun, stick_end) return (shadow_start, shadow_end) /// }}} {{{id=10| # this will make a 3D graph of the local normal, local north, and the apparent direction of the sun def plot_sun(data): var('u v') ground = parametric_plot3d((u, v, 0), (u, -1, 1), (v, -1, 1), opacity=0.5) all = default_plot() + ground for i, d in enumerate(data): az = r_from_d(d['az']) # radians from degrees el = r_from_d(d['el']) color_shift = 1 - float(i)/len(data) # start at 1, go to 0 by the end opac = 1 if el > 0 else .2 sun_line = line3d([(0,0,0), (0,1,0)], arrow_head=True, color=(1, color_shift, 0), opacity=opac) # start with a line going straight north sun_line = sun_line.rotateX(-el).rotateZ(az) # print str(d['datetime']), d['az'], d['el'] all += sun_line all.show(aspect_ratio=[1,1,1], frame=False) /// }}} {{{id=37| def default_plot(): Z = line3d([(0,0,0), (0,0,1)], arrow_head=True, color='black') north = line3d([(0,0,0), (0,1,0)], arrow_head=True, color='green') # north is in +y direction, normal is in +z direction (east is in +x direction) Z_text = text3d("Z", (0,0,1.1)) north_text = text3d("North", (0, 1.1, 0)) return Z + north + Z_text + north_text def show_projection(stick_base, stick_tip, normal, anchor, data, show_stick=True, normalize=False, points=False): # Generate 'base' plot g = default_plot() if show_stick: g += line3d([stick_base.list(), stick_tip.list()], color='red', thickness=2) # Plots shadows for each moment in time specified by data for i, d in enumerate(data): # Convert Horizons az. and el. to radians az = r_from_d(d['az']) el = r_from_d(d['el']) sun = sun_vector(az, el) result = get_stick_endpoint_shadows(normal, anchor, sun, stick_base, stick_tip) if result == None or normal.dot_product(sun) <= 0: continue shadow_start, shadow_end = result if normalize: if shadow_start.norm() > 0: shadow_start = shadow_start / shadow_start.norm() if shadow_end.norm() > 0: shadow_end = shadow_end / shadow_end.norm() color_shift = .8 - float(i)/len(data)*.8 # start at .8, go to 0 by the end color = (color_shift, color_shift, color_shift) if points: g += point3d(shadow_start.list(), color=color) g += point3d(shadow_end.list(), color=color) g += line3d([shadow_start.list(), shadow_end.list()], color=color) g.show(aspect_ratio=[1,1,1]) /// }}} {{{id=53| %python # Initializing parameters for casting shadows stick_base = vector([0,0,0]) stick_tip = stick_vector(0, pi/2, 1, stick_base) normal = vector([0,0,1]) anchor = vector([0,0,0]) start = datetime(2010, 5, 25, 12, 00, 1) # y m d h min s end = datetime(2011, 5, 25, 12, 00, 0) step = '5d' /// }}}

Sun at noon every 2 days for 1 year starting May 25, 2010

Analemma

{{{id=46| data = get_angle_data(start, end, step, utc_offset=int(-8)) plot_sun(data) /// No location specified. Defaulting to (lat, long) = (47, 238) Data received }}}

Projection of the analemma in Seattle from a stick pointing straight up

{{{id=55| show_projection(stick_base, stick_tip, normal, anchor, data, normalize=False, points=True) /// }}}

Sun data from Horizons - today, every 10 minutes from 7 AM to 5 PM

{{{id=57| start = datetime(2010, 6, 2, 7, 00, 1) # y m d h min s end = datetime(2010, 6, 2, 17, 00, 0) step = '10m' data = get_angle_data(start, end, step, utc_offset=int(-8)) plot_sun(data) /// Data received }}} {{{id=58| stick_base = vector([0,-.5,1]) stick_tip = vector([0,.5,1]) #stick_vector(0, pi/2, 1, stick_base) show_projection(stick_base, stick_tip, normal, anchor, data, normalize=False, points=False) /// }}} {{{id=59| stick_base = vector([-.5,0,1]) stick_tip = vector([.5,0,1]) #stick_vector(0, pi/2, 1, stick_base) show_projection(stick_base, stick_tip, normal, anchor, data, normalize=True, points=False) /// }}} {{{id=60| normal = vector([0,-1,0]) anchor = vector([0,0,0]) stick_base = vector([0,0,0]) stick_tip = vector([0,-1,0]) #stick_vector(0, pi/2, 1, stick_base) show_projection(stick_base, stick_tip, normal, anchor, data, normalize=True, points=False) /// }}}

Model of the sundial on the Physics and Astronomy building today

{{{id=61| normal = vector(polar_to_cartesian(7*pi/6, 0)) anchor = vector([0,0,0]) stick_base = vector([0,0,0]) stick_tip = stick_vector(pi, -pi/4, 1, stick_base) show_projection(stick_base, stick_tip, normal, anchor, data, normalize=False, points=False) /// }}}

Model of the sundial on the Physics and Astronomy building on 3/20/2010 (spring equinox)

{{{id=62| start = datetime(2010, 3, 20, 10, 00, 1) # y m d h min s end = datetime(2010, 3, 20, 18, 00, 0) step = '1m' data = get_angle_data(start, end, step, utc_offset=int(-8)) show_projection(stick_base, stick_tip, normal, anchor, data, normalize=False, points=False) /// Data received }}}

Equatorial sundial - today from 6 AM to 6 PM every hour

{{{id=67| normal = vector(polar_to_cartesian(0, r_from_d(47))) anchor = vector([0,0,0]) stick_base = vector(polar_to_cartesian(0, r_from_d(47))) stick_tip = stick_vector(0, r_from_d(47), 1, stick_base) start = datetime(2010, 6, 2, 6, 00, 1) # y m d h min s end = datetime(2010, 6, 2, 18, 00, 2) step = '1h' data = get_angle_data(start, end, step, utc_offset=int(-8)) show_projection(stick_base, stick_tip, normal, anchor, data, normalize=False, points=True) /// Data received }}}

Horizontal sundial today

{{{id=70| normal = vector((0,0,1)) anchor = vector((0,0,0)) stick_base = vector((0,0,0)) stick_tip = stick_vector(0, r_from_d(47), 1, stick_base) show_projection(stick_base, stick_tip, normal, anchor, data, normalize=True, points=True) /// }}}

Vertical sundial today

{{{id=72| normal = vector([0,-1,0]) anchor = vector(polar_to_cartesian(0, r_from_d(47))) show_projection(stick_base, stick_tip, normal, anchor, data, normalize=False, points=False) /// }}}

Polar sundial today

{{{id=71| normal = vector(polar_to_cartesian(0, r_from_d(47) + pi/2)) anchor = vector((0,0,-1)) show_projection(stick_base, stick_tip, normal, anchor, data, normalize=False, points=True) /// }}} {{{id=77| DAYS_IN_YEAR = 365.25636 SECONDS_IN_YEAR = DAYS_IN_YEAR*24*60*60 def earth_location(t): ''' t is the time past since perihelion in days returns an (r, theta) pair of the earth's position relative to the sun (add pi to get sun relative to earth) r isn't required for the sundial, but it's used if we want to plot the path of the earth in polar coords this follows the steps found here: http://en.wikipedia.org/wiki/Kepler's_laws_of_planetary_motion#Position_as_a_function_of_time ''' var('E th') a = 1.4960e8 # semi-major axis in km - half the measurement of the orbit at it's widest part ec = 0.0167 # eccentricity - a really accurate value of this would be nice M = 2*pi*t/DAYS_IN_YEAR # mean anomaly ecc_anomaly = find_root(M == E-ec*sin(E), 0, 2*pi) # eccentric anomaly true_anomaly = 2*arctan(sqrt((1+ec)/(1-ec))*tan(E/2)).subs(E=ecc_anomaly) r = (a*(1-ec^2)/(1+ec*cos(th))).subs(th=true_anomaly) return (r, true_anomaly) def graph_earth(t): var('x y') r, th = earth_location(t) # neat python trick to assign variables to entries in a tuple # get x and y values to plot on a graph x = float((r*cos(th)).subs(th=th)) y = float((r*sin(th)).subs(th=th)) g = point((x, y), color='green') + point((0,0), color='orange') # adding two points together makes them part of the same plot g.set_axes_range(-1.4960e8, 1.4960e8, -1.4960e8, 1.4960e8) # fix the axes so we don't zoom in and out g.set_aspect_ratio(1) return g /// }}} {{{id=19| # this will make an animation of the orbit, 1 frame per earth day # it takes a couple minutes to render, though animate([graph_earth(t) for t in xrange(0, 365)]).show(delay=1) /// }}} {{{id=20| DAYS_IN_YEAR = 365.242199 # shows the earth's position given a day of the year @interact def show_earth(t=(0,365)): print 'percentage of year:' print t/DAYS_IN_YEAR*100 graph_earth(t).show(aspect_ratio=1) /// }}} {{{id=21| # shows the earth's position at certain points in the year # you can see that 1/4 of the way through the year, the earth has actually moves more than 90 degrees! show_earth(DAYS_IN_YEAR/4) show_earth(DAYS_IN_YEAR/2) show_earth(DAYS_IN_YEAR*3/4) /// percentage of year: 25.0000000000000

percentage of year: 50.0000000000000

percentage of year: 75.0000000000000

}}} {{{id=80| /// }}} {{{id=82| /// }}} {{{id=83| /// }}}