```import math
import rhinoscriptsyntax as rs
import solarGeom as sg

#Midterm project sample code

#1. Capturing the buiding geometry
# This example takes an elliptical plan for each floor and generates a set of
# points where each floor have a point distribution that is a rotated version
# of the floor before.  This is captured as a double list, with the first index
# cycling over the number of floors, and second index capturing the points
# on each floor

def getPointsTower(numFloors, numPoints):
interval = 2*math.pi/numPoints
rot = interval/7
height = 200/numFloors
pOut = []
for floor in range(numFloors):
pOut.append([])
for i in range(numPoints):
t = i*2*math.pi/numPoints+floor*rot
x = 32.5*math.cos(t)
y = 55*math.sin(t)
z = floor*height
pOut[floor].append([x,y,z])
return pOut

tower = getPointsTower(30, 10)

#for floor in range(len(tower)):

# 2. Performing a solar incidence analysis on the geometry
# For now, let's just try to do this for a fixed hour of a given day
latitude = 40.75
longitude = -73.5
timezone = -5.0
day = sg.calc_dayOfYear("3/23")
hour = sg.calc_hourDecimal("13:00")
sun = sg.calc_sunVector(latitude, longitude, timezone, day, hour)

# Additionally, we will use the placeholder function which returns a normal intensity of radiation
#for an input of day in the year and hour of the day
# For your project you will be using the data from the file dirNormals.txt
return 1000*max(0, -math.cos(hour*math.pi/12) + 1/8.0 - (day-171)**2/115600)

intensity_n = directNormalIrad(day, hour)
print ("for hour "+str(hour)
+", directbeamintensity = "+str(intensity_n)
)

#3.  Designing fins based on the analysis -- making decisions about the
# number, shape and location
# As a preliminary step, we will just represent the simple analysis that we
# did in part 2, and generate a vector that is normal to each plane proportional
# to the beam intensity on the plane.

#Useful to visualize a vector - can get this function in the Rhino-Python primer
if base_point==None:
base_point = [0,0,0]
tip_point = rs.PointAdd(base_point, vecdir)
line = rs.AddLine(base_point, tip_point)
if line:
return rs.CurveArrows(line, 2)

def visualizeBeamIntensity(towerIn, sunIn, intensityIn):
sunUnit = rs.VectorUnitize(sunIn)
for floor in range(len(towerIn)):
numPoints = len(towerIn[floor])
for i in range(numPoints):
p1 = towerIn[floor][i]
p2 = towerIn[floor][(i+1)%numPoints]
v1 = rs.VectorSubtract(p2, p1)
v2 = [0,0,1]
n = rs.VectorCrossProduct(v1, v2)
if rs.VectorLength(n) > rs.UnitAbsoluteTolerance(10**(-3), True):
nUnit = rs.VectorUnitize(n)
cosTheta = rs.VectorDotProduct(nUnit, sunUnit)
if cosTheta > 0:
intensityIs = intensityIn*cosTheta
vBeam = rs.VectorScale(nUnit, intensityIs/10)

visualizeBeamIntensity(tower, sun, intensity_n)

# or I can visualize the envelope of sun coverage by visualizing the beam intensity on the tower
# throughout the "working" day
hourStart = sg.calc_hourDecimal("7:30")
hourEnd = sg.calc_hourDecimal("18:00")
R = 20
N = 20 #number of sample points between starting time and end time
for i in range(N+1):
h = hourStart+i*(hourEnd-hourStart)/N
sunNow = sg.calc_sunVector(latitude, longitude, timezone, day, h)
if sunNow[2]>0:
I_n = directNormalIrad(day, h)
visualizeBeamIntensity(tower, sunNow, I_n)```