# Name :          Simple Solar Study
# Description :   Calculates amount of direct sunlight hitting a target surface as a percentage
#                 of total sunlight available under open sky. Minimum analysis time is 24 hours, 
#                 maximum is the whole year. Solar angle resolution is per hour; target surface resolution
#                 determined by size and number of faces to analyze. Results are visually encoded by
#                 shading the faces and by adding daylight factor label to faces. Results can also
#                 be saved to a comma-separated-value (text) file
# Author :        Brian Monwai
# Usage :         1. Geo-reference the model
#                 2. Create single group containing faces to analyze.
#                 3. a - Faces that provide shade should be placed on their own layer
#                    b - Glass/partial transmittance faces should be on their own layer (adv mode)
#                 4. Select the group containing faces to analyze
#                 5. Run Simple Solar study from Plugins menu
#						6. Select start/end dates (all dates are inclusive), time zone, and any other options and press OK!
#
# History:        Rev  | Date       | Notes
#                 -----|------------|-----------------------------------
#                 0.98 | 2013-08-11 | Add ability to read weather data from TMY3 file 
#                      |            | Bug: in fine/tmy3 mode, adjust for surface tilt based on beam strength on horiz. surface, not to normal-to-beam surface
#                 0.97 | 2013-08-07 | Bug: remove Numeric class overloads, make help file search relative
#                 0.96 | 2013-08-06 | Add 'help' and 'about' documentation
#                      |            | Add no-gradient option
#                      |            | Add runtime to output msgbox
#                 0.95 | 2013-07-14 | Add undo, update status bar messages, add small circle at sunrise position
#                      |            | Debug parms to name solar path layer
#                      |            | Bug: fix translation from solar azimuth to x,y,z coords, get Bird/Hulstrum (Calculation mode = "fine") working
#                      |            | New feature: support tilted surfaces, output in W/m^2
#                 0.94 | 2013-06-20 | New feature: support for glass/transmittance layer (advanced only)
#                      |            | Enhancement: Allow multiple occluding layers (advanced only)
#                      |            | Add debugParms: messageLevel, LabelLayerName, lat/long, etc.
#                      |            | Move sun path and save results to advanced mode
#                 0.93 | 2013-03-17 | Restructure code to reduce beam intensity based on sun elevation
#                      |            | Add solar insolation (Bird and Hulstrom) calculations
#                      |            | Output analysis face entity ID and model units to csv file
#                      |            | Trace sun path instead of drawing rays from sun
#                 0.92 | 2013-02-24 | Use a point from selection as origin instead of (0,0,0)
#                      |            | Ignore hidden faces for occlusion test
#                 0.91 | 2013-02-23 | Use status bar for progress
#                      |            | Check for nil faces intersection
#                      |            | Fix leap year calc every 400 yrs
#                 0.90 | 2013-02-20 | First "public" release
#
#

=begin
Simple Solar Study
	Copyright (c) 2013 Brian Monwai
	
	THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
	IMPLIED. THE AUTHOR IS NOT LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
	LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
	OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
	THE SOFTWARE.

	Sources:
		solrad_ver16.xls (Greg Pelletier, Washington State Dept. of Ecology), who converted NOAA's web based solar calculator from Javasript to Excel macro)
		NOAA: http://www.esrl.noaa.gov/gmd/grad/solcalc/   see Javascript source, as well as their own excel versions here: http://www.esrl.noaa.gov/gmd/grad/solcalc/calcdetails.html
		Bird and Hulstrum (available here: )
		Meeus 
		solpos (C code from National Renewable Energy Laboratory, 25 Mar 1998)
		versions before 0.93 based on Sunposition v1.2.1, 2009, by Gabriel Miller

		Face centroid:
		thanks to Fredo http://sketchucation.com/forums/viewtopic.php?f=323&t=13378&p=430718&hilit=centroid#p430718
		and Wikipedia http://en.wikipedia.org/wiki/Centroid#Centroid_of_polygon
		
=end

require 'sketchup.rb'

module BM_SimpleSolarStudy

	# ---------------------
	# CONSTANTS
	# ---------------------
	VERSION = '0.98'
	
	unless file_loaded?(__FILE__)
		pm = UI.menu('Plugins').add_submenu('Simple solar study')
		if (pm)
			pm.add_item('Basic mode') { sss_main_proc(true) }   # 0.93
			pm.add_item('Advanced mode') { sss_main_proc(false) }     
			pm.add_item('Help') { sss_help_proc }
			pm.add_item('About') { sss_about_proc }
		end
	end
	file_loaded(__FILE__)


# ----------------------------------------------------------------------
# Bounding box for optimizing sun ray intersections
# ----------------------------------------------------------------------
class SSS_BoundingBox   # 1.10

	def initialize
		@bb = Geom::BoundingBox.new
	end # initialize
	
	# ---------------------------------
	# Public functions
	# ---------------------------------	
	def add(aface)
		@bb.add(aface.bounds)		
	end
	
	def isValid()
		return @bb.valid?
	end
	
	def calcBBPlanes()
		vx = Geom::Vector3d.new 1,0,0
		vy = Geom::Vector3d.new 0,1,0
		vz = Geom::Vector3d.new 0,0,1
		
		@bottom_plane = [Geom::Point3d.new(@bb.corner(0)), vz]    # specify with point and normal to plane
		@top_plane =    [Geom::Point3d.new(@bb.corner(4)), vz]
		@front_plane =  [Geom::Point3d.new(@bb.corner(0)), vy]
		@back_plane =   [Geom::Point3d.new(@bb.corner(2)), vy]
		@left_plane =   [Geom::Point3d.new(@bb.corner(2)), vx]
		@right_plane =  [Geom::Point3d.new(@bb.corner(1)), vx]
	end
	
	def inShadow(sunPt, analysisPt)
	
		startpos = analysisPt
		endpos = sunPt
		vector = endpos-startpos	
	
		line = [ startpos, vector]
	
		intpt_bot = Geom.intersect_line_plane(line, @bottom_plane)   
		return true if !intpt_bot.nil? && @bb.contains?(intpt_bot)		
		
		intpt_top = Geom.intersect_line_plane(line, @top_plane)   
		return true if !intpt_top.nil? && @bb.contains?(intpt_top)		
		
		intpt_frn = Geom.intersect_line_plane(line, @front_plane)   
		return true if !intpt_frn.nil? && @bb.contains?(intpt_frn)		
		
		intpt_bac = Geom.intersect_line_plane(line, @back_plane) 
		return true if !intpt_bac.nil? && @bb.contains?(intpt_bac)				
		
		intpt_lft = Geom.intersect_line_plane(line, @left_plane) 
		return true if !intpt_lft.nil? && @bb.contains?(intpt_lft)				
				
		intpt_rgt = Geom.intersect_line_plane(line, @right_plane)
		return true if !intpt_rgt.nil? && @bb.contains?(intpt_rgt)						
		
		#if intpt_bot.nil? && intpt_top.nil? && intpt_frn.nil? && intpt_bac.nil? && intpt_lft.nil? && intpt_rgt.nil?
		#	return false  # no intersection
		#end
	   #
		#if @bb.contains?(intpt_bot) || @bb.contains?(intpt_top) || @bb.contains?(intpt_frn) || @bb.contains?(intpt_bac) || @bb.contains?(intpt_lft) || @bb.contains?(intpt_rgt)
		#	return true  # intersection within bounding box
		#end
		
		return false
	end
	
	def drawBBPlanes(entities)
		entities.add_edges Geom.intersect_plane_plane(@bottom_plane, @front_plane)
		entities.add_edges Geom.intersect_plane_plane(@top_plane, @front_plane)
		entities.add_edges Geom.intersect_plane_plane(@bottom_plane, @back_plane)
		entities.add_edges Geom.intersect_plane_plane(@top_plane, @back_plane)
	end
	
	def drawBBFaces(entities)
		@frontFace  = entities.add_face @bb.corner(0), @bb.corner(1), @bb.corner(4), @bb.corner(5)   # front
		@backFace   = entities.add_face @bb.corner(2), @bb.corner(3), @bb.corner(6), @bb.corner(7)   # back
		@bottomFace = entities.add_face @bb.corner(0), @bb.corner(1), @bb.corner(2), @bb.corner(3)   # bottom
		@topFace    = entities.add_face @bb.corner(4), @bb.corner(5), @bb.corner(6), @bb.corner(7)   # top
		@leftFace   = entities.add_face @bb.corner(0), @bb.corner(2), @bb.corner(6), @bb.corner(4)   # left
		@rightFace  = entities.add_face @bb.corner(1), @bb.corner(3), @bb.corner(7), @bb.corner(5)   # right
	end

	def drawBBEdges(entities)
		for i in 0..7
			for j in 0..7
				if i != j
					entities.add_edges @bb.corner(i), @bb.corner(j)
				end
			end
		end
	end
	
	# ---------------------------------
	# Private functions
	# ---------------------------------	
	private
	# nothing private!
	
end  # class

# ------------------------------------------------------------------------------------
# sunData class
# 
#   This class stores the sun position at a given point in time, as calculated in the 
#   sunData::calcSun function. Public class functions are used to find azimuth/elevation/etc.
#
# ------------------------------------------------------------------------------------
class SSS_sunData   # 0.93

	def initialize
		@pressureMb = 1013 # barometric pressure in millibars (1013 = sea level)
		@ozoneCm = 0.35    # ozone thickness of atmosphere in cm (typ 0.05 to 0.4)
		@waterCm = 4.0     # water vapor thickness of atmosphere in cm (typ 0.01 to 6.5)
		@aod500 = 0.35     # aerosol optical depth at 500nm (typ 0.02 to 0.5, values > 0.5 represent clouds and volcanic dust, etc.)
		@aod380 = 0.35     # aerosol optical depth at 380nm (typ 0.1 to 0.5)
		@fscatter = 0.84   # proportional (cosine) foward scattering of incoming radiation (typ 0.84.  1.0 for all outgoing, 0.0 for all backward scattering )
		@albedo = 0.20     # surface albedo or ground reflectance (typ 0.2 for land, 0.25 for vegetation, 0.9 for snow)
		@SOLAR_CONSTANT = 1367   # watts/m^2  NREL		
		@sundist = 10000   # for positioning sun
		
	end # initialize

	# ---------------------------------
	# Public functions
	# ---------------------------------	
	def isLeapYear(yr) 	
		return ((yr % 4 == 0 && yr % 100 != 0) || yr % 400 == 0)
	end
		
	# ------------------------------------------------------------------------------------
	# Calculate sun position
	#
   # INPUTS	
	#   lat: decimal latitude (pos = north)
	#   long: decimal longitude (neg = west)
	#   year: 4 digit year
	#   month: 1 = january
	#   day: 1 - 31
	#   hour: 0 = midnight, 23 = 11:00pm
	#   timezone: -12 to +12
	#
	# RETURN										
	#   Sun position as Point3d. Use public class functions to get azimuth/elevation/etc.
	# ------------------------------------------------------------------------------------
	def calcSun(lat, long, year, month, day, hour, timezone)
	
		# convert calendar date to Julian Day
		jd = calcJD(year, month, day)
		
		# convert to Julian century
		timenow = hour - timezone
		jc = calcTimeJulianCent(jd + timenow/24.0)
		
		# UI.messagebox(sprintf("Julian day=%.6f century=%.6f", jd, jc))  # DEBUG
		
		# Limit extreme N/S latitudes
		if (lat < -89.8)
			lat = -89.8   
		elsif (lat > 89.8)
			lat = 89.8
		end
		
		# Geometric mean longitude (degrees)		
		@gmls = 280.46646 + jc * (36000.76983 + 0.0003032 * jc)
		@gmls = @gmls % 360

		# Geometric mean anomaly (degrees)		
		@gmas = 357.52911 + jc * (35999.05029 - 0.0001537 * jc)
		
		# Eccentricity of earth's orbit (unitless)
		@ecc = 0.016708634 - jc * (0.000042037 + 0.0000001267 * jc)
		
		# Equation of center for the sun (degrees)
		mrad = @gmas.degrees		
		@eqctr = Math.sin(mrad) * (1.914602 - jc * (0.004817 + 0.000014 * jc)) + Math.sin(2 * mrad) * (0.019993 - 0.000101 * jc) + Math.sin(3 * mrad) * 0.000289
		
		# True longitude of the sun (degrees)
		@sunLong = @gmls + @eqctr
		
		# True anomaly of the sun (degrees)
		@sunAnom = @gmas + @eqctr
		
		# Distance to the sun (AUs)
		@sunRadius = (1.000001018 * (1 - @ecc ** 2)) / (1 + @ecc * Math.cos(@sunAnom.degrees))
		
		# Apparent longitude of the sun (degrees)
		omega = 125.04 - 1934.136 * jc
		@sunAppLong = @sunLong - 0.00569 - 0.00478 * Math.sin(omega.degrees)
		
		# Mean obliquity of the ecliptic (degrees)
		seconds = 21.448 - jc*(46.8150 + jc*(0.00059 - jc*(0.001813)))
		@meanOblique = 23.0 + (26.0 + (seconds/60.0))/60.0
		
		# Corrected obliquity of the ecliptic (degrees)	
		@corrOblique = @meanOblique + 0.00256 * Math.cos(omega.degrees)

		# Right ascension of the sun (degrees) 
		@rtAscension = (Math.atan2(Math.cos(@corrOblique.degrees) * Math.sin(@sunAppLong.degrees), Math.cos(@sunAppLong.degrees))).radians
		
		# Declination of the sun (degrees)
		@declin = (Math.asin(Math.sin(@corrOblique.degrees) * Math.sin(@sunAppLong.degrees))).radians
		
		# Equation of time (difference between true solar time and mean solar time) (minutes)
		y = Math.tan(@corrOblique.degrees / 2.0)
		y = y * y
		@eTime = (y * Math.sin(2 * @gmls.degrees) - 2 * @ecc * Math.sin(@gmas.degrees) + 4 * @ecc * y * Math.sin(@gmas.degrees) * Math.cos(2 * @gmls.degrees) - 0.5 * y * y * Math.sin(4 * @gmls.degrees) - 1.25 * @ecc * @ecc * Math.sin(2 * @gmas.degrees)).radians * 4.0
		
		# True solar time (minutes)
		@solarTime = (hour/24.0 * 1440 + @eTime + 4.0 * long - 60.0 * timezone) % 1440
		
		# Hour angle (degrees)
		if (@solarTime / 4) < 0
			@hrAngle = (@solarTime / 4) + 180
		else 
			@hrAngle = (@solarTime / 4) - 180
		end
		
		# Solar zenith (degrees)
		@zenith = (Math.acos(Math.sin(lat.degrees)*Math.sin(@declin.degrees) + Math.cos(lat.degrees)*Math.cos(@declin.degrees)*Math.cos(@hrAngle.degrees))).radians
		
		# Solar elevation (angle above horizon) (degrees)
		@elev = 90 - @zenith
		
		# Adjust elevation for atmospheric refraction
		if (@elev > 85.0)
         refractionCorrection = 0.0
      else
         te = Math.tan(@elev.degrees)
         if (@elev > 5.0)
            refractionCorrection = 58.1 / te - 0.07 / (te ** 3) + 0.000086 / (te ** 5)
			elsif (@elev > -0.575) 
            refractionCorrection = 1735.0 + @elev * (-518.2 + @elev * (103.4 + @elev * (-12.79 + @elev * 0.711)))
         else
            refractionCorrection = -20.774 / te
         end
         refractionCorrection = refractionCorrection / 3600
      end
		
		@elevCorrected = @elev + refractionCorrection		
		
		# Azimuth (degrees clockwise from north)
		@azimuth = Math.acos(((Math.sin(lat.degrees) * Math.cos(@zenith.degrees)) - Math.sin(@declin.degrees)) / (Math.cos(lat.degrees) * Math.sin(@zenith.degrees)))
		if @hrAngle > 0
			@azimuth = (@azimuth.radians + 180) % 360			
		else				
			@azimuth = (540 - @azimuth.radians) % 360
		end
				
		
		# -----------------------------
		# translate to x,y,z	
		# -----------------------------
		north = Sketchup.active_model.shadow_info["NorthAngle"]
		if north.nil?
			north = 0.0
		end 
		
		z = @sundist * Math.sin(@elevCorrected.degrees)
		x = @sundist * Math.cos(@elevCorrected.degrees) * Math.sin((@azimuth + north).degrees)   # 0.95 switch sin/cos since north is generally in y dir
		y = @sundist * Math.cos(@elevCorrected.degrees) * Math.cos((@azimuth + north).degrees)
		
		@sunpt = Geom::Point3d.new(x, y, z)
				
		return @sunpt
		
	end #calcSun
	
	
	#-------------------------------------------------------------------------------------
	# Calculate solar radiation using Bird-Hulstrom method. Run after calcSun.
   # Unless otherwise noted, this function is converted from the solrad_ver16.xls macro by
	# Greg Pelletier, Washington State Dept of Ecology
	#
   # INPUTS
	#    year: 4 digit year								
	#    month: January = 1								
	#    day: 1 - 31									
   #
	# RETURN										
	#   Global radiation on a horizontal surface in Watts/m^2. 
	#   Use public class functions to get component parts of insolation
	#--------------------------------------------------------------------------------------
	def calcSolarRadiation(year, month, day)
		
		# ETR - extraterrestrial (outside the earth's atmosphere) direct beam intensity
		# from BIRD_08_16_2012.xls by Daryl R. Myers, National Renewable Energy Laboratory
		doy = calcDayOfYear(month, day, year)
		@etr = @SOLAR_CONSTANT * (1.00011 + 0.034221 * Math.cos(2 * Math::PI * (doy-1)/365) + 0.00128 * Math.sin(2 * Math::PI * (doy-1)/365) + 0.000719 * Math.cos(2 * (2 * Math::PI * (doy-1)/365)) + 0.000077 * Math.sin(2 * (2 * Math::PI * (doy-1)/365)) )

		tauA = 0.2758 * @aod380 + 0.35 * @aod500
		
		if @zenith < 89 

			@airMass = 1 / (Math.cos(@zenith.degrees) + 0.15 / (93.885 - @zenith) ** 1.253)    # gp Rob Annear's review of literature changes 1.25 to 1.253
 
			#Greg Pelletier 13-Aug-2012 correction based on comment from Mike Mickeslon and double check original Bird model eqn
			@tRayleigh = Math::E ** (-0.0903 * (@pressureMb * @airMass / 1013) ** 0.84 * (1 + @pressureMb * @airMass / 1013 - (@pressureMb * @airMass / 1013) ** 1.01)) 
			@tOzone = 1 - 0.1611 * (@ozoneCm * @airMass) * (1 + 139.48 * (@ozoneCm * @airMass)) ** -0.3035 - 0.002715 * (@ozoneCm * @airMass) / (1 + 0.044 * (@ozoneCm * @airMass) + 0.0003 * (@ozoneCm * @airMass) ** 2)       #gp Rob changed -0.3034 to -0.3035
			@tGases = Math::E ** (-0.0127 * (@airMass * @pressureMb / 1013) ** 0.26)         #gp correction per NREL
			@tWater = 1 - 2.4959 * @airMass * @waterCm / ((1 + 79.034 * @waterCm * @airMass) ** 0.6828 + 6.385 * @waterCm * @airMass)
			@tAerosol = Math::E ** (-(tauA ** 0.873) * (1 + tauA - tauA ** 0.7088) * @airMass ** 0.918)
			taa = 1 - 0.1 * (1 - @airMass + @airMass ** 1.06) * (1 - @tAerosol)
			rs = 0.0685 + (1 - @fscatter) * (1 - @tAerosol / taa)
			id = 0.9662 * @etr * @tAerosol * @tWater * @tGases * @tOzone * @tRayleigh
			ias = @etr * Math.cos(@zenith.degrees) * 0.79 * @tOzone * @tGases * @tWater * taa * ((0.5 * (1 - @tRayleigh) + @fscatter * (1 - (@tAerosol / taa))) / (1 - @airMass + (@airMass) ** 1.02))   #gp Rob changed (0.5 to ((0.5 and TAA)) to TAA)))

		else
			@airMass = 0
			@tRayleigh = 0
			@tOzone = 0
			@tGases = 0
			@tWater = 0
			@tAerosol = 0
			taa = 0
			rs = 0
			id = 0
			ias = 0
		end

		if @zenith < 90
			idnH = id * Math.cos(@zenith.degrees)
		else
			idnH = 0
		end
		
		if @airMass > 0
			gh = (idnH + ias) / (1 - @albedo * rs)
		else
			gh = 0	
		end

		@directBeam = id    # on surface normal to beam
		@directHoriz = idnH
		@globalHoriz = gh
		@diffuseHoriz = gh - idnH
				
		return @globalHoriz
		
	end #calcSolarRadiation
	
	
	# **********************************************************************
	#  Public functions to set private class data  # 0.94
	#
		
	def SetPressure(v) @pressureMb = v end
	def SetOzone(v) @ozoneCm = v end
	def SetWater(v) @waterCm = v end
	def SetAerosol_500(v) @aod500 = v end
	def SetAerosol_380(v) @aod380 = v end
	def SetFwdScatter(v) @fscatter = v end
	def SetReflectance(v) @albedo = v end
	def setSunDistance(v) @sundist = v end
	
	
	# **********************************************************************
	#  Public functions to read data from class
	#
	
	# sun position info
	def getSunPosition() return @sunpt end
	def getAzimuth() return @azimuth end
	def getElevAngleNonRefracted() return @elev end
	def getElevAngle() return @elevCorrected end
	def getZenith() return @zenith end
	def getHrAngle() return @hrAngle end
	def getSolarTime() return @solarTime end
	def getEqnTime() return @eTime end
	def getDeclination() return @declin end
	def getRadius() return @sunRadius end
	def getEccentricity() return @ecc end
	def getGeometricMeanAnomaly() return @gmas end
	def getGeometricMeanLongitude() return @gmls end
	
	# solar insolation info
	def getDirectBeam() return @directBeam end
	def getDirectHoriz() return @directHoriz end
	def getGlobalHoriz() return @globalHoriz end
	def getDiffuseHoriz() return @diffuseHoriz end
	
	
	# ----------------------------------------------------------------------
	#  Calculate Julian day
	#
	# INPUTS
	#    year: 4 digit year								
	#    month: January = 1								
	#    day: 1 - 31									
   #
	# RETURN										
	#    Julian day
	# NOTES											
	#    Number is returned for start of day.  Fractional days (hours) are accounted for elsewhere
	# ----------------------------------------------------------------------
	def calcJD(year, month, day)	
		if (month <= 2)
			year -= 1
			month += 12
		end
		
		a = (year/100).to_int
		b = 2 - a + (a/4).to_int

		jd = (365.25*(year + 4716)).to_int + (30.6001*(month+1)).to_int + day + b - 1524.5
		return jd
	end
	
	# -------------------------------------------------------------------------
	# Convert Julian Day to centuries since J2000.0
	#			
	# INPUT
	#   jd: Julian Day 	
   #
	# RETURN										
	#   julian century
	# -------------------------------------------------------------------------
	def calcTimeJulianCent(jd)	
		t = (jd - 2451545.0)/36525.0
		return t
	end
	
	# -------------------------------------------------------------------------
	#  Find numerical day-of-year from month, day and leap year info  
	#
	#  INPUTS
	#    month: January = 1							
	#    day  : 1 - 31								
   #    year : 4 digit year
	#
	#  RETURN									
	#    The numerical day of year				
	# -------------------------------------------------------------------------
	def calcDayOfYear(mn, dy, year)
 	
		if isLeapYear(year)
			k = 1
		else
			k = 2
		end
		doy = ((275 * mn)/9).to_int - k * ((mn + 9)/12).to_int + dy -30;
		return doy
	end
	
	
	# ---------------------------------
	# All subsequent functions are private 
	# ---------------------------------
   private

	# nothing private!
	
end # class SSS_sunData    # 0.93
	
######################################################################

# ------------------------------------------------------------------------------------
# weatherFileTMY3 class
# 
#   This class reads the entire weather file and saves the solar radiation in an array
#   Note: Data for leap years is used in TMY calculations, but TMY3 file format does not list Feb 29
#
# ------------------------------------------------------------------------------------
class SSS_weatherFileTMY3   # 0.98

	def initialize
		@header = ""   # station code, name, state, timezone, latitude, longitude, elev
		@directBeam = Array.new
		@diffuseHoriz = Array.new
		@directHoriz = Array.new	
		@errorStr = ""
	end # initialize

	# ---------------------------------
	# Public functions
	# ---------------------------------	
	
	# read weather file and store data in array
	# return true if successful, false if not (in which case, 
	# user should read errorStr
	def readFile(fh)		
	
		@header = fh.gets   # first line: station ID
		lineArr = @header.split(',')
		if lineArr.length != 7 
			@errorStr = "Problem with weather file header"
			return false
		end
		
		fh.gets             # second line: column headers
		
		while line = fh.gets  # get data
			lineArr = line.split(',')
			if lineArr.length != 68
				@errorStr = "Problem with weather file data at line " + fh.lineno.to_s
				return false
			end
			wfDate = lineArr[0]
			wfTime = lineArr[1]
			wfDirectNorm = lineArr[7].to_i
			wfDiffuseHoriz = lineArr[10].to_i
			wfDirectHoriz = lineArr[4].to_i - wfDiffuseHoriz
			
			@directBeam.push(wfDirectNorm)
			@diffuseHoriz.push(wfDiffuseHoriz)
			@directHoriz.push(wfDirectHoriz)			
		end
		
		if fh.lineno != 8762
			@errorStr = "Only read " + fh.lineno.to_s + " lines in weather file instead of 8762"
			return false
		end
		
		return true
	end
	
	def getError()
		return @errorStr
	end
	
	def getHeader()  # station code, name, state, timezone, latitude, longitude, elev
		return @header
	end
	
	def getDirectBeam(hour)	
		return @directBeam[hour]
	end
	
	def getDiffuseHoriz(hour) 
		return @diffuseHoriz[hour]	
	end
	
	def getDirectHoriz(hour)
		return @directHoriz[hour]	
	end

	# ---------------------------------
	# Private functions
	# ---------------------------------	
	private
	
	# nothing private!
	
end   # class SSS_weatherFileTMY3   # 0.98
	
######################################################################
	
	
# -----------------------------------------------------------
# Adjust solar intensity for tilted surfaces  # 0.95
#
# INPUTS: solarZenith = solar angle from vertical (z-axis), in radians
#         solarAzimuth = solar angle from north, in radians
#         directBeam = direct beam intensity on horizontal surface
#         sunpt = Point3d sun position
#
# NOTES: If analysis surfaces are not facing upwards, they will be modified to do so!
#        If angle between normal and sun is > 90, result will be 0
#
# -----------------------------------------------------------
def self.adjustForTilt(directBeam, sunpt, solarZenith, solarAzimuth, surface)

	surfaceNorm = surface.normal
	
	if surfaceNorm.z < 0   # sketchup face is pointing in negative direction	
		status = surface.reverse!    # reverse the face and get new normal
		surfaceNorm = surface.normal
	end
	
	#x_vec = Geom::Vector3d.new(1,0,0)
	#y_vec = Geom::Vector3d.new(0,1,0)
	z_vec = Geom::Vector3d.new(0,0,1)
	sun_vec = Geom::Vector3d.new(sunpt.x, sunpt.y, sunpt.z)

	# ----------------------------------
	# Check for tilt of the surface (from horiz)
	# ----------------------------------
	
	surfaceTilt = (surfaceNorm.angle_between(z_vec)).radians      # assert: angle_between always returns between 0 and 180

	if surfaceTilt < 0.2     
	
		# no signficant tilt. just adjust for sun position		
		return directBeam * Math.cos(solarZenith)  # sun not as strong near horizons	

	else 
	
				
		if (surfaceNorm.angle_between(sun_vec)).radians > 90     # face is pointing away from sun
			return 0
		end
		
		# get surface aspect (angle from north)
		northAngle = Sketchup.active_model.shadow_info["NorthAngle"]
		if northAngle.nil?			
			northAngle = 0.0
			north_vec = Geom::Vector3d.new(0,1,0)
		else
			north_ptx = Math.cos((90-northAngle).degrees)
			north_pty = Math.cos(northAngle.degrees)
			north_vec = Geom::Vector3d.new(north_ptx, north_pty, 0)
		end	

		#==debug
		#mod = Sketchup.active_model # Open model
		#ent = mod.entities # All entities in model
		#pto = Geom::Point3d.new(0,0,0)
		#ptn = Geom::Point3d.new(north_vec.x, north_vec.y, 0)
		#ent.add_edges pto, ptn
		#==debug
	
		# project normal onto xy plane and get angle from north vector
		surfaceNormProj = surfaceNorm
		surfaceNormProj.z = 0	
		surfaceAspect = surfaceNormProj.angle_between(north_vec)   # in radians

		#==debug
		#puts "Surface Aspect = " + surfaceAspect.radians.to_s
		#ptsa = Geom::Point3d.new(surfaceNormProj.x, surfaceNormProj.y, 0)
		#ent.add_edges pto, ptsa
		#==debug
		
		# from Martin Rymes, NREL, solpos.c
		ca = Math.cos(solarAzimuth)
		cp = Math.cos(surfaceAspect)
		ct = Math.cos(surfaceTilt.degrees)
		cz = Math.cos(solarZenith)
		sa = Math.sin(solarAzimuth)
		sp = Math.sin(surfaceAspect)
		st = Math.sin(surfaceTilt.degrees)
		sz = Math.sin(solarZenith)
		
		tiltFactor = cz * ct + sz * st * ( ca * cp + sa * sp)
		
		if tiltFactor > 0.0
			return directBeam * tiltFactor
		else 
			return 0
		end
	end
	
end

# -----------------------------------------------------------
# Adjust solar intensity for sun zenith  # 0.98
# (For horizontal surfaces only)
#
# INPUTS: solarZenith = solar angle from vertical (z-axis), in radians
#         directBeam = direct beam solar intensity normal to surface
#
# -----------------------------------------------------------
def self.adjustForHorizon(directBeam, solarZenith)
	return directBeam * Math.cos(solarZenith)  # sun not as strong near horizons		
end


# ---------------------------------	
# open text file for writing
# returns file handle
# set mode to 0 for results file (writing)
#             1 for solar calcs file (writing)
#             2 for irradiance calcs file (writing) # 0.95
# ---------------------------------
def self.createFile(mode = 0)
	mod = Sketchup.active_model
	if mod.path != ""
		dirname = File.dirname(mod.path)
		basename = File.basename(mod.path, ".skp")	  # from testmodel.skp returns testmodel
	else
		if (/darwin/ =~ RUBY_PLATFORM) != nil   # mac os
			dirname = "~/Documents"
			basename = "untitled"
		else
			#dirname = File.expand_path('~')  # this doesn't work
			dirname = ""
			basename = "untitled"
		end
	end

	if 0 == mode
		model_name = basename + "-solardata.csv"			
		actualFile = UI.savepanel("Save results as comma-separated-value file", dirname, model_name)
	elsif 1 == mode
		t = Time.new
		model_name = basename + "-solarcalcs" + sprintf("%4d%02d%02d%02d%02d.csv", t.year, t.month, t.day, t.hour, t.min)
		actualFile = UI.savepanel("Save solar position calculation file", dirname, model_name)
	elsif 2 == mode
		t = Time.new
		model_name = basename + "-irradiancecalcs" + sprintf("%4d%02d%02d%02d%02d.csv", t.year, t.month, t.day, t.hour, t.min)
		actualFile = UI.savepanel("Save irradiance calculation file", dirname, model_name)		
	end
	
	return false if !actualFile
	fh = File.open(actualFile, 'w')
	return fh
end 

# ---------------------------------	
# open text file for reading
# returns file handle
# set mode to 0 for tmy3 weather file
# ---------------------------------
def self.openFile(mode = 0)
	mod = Sketchup.active_model
	if mod.path != ""
		dirname = File.dirname(mod.path)
		basename = File.basename(mod.path, ".skp")	  # from testmodel.skp returns testmodel
	else
		if (/darwin/ =~ RUBY_PLATFORM) != nil   # mac os
			dirname = "~/Documents"
			basename = "untitled"
		else
			#dirname = File.expand_path('~')  # this doesn't work
			dirname = ""
			basename = "untitled"
		end
	end

	if 0 == mode
		actualFile = UI.openpanel("Open comma-separated TMY3 weather file")			
	end
	
	return false if !actualFile
	fh = File.open(actualFile, 'r')
	return fh
end 


def self.closeFile(fh)
   fh.close
end

def self.writeFile(fh, str)
    fh.write(str)
end	
	
# returns float
def self.remapVal(val, current_lower, current_upper, new_lower, new_upper)
	if (val < current_lower || val > current_upper)
		return nil
	elsif (current_upper < current_lower || new_upper < new_lower)
		return nil
	else
    range1 = current_upper - current_lower
    range2 = new_upper - new_lower
    scale = (val.to_f - current_lower.to_f) / range1.to_f
    newval = (scale * range2.to_f) + new_lower
    return newval.to_f
	end
end

# 0.00 = black, 1.00 = white
def self.percentgray2rgb(val)
	r = remapVal(val, 0, 1, 0, 255)
	r = r.to_i
	rgb = r << 16 | r << 8 | r    # bitwise shift and bitwise OR
	return rgb
end

# saturation and brightness fixed at 100%
# see http://www.cs.rit.edu/~ncs/color/t_convert.html
def self.hue2rgb(hue)

	#void HSVtoRGB( float *r, float *g, float *b, float h, float s, float v )
	#int i;
	#float f, p, q, t;
	
	rgb = [0,0,0]
	
	v = 1.0
	s = 1.0
	h = hue.to_f / 60.0  # sector 0 to 5
	i = h.floor          
	f = h - i;			   # factorial part of h
	p = v * ( 1 - s );
	q = v * ( 1 - s * f );
	t = v * ( 1 - s * ( 1 - f ) );
	
	case i 
		when 0 
			rgb = [v,t,p]
		when 1 
			rgb = [q,v,p]
		when 2 
			rgb = [p,v,t]
		when 3 
			rgb = [p,q,v]
		when 4 
			rgb = [t,p,v]
		else   
			rgb = [v,p,q]
	end
	
	rgb[0] = remapVal(rgb[0], 0, 1, 0, 255)
	rgb[1] = remapVal(rgb[1], 0, 1, 0, 255)
	rgb[2] = remapVal(rgb[2], 0, 1, 0, 255)
	#puts rgb.inspect # DEBUG
	rgbnum = rgb[0].to_i << 16 | rgb[1].to_i << 8 | rgb[2].to_i    # bitwise shift and bitwise OR
	return rgbnum
	
end

# http://en.wikipedia.org/wiki/HSL_and_HSV
# h [0..360]  s [0..1]  v [0..1]
def self.hsv2rgb(hue, sat, val)
	
	rgb = [0,0,0]
	
	v = val.to_f
	s = sat.to_f
		
	c = v * s
	h = hue.to_f / 60.0  # sector 0 to 5
	tmp = (h % 2) - 1
	x = c * (1 - tmp.abs)
	
	if 0 <= h && h < 1
		rgb = [c,x,0]
	elsif 1 <= h && h < 2
		rgb = [x,c,0]
	elsif 2 <= h && h < 3
		rgb = [0,c,x]
	elsif 3 <= h && h < 4
		rgb = [0,x,c]
	elsif 4 <= h && h < 5
		rgb = [x,0,c]
	elsif 5 <= h && h < 6
		rgb = [c,0,x]
	else
		rgb = [0,0,0]	
	end
	
	m = v - c
	
	rgb[0] = rgb[0] + m
	rgb[1] = rgb[1] + m
	rgb[2] = rgb[2] + m
	
	rgb[0] = remapVal(rgb[0], 0, 1, 0, 255)
	rgb[1] = remapVal(rgb[1], 0, 1, 0, 255)
	rgb[2] = remapVal(rgb[2], 0, 1, 0, 255)
	#puts rgb.inspect # DEBUG
	rgbnum = rgb[0].to_i << 16 | rgb[1].to_i << 8 | rgb[2].to_i    # bitwise shift and bitwise OR
	return rgbnum
	
end


# ---------------------------------------------------------
# Given a face, find its centroid 
# Works for faces parallel to the ground plane. 
#
# Returns a Point3D
#
# see Fredo http://sketchucation.com/forums/viewtopic.php?f=323&t=13378&p=430718&hilit=centroid#p430718
# and Wikipedia http://en.wikipedia.org/wiki/Centroid#Centroid_of_polygon
#
#----------------------------------------------------------
def self.getCentroid(face)
   v = face.outer_loop.vertices	
   area = 0.0
	vdxmax1 = v.length - 1
	vdxmax2 = v.length - 2
	cx = 0
	cy = 0
	
	for i in 0 .. vdxmax2
		# calculate the area
		pt1 = v[i].position
		pt2 = v[i+1].position
		area += (pt1.x * pt2.y) - (pt2.x * pt1.y)
		cx += (pt1.x + pt2.x) * (pt1.x * pt2.y - pt2.x * pt1.y)
		cy += (pt1.y + pt2.y) * (pt1.x * pt2.y - pt2.x * pt1.y)
	end		
	# wrap case: x-sub-n = x-sub-0
	pt1 = v[vdxmax1].position
	pt2 = v[0].position
	area += (pt1.x * pt2.y) - (pt2.x * pt1.y)
	cx += (pt1.x + pt2.x) * (pt1.x * pt2.y - pt2.x * pt1.y)
	cy += (pt1.y + pt2.y) * (pt1.x * pt2.y - pt2.x * pt1.y)
	
	area = area / 2.0
	cx = cx / (6.0 * area)
	cy = cy / (6.0 * area)
	
   p = Geom::Point3d.new(cx, cy, 0)   
	#puts "face area: " + area.to_s + " centroid: " + p.to_s  DEBUG
	return p
	
end

# ----------------------------------------------------------------------
# Check if the face blocks sun Point from analysis Point
#
# RETURN
#   true if in shadow
# ----------------------------------------------------------------------
def self.inShadow(sunPt, analysisPt, face, modelEntities, debugParms)

	startpos = analysisPt
	endpos = sunPt
	vector = endpos-startpos
	plane = face.plane
	
	line = [ startpos, vector]
	
	intpt = Geom.intersect_line_plane(line, plane)   # locate intersection
	if intpt.nil?   # 0.91
		# DEBUG -- intersection should never be nil
		#modelEntities.add_edges startpos, endpos 
		#face.material = 0xff0000
		return false  # no intersection
	end
	
	#DEBUG
	#modelEntities.add_cpoint intpt
	#modelEntities.add_edges startpos, endpos 
	
	# test if the intersection pt is on the face
	result = face.classify_point(intpt)

	if (debugParms['MessageLevel'].to_i & 0x4) > 0
		puts result.to_s   # 0=unknown 1=on face 2=on vertex 4=on edge 16=outside 32=not on plane
	end
	
	# ideally use centroid analysis pt, so we only check for 16
	#if result == 16 || result == 2 || result == 4  # no intersection, so face did not block sun
	if result == 16  # no intersection, so face did not block sun
		return false
	else
		return true    # this face blocks the sun
	end
	
end
 
def self.sss_help_proc()
	help_file = File.dirname(__FILE__) + "/Simple_Solar_Study/help.html"	
	if (File.exist?(help_file))	
		UI.openURL "file://" + help_file
	else
		UI.messagebox("Error: " + help_file + " not found")
	end
end  # end help 
 
def self.sss_about_proc()
	aboutstr = "Simple Solar Study, v" + VERSION.to_s + "\nBy Brian Monwai\n\n"
   aboutstr = aboutstr + "Code for the calculating the Julian day, the Julian century, and an implementation of the Bird and Hulstrum clear sky radiation model are translated from solrad_ver16.xls with permission from Greg Pelletier, Washington State Dept of Ecology.\n"
	aboutstr = aboutstr + "\n"
	aboutstr = aboutstr + "Other calculations and algorithms are based on Excel formulas in BIRD_08_16_2012.xls by Daryl Myers, National Renewable Energy Laboratory, "
	aboutstr = aboutstr + "Javascript code in the online NOAA Solar Calculator at http://www.esrl.noaa.gov/gmd/grad/solcalc/, and C code from solpos.c (1998) by Martin Rymes, National Renewable Energy Laboratory. "
	aboutstr = aboutstr + "Versions of SimpleSolarStudy before v0.93 were based on the SketchUp plugin Sunposition, v1.2.1 by Gabriel Miller, 2009.\n"
	aboutstr = aboutstr + "\n"
	aboutstr = aboutstr + "DISCLAIMER. None of the authors or their respective organizations are responsible "
	aboutstr = aboutstr + "for any damages, including incidental or consequential damages, arising from use "
	aboutstr = aboutstr + "or misuse of this software, or from results achieved or conclusions drawn by others.\n"
	aboutstr = aboutstr + "\n"
	aboutstr = aboutstr + "NO WARRANTIES. This software and any related documentation is provided "
	aboutstr = aboutstr + "as is without warranty of any kind, either express or implied, including, without "
	aboutstr = aboutstr + "limitation, the implied warranties or merchantability, fitness for a particular purpose, "
	aboutstr = aboutstr + "or noninfringement. The entire risk arising out of use or performance of the software "
	aboutstr = aboutstr + "remains with the user.\n"
	aboutstr = aboutstr + "\n"
	aboutstr = aboutstr + "NO LIABILITY FOR DAMAGES. In no event shall the Washington State Department of Ecology, "
	aboutstr = aboutstr + "NREL, NOAA, the developers of the software, or their suppliers be liable for any damages "
	aboutstr = aboutstr + "whatsoever (including, without limitation, damages for loss of business profits, "
	aboutstr = aboutstr + "business interruption, loss of business information, or any other pecuniary loss) "
	aboutstr = aboutstr + "arising out of the use of or inability to use this product, even if the "
	aboutstr = aboutstr + "Washington State Department of Ecology, NREL, NOAA, the software developers or their suppliers "
	aboutstr = aboutstr + "have been advised of the possibility of such damages.\n"
	UI.messagebox(aboutstr)
end  # end about

def self.sss_main_proc(basic_mode = true)

mod = Sketchup.active_model # Open model
modelEnts = mod.active_entities
sel = mod.selection # Current selection


# DEBUG  0.93
#mysun = SSS_sunData.new
#mysun.calcSun(40, -105, 2012, 6, 21, 1, -6)
#UI.messagebox(sprintf("gmls=%.6f gmas=%.6f ecc=%.6f rad=%.6f declin=%.6f eqntime=%.6f solartime=%.6f hrangle=%.6f zenith=%.6f elevangle=%.6f azimuth=%.6f", mysun.getGeometricMeanLongitude()	, mysun.getGeometricMeanAnomaly(), mysun.getEccentricity(), mysun.getRadius(), mysun.getDeclination(), mysun.getEqnTime(), mysun.getSolarTime(), mysun.getHrAngle(), mysun.getZenith(), mysun.getElevAngle(), mysun.getAzimuth() ))		
#return

# ------------------------
# Make sure selection is a Group
# ------------------------
if (sel.length != 1)
    UI.messagebox("Please select a single Group of faces to analyze")    
    return 
end 

# Look at all of the entities in the selection.
grp = sel[0]
trans = grp.transformation  # get translation of group relative to world coords  

if !grp.kind_of?(Sketchup::Group)
    UI.messagebox("The selection must be a Group")
    return 
end

# --------------------------
# Get user input
# --------------------------
myLayers = mod.layers
myLayersIdx = myLayers.count - 1 
myLayerNamesArr = []
myLayerNamesHash = Hash.new
myLayerNames = ""

glassLayerNamesHash = Hash.new  # 0.94
glassLayerNamesArr = Array.new

for i in 0 .. myLayersIdx
    myLayerNamesArr[i] = myLayers[i].name	 
	 myLayerNamesHash[myLayers[i].name] = true
    if (i == 0)
       myLayerNames = myLayers[i].name
    else
       myLayerNames = myLayerNames + "|" + myLayers[i].name
    end
end

t = Time.now	
y = t.gmtime.year
monthNames = { "Jan" => 0, "Feb" => 1, "Mar" => 2, "Apr" => 3, "May" => 4, "Jun" => 5, "Jul" => 6, "Aug" => 7, "Sep" => 8, "Oct" => 9, "Nov" => 10, "Dec" => 11}

# extra parms
debugParms = Hash.new
debugParms['SaveSolarCalcs'] = "No"
debugParms['SaveIrradianceCalcs'] = "No"
# decide which messages are output to ruby console
debugParms['MessageLevel'] = "1"   
debugParms['LabelLayerName'] = "Solar study labels"  # 0.94
debugParms['SunpathLayerName'] = "Solar study sunpath"  # 0.95

if (basic_mode)

	prompts1 = ["Layer with shade structures: ", 
				  "Time zone", 
				  "Start month", "Start day", 
				  "End month", "End day", 
				  "Year", 
				  "Calculation mode", 
				  "Label results?", 
				  "Gradient"]
				  
	defaults1 = [myLayerNamesArr[0], 
	            "-6",
					"Jan", "1", 
					"Dec", "31", 
					y.to_s, 
					"Medium", 
					"No", 
					"Color ramp 2"]
					
	lists1 = [myLayerNames, 
				"-12|-11|-10|-9|-8|-7|-6|-5|-4|-3|-2|-1|0|1|2|3|4|5|6|7|8|9|10|11|12", 
				"Jan|Feb|Mar|Apr|May|Jun|Jul|Aug|Sep|Oct|Nov|Dec", "", 
				"Jan|Feb|Mar|Apr|May|Jun|Jul|Aug|Sep|Oct|Nov|Dec", "", 
				"", 
				"Coarse|Medium|Fine|TMY3", 
				"Yes|No", 
				"None|Color ramp 2|Color ramp 1|Grayscale"]				

	userinput = UI.inputbox prompts1, defaults1, lists1, "Simple Solar Study - Basic Mode"
	if !userinput || userinput.nil?
		return 
	end

	ocLayerName = userinput[0]
	userTimeZone = userinput[1].to_i
	userStartMonth = monthNames[userinput[2]]
	userStartDate = userinput[3].to_i
	userEndMonth = monthNames[userinput[4]]
	userEndDate = userinput[5].to_i
	userYear = userinput[6].to_i
	calculationMethod = userinput[7]
	drawPath = "No"   # 0.94
	labelResults = userinput[8]
	outputMode = "Daylight factor"
	gradientType = userinput[9]
	saveResults = "No"   # 0.94
	debugString = ""
	
else  # advanced mode
	prompts2 = ["Layer(s) with shade structures: ",   # comma-separated  # 0.94
				  "Layer(s) of glazing: ",               # comma-separated  # 0.94
				  "Glazing transmittance [0-100]: ",     # comma-separated, to correspond with above layers # 0.94
				  "Time zone", 
				  "Start month", "Start day", 
				  "End month", "End day", 
				  "Year", 
				  "Calculation mode",
				  "Trace sun path?", 
				  "Label results?", 
				  "Result units",
				  "Gradient", 				  
				  "Save results to file?",
				  "Advanced parameters"]
	defaults2 = ["",                    # shade
	            "",                     # glass
					"70",                   # transmittance
	            "-6",                   # tz
					"Jan", "1",             # start
					"Dec", "31",            # end
					y.to_s,                 # year
					"Medium",               # calc mode
					"No",                   # trace path
					"No",                   # label
					"Daylight factor",      # output mode
					"Color ramp 2",         # gradient
					"No",                   # save to file
					""]                     # debug parms
	lists2 = ["",   
	          "",   
				 "",
				"-12|-11|-10|-9|-8|-7|-6|-5|-4|-3|-2|-1|0|1|2|3|4|5|6|7|8|9|10|11|12", 
				"Jan|Feb|Mar|Apr|May|Jun|Jul|Aug|Sep|Oct|Nov|Dec", "", 
				"Jan|Feb|Mar|Apr|May|Jun|Jul|Aug|Sep|Oct|Nov|Dec", "", 
				"", 
				"Coarse|Medium|Fine|TMY3", 
				"Yes|No",				
				"Yes|No", 
				"Daylight factor|W/m^2 (calc mode=fine)",
				"None|Color ramp 2|Color ramp 1|Grayscale", 
				"Yes|No",
				""]

	userinput = UI.inputbox prompts2, defaults2, lists2, "Simple Solar Study - Advanced Mode"
	if !userinput || userinput.nil?
		return 
	end

	ocLayerName = userinput[0]
	glassLayerNamesArr = userinput[1].split(',')  # 0.94
	glassTransmittanceArr = userinput[2].split(',') # 0.94
	userTimeZone = userinput[3].to_i
	userStartMonth = monthNames[userinput[4]]
	userStartDate = userinput[5].to_i
	userEndMonth = monthNames[userinput[6]]
	userEndDate = userinput[7].to_i
	userYear = userinput[8].to_i
	calculationMethod = userinput[9]
	drawPath = userinput[10]
	labelResults = userinput[11]
	outputMode = userinput[12]  # 0.95
	gradientType = userinput[13]	
	saveResults = userinput[14]
	debugString = userinput[15]   # 0.94

	debugArr = Array.new
	debugArr = debugString.gsub(/\s+/, '').split(',')  # remove all spaces and then split on commas
	debugArr.each { |a|
		pair = Array.new
		pair = a.split('=')  # split on equal sign
		debugParms[pair[0]] = pair[1]
	}
	
	
end  # UI: basic or advanced mode

# -----------------------------------------
# gather all potentially occluding faces
# -----------------------------------------
ocFaces = []
glassFaces = []
glassFacesVt = []

# Allow multiple shade layers and glass layers     # 0.94
ocLayerNamesHash = Hash.new
ocLayerNamesArr = ocLayerName.split(',')
len = ocLayerNamesArr.length - 1
for i in 0 .. len
	a = ocLayerNamesArr[i]
	if !myLayerNamesHash.key?(a)  
		UI.messagebox("Shade layer " + a + " does not exist")
		return
	end
	ocLayerNamesHash[a] = true 
end

# check for valid layers
len = glassLayerNamesArr.length - 1
for i in 0 .. len 
	a = glassLayerNamesArr[i] 
	if !myLayerNamesHash.key?(a)  
		UI.messagebox("Glazing layer " + a + " does not exist")
		return
	end
	# hash key = Layer Name, value = transmittance percent
	if glassTransmittanceArr.length == glassLayerNamesArr.length
		glassLayerNamesHash[a] = glassTransmittanceArr[i].to_f / 100.0  # map transmittance values to each layer
	else
		glassLayerNamesHash[a] = glassTransmittanceArr[0].to_f / 100.0  # ELSE apply first transmittance value to all layers
	end
end

# loop through faces and add to arrays
# and bounding box  # 1.10
#ocBBox = SSS_BoundingBox.new
#glassBBox = SSS_BoundingBox.new
modelEnts.each { |ent|
    if ent.kind_of?(Sketchup::Face)
		if ocLayerNamesHash.key?(ent.layer.name)   # if face's layer is in the hash, push face to array   # 0.94
         ocFaces.push(ent)
			#ocBBox.add(ent)
		elsif glassLayerNamesHash.key?(ent.layer.name)   
		   glassFaces.push(ent)
			glassFacesVt.push(glassLayerNamesHash[ent.layer.name])   # add visual transmittance
			#glassBBox.add(ent)
		end
    end
}
#ocBBox.calcBBPlanes()
#ocBBox.drawBBEdges(modelEnts)
#ocBBox.drawBBFaces(modelEnts)
#glassBBox.calcBBPlanes()
#glassBBox.drawBBEdges(modelEnts)




# -------------------------------
# Status bar as progress bar
# -------------------------------
startTime = Time.new   # 0.96
puts "Starting Simple Solar Study"
Sketchup.set_status_text "Calculating solar positions", SB_PROMPT    # 0.91
mod.start_operation("SimpleSolarStudy")  # 0.95 for undo

# -------------------------
# Create annual sun position
# -------------------------

# get latitude and longitude of point in the group
if !(mod.georeferenced? || (debugParms.key?("Lat") && debugParms.key?("Long")))   # 0.94 georeference OR provide lat&long
	UI.messagebox('This model needs to be georeferenced.')
	return
end

#        longitude (neg = west)
#           |         latitude (pos = north)
#           |          |
#           |          |
myLoc = [-122.332071, 47.606202, 0]
origin = Geom::Point3d.new(0,0,0)   
ents = grp.entities
ents.each { |ent| 
    if ent.kind_of?(Sketchup::Face)        
		  origin = ent.vertices[0].position    # grab the first vertex of the first face instead of [0,0,0]  # 0.92
		  origin.transform! trans              # convert point from group coords to world coords
		  break
    end
}
myLoc = mod.point_to_latlong(origin)  # expects sketchup 6.0+ to return point3D   
latitude = myLoc[1].to_f
longitude = myLoc[0].to_f  

# allow advanced user override 
if debugParms.key?("Lat")  # 0.94
	latitude = debugParms['Lat'].to_f
end
if debugParms.key?("Long")
	longitude = debugParms['Long'].to_f
end


# DEBUG
if (debugParms['MessageLevel'].to_i & 0x1) > 0
	puts "Origin=" + origin.inspect + " myLoc=" + myLoc.inspect + " latitude=" + latitude.to_s + " longitude=" + longitude.to_s
end


daysInMonthNotLeap = [31,28,31,30,31,30,31,31,30,31,30,31]
daysInMonthLeap   =  [31,29,31,30,31,30,31,31,30,31,30,31]
if (userYear % 4 == 0 && userYear % 100 != 0) || userYear % 400 == 0    # 0.91
  daysInMonth = daysInMonthLeap
else
  daysInMonth = daysInMonthNotLeap
end

# simple date checks
if (userStartMonth > userEndMonth) || (userStartMonth == userEndMonth && userStartDate > userEndDate)
	UI.messagebox('End date cannot be earlier than start date')
	return 
elsif (userStartDate > daysInMonth[userStartMonth]) || (userEndDate > daysInMonth[userEndMonth])
	UI.messagebox('There are not that many days in the month')
	return 	
end 

# ------------------------------
# Initialize data structures
# ------------------------------
mysun = SSS_sunData.new   # 0.93
mytmy3 = SSS_weatherFileTMY3.new # 0.98
if calculationMethod == "TMY3" 
	if mysun.isLeapYear(userYear) 
		UI.messagebox("TMY3 cannot be used on a leap year")
		return
	end
	wfh = openFile(0)
	if wfh
		fileok = mytmy3.readFile(wfh) if wfh   # suck file into array
		if (!fileok)
			UI.messagebox(mytmy3.getError())
			return
		else
			puts mytmy3.getHeader()
		end
	else
		UI.messagebox("Could not open weather file")
		return
	end
end

# 0.94
if debugParms.key?("Pressure")
	mysun.SetPressure(debugParms["Pressure"].to_f)   
end
if debugParms.key?("Ozone")
	mysun.SetOzone(debugParms["Ozone"].to_f)   
end
if debugParms.key?("Water")
	mysun.SetWater(debugParms["Water"].to_f)   
end
if debugParms.key?("Aerosol-500")
	mysun.SetAerosol_500(debugParms["Aerosol-500"].to_f)   
end
if debugParms.key?("Aerosol-380")
	mysun.SetAerosol_380(debugParms["Aerosol-380"].to_f)   
end
if debugParms.key?("FwdScatter")
	mysun.SetFwdScatter(debugParms["FwdScatter"].to_f)   
end
if debugParms.key?("Reflectance")
	mysun.SetReflectance(debugParms["Reflectance"].to_f)   
end
if debugParms.key?("SunDist")
	mysun.setSunDistance(debugParms["SunDist"].to_i) 
end

sunspots = Array.new      # sun position
sunstrength = Array.new   # solar intensity 0.93
sunzenith = Array.new     # keep track for tilt calculations 0.95
sunazimuth = Array.new    # keep track for tilt calculations 0.95
sundiffuse_hz = Array.new 
sundirect_hz = Array.new
clearsundiffuse_hz = Array.new  # for max daylight in TMY3 mode  # 0.98
clearsundirect_hz = Array.new
#displayRays = Array.new  # for drawing ray from sun to origin. commented out 0.93
timesteps = 0
prevspt = Geom::Point3d.new()

if debugParms['SaveSolarCalcs'] == "Yes"
	cfh = createFile(1)
	if calculationMethod == "Fine" || calculationMethod == "TMY3"  # Bird model - solar irradiance
		writeFile(cfh, "Date,Time,Declination,EqnTime,HourAngle,Elev(deg),Azimuth(deg),DirectBeamNormal(W/m^2),DirectHoriz(W/m^2),DiffuseHoriz(W/m^2),GlobalHoriz(W/m^2)\n") if cfh
	else
		writeFile(cfh, "Date,Time,Declination,EqnTime,HourAngle,Elev(deg),Azimuth(deg),Sunstrength\n") if cfh
	end
end

if drawPath == "Yes"   # 0.95
	sp_layer = myLayers.add debugParms['SunpathLayerName']
	mod.active_layer = sp_layer
end

# -----------------------------
# Loop through months
# -----------------------------
for mdx in userStartMonth .. userEndMonth do
    m = mdx + 1
    maxDays = daysInMonth[mdx]

    # -----------------------------	 
	 # Loop through days
	 # -----------------------------
    for day in 1 .. maxDays do
	 
		if (mdx == userStartMonth && day < userStartDate) || (mdx == userEndMonth && day > userEndDate)
			next
		end
	
		sunrise = true
		for step in 0 .. 23
			timesteps += 1
			
			# ---------------------------------
			# calculate sun position
			# ---------------------------------
			spt = mysun.calcSun(latitude, longitude, userYear, m, day, step, userTimeZone)   # 0.93
			elev = mysun.getElevAngle()
			# Date,Time,Declination,EqnTime,HourAngle,Elev,Azimuth,Sunstrength
			writeFile(cfh, sprintf("%04d-%02d-%02d,%02d:00,%.4f,%.4f,%.4f,%.4f,%.4f", userYear, m, day, step, mysun.getDeclination(), mysun.getEqnTime(), mysun.getHrAngle(), mysun.getElevAngle(), mysun.getAzimuth())) if cfh

			# ---------------------------------
			# calculate solar intensity and establish 
			# sunshineStr for printing 
			# ---------------------------------
			if calculationMethod == "Medium"  
				sunshine = Math.sin(elev.degrees)  # sun not as strong near horizons
				sunshineStr = sprintf(",%.4f\n", sunshine)
			elsif calculationMethod == "Fine"   # Bird model/solar irradiance
				mysun.calcSolarRadiation(userYear, m, day)  				
				sunshine = mysun.getDirectBeam()   # sun intensity normal to beam
				diffuse_hz = mysun.getDiffuseHoriz() 
				direct_hz = mysun.getDirectHoriz()
				sunshineStr = sprintf(",%.4f,%.4f,%.4f,%.4f\n", sunshine, direct_hz, diffuse_hz, mysun.getGlobalHoriz())
			elsif calculationMethod == "TMY3"   # based on observed data, read from weather file   # 0.98
				mysun.calcSolarRadiation(userYear, m, day)  # need this for clear day max light			
				cleardiffuse_hz = mysun.getDiffuseHoriz()   # need this for clear day max light			
				cleardirect_hz = mysun.getDirectHoriz()     # need this for clear day max light			
				sunshine = mytmy3.getDirectBeam(timesteps - 1)
				diffuse_hz = mytmy3.getDiffuseHoriz(timesteps - 1) 
				direct_hz = mytmy3.getDirectHoriz(timesteps - 1)
				sunshineStr = sprintf(",%.4f,%.4f,%.4f,%.4f\n", sunshine, direct_hz, diffuse_hz, direct_hz + diffuse_hz)								
			elsif calculationMethod == "Coarse"
				sunshineStr = ",0\n"
			end
			
			# ------------------------------------
			# Add to array only if we're above horizon
			# ------------------------------------
			if (!spt.nil? && elev > 0.0 && elev < 180.0)
			
				sunspots.push(spt) # add location			
			
				if drawPath == "Yes" && day == 1 && m < 7  # trace path for first day of every month, Jan - Jun     # 0.93
					if sunrise
						prevspt = spt
						sunrise = false
						n = Geom::Vector3d.new(0,0,1)
						modelEnts.add_circle spt, n, 20
					else
						modelEnts.add_edges prevspt, spt    # connect the suns together
						prevspt = spt
					end
				end			
				# draw edges?
				#if day == 1 && step % 4 == 0  # we're only going to display sun rays on first of the month, every four hours
				#	displayRays.push(1)
				#else
			   #	displayRays.push(0)
				#end
				
				# solar strength?   # 0.93
				if calculationMethod == "Medium"
					sunstrength.push(sunshine)				
				elsif	calculationMethod == "Fine" 
					sunzenith.push(mysun.getZenith())
					sunazimuth.push(mysun.getAzimuth())
					sundiffuse_hz.push(diffuse_hz)
					sundirect_hz.push(direct_hz)
					sunstrength.push(sunshine)		
				elsif calculationMethod == "TMY3"  # 0.98
					sunzenith.push(mysun.getZenith())
					sunazimuth.push(mysun.getAzimuth())
					clearsundiffuse_hz.push(cleardiffuse_hz)  # need this for clear day max light			
					clearsundirect_hz.push(cleardirect_hz)    # need this for clear day max light			
					sundiffuse_hz.push(diffuse_hz)
					sundirect_hz.push(direct_hz)
					sunstrength.push(sunshine)							
				elsif calculationMethod == "Coarse"  
					sunshineStr = ",1\n"
					sunstrength.push(1.0)							
				end
						
				
			end  # end if elev > 0
			
			writeFile(cfh, sunshineStr) if cfh
			
		 end # hours
   end #days
end #months

closeFile(cfh) if cfh
	

# -------------------------
# Now start analysis
# -------------------------
if debugParms['SaveIrradianceCalcs'] == "Yes" && (calculationMethod == "Fine" || calculationMethod == "TMY3")
	ifh = createFile(2)
	writeFile(ifh, "face-index,sunpos-index,zenith,azimuth,directBeamNormal,diffuseHoriz,directHoriz,vt,adjustedDirectHz,result,cumResult\n") if ifh
end

# add analysis faces to array
facecount = 0
faces = []
centroids = []

ents = grp.entities
ents.each { |ent|
    if ent.kind_of?(Sketchup::Face)        
        facecount = facecount + 1
        faces.push(ent)
    end
}

facesMaxIdx = faces.length - 1
sunspotsMaxIdx = sunspots.length - 1
ocFacesMaxIdx = ocFaces.length - 1
glassFacesMaxIdx = glassFaces.length - 1
results = Array.new(faces.length) { |i| 0 }  # store the results of sun strikes in new array. index will match analysis faces

if (debugParms['MessageLevel'].to_i & 0x1 ) > 0
	puts "Analysis pts=" + faces.length.to_s + " timesteps=" + timesteps.to_s + " sunpositions=" + sunspots.length.to_s + " occluding faces=" + ocFaces.length.to_s + " glazing faces=" + glassFaces.length.to_s   # 0.94
end


# -------------------------
# loop through analysis points, calculating incident light
# -------------------------
for ii in 0 .. facesMaxIdx do
    
	 # Progress bar	 
	 ii_plus = ii + 1
	 Sketchup.set_status_text "Measuring irradiance on face " + ii_plus.to_s + " of " + facecount.to_s, SB_PROMPT    # 0.91
	 
    # Find the centroid of the face
	 p = getCentroid(faces[ii])
	 centroids[ii] = p   # save the centroids in group coords 
    p2 = p.transform trans  # convert point from group coords to world coords 
     
	 # loop through sun positions
    for jj in 0 .. sunspotsMaxIdx do		 
		lightFactor = 1.00   # 0.94

		#==debug
		#if drawEdges == "Yes" && displayRays[jj] == 1    # option to draw edge between analysis pt and sun pt
		#    modelEnts.add_edges p2, sunspots[jj]		
		#end
		#==debug
				
		#if ocFacesMaxIdx > 16    # check bounding box first  # 1.10
			#ocBBox.inShadow(sunspots[jj], p2)
		#end
		
		# loop through occluding faces (only if face is not "hidden" in sketchup # 0.92)
		for kk in 0 .. ocFacesMaxIdx do			 			 
			 if ocFaces[kk].hidden? == false && inShadow(sunspots[jj], p2, ocFaces[kk], modelEnts, debugParms)   # don't count if sun is blocked
				lightFactor = 0.00
				break  # sun is blocked; quit loop early
			end
		end 
		
		# loop through glass faces
		if lightFactor > 0
			for mm in 0 .. glassFacesMaxIdx do
				if glassFaces[mm].hidden? == false && inShadow(sunspots[jj], p2, glassFaces[mm], modelEnts, debugParms)   # don't count if sun is blocked
					lightFactor = lightFactor * glassFacesVt[mm]    # glassFacesVt is between 0.0 and 1.0 -- allows for successive panes of glass
				end
			end
		end
		
		if calculationMethod == "Fine" || calculationMethod == "TMY3"  # 0.95
			# account for tilt, glass, and diffuse (assume diffuse is same for any angle)
			adjustedIntensity = (adjustForTilt(sundirect_hz[jj], sunspots[jj], sunzenith[jj].degrees, sunazimuth[jj].degrees, faces[ii]))  # 0.98 use direct hz instead of normal
			result = sundiffuse_hz[jj] + (adjustedIntensity * lightFactor)
			results[ii] += result
			writeFile(ifh, sprintf("%d,%d,%.1f,%.1f,%.1f,%.1f,%.1f,%.2f,%.1f,%.1f,%.1f\n", ii, jj, sunzenith[jj], sunazimuth[jj], sunstrength[jj], sundiffuse_hz[jj], sundirect_hz[jj], lightFactor, adjustedIntensity, result, results[ii])) if ifh
		else 
			results[ii] += sunstrength[jj] * lightFactor        # increment sun strikes if no shadow  # 0.93   # factor in glass # 0.94
		end
				
	 end
	 
end
closeFile(ifh) if ifh


# -------------------------
# now color faces according to proportion of max daylight
# -------------------------
Sketchup.set_status_text "Gathering results for " + facecount.to_s + " faces", SB_PROMPT 

# get max daylight value
maxDaylight = 0.0
if calculationMethod == "Medium"
	sunstrength.each { |v| maxDaylight += v.to_f } 
elsif	calculationMethod == "Fine" 
	for jj in 0 .. sunspotsMaxIdx do		 
		maxDaylight += adjustForHorizon(sundirect_hz[jj], sunzenith[jj].degrees) + sundiffuse_hz[jj]   
	end
elsif	calculationMethod == "TMY3"   # 0.98  get daylight factor relative to clear sky
	for jj in 0 .. sunspotsMaxIdx do		 
		maxDaylight += adjustForHorizon(clearsundirect_hz[jj], sunzenith[jj].degrees) + clearsundiffuse_hz[jj]   
	end
elsif calculationMethod == "Coarse"
	maxDaylight = sunstrength.length.to_f
end

maxDaylightFactor = 0.0
minDaylightFactor = 1.0
for ii in 0 .. facesMaxIdx do
		 
	daylightFactor = results[ii].to_f / maxDaylight
	
	# find max and min
	if daylightFactor > maxDaylightFactor    
		maxDaylightFactor = daylightFactor 
	end
	if daylightFactor < minDaylightFactor    
		minDaylightFactor = daylightFactor 
	end
	
	# color the front and back faces of the analyzed surface
	if gradientType == "Color ramp 1"
		huescale = remapVal(1 - daylightFactor, 0, 1, 240, 360)  # gradient from blue (not much sun) to red (lots)
		if (debugParms['MessageLevel'].to_i & 0x2) > 0
			puts "daylight factor: " + sprintf("%.2f", daylightFactor) + "; Hue: " + sprintf("%.2f", huescale) # DEBUG
		end
		rgb = hue2rgb(huescale)
	elsif gradientType == "Color ramp 2"
		huescale = remapVal(daylightFactor, 0, 1, 120, 240)  # gradient from red (sunny) to green (not much)
		if (debugParms['MessageLevel'].to_i & 0x2) > 0
			puts "daylight factor: " + sprintf("%.2f", daylightFactor) + "; Hue: " + sprintf("%.2f", huescale) # DEBUG
		end
		rgb = hsv2rgb(huescale, 1, 1)
	else 
		rgb = percentgray2rgb(daylightFactor)		
	end
	#puts "0x" + rgb.to_s(16)  # DEBUG
	
	if gradientType != "None"    # 0.96
		faces[ii].material = rgb
		faces[ii].back_material = rgb	
	end
	
	# add text
	if labelResults == "Yes"
		df_layer = myLayers.add debugParms['LabelLayerName']
		mod.active_layer = df_layer
		p = centroids[ii]
		p2 = p.transform trans  # convert point from group coords to world coords 
		if outputMode == "Daylight factor"
			modelEnts.add_text sprintf("%.2f", daylightFactor), p2
		else
			modelEnts.add_text sprintf("%d", results[ii]), p2
		end
	end
	# add normals
	if debugParms['ShowNormals'] == "Yes"
		df_layer = myLayers.add debugParms['LabelLayerName']
		mod.active_layer = df_layer
		p = centroids[ii]
		p2 = p.transform trans  # convert point from group coords to world coords 
		surfaceNorm = faces[ii].normal
		surfaceNorm.length = surfaceNorm.length * daylightFactor * 10
		surfaceNormPt = Geom::Point3d.new(p.x + surfaceNorm.x, p.y + surfaceNorm.y, p.z + surfaceNorm.z)
		surfaceNormPt.transform! trans  # convert point from group coords to world coords 
		modelEnts.add_edges p2, surfaceNormPt
	end
	
end

# ----------------------
# Results
# ----------------------
outstr_bas = "Faces analyzed: " + facecount.to_s + "\nTotal hours: " + timesteps.to_s + "\nHours of sunlight: " + sunspots.length.to_s + "\nShading faces: " + ocFaces.length.to_s +                                                                                                                               "\nMax daylight factor: " + sprintf("%.2f", maxDaylightFactor) + "\nMin daylight factor: " + sprintf("%.2f", minDaylightFactor)
outstr_adv = "Faces analyzed: " + facecount.to_s + "\nTotal hours: " + timesteps.to_s + "\nHours of sunlight: " + sunspots.length.to_s + "\nShading faces: " + ocFaces.length.to_s + "\nFaces with transmittance: " + glassFaces.length.to_s + "\nMax irradiance (single face): " + sprintf("%.2f", maxDaylight) + "\nMax daylight factor: " + sprintf("%.2f", maxDaylightFactor) + "\nMin daylight factor: " + sprintf("%.2f", minDaylightFactor)
outstr_fh  = "Faces analyzed: " + facecount.to_s + " #Total hours: " + timesteps.to_s + " #Hours of sunlight: " + sunspots.length.to_s + " #Shading faces: " + ocFaces.length.to_s + " #Faces with transmittance: " + glassFaces.length.to_s + "\nMax irradiance (single face): " + sprintf("%.2f", maxDaylight) + " #Max daylight factor: " + sprintf("%.2f", maxDaylightFactor) + " #Min daylight factor: " + sprintf("%.2f", minDaylightFactor)
modelUnitAbbrev = ["inches", "feet", "mm", "cm", "m"]   # 0.93
modelUnits = modelUnitAbbrev[Sketchup.active_model.options["UnitsOptions"]["LengthUnit"]]

if saveResults == "Yes"
	fh = createFile(0)
	if (fh)
		writeFile(fh, "#Simple Solar Study " + sprintf("%i%02d%02d-%i%02d%02d", userYear, userStartMonth + 1, userStartDate, userYear, userEndMonth + 1, userEndDate) + "\n")
		writeFile(fh, "#Latitude=" + latitude.to_s + " #Longitude=" + longitude.to_s + "\n")
		writeFile(fh, "#" + outstr_fh + "\n")
		if outputMode == "Daylight factor"
			writeFile(fh, "index,faceID,centroid-x(" + modelUnits + "),centroid-y(" + modelUnits + "),daylightFactor\n")   # 0.93
		else
			writeFile(fh, "index,faceID,centroid-x(" + modelUnits + "),centroid-y(" + modelUnits + "),Irradiance\n")   # 0.95			
		end
		for ii in 0 .. facesMaxIdx do
			daylightFactor = results[ii].to_f / maxDaylight
			p = centroids[ii]
			p2 = p.transform trans  # convert point from group coords to world coords
			if outputMode == "Daylight factor"
				writeFile(fh, sprintf("%i,%s,%d,%d,%0.4f\n", ii, faces[ii].entityID.to_s, p2.x, p2.y, daylightFactor))  # 0.93
			else
				writeFile(fh, sprintf("%i,%s,%d,%d,%d\n", ii, faces[ii].entityID.to_s, p2.x, p2.y, results[ii]))  # 0.95				
			end
		end
		closeFile(fh)
	end
end

mod.commit_operation  # 0.95  for undo

timeDiff = Time.new - startTime  # 0.96
timeStr = "Simple Solar Study complete in " + timeDiff.to_s + " seconds"
puts timeStr
if basic_mode
	UI.messagebox(outstr_bas + "\n" + timeStr)
else
	UI.messagebox(outstr_adv + "\n" + timeStr)
end

end # main_proc
	
end #module