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)):
#	rs.AddPoints(tower[floor])
 
# 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
def directNormalIrad(day, hour):
	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
def AddVector(vecdir, base_point):
	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)
					AddVector(vBeam, p1)
 
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)