{{{id=34|
%python
from sunmodel import *
import sys

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 null 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, 13, 00, 1) # y m d h min s
end = datetime(2011, 5, 25, 13, 00, 0)
step = '5d'
///
}}}

{{{id=65|
%python
from sage.plot.plot3d.shapes2 import frame3d

def build_grid(x_values, y_values, x_res, y_res):
    '''
    returns a list of points, given min/max values and a resolution (# of spaces) for x and y
    first two arguments are (min,max) for x and y, next two are resolutions for x, y
    '''
    x_spacing = (x_values[1] - x_values[0]) / x_res
    y_spacing = (y_values[1] - y_values[0]) / y_res
    data_points = []
    for j in range(y_res+1):
        for i in range(x_res+1):
            data_points.append((x_values[0] + x_spacing*i, y_values[0] + y_spacing*j, 0))
    return data_points
    
def plot_wall(wall):
    min = wall['min']
    max = wall['max']
    
    points = [min]
    points.append((min[0], min[1], max[2]))
    points.append(max)
    points.append((max[0], max[1], min[2]))
    return polygon3d(points, color='blue', opacity=.5)

def get_point_illumination(p, sun_vector, walls):
    'returns the illumination value for this point, given a sun vector and walls'
    for w in walls:
        intersect = get_shadow_vector(vector(w["normal"]), vector(w["min"]), sun_vector, vector(p)).list()
        for i in range(3):
            # if the intersection is outside the wall along this axis, the intersection can't be within the wall
            if intersect[i] > w['max'][i] or intersect[i] < w['min'][i]:
                break # move on to the next wall
            if i == 2:
                # the plane intersection is within the boundaries on all axes, since we haven't hit the BREAK yet
                # this means that the intersection is within this wall, and the point is not illuminated
                return 0
    return 1

def get_total_point_illumination(p, sun_vectors, walls):
    'returns the total point illumination given a point, a list of sun vectors, and a list of walls'
    value = 0
    for v in sun_vectors:
        value += get_point_illumination(p, v, walls)
    return value

def get_illumination_map(points, sun_vectors, walls):
    'returns a list of illumination values for every point in points, given sun vectors and a list of walls'
    return [get_total_point_illumination(p, sun_vectors, walls) for p in points]
    
def plot_illumination(points, sun_vectors, walls, labels=False, point_size=10):
    'show a plot of illumination values for the given list of points, sun vectors, and walls'
    g = text3d('UW Farm', (-10,-10,0))
    for w in walls:
        g += plot_wall(w)
       
    map = get_illumination_map(points, sun_vectors, walls) # list of illumination values for every point in points
    min_val = map[0]
    max_val = map[0]
    for val in map:
        min_val = min(val, min_val)
        max_val = max(val, max_val)
        
    for i in xrange(len(points)):
        val = map[i]
        # red-yellow is harder to see.  we'll use white-black
        # green is 0 for max_val (makes red) and 1 for min_val (makes yellow)
        #green = 1 - float(val - min_val)/(max_val - min_val)
        
        # 0 for min, 1 for max
        color = float(val - min_val)/(max_val - min_val)
        
        # plot a 3d point at the i'th point's location with the color determined above
        g += point3d(points[i], color=(color, color, color), size=point_size)
        
        if labels:
            g += text3d(str(val), (points[i][0], points[i][1], 10)) # text is 1 unit above point
     
    g.show(aspect_ratio=1)
    
def get_sun_vectors(data):
    #                                      v North is 37.5 degrees off from +y in our plot
    return [sun_vector(r_from_d(d['az'] - 37.5), r_from_d(d['el'])) for d in data]
///
}}}

<p>Make the walls and points</p>

{{{id=73|
%python
walls = [\
    # foege
    {"normal": (1,0,0), "min":(137,43,0), "max":(137,318,60)},\
    {"normal": (0,1,0), "min":(79,318,0), "max":(137,318,60)},\
    {"normal": (1,0,0), "min":(79,318,0), "max":(79,369,60)},\
    {"normal": (0,-1,0), "min":(79,369,0), "max":(137,369,80)},\
    {"normal": (1,0,0), "min":(137,369,0), "max":(137,586,80)},\
    # ocean
    {"normal": (-1,0,0), "min":(251,0,0), "max":(251,190,70)},\
    {"normal": (0,1,0), "min":(251,190,0), "max":(353,190,70)},\
    # hitchcock
    {"normal": (0,-1,0), "min":(354,367,0), "max":(467,367,60)},\
    {"normal": (-1,0,0), "min":(354,367,0), "max":(354,540,60)},\
    # east building
    {"normal": (-1,0,0), "min":(419,127,0), "max":(419,255,150)},\
    {"normal": (0,1,0), "min":(419,255,0), "max":(531,255,150)},\
    {"normal": (-1,0,0), "min":(531,255,0), "max":(531,438,150)}\
]

points = build_grid((140,250), (0,180), 6, 8) + build_grid((140,350), (200,580), 10, 20)
///
}}}

<p>Shadow of a single sun vector</p>

{{{id=75|
data = [{'az': 127, 'el':35}]
sun_vectors = get_sun_vectors(data)
plot_sun(data)
///
}}}

{{{id=74|
plot_illumination(points, sun_vectors, walls, labels=False)
///
}}}

<p>Overlay of illumination and campus map</p>
<p><img src="http://students.washington.edu/bigmack/math480/overlay.jpg" alt="" /></p>

<p>Illumination during summer solstice</p>

{{{id=76|
start = datetime(2010, 6, 21, 6, 00, 1) # y m d h min s
end = datetime(2010, 6, 21, 20, 00, 0)
step = '20m'
data = get_angle_data(start, end, step, utc_offset=int(-7))
sun_vectors = get_sun_vectors(data)
plot_sun(data)
///
No location specified.  Defaulting to (lat, long) = (47, 238)
Data received
}}}

{{{id=63|
plot_illumination(points, sun_vectors, walls, labels=False)
///
}}}

<p>Illumination during winter solstice</p>

{{{id=77|
start = datetime(2010, 12, 21, 6, 00, 1) # y m d h min s
end = datetime(2010, 12, 21, 20, 00, 0)
step = '20m'
data = get_angle_data(start, end, step, utc_offset=int(-7))
sun_vectors = get_sun_vectors(data)
plot_sun(data)
///
Data received
}}}

{{{id=66|
plot_illumination(points, sun_vectors, walls, labels=False)
///
}}}

<p>Illumination during spring equinox</p>

{{{id=71|
start = datetime(2010, 3, 20, 6, 00, 1) # y m d h min s
end = datetime(2010, 3, 20, 20, 00, 0)
step = '20m'
data = get_angle_data(start, end, step, utc_offset=int(-7))
sun_vectors = get_sun_vectors(data)
plot_sun(data)
///
Data received
}}}

{{{id=80|
plot_illumination(points, sun_vectors, walls, labels=False)
///
}}}

{{{id=72|

///
}}}

{{{id=85|

///
}}}