#!/usr/bin/python2

# TODO: Have it sap air from the greenhouse at temperature proportional to height

# Gmodel: A greenhouse simulator
#
# NOTE: This simulator was initially written with my greenhouse in mind.  
# Thus, it may not encompass the designs of *everyone*'s greenhouse.  Also
# note that there are many factors that are difficult to model (such as IR)
# and hard to get good data for (such as the particular sunlight absorption
# details of a given greenhouse).  Be warned!  :)
#
# Apart from these, bombs away with criticisms, snide remarks, and most
# constructively, feature and buxfix patches.  :)  Mail me at 
# meme@daughtersoftiresias.org.  In all patches, be sure to strip out any
# configuration changes you've made that you don't have a good reason to
# make part of the defaults.
#
# TODOS:
#
# * Our current water model assumes that water is in cylindrical barrels,
#   touching the ground, only on the south side, painted black, covering the
#   entire side, absorbs all heat instantly, has a surface temperature
#   equivalent to the water on all sides.  Can we reduce these assumptions?
# * Our IR model could use a good analysis.  :P  The big problem with it
#   is that, while for a fixed number of bodies with known geometry, and a
#   lot of raycasting, you can figure this problem out.  For the situation of
#   a greenhouse (with shelves, pots, etc) and an outside environment, it's 
#   essentially unmodelable.  Thus, I have to make some pretty major
#   assumptions, and I hope that it doesn't throw the results off too much.
# * We're assuming that the indoor ground never gets sunlight due to the
#   barrels on the south side.  Remove this assumption. 
# * TempAtHeight needs to be flushed out.  In fact, our infiltration model
#   could use some work, such as removing the flow fudge factor.
# * The various constants that are treated as unchanging with temperature,
#   but actually do change with temperature in the real-world.

# // I N C L U D E S ////////////////////////////////////////////////////////////

import sys, random
from math import *

import climate_data

# // D E F I N E S //////////////////////////////////////////////////////////////

# Greenhouse defines
VerticalWallHeight =	1.2	# Meters
PeakHeight = 		2.0	# Meters
Length = 		2.0	# Meters
Width = 		1.3	# Meters
WestClearPercent = 	0.3	# %/100
EastClearPercent =	0.3	# %/100
WaterDiameter = 	0.28	# Meters
WaterHeight = 		0.29	# Meters
WaterContainers =	7	# Count
OpaqueInsul_R = 	2.0	# K*m^2/W -- 0.1761 * American R-values
ClearInsul_R =	 	1.43	# K*m^2/W -- 0.1761 * American R-values
WaterInsul_R = 		0.6	# K*m^2/W -- 0.1761 * American R-values
GroundCoverInsul_R = 	0.2	# K*m^2/W -- 0.1761 * American R-values
OpaqueInsulIRReflect = 	0.97	# %/100
ClearInsulIRAbsorb = 	0.2	# %/100	(a complete guess)
GroundReflection = 	0.05	# %/100	-- Reflectiveness of the ground south of the greenhouse
SnowReflection = 	0.9	# %/100
SnowPatchinessStart = 	0.03	# m -- The snow thickness for which patches start to show through
PercentAbsorbed = 	0.7	# %/100 -- Light absorbed in the greenhouse.
TransmissionLoss = 	0.75	# %/100 -- Light lost passing through the glazing.
WaterRadiationInside = 	0.75	# %/100 -- Radiation from water that is ultimately redirected *into* the greenhouse, not out of it.
GlazingReabsorb	= 	0.2	# %/100 -- Energy returned to inside through glazing IR/sunlight absorption.  Rarely about 50%.
GlazingReabsorbWater = 	0.1	# %/100	-- Percent of IR energy returned by the glazing that goes into the water.
GlazingReabsorbSoil =   0.2	# %/100 -- Percent of IR energy returned by the glazing that goes into the soil.
GreenhouseBlackbody =	0.7	# %/100 of the energy that a perfect blackbody would radiate that the greenhouse radiates.
WaterBlackbody = 	0.9	# %/100 of the energy that a perfect blackbody would radiate that the water radiates.
SoilBlackbody = 	0.8	# %/100 of the energy that a perfect blackbody would radiate that the soil radiates.
OutsideBlackbody = 	0.8	# %/100 of the energy that a perfect blackbody would radiate that the outside radiates.
HeaterWatts = 		2637.0	# Watts -- BTU/h * 0.293
HeaterPilotWatts = 	20.0	# Watts -- BTU/h * 0.293
HeaterThermostatOn = 	8.0	# Celcius
HeaterThermostatOff = 	12.0	# Celcius
HeaterEfficiency = 	0.9	# Percent of the heat from the heater that makes it into the inside air.
FuelEnergyDensity = 	5.0e7	# J/kg -- Propane = 5.0e7.
FuelQuantity = 		11.5	# Kilograms -- Pounds / 2.2046
InfiltrationArea = 	0.025	# m^2.
AvgHoleMinDiameter = 	0.07	# m -- give a bit of extra if holes are greatly elongated.
AvgHoleLength = 	0.8	# m
AvgBottomHoleHeight = 	0.3	# m
AvgTopHoleHeight = 	0.32	# m
VentThermostatOpen = 	41.0	# Celcius
VentThermostatClose = 	36.0	# Celcius
VentArea = 		0.03	# m^2
VentFanPower = 		0.02	# m^3/s
VentHeight = 		2.0	# m
SoilSpecificHeat = 	1100.0	# J/kg*K (Dry soil ~= 800; wet soil ~= 1500)
SoilDensity = 		1800.0	# kg/m^3 (Organic, silty clay ~= 1040; gravel and sand ~= 2200.  Typical ~= 2000)
SoilInsul_R = 		3.0	# K*m^2/(W*m) -- 0.1761 * American R-values per meter
SoilConductivity = 	0.9	# W/(m*C) -- 0.17 to 1.13
DayTemperatureCurve = [		# %/100 of daily high/low range.  This curve is most definately an approximation.
0.05,
0.0,	#It's coldest right *after* dawn.
0.15,
0.35,
0.5,
0.65,
0.75,
0.83,
0.9,
0.94,
0.98,
1.0,	#It's warmest *after* midday.
0.94,
0.87,
0.75,
0.65
]

NightTemperatureCurve = [	# %/100 of daily high/low range
0.65,
0.5,
0.4,
0.32,
0.25,
0.2,
0.15,
0.10,
0.05
]

ClimateData = climate_data.GetClimateData()

#for i in range(len(ClimateData)):	# Set fixed highs/lows
#  ClimateData[i]=(293.8, 280.4, 293.8, 280.4)

# Run defines
InsideTemp =	 	25.0	# Celcius
WaterTemp =	 	30.0	# Celcius
Latitude = 		40.5	# Degrees
Date = 			-20.0	# Days from the winter solstice (~Dec 22nd)
TemperatureExpPower = 	3.0	# The higher the number, the rarer the temperature extremes will be.
DayTempSmoothing = 	0.5	# %/100 -- how much weight to give to the previous day's temperature when looking at the current day's temp.
CloudTempchangeDependence = 18.0 # The higher the number, the less clouds track temperature changes.
CloudRandomPower = 	4.0	# The higher the number, the less randomness there is to cloudy days.
AverageCloudCover = 	0.3	# %/100 of the sky that's obscured by clouds on average.
AverageWindspeed = 	5.0	# m/s -- mph * 0.44704
WindExpPower = 		1.0	# The higher the number, the rarer the wind will be extreme.
WindAverageFactor = 	0.8	# How much to weigh how windy it was when determining how windy it is now.
WindUpdateFrequency = 	3629	# How many seconds between windspeed changes.
PrecipRandomPower = 	2.5	# Polarizes the random component of precipitation selection
PrecipChance = 		0.05	# Odds of being in a storm (as opposed to not being in a storm)
AvgStormLength = 	6000	# Average length of a storm
MaxStormStrength = 	0.00002	# Maximum precipitation in a storm, measured in meters of snow per second.
SnowSublimationScalar =	8e-10	# A scalar that helps determine how quickly sublimation of snow from sunlight occurs.
MeltRateFactor =	4e-7	# A scalar used to determine how quickly snow melts relative to the above-freezing temperature difference.
ObstructionStartAngle = 7.0	# Degrees -- The sun altitude in which it starts to be blocked by horizon features
ObstructionEndAngle = 	0.0	# Degrees -- The sun altitude in which the sun is completely blocked by horizon features
MinSolarIrradiance =    250.0   # W/m^2
MaxSolarIrradiance = 	1000.0	# W/m^2
TimeIncrement = 	30	# Seconds
CurTime =		43200	# Seconds
SoilStartingTemp = 	20.0	# Celcius
SoilDepth = 		2.0	# Meters.  The depth of the soil to calculate for.
SoilDepthIntervals = 	0.1	# Meters.
#SoilDiffusivity = 	1.6e-5	# m^2/s
#SoilHeatTransferCoeff = 40.0	# W/(m^2*C)
VentOpen = 		False	# True or false
HeaterOn = 		False	# True or false
DebugInside = 		False	# Boolean
DebugWater = 		False	# Boolean
DebugSoil = 		False	# Boolean
DebugWatts = 		False	# Boolean
DebugExtras = 		True	# Boolean

# Constants
Pi = 			3.14159265358979323
KelvinConvert = 	273.15	# Degrees Kelvin for the freezing point of water
StefanBoltzmannConst = 	5.6704e-8 # W*m^-2*K^-4
SpecHeatCapacityAir = 	1230.0	# J/(K*m^3)
SpecHeatCapacityWater = 4185500.0 # J/(K*m^3)
SpecHeatCapacityOther = 70000.0	# J/(K*m^3) -- Extra for the frame, pots, plants, etc.  We'll assume that these are temperature-equalized with the air.
AirViscosity = 		17.4e-6	# Poiseuille (Ns/m^2) (10 Poise = 1 Poiseuille)
AirPressure = 		101325.0 # Pascals at STP
AirDensity = 		1.2	# kg/m^3
AxialTilt = 		23.45	# Degrees
GravityAccel = 		9.81	# Meters per second
MinConvectionRate = 	0.0003	# m^3/(K*m^2*s) -- A complete and utter guess.  I don't want to have to get into things like Reynolds numbers, so this is a fudge factor to account for the mixing of two stationary bodies of air.
AirInfluxResistance = 	0.015	# A fudge factor to the resistance of air influx through vents and holes.  A lower number means more rapid airflow.

# Calculated defines
InsideTemp += KelvinConvert
WaterTemp += KelvinConvert
HeaterThermostatOn += KelvinConvert
HeaterThermostatOff += KelvinConvert
VentThermostatOpen += KelvinConvert
VentThermostatClose += KelvinConvert
ObstructionStartAngle *= Pi / 180
ObstructionEndAngle *= Pi / 180
Latitude *= Pi / 180
AxialTilt *= Pi / 180
HeightDiff = PeakHeight - VerticalWallHeight
Diagonal = sqrt(HeightDiff * HeightDiff + Width * Width / 4)
NorthVerticalWallArea = VerticalWallHeight * Length
NorthPeakArea = Diagonal * Length
NorthArea = NorthVerticalWallArea + NorthPeakArea
NorthPeakAngle = atan((Width / 2) / HeightDiff)
SouthVerticalWallArea = VerticalWallHeight * Length
SouthPeakArea = Diagonal * Length
SouthArea = SouthVerticalWallArea + SouthPeakArea
SouthPeakAngle = atan((Width / 2) / HeightDiff)
WestArea = VerticalWallHeight * Width + HeightDiff * Width / 2
EastArea = VerticalWallHeight * Width + HeightDiff * Width / 2
WestClearArea = WestArea * WestClearPercent
WestOpaqueArea = WestArea * (1.0 - WestClearPercent)
EastClearArea = EastArea * EastClearPercent
EastOpaqueArea = EastArea * (1.0 - EastClearPercent)
ClearArea = SouthArea + EastClearArea + WestClearArea
OpaqueArea = NorthArea + EastOpaqueArea + WestOpaqueArea
Area = OpaqueArea + ClearArea
GreenhouseVolume = Width * Length * (VerticalWallHeight + PeakHeight / 2)
GreenhouseHeatCapacity = (SpecHeatCapacityAir + SpecHeatCapacityOther) * GreenhouseVolume
WaterTopArea = WaterContainers * WaterDiameter * WaterDiameter / 4 * Pi
WaterSideArea = WaterContainers * WaterDiameter * WaterHeight * Pi
WaterCrossSection = Width * WaterHeight
WaterArea = WaterTopArea + WaterSideArea
WaterRadiationArea = WaterTopArea + 2 * WaterHeight * Length + WaterDiameter * WaterHeight * Pi
WaterVolume = WaterTopArea * WaterHeight
WaterHeatCapacity = SpecHeatCapacityWater * WaterVolume
SoilInsul_R = SoilDepthIntervals / SoilConductivity
SoilHeatCapacity = Width * Length * SoilSpecificHeat * SoilDensity * SoilDepthIntervals
SoilTemp = [SoilStartingTemp + KelvinConvert] * int(SoilDepth / SoilDepthIntervals)
GlazingReabsorbGreenhouse = 1.0 - GlazingReabsorbWater - GlazingReabsorbSoil
WaterRadiationOutside = 1.0 - WaterRadiationInside
WaterRadiationOutsideArea = WaterRadiationArea * WaterRadiationOutside
if WaterRadiationArea * WaterRadiationOutside > ClearArea:
  WaterRadiationOutsideArea = ClearArea
WaterRadiationInsideArea = WaterRadiationArea - WaterRadiationOutsideArea
WaterRadiationSoilArea = WaterArea

# // F U N C T I O N S //////////////////////////////////////////////////////////

def AreaD(d):
  return d * d / 4 * Pi
  
def AreaR(r):
  return r * r * Pi

def ComputeDeclination(AxialTilt, Date):
  return AxialTilt * -cos((Date / 356.24) * (2 * Pi))
  
def ComputeSunrise(Latitude, Declination):
  return 43200.0 - 86400.0 * acos(-sin(Latitude) * sin(Declination) / (cos(Latitude) * cos(Declination))) / (2 * Pi)

def ComputeSunset(Latitude, Declination):	
  return 86400.0 * acos(-sin(Latitude) * sin(Declination) / (cos(Latitude) * cos(Declination))) / (2 * Pi) + 43200.0
  
def Yesterday(Date):
  Date -= 1
  if Date < -183:
    Date += 365
  return Date
  
def Tomorrow(Date):
  Date += 1
  if Date > 182:
    Date -= 365
  return Date

def NormalizeDate(Date):
  newdate = Date + 9
  if newdate < 0:
    newdate += 365
  return int(newdate)

def GetHigh(Date):
  Date = NormalizeDate(Date)
  Rnd=pow(random.random(), TemperatureExpPower)
  if random.random() > 0.5:
    return ClimateData[Date][0] * (1.0 - Rnd) + ClimateData[Date][2] * Rnd
  else:
    Opposite = pow(ClimateData[Date][0], 2) / ClimateData[Date][2]
    return ClimateData[Date][0] * (1.0 - Rnd) + Opposite * Rnd

def GetLow(Date):
  Date = NormalizeDate(Date)
  Rnd=pow(random.random(), TemperatureExpPower)
  if random.random() > 0.5:
    return ClimateData[Date][1] * (1.0 - Rnd) + ClimateData[Date][3] * Rnd
  else:
    Opposite = pow(ClimateData[Date][1], 2) / ClimateData[Date][3]
    return ClimateData[Date][1] * (1.0 - Rnd) + Opposite * Rnd
    
def GetCloudiness(YesterdayHigh, TodayHigh, CloudTempchangeDependence, CloudRandomPower, AverageCloudCover):
  TempDependentCloudiness = 0.5 + ((YesterdayHigh - TodayHigh) / CloudTempchangeDependence)
  RandomCloudiness = pow(random.random(), CloudRandomPower)
  if random.random() > 0.5:
    RandomCloudiness = -RandomCloudiness
  Cloudiness = TempDependentCloudiness + RandomCloudiness
  Cloudiness += AverageCloudCover - 0.5
  if Cloudiness > 1:
    Cloudiness = 1
  elif Cloudiness < 0:
    Cloudiness = 0
  return Cloudiness

def ComputeOutsideTemp(CurTime, YesterdayHigh, YesterdayLow, TodayHigh, TodayLow, Sunrise, Sunset, DayTemperatureCurve, NightTemperatureCurve):
  Percent2 = 0.0
  PastPeak = False
  PastMin = False
  if CurTime > Sunrise and CurTime < Sunset:
    PercentOfDay = (CurTime - Sunrise) / (Sunset - Sunrise)
    Lower = int(PercentOfDay * len(DayTemperatureCurve))
    Percent = PercentOfDay * len(DayTemperatureCurve) - Lower
    Upper = Lower + 1
    if Upper == len(DayTemperatureCurve):	#Special case: we're precisely at sunset
      Percent2 = DayTemperatureCurve[-1]
    else:
      Percent2 = DayTemperatureCurve[Lower] * (1.0 - Percent) + DayTemperatureCurve[Upper] * Percent

    Min = 0
    for i in range(len(DayTemperatureCurve) - 1):
      if DayTemperatureCurve[i + 1] > DayTemperatureCurve[Min]:
        break
      Min = i
    Min += 1
    if Lower >= Min:
      PastMin = True

    if PastMin == True:
      Peak = Min
      for i in range(len(DayTemperatureCurve) - Min - 1):
        if DayTemperatureCurve[Peak + 1] < DayTemperatureCurve[Peak]:
          break
        Peak += 1
      Peak += 1
      if Lower >= Peak:
        PastPeak = True

  else:
    if CurTime < Sunrise:
      CurTime += 86400.0
    else:
      PastMin = True
      PastPeak = True
    Sunrise += 86400.0
    PercentOfNight = (CurTime - Sunset) / (Sunrise - Sunset)
    Lower = int(PercentOfNight * len(NightTemperatureCurve))
    Percent = PercentOfNight * len(NightTemperatureCurve) - Lower
    Upper = Lower + 1
    if Upper == len(NightTemperatureCurve):	#Special case: we're precisely at sunrise
      Percent2 = NightTemperatureCurve[-1]
    else:
      Percent2 = NightTemperatureCurve[Lower] * (1.0 - Percent) + NightTemperatureCurve[Upper] * Percent

  if PastPeak:
    return Percent2 * (TodayHigh - TodayLow) + TodayLow
  elif PastMin:
    return Percent2 * (TodayHigh - YesterdayLow) + YesterdayLow
  else:
    return Percent2 * (YesterdayHigh - YesterdayLow) + YesterdayLow

def ComputeConduction(Temp1, Temp2, Area, Insul_R):
  return (Temp1 - Temp2) * Area / Insul_R

def ComputeRadiation(Temp1, Temp2, Area, Blackbody1, Blackbody2):
  return StefanBoltzmannConst * Area * (pow(Temp1, 4) * Blackbody1 - pow(Temp2, 4) * Blackbody2)

def ComputeGreenhouseRadiation(InsideTemp, OutsideTemp, ClearArea, OpaqueArea, OpaqueInsulIRReflect, Blackbody1, Blackbody2):
  return StefanBoltzmannConst * (ClearArea + OpaqueArea * (1.0 - OpaqueInsulIRReflect)) * (pow(InsideTemp, 4) * Blackbody1 - pow(OutsideTemp, 4) * Blackbody2)

def ComputeSunAltitude(CurTime, Latitude, Declination):
  TimeAngle = (CurTime + 43200) / 86400.0 * (2 * Pi)
  return asin(sin(Latitude) * sin(Declination) + cos(Latitude) * cos(Declination) * cos(TimeAngle))

def ComputeSunAzimuth(CurTime, Declination, SunAltitude):
  #Noon is zero, negative is to the west, positive is to the east.
  TimeAngle = (CurTime + 43200) / 86400.0 * (2 * Pi)
  Azimuth = asin(-cos(Declination) * sin(TimeAngle) / cos(SunAltitude))
  #Note: This Azimuth isn't corrected for asin phase errors; however, they only occur at night, so it's not an issue.
  return Azimuth

def ComputeSolarObstructMultiplier(SunAltitude, ObstructionStartAngle, ObstructionEndAngle):
  if SunAltitude >= ObstructionStartAngle:
    return 1.0
  elif SunAltitude <= ObstructionEndAngle:
    return 0.0
  else:
    return (SunAltitude - ObstructionEndAngle) / (ObstructionStartAngle - ObstructionEndAngle)

def ComputeCurSolarIrradiance(MinSolarIrradiance, MaxSolarIrradiance, SolarObstructMultiplier, CurGroundReflection, Cloudiness):
  CloudReducedIrradiance = MinSolarIrradiance + (1.0 - Cloudiness) * (MaxSolarIrradiance - MinSolarIrradiance)
  return CloudReducedIrradiance * SolarObstructMultiplier * (1.0 + CurGroundReflection)

def ComputeWaterIrradiance(WaterTopArea, WaterCrossSection, SunAltitude, SunAzimuth, CurSolarIrradiance):
  WaterTopSunPercent = cos(SunAzimuth) * cos(SunAltitude)
  WaterCrossSectionSunPercent = cos(SunAzimuth) * cos(SunAltitude - (Pi / 2))
  return CurSolarIrradiance * (WaterTopArea * WaterTopSunPercent + WaterCrossSection * WaterCrossSectionSunPercent)

def ComputeGreenhouseIrradiance(SouthVerticalWallArea, SouthPeakArea, SouthPeakAngle, WestClearArea, EastClearArea, SunAltitude, SunAzimuth, CurSolarIrradiance):
  SouthVerticalWallSunPercent = cos(SunAzimuth) * cos(SunAltitude)
  SouthPeakSunPercent = cos(SunAzimuth) * cos(SunAltitude - SouthPeakAngle)
  if SunAzimuth < 0:
    WestSunPercent = cos(SunAzimuth + (Pi / 2)) * cos(SunAltitude)
    return CurSolarIrradiance * (WestClearArea * WestSunPercent + SouthVerticalWallArea * SouthVerticalWallSunPercent + SouthPeakArea * SouthPeakSunPercent)
  else:
    EastSunPercent = cos(SunAzimuth - (Pi / 2)) * cos(SunAltitude)
    return CurSolarIrradiance * (EastClearArea * EastSunPercent + SouthVerticalWallArea * SouthVerticalWallSunPercent + SouthPeakArea * SouthPeakSunPercent)
    
def SolvePipeAirflow(PipeDiameter, PipeLength, FluidViscosity, Pressure1, Pressure2, AirInfluxResistance):
  R_Squared = PipeDiameter * PipeDiameter / 4
  FlowArea = R_Squared * Pi
  FlowRate = FlowArea * (Pressure1 - Pressure2) * R_Squared / (8 * FluidViscosity * Length + AirInfluxResistance)
  return FlowRate

def SolveHeightDependentFlowRate(InsideTemp, OutsideTemp, AirDensity, GravityAccel, MaxHeight, MinHeight, PipeDiameter, PipeLength, AirViscosity, WindPressureBoost, AirInfluxResistance):
  InsideDensity = AirDensity * KelvinConvert / InsideTemp
  OutsideDensity = AirDensity * KelvinConvert / OutsideTemp
  InsidePressure = InsideDensity * GravityAccel * (MaxHeight - MinHeight)
  OutsidePressure = OutsideDensity * GravityAccel * (MaxHeight - MinHeight)
  HeightDependentFlowRate = SolvePipeAirflow(PipeDiameter, PipeLength * 2.0, AirViscosity, OutsidePressure + WindPressureBoost, InsidePressure, AirInfluxResistance)
  return HeightDependentFlowRate


def SolveFlowRate(InsideTemp, OutsideTemp, AirDensity, GravityAccel, InfiltrationArea, AvgHoleMinDiameter, AvgHoleLength, AvgTopHoleHeight, AvgBottomHoleHeight, AirViscosity, WindPressureBoost, AirInfluxResistance, MinFlowRate, FanRate):
  FlowRate = SolveHeightDependentFlowRate(InsideTemp, OutsideTemp, AirDensity, GravityAccel, AvgTopHoleHeight, AvgBottomHoleHeight, AvgHoleMinDiameter, AvgHoleLength, AirViscosity, WindPressureBoost * InfiltrationArea, AirInfluxResistance)
  FlowRate *= (InfiltrationArea / AreaD(AvgHoleMinDiameter)) / 2.0  
  FlowRate += InfiltrationArea * MinFlowRate * (InsideTemp - OutsideTemp)
  FlowRate += FanRate
  return FlowRate

def GetWindspeed(AverageWindspeed, Windspeed, WindExpPower, WindAverageFactor):
  Rnd = pow(random.random(), WindExpPower)
  if random.random() < pow(0.5, WindExpPower):
    Rnd = 1.0 / (Rnd + 0.1)
  NewWind = Rnd * AverageWindspeed
  return Windspeed * WindAverageFactor + NewWind * (1.0 - WindAverageFactor)

def GetWindPressureBoost(Windspeed):
  return random.random() * Windspeed
  
def TempAtHeight(AvgTemp, Height, MaxHeight):
  return AvgTemp	# TODO: I don't know how to do this without a lot of trouble.
    
def GetPrecipitationRate(Cloudiness, Precipitation, PrecipRandomPower, PrecipChance, AvgStormLength, MaxStormStrength, TimeIncrement):
  if Precipitation == 0:
    if random.random() < PrecipChance * TimeIncrement / AvgStormLength:
      Precipitation = pow(random.random(), PrecipRandomPower) * MaxStormStrength
  else:
    if random.random() < (1.0 - PrecipChance) * TimeIncrement / AvgStormLength:
      Precipitation = 0;
  return Precipitation
  
def HandleSnowChange(OutsideTemp, PrecipitationRate, CurSolarIrradiance, Cloudiness, SunAltitude, SnowSublimationScalar, MeltRateFactor):
  SnowChange = 0.0
  if OutsideTemp > KelvinConvert:
    SnowChange -= MeltRateFactor * (OutsideTemp - KelvinConvert)
  if PrecipitationRate > 0:
    if OutsideTemp < KelvinConvert:
      SnowChange += PrecipitationRate
    else:
      SnowChange -= PrecipitationRate
  SnowChange -= SnowSublimationScalar * CurSolarIrradiance * (1.0 - Cloudiness) * sin(SunAltitude)
  return SnowChange

def GetGroundReflection(SnowAccumulation, GroundReflection, SnowReflection, SnowPatchinessStart):
  if SnowAccumulation > SnowPatchinessStart:
    return SnowReflection
  else:
    percent = SnowAccumulation / SnowPatchinessStart
    return SnowReflection * percent + GroundReflection * (1.0 - percent)

def GetDateStr(Date):
  newdate = NormalizeDate(Date)
  month = ""
  day = 0
  if newdate < 31:
    month = "January"
    day = newdate
  elif newdate < (31 + 28):
    month = "February"
    day = newdate - 31
  elif newdate < (31 + 28 + 31):
    month = "March"
    day = newdate - (31 + 28)
  elif newdate < (31 + 28 + 31 + 30):
    month = "April"
    day = newdate - (31 + 28 + 31)
  elif newdate < (31 + 28 + 31 + 30 + 31):
    month = "May"
    day = newdate - (31 + 28 + 31 + 30)
  elif newdate < (31 + 28 + 31 + 30 + 31 + 30):
    month = "June"
    day = newdate - (31 + 28 + 31 + 30 + 31)
  elif newdate < (31 + 28 + 31 + 30 + 31 + 30 + 31):
    month = "July"
    day = newdate - (31 + 28 + 31 + 30 + 31 + 30)
  elif newdate < (31 + 28 + 31 + 30 + 31 + 30 + 31 + 31):
    month = "August"
    day = newdate - (31 + 28 + 31 + 30 + 31 + 30 + 31)
  elif newdate < (31 + 28 + 31 + 30 + 31 + 30 + 31 + 31 + 30):
    month = "September"
    day = newdate - (31 + 28 + 31 + 30 + 31 + 30 + 31 + 31)
  elif newdate < (31 + 28 + 31 + 30 + 31 + 30 + 31 + 31 + 30 + 31):
    month = "October"
    day = newdate - (31 + 28 + 31 + 30 + 31 + 30 + 31 + 31 + 30)
  elif newdate < (31 + 28 + 31 + 30 + 31 + 30 + 31 + 31 + 30 + 31 + 30):
    month = "November"
    day = newdate - (31 + 28 + 31 + 30 + 31 + 30 + 31 + 31 + 30 + 31)
  else:
    month = "December"
    day = newdate - (31 + 28 + 31 + 30 + 31 + 30 + 31 + 31 + 30 + 31 + 30)
    
  day += 1	# Make it not be zero-indexed.
    
  return "%s %2d" % (month[:3], day)

def GetTimeStr(CurTime):
  CurTime = int(CurTime)
  suffix = ""
  if CurTime >= 43200:
    suffix = "PM"
    CurTime -= 43200
  else:
    suffix = "AM"
  hours = CurTime / 3600
  CurTime -= hours * 3600
  minutes = CurTime / 60
  CurTime -= minutes * 60
  seconds = CurTime
  if hours == 0:
    hours = 12
  return "%02d:%02d:%02d %s" % (hours, minutes, seconds, suffix)

# // M A I N ////////////////////////////////////////////////////////////////////

Declination = ComputeDeclination(AxialTilt, Date)
Sunrise = ComputeSunrise(Latitude, Declination)
Sunset = ComputeSunset(Latitude, Declination)
TodayHigh = GetHigh(Date)
TodayLow = GetLow(Date)
YesterdayHigh = GetHigh(Yesterday(Date))
YesterdayLow = GetLow(Yesterday(Date))
Cloudiness = GetCloudiness(YesterdayHigh, TodayHigh, CloudTempchangeDependence, CloudRandomPower, AverageCloudCover)
Windspeed = AverageWindspeed
WindPressureBoost = GetWindPressureBoost(Windspeed)
SnowAccumulation = 0.0
PrecipitationRate = 0.0
while True:
  CurTime += TimeIncrement
  if CurTime >= 86400:
    CurTime = 0
    Date = Tomorrow(Date)
    Declination = ComputeDeclination(AxialTilt, Date)
    Sunrise = ComputeSunrise(Latitude, Declination)
    Sunset = ComputeSunset(Latitude, Declination)
    YesterdayHigh = TodayHigh
    YesterdayLow = TodayLow
    TodayHigh = GetHigh(Date) * (1.0 - DayTempSmoothing) + YesterdayHigh * DayTempSmoothing
    TodayLow = GetLow(Date) * (1.0 - DayTempSmoothing) + YesterdayLow * DayTempSmoothing
    Cloudiness = GetCloudiness(YesterdayHigh, TodayHigh, CloudTempchangeDependence, CloudRandomPower, AverageCloudCover)

  OutsideTemp = ComputeOutsideTemp(CurTime, YesterdayHigh, YesterdayLow, TodayHigh, TodayLow, Sunrise, Sunset, DayTemperatureCurve, NightTemperatureCurve)
  
  OldInsideTemp=InsideTemp
  OldWaterTemp=WaterTemp
  OldSoilTemp=SoilTemp[0]

  WaterHeatConduction = ComputeConduction(WaterTemp, InsideTemp, WaterArea, WaterInsul_R)
  if DebugInside:
    print "Water ->Conduction-> Inside:\t", WaterHeatConduction
  InsideTemp += WaterHeatConduction * TimeIncrement / GreenhouseHeatCapacity
  if DebugWater:
    print "Inside ->Conduction-> Water:\t", -WaterHeatConduction
  WaterTemp -= WaterHeatConduction * TimeIncrement / WaterHeatCapacity

  GlazingHeatConduction = ComputeConduction(InsideTemp, OutsideTemp, OpaqueArea, OpaqueInsul_R)
  GlazingHeatConduction += ComputeConduction(InsideTemp, OutsideTemp, ClearArea, ClearInsul_R)
  if DebugInside:
    print "Outside ->Conduction-> Inside:\t", -GlazingHeatConduction
  InsideTemp -= GlazingHeatConduction * TimeIncrement / GreenhouseHeatCapacity
  
  WaterRadiation = ComputeRadiation(WaterTemp, InsideTemp, WaterRadiationInsideArea, WaterBlackbody, GreenhouseBlackbody)
  if DebugWater:
    print "Inside ->Radiation-> Water:\t", -WaterRadiation
  WaterTemp -= WaterRadiation * TimeIncrement / WaterHeatCapacity
  if DebugInside:
    print "Water ->Radiation-> Inside:\t", WaterRadiation
  InsideTemp += WaterRadiation * TimeIncrement / GreenhouseHeatCapacity

  WaterRadiation = ComputeRadiation(WaterTemp, OutsideTemp, WaterRadiationOutsideArea, WaterBlackbody, OutsideBlackbody)
  if DebugWater:
    print "Outside ->Radiation-> Water:\t", -WaterRadiation
  WaterTemp -= WaterRadiation * TimeIncrement / WaterHeatCapacity
  if DebugInside:
    print "Outside, Water ->RadiationReabsorb-> Inside:\t", WaterRadiation * GlazingReabsorb * GlazingReabsorbGreenhouse
  InsideTemp += WaterRadiation * GlazingReabsorb * GlazingReabsorbGreenhouse * TimeIncrement / GreenhouseHeatCapacity
  if DebugWater:
    print "Outside, Water ->RadiationReabsorb-> Water:\t", WaterRadiation * GlazingReabsorb * GlazingReabsorbWater
  WaterTemp += WaterRadiation * GlazingReabsorb * GlazingReabsorbWater * TimeIncrement / WaterHeatCapacity
  if DebugSoil:
    print "Outside, Water ->RadiationReabsorb-> Soil:\t", WaterRadiation * GlazingReabsorb * GlazingReabsorbSoil
  SoilTemp[0] += WaterRadiation * GlazingReabsorb * GlazingReabsorbSoil * TimeIncrement / SoilHeatCapacity

  GreenhouseRadiation = ComputeGreenhouseRadiation(InsideTemp, OutsideTemp, ClearArea - WaterCrossSection, OpaqueArea, OpaqueInsulIRReflect, GreenhouseBlackbody, OutsideBlackbody)
  if DebugInside:
    print "Outside ->Radiation-> Inside:\t", -GreenhouseRadiation
  InsideTemp -= GreenhouseRadiation * TimeIncrement / GreenhouseHeatCapacity
  if DebugInside:
    print "Outside, Inside ->RadiationReabsorb-> Inside:\t", GreenhouseRadiation * GlazingReabsorb * GlazingReabsorbGreenhouse
  InsideTemp += GreenhouseRadiation * GlazingReabsorb * GlazingReabsorbGreenhouse * TimeIncrement / GreenhouseHeatCapacity
  if DebugWater:
    print "Outside, Inside ->RadiationReabsorb-> Water:\t", GreenhouseRadiation * GlazingReabsorb * GlazingReabsorbWater
  WaterTemp += GreenhouseRadiation * GlazingReabsorb * GlazingReabsorbSoil * TimeIncrement / WaterHeatCapacity
  if DebugSoil:
    print "Outside, Inside ->RadiationReabsorb-> Soil:\t", GreenhouseRadiation * GlazingReabsorb * GlazingReabsorbSoil
  SoilTemp[0] += GreenhouseRadiation * GlazingReabsorb * GlazingReabsorbSoil * TimeIncrement / SoilHeatCapacity

  GreenhouseSoilConduction = ComputeConduction(SoilTemp[0], InsideTemp, Width * Length - WaterTopArea, SoilInsul_R / 2.0 + GroundCoverInsul_R)
  if DebugInside:
    print "Soil ->Conduction-> Inside:\t", GreenhouseSoilConduction
  InsideTemp += GreenhouseSoilConduction * TimeIncrement / GreenhouseHeatCapacity
  if DebugSoil:
    print "Inside ->Conduction-> Soil:\t", -GreenhouseSoilConduction
  SoilTemp[0] -= GreenhouseSoilConduction * TimeIncrement / SoilHeatCapacity

  WaterSoilConduction = ComputeConduction(SoilTemp[0], WaterTemp, WaterArea, SoilInsul_R / 2.0 + WaterInsul_R + GroundCoverInsul_R)
  if DebugWater:
    print "Soil ->Conduction-> Water:\t", WaterSoilConduction
  WaterTemp += WaterSoilConduction * TimeIncrement / WaterHeatCapacity
  if DebugSoil:
    print "Water ->Conduction-> Soil:\t", -WaterSoilConduction
  SoilTemp[0] -= WaterSoilConduction * TimeIncrement / SoilHeatCapacity
  
  GreenhouseSoilRadiation = ComputeRadiation(SoilTemp[0], InsideTemp, (Width * Length - WaterTopArea) * OpaqueArea / Area, SoilBlackbody, GreenhouseBlackbody)
  if DebugInside:
    print "Soil ->Radiation-> Inside:\t", GreenhouseSoilRadiation
  InsideTemp += GreenhouseSoilRadiation * TimeIncrement / GreenhouseHeatCapacity
  if DebugSoil:
    print "Inside ->Radiation-> Soil:\t", -GreenhouseSoilRadiation
  SoilTemp[0] -= GreenhouseSoilRadiation * TimeIncrement / SoilHeatCapacity

  OutsideSoilRadiation = ComputeRadiation(SoilTemp[0], OutsideTemp, (Width * Length - WaterTopArea) * (ClearArea - WaterCrossSection) / Area, SoilBlackbody, OutsideBlackbody)
  if DebugSoil:
    print "Outside ->Radiation-> Soil:\t", -OutsideSoilRadiation
  SoilTemp[0] -= GreenhouseSoilRadiation * TimeIncrement / SoilHeatCapacity
  if DebugInside:
    print "Outside, Soil ->RadiationReabsorb-> Inside:\t", OutsideSoilRadiation * GlazingReabsorb * GlazingReabsorbGreenhouse
  InsideTemp += GreenhouseSoilRadiation * GlazingReabsorb * GlazingReabsorbGreenhouse * TimeIncrement / GreenhouseHeatCapacity
  if DebugWater:
    print "Outside, Soil ->RadiationReabsorb-> Water:\t", OutsideSoilRadiation * GlazingReabsorb * GlazingReabsorbWater
  WaterTemp += GreenhouseSoilRadiation * GlazingReabsorb * GlazingReabsorbWater * TimeIncrement / WaterHeatCapacity
  if DebugSoil:
    print "Outside, Soil ->RadiationReabsorb-> Soil:\t", OutsideSoilRadiation * GlazingReabsorb * GlazingReabsorbSoil
  SoilTemp[0] += GreenhouseSoilRadiation * GlazingReabsorb * GlazingReabsorbSoil * TimeIncrement / SoilHeatCapacity

  WaterSoilRadiation = ComputeRadiation(SoilTemp[0], WaterTemp, WaterTopArea + WaterCrossSection, SoilBlackbody, WaterBlackbody)
  if DebugSoil:
    print "Water ->Radiation-> Soil:\t", -WaterSoilRadiation
  SoilTemp[0] -= WaterSoilRadiation * TimeIncrement / SoilHeatCapacity
  if DebugWater:
    print "Soil ->Radiation-> Water:\t", WaterSoilRadiation
  WaterTemp += WaterSoilRadiation * TimeIncrement / WaterHeatCapacity
  
  OrigSoilTemp = SoilTemp[:]
  for i in range(len(SoilTemp) - 1):
    InterlayerConduction = ComputeConduction(OrigSoilTemp[i], OrigSoilTemp[i+1], Width * Length, SoilInsul_R)
    SoilTemp[i+1] += InterlayerConduction * TimeIncrement / SoilHeatCapacity
    SoilTemp[i] -= InterlayerConduction * TimeIncrement / SoilHeatCapacity
    
    #For conduction to the outside, assume that the distance of soil that insulates is related to the diagonal to the surface.
    OutsideConduction = ComputeConduction(OrigSoilTemp[i], OutsideTemp, (Length + Width) * 2 * SoilDepthIntervals, (i + 0.5) * SoilDepthIntervals * sqrt(2.0) / SoilConductivity)
    SoilTemp[i] -= OutsideConduction * TimeIncrement / SoilHeatCapacity
    
  CurGroundReflection = GetGroundReflection(SnowAccumulation, GroundReflection, SnowReflection, SnowPatchinessStart)

  SunAltitude = ComputeSunAltitude(CurTime, Latitude, Declination)
  SunAzimuth = ComputeSunAzimuth(CurTime, Declination, SunAltitude)
  SolarObstructMultiplier = ComputeSolarObstructMultiplier(SunAltitude, ObstructionStartAngle, ObstructionEndAngle)
  CurSolarIrradiance = ComputeCurSolarIrradiance(MinSolarIrradiance, MaxSolarIrradiance, SolarObstructMultiplier, CurGroundReflection, Cloudiness)
  WaterIrradiance = ComputeWaterIrradiance(WaterTopArea, WaterCrossSection, SunAltitude, SunAzimuth, CurSolarIrradiance)
  if DebugWater:
    print "Sun ->Radiation-> Water:\t", WaterIrradiance * 1.0 * TransmissionLoss
  WaterTemp += (WaterIrradiance * 1.0 * TransmissionLoss) * TimeIncrement / WaterHeatCapacity
  if DebugInside:
    print "Sun, Water ->RadiationReabsorb-> Inside:\t", WaterIrradiance * (1.0 - 1.0 * TransmissionLoss) * GlazingReabsorb * GlazingReabsorbGreenhouse
  InsideTemp += WaterIrradiance * (1.0 - 1.0 * TransmissionLoss) * GlazingReabsorb * GlazingReabsorbGreenhouse * TimeIncrement / GreenhouseHeatCapacity
  if DebugWater:
    print "Sun, Water ->RadiationReabsorb-> Water:\t", WaterIrradiance * (1.0 - 1.0 * TransmissionLoss) * GlazingReabsorb * GlazingReabsorbWater
  WaterTemp += WaterIrradiance * (1.0 - 1.0 * TransmissionLoss) * GlazingReabsorb * GlazingReabsorbWater * TimeIncrement / GreenhouseHeatCapacity
  if DebugSoil:
    print "Sun, Water ->RadiationReabsorb-> Soil:\t", WaterIrradiance * (1.0 - 1.0 * TransmissionLoss) * GlazingReabsorb * GlazingReabsorbSoil
  SoilTemp[0] += WaterIrradiance * (1.0 - 1.0 * TransmissionLoss) * GlazingReabsorb * GlazingReabsorbSoil * TimeIncrement / GreenhouseHeatCapacity

  GreenhouseIrradiance = ComputeGreenhouseIrradiance(SouthVerticalWallArea, SouthPeakArea, SouthPeakAngle, WestClearArea - WaterCrossSection, EastClearArea, SunAltitude, SunAzimuth, CurSolarIrradiance)
  if DebugInside:
    print "Sun ->Radiation-> Inside:\t", GreenhouseIrradiance * PercentAbsorbed * TransmissionLoss
  InsideTemp += (GreenhouseIrradiance * PercentAbsorbed * TransmissionLoss) * TimeIncrement / GreenhouseHeatCapacity
  if DebugInside:
    print "Sun, Inside ->RadiationReabsorb-> Inside:\t", GreenhouseIrradiance * (1.0 - PercentAbsorbed * TransmissionLoss) * GlazingReabsorb * GlazingReabsorbGreenhouse
  InsideTemp += GreenhouseIrradiance * (1.0 - PercentAbsorbed * TransmissionLoss) * GlazingReabsorb * GlazingReabsorbGreenhouse * TimeIncrement / GreenhouseHeatCapacity
  if DebugWater:
    print "Sun, Inside ->RadiationReabsorb-> Water:\t", GreenhouseIrradiance * (1.0 - PercentAbsorbed * TransmissionLoss) * GlazingReabsorb * GlazingReabsorbWater
  WaterTemp += GreenhouseIrradiance * (1.0 - PercentAbsorbed * TransmissionLoss) * GlazingReabsorb * GlazingReabsorbWater * TimeIncrement / GreenhouseHeatCapacity
  if DebugSoil:
    print "Sun, Inside ->RadiationReabsorb-> Soil:\t", GreenhouseIrradiance * (1.0 - PercentAbsorbed * TransmissionLoss) * GlazingReabsorb * GlazingReabsorbSoil
  SoilTemp[0] += GreenhouseIrradiance * (1.0 - PercentAbsorbed * TransmissionLoss) * GlazingReabsorb * GlazingReabsorbSoil * TimeIncrement / GreenhouseHeatCapacity
  
  # Assume for now, with the low sun angles, that the floor of the greenhouse
  # never gets sunlight because it's blocked by the barrels on the south side.
  #Continue onward.
      
  if HeaterOn == False and InsideTemp <= HeaterThermostatOn and FuelQuantity > 0:
    HeaterOn = True
  elif HeaterOn == True and InsideTemp >= HeaterThermostatOff:
    HeaterOn = False
  
  if HeaterOn:
    if DebugInside:
      print "Heater ->On-> Inside:\t", (HeaterWatts + HeaterPilotWatts) * HeaterEfficiency * TimeIncrement / GreenhouseHeatCapacity
    InsideTemp += (HeaterWatts + HeaterPilotWatts) * HeaterEfficiency * TimeIncrement / GreenhouseHeatCapacity
    FuelQuantity -= (HeaterWatts + HeaterPilotWatts) * TimeIncrement / FuelEnergyDensity
    if FuelQuantity < 0:
      print "--== Heater fuel has run out! ==--"
  else:
    if DebugInside:
      print "Heater ->Pilot-> Inside:\t", HeaterPilotWatts * HeaterEfficiency * TimeIncrement / GreenhouseHeatCapacity
    InsideTemp += HeaterPilotWatts * HeaterEfficiency * TimeIncrement / GreenhouseHeatCapacity
    FuelQuantity -= HeaterPilotWatts * TimeIncrement / FuelEnergyDensity
    if FuelQuantity < 0:
      print "--== Heater fuel has run out! ==--"

  if CurTime % WindUpdateFrequency == 0:
    Windspeed = GetWindspeed(AverageWindspeed, Windspeed, WindExpPower, WindAverageFactor)
    WindPressureBoost = GetWindPressureBoost(Windspeed)
      
  PrecipitationRate = GetPrecipitationRate(Cloudiness, PrecipitationRate, PrecipRandomPower, PrecipChance, AvgStormLength, MaxStormStrength, TimeIncrement)
  
  SnowAccumulation += HandleSnowChange(OutsideTemp, PrecipitationRate, CurSolarIrradiance, Cloudiness, SunAltitude, SnowSublimationScalar, MeltRateFactor) * TimeIncrement
  if SnowAccumulation < 0:
    SnowAccumulation = 0.0

  AirFlow = 0
  TopInsideTemp = None
  BottomInsideTemp = None
  if VentOpen == False and InsideTemp >= VentThermostatOpen:
    VentOpen = True
  elif VentOpen == True and InsideTemp <= VentThermostatClose:
    VentOpen = False
    
  if VentOpen == True:
    CombinedArea = InfiltrationArea + VentArea
    NewAvgTopHeight = VentHeight * (VentArea / CombinedArea) + AvgTopHoleHeight * (InfiltrationArea / CombinedArea)
    NewAvgMinDiameter = 2.0 * sqrt(VentArea / Pi) * (VentArea / CombinedArea) + AvgHoleMinDiameter * (InfiltrationArea / CombinedArea)
    BottomInsideTemp = TempAtHeight(InsideTemp, AvgBottomHoleHeight, PeakHeight)
    TopInsideTemp = TempAtHeight(InsideTemp, NewAvgTopHeight, PeakHeight)
    AirFlow = SolveFlowRate(BottomInsideTemp, OutsideTemp, AirDensity, GravityAccel, CombinedArea, NewAvgMinDiameter, AvgHoleLength, NewAvgTopHeight, AvgBottomHoleHeight, AirViscosity, WindPressureBoost, AirInfluxResistance, MinConvectionRate, VentFanPower)
  else:
    BottomInsideTemp = TempAtHeight(InsideTemp, AvgBottomHoleHeight, PeakHeight)
    TopInsideTemp = TempAtHeight(InsideTemp, AvgTopHoleHeight, PeakHeight)
    AirFlow = SolveFlowRate(BottomInsideTemp, OutsideTemp, AirDensity, GravityAccel, InfiltrationArea, AvgHoleMinDiameter, AvgHoleLength, AvgTopHoleHeight, AvgBottomHoleHeight, AirViscosity, WindPressureBoost, AirInfluxResistance, MinConvectionRate, 0.0)

  HeatLoss = (TopInsideTemp - OutsideTemp) * AirFlow * SpecHeatCapacityAir
  if DebugInside:
    print "Outside ->Infiltration -> Inside:\t", -HeatLoss
  InsideTemp -= HeatLoss * TimeIncrement / GreenhouseHeatCapacity

  soilstr=""
  for temp in SoilTemp:
    soilstr+="%.0fC, " % (temp - KelvinConvert)
  outstr=GetDateStr(Date)+" "+GetTimeStr(CurTime)+": Outside=%.0fC, Inside=%.0fC, Water=%.0fC, Soil=" % (OutsideTemp-KelvinConvert, InsideTemp-KelvinConvert, WaterTemp-KelvinConvert)+soilstr[:-2]
  print outstr[0:80]
  if DebugWatts:
    InsideWatts=(InsideTemp - OldInsideTemp) / TimeIncrement * GreenhouseHeatCapacity
    WaterWatts=(WaterTemp - OldWaterTemp) / TimeIncrement * WaterHeatCapacity
    SoilWatts=(SoilTemp[0] - OldSoilTemp) / TimeIncrement * SoilHeatCapacity
    print "Inside %.1fW; Water %.1fW; Soil %.1fW" % (InsideWatts, WaterWatts, SoilWatts)
  if DebugExtras:
    HeaterStr = "OFF"
    if HeaterOn:
      HeaterStr = "ON"
    VentStr = "CLOSED"
    if VentOpen:
      VentStr = "OPEN"
    print "Vent %s; Heater %s; Fuel %.1fkg; Insolation %.1fW/m^2; Cloudiness %.0f%s; Wind %.1fm/s; Airflow %.3fm^3/s; Precipitation %.1fcm/h; SnowAccumulation %.1fcm; GroundReflection %.0f%s" % (VentStr, HeaterStr, FuelQuantity, CurSolarIrradiance, Cloudiness * 100, "%", Windspeed, AirFlow, PrecipitationRate * 360000, SnowAccumulation * 100, CurGroundReflection * 100, "%")