# -*- coding: utf-8 -*-
#######################################################
# author :Kensaku Okada [kensakuokada@email.arizona.edu]
# create date : 06 Nov 2016
# last edit date: 14 Dec 2016
#######################################################

##########import package files##########
from scipy import stats
import datetime
import sys
import os as os
import numpy as np
import ShadingCurtain
import matplotlib.pyplot as plt
import math
import CropElectricityYeildSimulatorConstant as constant
import Util
import OPVFilm
import Lettuce
import PlantGrowthModelE_J_VanHenten
import PlantGrowthModelS_Pearson1997
import SolarIrradianceMultiSpanRoof
from dateutil.relativedelta import *
import GreenhouseEnergyBalance as energyBalance

#######################################################

# def setSimulationSpecifications(simulatorClass):
#     '''
#     reference of the model:
#     :return:
#     '''

def calcOPVmoduleSolarIrradianceGHRoof(simulatorClass, roofDirectionNotation=constant.roofDirectionNotation):
    '''
    calculate the every kind of solar irradiance (W/m^2) to the OPV film. The tilt angel and direction angle  of OPV panes is define in CropElectricityYeildSimulatorConstant.py
    Reference of the model: A. yano et. at., 2009, "Electrical energy generated by photovoltaic modules mounted inside the roof of a north–south oriented greenhouse" and M. kacira et. al.,  ""
    Reference URL: "https://www.actahort.org/books/1037/1037_9.htm" and "https://www.sciencedirect.com/science/article/pii/S0960148104000060"
    '''

    year = simulatorClass.getYear()
    month = simulatorClass.getMonth()
    day = simulatorClass.getDay()
    hour = simulatorClass.getHour()
    hourlyHorizontalDiffuseOuterSolarIrradiance = simulatorClass.getImportedHourlyHorizontalDiffuseSolarRadiation()
    hourlyHorizontalDirectOuterSolarIrradiance = simulatorClass.getImportedHourlyHorizontalDirectSolarRadiation()

    # [rad] symbol: delta
    hourlyDeclinationAngle = OPVFilm.calcDeclinationAngle(year, month, day)
    # print "hourlyDecl	inationAngle:{}".format(np.degrees(hourlyDeclinationAngle))
    # print "hourlyDeclinationAngle:{}".format(hourlyDeclinationAngle)

    # [rad] symbol: omega
    hourlySolarHourAngle = OPVFilm.getSolarHourAngleKacira2003(hour)
    # print ("hourlySolarHourAngle by Kacira2003:{}".format(np.degrees(hourlySolarHourAngle)))

    # # [rad] symbol: omega
    # hourlySolarHourAngle = OPVFilm.getSolarHourAngleYano2009(hour)
    # print ("hourlySolarHourAngle by Yano2009:{}".format(np.degrees(hourlySolarHourAngle)))

    # [rad] symbol: alpha. elevation angle = altitude angle
    hourlySolarAltitudeAngle = OPVFilm.calcSolarAltitudeAngle(hourlyDeclinationAngle, hourlySolarHourAngle)
    # print "np.degrees(hourlySolarAltitudeAngle):{}".format(np.degrees(hourlySolarAltitudeAngle))
    # print "hourlySolarAltitudeAngle:{}".format(hourlySolarAltitudeAngle)

    # set the solar altitude angle, which is necessary to calculate the solar radiation through multispan roof
    simulatorClass.hourlySolarAltitudeAngle = hourlySolarAltitudeAngle

    # [rad] symbol: beta. azimuth angle
    hourlySolarAzimuthAngle = OPVFilm.calcSolarAzimuthAngle(hourlyDeclinationAngle, hourlySolarAltitudeAngle, hourlySolarHourAngle)
    # print "hourlySolarAzimuthAngle:{}".format(hourlySolarAzimuthAngle)

    # set the solar azimuth angle, which is necessary to calculate the solar radiation through multispan roof
    simulatorClass.hourlySolarAzimuthAngle = hourlySolarAzimuthAngle

    # used only in Kacira 2003
    # [rad] symbol: theta_z
    hourlyZenithAngle = math.radians(90.0) - hourlySolarAltitudeAngle
    # print "math.radians(90.0):{}".format(math.radians(90.0))
    # print "hourlyZenithAngle:{}".format(hourlyZenithAngle)

    # if the direction of greenhouse is north-south and the roof tilt direction is east-west
    if roofDirectionNotation == "EastWestDirectionRoof":
        # [rad] symbol: phi_p
        # module azimuth angle (yano 2009) == surface azimuth angle (kacira 2003)
        # if the OPV module facing east
        hourlyModuleAzimuthAngleEast = math.radians(-90.0)
        # hourlyModuleAzimuthAngleEast = math.radians(180.0)
        # if the OPV module facing west
        hourlyModuleAzimuthAngleWest = math.radians(90.0)
        # hourlyModuleAzimuthAngleWest = math.radians(0.0)
    # if the direction of greenhouse is east-west and the roof tilt direction is north-south
    elif roofDirectionNotation == "NorthSouthDirectionRoof":
        hourlyModuleAzimuthAngleNorth = math.radians(180.0)
        # if the OPV module facing west
        hourlyModuleAzimuthAngleSouth = math.radians(0.0)

    # set the module azimuth angle, which is necessary to calculate the solar radiation through multispan roof
    simulatorClass.hourlyModuleAzimuthAngleEast = hourlyModuleAzimuthAngleEast
    simulatorClass.hourlyModuleAzimuthAngleWest = hourlyModuleAzimuthAngleWest

    # this computation is necessary to calculate the horizontal incidence angle for horizontal direct solar irradiance. This data is used at getDirectHorizontalSolarRadiation function
    hourlyModuleAzimuthAngleSouth = math.radians(0.0)
    hourlyHorizontalSolarIncidenceAngle = OPVFilm.calcSolarIncidenceAngleYano2009(hourlySolarAltitudeAngle, hourlySolarAzimuthAngle, hourlyModuleAzimuthAngleSouth, 0)
    # print "hourlyHorizontalSolarIncidenceAngle:{}".format(hourlyHorizontalSolarIncidenceAngle)

    if roofDirectionNotation == "EastWestDirectionRoof":
        #The incident angle of the beam sunlight on the module surface. [rad] symbol: theta_I
        hourlySolarIncidenceAngleEastDirection = OPVFilm.calcSolarIncidenceAngleYano2009(hourlySolarAltitudeAngle, hourlySolarAzimuthAngle, hourlyModuleAzimuthAngleEast)
        hourlySolarIncidenceAngleWestDirection = OPVFilm.calcSolarIncidenceAngleYano2009(hourlySolarAltitudeAngle, hourlySolarAzimuthAngle, hourlyModuleAzimuthAngleWest)
        # print("hourlySolarIncidenceAngleEastDirection:{}".format(hourlySolarIncidenceAngleEastDirection))
        # print("hourlySolarIncidenceAngleWestDirection:{}".format(hourlySolarIncidenceAngleWestDirection))
    # if the direction of greenhouse is east-west and the roof tilt direction is north-south
    elif roofDirectionNotation == "NorthSouthDirectionRoof":
        # The suitability of the output value is not examined because our greenhouse was "EastWestDirectionRoof" (= north-south direction greenhouse)
        hourlySolarIncidenceAngleEastDirection = OPVFilm.calcSolarIncidenceAngleYano2009(hourlySolarAltitudeAngle, hourlySolarAzimuthAngle, hourlyModuleAzimuthAngleNorth)
        hourlySolarIncidenceAngleWestDirection = OPVFilm.calcSolarIncidenceAngleYano2009(hourlySolarAltitudeAngle, hourlySolarAzimuthAngle, hourlyModuleAzimuthAngleSouth)
    # print ("hourlySolarIncidenceAngleEastDirection:{}".format(hourlySolarIncidenceAngleEastDirection))
    # print ("hourlySolarIncidenceAngleWestDirection:{}".format(hourlySolarIncidenceAngleWestDirection))

    # set the incidence angle
    simulatorClass.hourlySolarIncidenceAngleEastDirection = hourlySolarIncidenceAngleEastDirection
    simulatorClass.hourlySolarIncidenceAngleWestDirection = hourlySolarIncidenceAngleWestDirection

    # np.set_printoptions(threshold=np.inf)
    # print("hourlySolarIncidenceAngleEastDirection:{}".format(hourlySolarIncidenceAngleEastDirection))
    # print("hourlySolarIncidenceAngleWestDirection:{}".format(hourlySolarIncidenceAngleWestDirection))
    # np.set_printoptions(threshold=1000)

    # estimated horizontal solar irradiances [W m^-2]. these values are used only when estimating solar radiations.
    # symbol: I_DH.
    directHorizontalSolarRadiation = OPVFilm.getDirectHorizontalSolarRadiation(hourlySolarAltitudeAngle, hourlyHorizontalSolarIncidenceAngle)
    # print "directHorizontalSolarRadiation:{}".format(directHorizontalSolarRadiation)
    # set the data. this is necessary when estimating the solar irradiance under multi-span greenhouse with estimated data
    simulatorClass.directHorizontalSolarRadiation = directHorizontalSolarRadiation
    # symbol: I_S
    diffuseHorizontalSolarRadiation = OPVFilm.getDiffuseHorizontalSolarRadiation(hourlySolarAltitudeAngle, hourlyHorizontalSolarIncidenceAngle)
    simulatorClass.diffuseHorizontalSolarRadiation = diffuseHorizontalSolarRadiation
    # print "diffuseHorizontalSolarRadiation:{}".format(diffuseHorizontalSolarRadiation)
    # symbol: I_HT
    totalHorizontalSolarRadiation = directHorizontalSolarRadiation + diffuseHorizontalSolarRadiation
    simulatorClass.totalHorizontalSolarRadiation = totalHorizontalSolarRadiation

    # print "totalHorizontalSolarRadiation:{}".format(totalHorizontalSolarRadiation)

    # tilted surface  solar radiation [W m^-2], real / estimated value branch is calculated in this functions
    # symbol: I_TD (= H_b at Kacira 2004). direct beam radiation on the tilted surface
    # print ("call getDirectTitledSolarRadiation for east direction OPV")
    directTiltedSolarRadiationEastDirection = OPVFilm.getDirectTitledSolarRadiation(simulatorClass, hourlySolarAltitudeAngle, hourlySolarIncidenceAngleEastDirection, \
                                                                                    hourlyHorizontalDirectOuterSolarIrradiance)
    # print ("call getDirectTitledSolarRadiation for west direction OPV")
    directTiltedSolarRadiationWestDirection = OPVFilm.getDirectTitledSolarRadiation(simulatorClass, hourlySolarAltitudeAngle, hourlySolarIncidenceAngleWestDirection, \
                                                                                    hourlyHorizontalDirectOuterSolarIrradiance)

    # print("directTiltedSolarRadiationEastDirection:{}".format(directTiltedSolarRadiationEastDirection))
    # print("directTiltedSolarRadiationWestDirection:{}".format(directTiltedSolarRadiationWestDirection))


    # symbol: I_TS  (= H_d_p at Kacira 2004). diffused radiation on the tilted surface.
    diffuseTiltedSolarRadiation = OPVFilm.getDiffuseTitledSolarRadiation(simulatorClass, hourlySolarAltitudeAngle, diffuseHorizontalSolarRadiation, \
                                                                         hourlyHorizontalDiffuseOuterSolarIrradiance)
    # print "diffuseTiltedSolarRadiation:{}".format(diffuseTiltedSolarRadiation)

    # symbol: I_Trho (= H_gr at Kacira 2004) (albedo radiation = reflectance from the ground)
    albedoTiltedSolarRadiation = OPVFilm.getAlbedoTitledSolarRadiation(simulatorClass, hourlySolarAltitudeAngle, totalHorizontalSolarRadiation, \
                                                                       hourlyHorizontalDirectOuterSolarIrradiance+hourlyHorizontalDiffuseOuterSolarIrradiance)

    return directTiltedSolarRadiationEastDirection, directTiltedSolarRadiationWestDirection, diffuseTiltedSolarRadiation, albedoTiltedSolarRadiation

    # ####################################################################################################
    # # Stop execution here...
    # sys.exit()
    # # Move the above line to different parts of the assignment as you implement more of the functionality.
    # ####################################################################################################

def getDailyElectricityYieldperArea(simulatorClass,hourlyOPVTemperature, directSolarRadiationToOPV, diffuseSolarRadiationToOPV,groundReflectedSolarradiationToOPV):
    '''
    calculate the daily electricity yield per area (m^2).
    :param hourlyOPVTemperature: [celsius]
    :param directSolarRadiationToOPV: [W/m^2]
    :param diffuseSolarRadiationToOPV: [W/m^2]
    :param groundReflectedSolarradiationToOPV:[W/m^2]
    :return:
    '''
    # print "total solar irradiance:{}".format(directSolarRadiationToOPV+diffuseSolarRadiationToOPV+groundReflectedSolarradiationToOPV)

    # [W/m^2] == [J/s/m^2] -> [J/m^2] per day
    dailyJopvoutperArea = OPVFilm.calcOPVElectricEnergyperArea(simulatorClass, hourlyOPVTemperature, directSolarRadiationToOPV+diffuseSolarRadiationToOPV+groundReflectedSolarradiationToOPV)
    # print "dailyJopvout:{}".format(dailyJopvout)

    return dailyJopvoutperArea


def setDirectSolarIrradianceThroughMultiSpanRoof(simulatorClass):
    '''
    calculate the solar irradiance to multi-span roof.
    this calculates the solar irradiance to the single span.
    the variables names follow the symbols in the reference.
    Reference of the model: T. Soriano. et al, 2004, "A Study of Direct Solar Radiation Transmission in Asymmetrical Multi-span Greenhouses using Scale Models and Simulation Models"
    source: https://www.sciencedirect.com/science/article/pii/S1537511004000455

    '''
    # get the direct solar radiation [W/m^2]. these values are not directly used to calculate the transmittance but just used to check the existance of solar irradicne at each hour.
    directSolarRadiationToOPVEastFacingRoof = simulatorClass.getDirectSolarRadiationToOPVEastDirection()
    directSolarRadiationToOPVWestFacingRoof = simulatorClass.getDirectSolarRadiationToOPVWestDirection()
    # print("directSolarRadiationToOPVEastFacingRoof: {}".format(directSolarRadiationToOPVEastFacingRoof))
    # print("directSolarRadiationToOPVWestFacingRoof: {}".format(directSolarRadiationToOPVWestFacingRoof))

    # module azimuth of each roof facing the opposite direction [rad], which is a scalar value
    hourlyModuleAzimuthAngleEast = simulatorClass.hourlyModuleAzimuthAngleEast
    hourlyModuleAzimuthAngleWest = simulatorClass.hourlyModuleAzimuthAngleWest
    # print("hourlyModuleAzimuthAngleEast: {}".format(hourlyModuleAzimuthAngleEast))
    # print("hourlyModuleAzimuthAngleWest: {}".format(hourlyModuleAzimuthAngleWest))

    # angle between the incident ray and the horizontal axis perpendicular to the greenhouse span. This angle is symbolized with E in the reference paper [rad]
    EPerpendicularEastOrNorthFacingRoof = SolarIrradianceMultiSpanRoof.getAngleBetweenIncientRayAndHorizontalAxisPerpendicularToGHSpan(simulatorClass, hourlyModuleAzimuthAngleEast)
    EPerpendicularWestOrSouthFacingRoof = SolarIrradianceMultiSpanRoof.getAngleBetweenIncientRayAndHorizontalAxisPerpendicularToGHSpan(simulatorClass, hourlyModuleAzimuthAngleWest)

    # np.set_printoptions(threshold=np.inf)
    # print("EPerpendicularEastOrNorthFacingRoof: {}".format(EPerpendicularEastOrNorthFacingRoof))
    # print("EPerpendicularWestOrSouthFacingRoof: {}".format(EPerpendicularWestOrSouthFacingRoof))
    # np.set_printoptions(threshold=1000)

    # # Referring to Soriano. et al, (2004), it was found that we can get the direct solar irradiance to horizontal surface inside multi-span greenhouse just by
    # # multiplying the outer solar irradiance to horizontal surface with
    # # angle between the incident ray and the horizontal axis perpendicular to the greenhouse span. This angle is nit symbolized in the reference paper.
    # # the following angles should be same in our case, but both were separately calculated for program expandability
    # EParallelEastOrNorthFacingRoof = SolarIrradianceMultiSpanRoof.getAngleBetweenIncientRayAndHorizontalAxisParallelToGHSpan(simulatorClass, hourlyModuleAzimuthAngleEast)
    # EParallelWestOrSouthFacingRoof = SolarIrradianceMultiSpanRoof.getAngleBetweenIncientRayAndHorizontalAxisParallelToGHSpan(simulatorClass, hourlyModuleAzimuthAngleWest)
    # # print("EParallelEastOrNorthFacingRoof: {}".format(EParallelEastOrNorthFacingRoof))
    # # print("EParallelWestOrSouthFacingRoof: {}".format(EParallelWestOrSouthFacingRoof))
    #
    # # get the T_mat for parallel irradiance
    # T_matForParallelIrrEastOrNorthFacingRoof = SolarIrradianceMultiSpanRoof.getTransmittanceForParallelIrrThroughMultiSpanRoof(simulatorClass,EParallelEastOrNorthFacingRoof)
    # T_matForParallelIrrWestOrSouthFacingRoof = SolarIrradianceMultiSpanRoof.getTransmittanceForParallelIrrThroughMultiSpanRoof(simulatorClass,EParallelWestOrSouthFacingRoof)
    # print("T_matForParallelIrrEastOrNorthFacingRoof: {}".format(T_matForParallelIrrEastOrNorthFacingRoof))
    # print("T_matForParallelIrrWestOrSouthFacingRoof: {}".format(T_matForParallelIrrWestOrSouthFacingRoof))

    # get the incidence angles
    hourlySolarIncidenceAngleEastDirection = simulatorClass.hourlySolarIncidenceAngleEastDirection
    hourlySolarIncidenceAngleWestDirection = simulatorClass.hourlySolarIncidenceAngleWestDirection

    # get the direct solar irradiance on each axis
    directSolarIrradiancePerpendicularToOPVEastDirection = directSolarRadiationToOPVEastFacingRoof * np.cos(hourlySolarIncidenceAngleEastDirection)
    directSolarIrradianceParallelToOPVEastDirection = directSolarRadiationToOPVEastFacingRoof * np.sin(hourlySolarIncidenceAngleEastDirection)
    directSolarIrradiancePerpendicularToOPVWestDirection = directSolarRadiationToOPVWestFacingRoof * np.cos(hourlySolarIncidenceAngleWestDirection)
    directSolarIrradianceParallelToOPVWestDirection = directSolarRadiationToOPVWestFacingRoof * np.sin(hourlySolarIncidenceAngleWestDirection)

    # np.set_printoptions(threshold=np.inf)
    # print("directSolarIrradiancePerpendicularToOPVEastDirection: {}".format(directSolarIrradiancePerpendicularToOPVEastDirection))
    # print("directSolarIrradianceParallelToOPVEastDirection: {}".format(directSolarIrradianceParallelToOPVEastDirection))
    # print("directSolarIrradiancePerpendicularToOPVWestDirection: {}".format(directSolarIrradiancePerpendicularToOPVWestDirection))
    # print("directSolarIrradianceParallelToOPVWestDirection: {}".format(directSolarIrradianceParallelToOPVWestDirection))
    # np.set_printoptions(threshold=1000)

    # the the T_mat for parpendicular irradiance
    # print("getTransmittanceForPerpendicularIrrThroughMultiSpanRoofFacingEastOrNorth start")

    # to avoide the error "RuntimeError: maximum recursion depth exceeded", the maximum recursion limitation is increased.
    # sys.setrecursionlimit(constant.mMax)
    # print("sys.getrecursionlimit():{}".format(sys.getrecursionlimit()))
    T_matForPerpendicularIrrEastOrNorthFacingRoof = SolarIrradianceMultiSpanRoof.getTransmittanceForPerpendicularIrrThroughMultiSpanRoofFacingEastOrNorth(\
                                                                simulatorClass, directSolarIrradiancePerpendicularToOPVEastDirection, EPerpendicularEastOrNorthFacingRoof)
    # print("getTransmittanceForPerpendicularIrrThroughMultiSpanRoofFacingWestOrSouth start: {}")
    T_matForPerpendicularIrrWestOrSouthFacingRoof = SolarIrradianceMultiSpanRoof.getTransmittanceForPerpendicularIrrThroughMultiSpanRoofFacingWestOrSouth(\
                                                                simulatorClass, directSolarIrradiancePerpendicularToOPVWestDirection, EPerpendicularWestOrSouthFacingRoof)
    # roll back the recursive limitation setting. The default number should be changed according to each local env.
    # sys.setrecursionlimit(constant.defaultIterationLimit)

    # set the data
    simulatorClass.T_matForPerpendicularIrrEastOrNorthFacingRoof = T_matForPerpendicularIrrEastOrNorthFacingRoof
    simulatorClass.T_matForPerpendicularIrrWestOrSouthFacingRoof = T_matForPerpendicularIrrWestOrSouthFacingRoof

    # the overall transmittance of multispanroof. The solar irradiance inside the greenhouse can be derived only by multiplying this with the outer solar irradiance for horizontal surface
    integratedT_mat = SolarIrradianceMultiSpanRoof.getIntegratedT_matFromBothRoofs(T_matForPerpendicularIrrEastOrNorthFacingRoof, T_matForPerpendicularIrrWestOrSouthFacingRoof)
    # set the data
    simulatorClass.integratedT_mat = integratedT_mat

    # np.set_printoptions(threshold=np.inf)
    # print("T_matForPerpendicularIrrEastOrNorthFacingRoof: {}".format(T_matForPerpendicularIrrEastOrNorthFacingRoof))
    # print("T_matForPerpendicularIrrWestOrSouthFacingRoof: {}".format(T_matForPerpendicularIrrWestOrSouthFacingRoof))
    # print("integratedT_mat:{}".format(integratedT_mat))
    # np.set_printoptions(threshold=1000)

    # get the solar irradiance inside
    if constant.ifUseOnlyRealData == True:
      hourlyDirectSolarRadiationAfterMultiSpanRoof = integratedT_mat * simulatorClass.hourlyHorizontalDirectOuterSolarIrradiance
    # this uses the estimated direct solar irradiance (without real data)
    else:
      hourlyDirectSolarRadiationAfterMultiSpanRoof = integratedT_mat * simulatorClass.directHorizontalSolarRadiation

    # set the solar irradiance [W/m^2]
    simulatorClass.setHourlyDirectSolarRadiationAfterMultiSpanRoof(hourlyDirectSolarRadiationAfterMultiSpanRoof)

    # unit change of the imported outer solar radiation: [W m^-2] -> [umol m^-2 s^-1] == PPFD
    hourlyDirectPPFDTAfterMultiSpanRoof = Util.convertFromWattperSecSquareMeterToPPFD(hourlyDirectSolarRadiationAfterMultiSpanRoof)
    # set the solar irradiance [umol m^-2 s^-1] == PPFD
    simulatorClass.setHourlyDirectPPFDAfterMultiSpanRoof(hourlyDirectPPFDTAfterMultiSpanRoof)

    # convert the unit into PPFD snd DLI
    directDLIAfterMultiSpanRoof = Util.convertFromHourlyPPFDWholeDayToDLI(hourlyDirectPPFDTAfterMultiSpanRoof)
    simulatorClass.setHourlyDirectPPFDAfterMultiSpanRoof(directDLIAfterMultiSpanRoof)


def setSolarIrradianceToPlants(simulatorClass):
    '''
    # calculate the light intensity to plants after penetrating the roof, considering the sidewall material transmittance, shading curtain, and the greenhouse structure shading
    it was assumed ground reflectance does not significantly affect the solar irradiance to plants
    '''

    directSolarIrradianceBeforeShadingCurtain = OPVFilm.getDirectSolarIrradianceBeforeShadingCurtain(simulatorClass)
    # set the data to the object
    simulatorClass.directSolarIrradianceBeforeShadingCurtain = directSolarIrradianceBeforeShadingCurtain
    diffuseSolarIrradianceBeforeShadingCurtain = OPVFilm.getDiffuseSolarIrradianceBeforeShadingCurtain(simulatorClass)
    # set the data to the object
    simulatorClass.diffuseSolarIrradianceBeforeShadingCurtain = diffuseSolarIrradianceBeforeShadingCurtain

    # get the shading curtain transmittance
    ShadingCurtain.getHourlyShadingCurtainDeploymentPatternChangingEachMonthMain(simulatorClass)
    # #############command to print out all array data
    # np.set_printoptions(threshold=np.inf)
    # print("simulatorClass.transmittanceThroughShadingCurtainChangingEachMonth:{}".format(simulatorClass.transmittanceThroughShadingCurtainChangingEachMonth))
    # np.set_printoptions(threshold=1000)
    # #############

    # calculate the light intensity to plants [W m-2]
    directSolarIrradianceToPlants = OPVFilm.getDirectSolarIrradianceToPlants(simulatorClass, directSolarIrradianceBeforeShadingCurtain)
    diffuseSolarIrradianceToPlants = OPVFilm.getDiffuseSolarIrradianceToPlants(simulatorClass, diffuseSolarIrradianceBeforeShadingCurtain)

    # # change this part if you can to see how plant fresh weight changes with the change of solar irradiance to plants
    # directSolarIrradianceToPlants = directSolarIrradianceToPlants * 2.0
    # diffuseSolarIrradianceToPlants = diffuseSolarIrradianceToPlants * 2.0

    # set the data to the object
    simulatorClass.directSolarIrradianceToPlants = directSolarIrradianceToPlants
    simulatorClass.diffuseSolarIrradianceToPlants = diffuseSolarIrradianceToPlants

    # #############command to print out all array data
    # np.set_printoptions(threshold=np.inf)
    # print("directSolarIrradianceToPlants:{}".format(directSolarIrradianceToPlants))
    # print("diffuseSolarIrradianceToPlants:{}".format(diffuseSolarIrradianceToPlants))
    # np.set_printoptions(threshold=1000)
    # #############

    # unit change of the imported outer solar radiation: [W m^-2] -> [umol m^-2 s^-1] == PPFD
    directPPFDToPlants = Util.convertFromWattperSecSquareMeterToPPFD(directSolarIrradianceToPlants)
    diffusePPFDToPlants = Util.convertFromWattperSecSquareMeterToPPFD(diffuseSolarIrradianceToPlants)
    # set the solar irradiance [umol m^-2 s^-1] == PPFD
    simulatorClass.directPPFDToPlants = directPPFDToPlants
    simulatorClass.diffusePPFDToPlants = diffusePPFDToPlants

    # convert the unit into PPFD snd DLI
    directDLIToPlants = Util.convertFromHourlyPPFDWholeDayToDLI(directPPFDToPlants)
    diffuseDLIToPlants = Util.convertFromHourlyPPFDWholeDayToDLI(diffusePPFDToPlants)
    simulatorClass.directDLIToPlants = directDLIToPlants
    simulatorClass.diffuseDLIToPlants = diffuseDLIToPlants
    # #############command to print out all array data
    # np.set_printoptions(threshold=np.inf)
    # print("directDLIToPlants:{}".format(directDLIToPlants))
    # print("diffuseDLIToPlants:{}".format(diffuseDLIToPlants))
    # np.set_printoptions(threshold=1000)
    # #############

# def setThermalTimeToPlants(simulatorClass):
#   '''
#       calc/set the thermal time to the object: average Celsius temperature per day * days [Celsius d]
#       On the model, since it was assumed the temperature in the greenhouse is maintained at the set point by cooling system (pad and fan system), this function is not used.
#   '''
#   importedHourlyAirTemperature = simulatorClass.getImportedHourlyAirTemperature()
#
#   # TODO assume the greenhouse temperature from the outer air temperature
#   airTemperatureInGreenhouse = importedHourlyAirTemperature
#   simulatorClass.setThermalTimeToPlants(airTemperatureInGreenhouse)


def getPlantYieldSimulation(simulatorClass):
    '''
    calculate the daily plant yield

    :param cultivationDaysperHarvest: [days / harvest]
    :param OPVAreaCoverageRatio: [-] range(0-1)
    :param directPPFDToOPV: hourly average [umol m^-2 s^-1] == PPFD
    :param diffusePPFDToOPV: hourly average [umol m^-2 s^-1] == PPFD
    :param groundReflectedPPFDToOPV: hourly average [umol m^-2 s^-1] == PPFD
    :param hasShadingCurtain: Boolean
    :param ShadingCurtainDeployPPFD: float [umol m^-2 s^-1] == PPFD
    :param cropElectricityYieldSimulator1: object
    :return:
    '''
    plantGrowthModel = simulatorClass.getPlantGrowthModel()

    # get cultivation days per harvest. this may not be used in some plant growth models
    cultivationDaysperHarvest = simulatorClass.getCultivationDaysperHarvest()
    # OPVAreaCoverageRatio = simulatorClass.getOPVAreaCoverageRatio()
    # hasShadingCurtain = simulatorClass.getIfHasShadingCurtain()
    # ShadingCurtainDeployPPFD = simulatorClass.getShadingCurtainDeployPPFD()

    # # This unit conversion was done at getSolarIrradianceToPlants
    # # calculate the light intensity to plants
    # # hourly average PPFD [umol m^-2 s^-1]
    # hourlyInnerPPFDToPlants = OPVFilm.calcHourlyInnerLightIntensityPPFD(directPPFDToMultiSpanRoof + diffusePPFDToMultiSpanRoof + groundReflectedPPFDToMultiSpanRoof, \
    #     OPVAreaCoverageRatio, constant.OPVPARTransmittance, hasShadingCurtain,ShadingCurtainDeployPPFD, simulatorClass)

    # np.set_printoptions(threshold=np.inf)
    # print "OPVAreaCoverageRatio:{}, directPPFDToOPV+diffusePPFDToOPV+groundReflectedPPFDToOPV:{}".format(OPVAreaCoverageRatio, directPPFDToOPV+diffusePPFDToOPV+groundReflectedPPFDToOPV)
    # np.set_printoptions(threshold=1000)

    # calculate the daily increase of unit fresh weight
    # this model considers only solar irradiance, and so this will not be so practical
    # the simulated cultivar is butter head lettuce
    if plantGrowthModel == constant.A_J_Both_Modified_TaylorExpantionWithFluctuatingDLI:
        #unit  [g/head]
        shootFreshMassList, unitDailyFreshWeightIncrease, accumulatedUnitDailyFreshWeightIncrease, unitHarvestedFreshWeight = \
            Lettuce.calcUnitDailyFreshWeightBoth2003TaylorExpantionWithVaryingDLI(simulatorClass.directPPFDToPlants + simulatorClass.diffusePPFDToPlants, cultivationDaysperHarvest, simulatorClass)
        # print "shootFreshMassList.shape:{}".format(shootFreshMassList.shape)

    # the simulated cultivar is Berlo and Norden
    elif plantGrowthModel == constant.E_J_VanHenten1994:
        # unit  [g/head]
        shootFreshMassList, \
        unitDailyFreshWeightIncrease, \
        accumulatedUnitDailyFreshWeightIncrease, \
        unitHarvestedFreshWeight = \
            PlantGrowthModelE_J_VanHenten.calcUnitDailyFreshWeightE_J_VanHenten1994(simulatorClass)
        # print("shootFreshMassList.shape[0]:{}".format(shootFreshMassList.shape[0]))
        # print("unitDailyFreshWeightIncrease.shape[0]:{}".format(unitDailyFreshWeightIncrease.shape[0]))
        # print("accumulatedUnitDailyFreshWeightIncrease.shape[0]:{}".format(accumulatedUnitDailyFreshWeightIncrease.shape[0]))
        # print("unitHarvestedFreshWeight.shape[0]:{}".format(unitHarvestedFreshWeight.shape[0]))

        # set the data
        simulatorClass.shootFreshMassList = shootFreshMassList

        # Be careful! this model returns hourly weight, not daily weight. so convert the hourly value into daily value.
        dailyShootFreshMassList = shootFreshMassList[23::constant.hourperDay]
        # print("dailyShootFreshMassList:{}".format(dailyShootFreshMassList))

        # dailyUnitDailyFreshWeightIncrease = np.array(sum[ unitDailyFreshWeightIncrease[constant.hourperDay*(i-1):constant.hourperDay*i]] \
        #                                              for i in range (0, unitDailyFreshWeightIncrease.shape[0]/constant.hourperDay ))
        dailyUnitDailyFreshWeightIncrease = Lettuce.getFreshWeightIncrease(dailyShootFreshMassList)
        dailyAccumulatedUnitDailyFreshWeightIncrease = Lettuce.getAccumulatedFreshWeightIncrease(dailyShootFreshMassList)
        dailyUnitHarvestedFreshWeight = Lettuce.getHarvestedFreshWeight(dailyShootFreshMassList)
        # print("dailyUnitDailyFreshWeightIncrease.shape:{}".format(dailyUnitDailyFreshWeightIncrease.shape))
        # print("dailyAccumulatedUnitDailyFreshWeightIncrease.shape:{}".format(dailyAccumulatedUnitDailyFreshWeightIncrease.shape))
        # print("dailyUnitHarvestedFreshWeight.shape:{}".format(dailyUnitHarvestedFreshWeight.shape))
        # print("dailyUnitHarvestedFreshWeight:{}".format(dailyUnitHarvestedFreshWeight))

    # this model was coded, but the result was not better than constant.E_J_VanHenten1994
    elif plantGrowthModel == constant.S_Pearson1997:
        # unit  [g/head]
        dailyShootFreshMassList, \
        dailyUnitDailyFreshWeightIncrease, \
        dailyAccumulatedUnitDailyFreshWeightIncrease, \
        dailyUnitHarvestedFreshWeight = \
          PlantGrowthModelS_Pearson1997.calcUnitDailyFreshWeightS_Pearson1997(simulatorClass)

    else:
      print ("no valid model name is assigned. Stop the simulation. Please choose a registered one")
      ####################################################################################################
      # Stop execution here...
      sys.exit()
      # Move the above line to different parts of the assignment as you implement more of the functionality.
      ####################################################################################################

    # set the values to the object
    simulatorClass.dailyShootFreshMass = dailyShootFreshMassList
    simulatorClass.dailyUnitDailyFreshWeightIncrease = dailyUnitDailyFreshWeightIncrease
    simulatorClass.dailyAccumulatedUnitDailyFreshWeightIncrease = dailyAccumulatedUnitDailyFreshWeightIncrease
    simulatorClass.dailyUnitHarvestedFreshWeight = dailyUnitHarvestedFreshWeight

    return dailyShootFreshMassList, dailyUnitDailyFreshWeightIncrease, dailyAccumulatedUnitDailyFreshWeightIncrease, dailyUnitHarvestedFreshWeight


def getTotalDLIToPlants(OPVAreaCoverageRatio, directPPFDToOPV, diffusePPFDToOPV, groundReflectedPPFDToOPV, hasShadingCurtain, ShadingCurtainDeployPPFD, \
                        cropElectricityYieldSimulator1):
    '''
    the daily light integral to plants for the given simulation period.
    :param OPVAreaCoverageRatio:
    :param directPPFDToOPV:
    :param diffusePPFDToOPV:
    :param groundReflectedPPFDToOPV:
    :param hasShadingCurtain:
    :param ShadingCurtainDeployPPFD:
    :param cropElectricityYieldSimulator1: instance
    :return:
    '''

    # calculate the light intensity to plants
    # hourly average PPFD [umol m^-2 s^-1]
    hourlyInnerPPFDToPlants = OPVFilm.calcHourlyInnerLightIntensityPPFD(directPPFDToOPV+diffusePPFDToOPV+groundReflectedPPFDToOPV, \
        OPVAreaCoverageRatio, constant.OPVPARTransmittance, hasShadingCurtain,ShadingCurtainDeployPPFD, cropElectricityYieldSimulator1)

    # convert PPFD to DLI
    innerDLIToPlants = Util.convertFromHourlyPPFDWholeDayToDLI(hourlyInnerPPFDToPlants)
    # print "innerDLIToPlants:{}".format(innerDLIToPlants)

    return innerDLIToPlants


def penalizeDailyHarvestedFreshWeightPerHead(dailyHarvestedFreshWeightPerHead, simulatorClass):
    '''
    the function was made based on the data of plant fresh weights for 400 600, and 800 PPFD (umiol m^-2, s-1) in the source below:
    Table 1 at "Effects of different light intensities on anti-oxidative enzyme activity, quality and biomass in lettuce, Weiguo Fu, Pingping Li, Yanyou Wu, Juanjuan Tang"

    The parameters were derived with the solver of Excel 2007, the process is written in "penalizePlantYieldBySolarRadiation.xlsx"

    Table 1. Quality and biomass of above-ground part of lettuce under different light intensity treatments
    Light intensity (μmol m-2 s-1)        Biomass of above-ground part (g plant-1, FW)
    100                                   127.98 ± 8.32
    200                                   145.65 ± 7.53
    400                                   158.45 ± 6.21
    600                                   162.89 ± 7.13
    800                                   135.56 ± 5.76
    '''

    penalizedUnitDailyHarvestedFreshWeight = np.zeros(dailyHarvestedFreshWeightPerHead.shape[0])

    # the DLI including both direct and diffuse to plants
    totalDLItoPlants = simulatorClass.totalDLItoPlants
    # print("totalDLItoPlants:{}".format(totalDLItoPlants))

    # get the average DLI of each cultivation cycle.
    if  simulatorClass.getPlantGrowthModel()  == constant.A_J_Both_Modified_TaylorExpantionWithFluctuatingDLI:
      averageDLIonEachCycle = simulatorClass.getAverageDLIonEachCycle()
    else:
      averageDLIonEachCycle = np.zeros(simulatorClass.totalDLItoPlants.shape[0])

      nextCultivationStartDay = 0
      for i in range(0, simulatorClass.totalDLItoPlants.shape[0]):

        # if the date is not the harvest date, then skip.
        if dailyHarvestedFreshWeightPerHead[i] == 0.0:
            continue
        # Right row, E_J_VanHenten1994 is assumed
        else:
          # calc the average DLI during each cultivation cycle
          averageDLIonEachCycle[i] = np.mean(totalDLItoPlants[nextCultivationStartDay:i + 1])
          print("i:{}, averageDLIonEachCycle:{}".format(i, averageDLIonEachCycle[i]))
          # update lastHarvestDay
          # It was assumed to take 3 days to the next cultivation cycle assuming "transplanting  shock  prevented growth during the first 48 h", and it takes one day for preparation.
          nextCultivationStartDay = i + 3

      # print("averageDLIonEachCycle:{}".format(averageDLIonEachCycle))

    # parameters, which is derived from the data of the reference
    photoPriod = {"hour":14.0}
    optimumLightIntensityDLI = {"mol m-2 d-1": 26.61516313}
    # maximumYieldFW = {"g unit-1": 164.9777479}
    maximumYieldFW = getPenalizedUnitFreshWeight(optimumLightIntensityDLI["mol m-2 d-1"])
    # print("maximumYieldFW:{}".format(maximumYieldFW))
    # convert PPFD to DLI
    # optimumLightIntensityPPFD = {"umol m-2 s-1": 524.1249999}
    # optimumLightIntensityDLI = {"mol m-2 d-1": optimumLightIntensityPPFD["umol m-2 s-1"] * constant.secondperMinute * constant.minuteperHour * photoPriod["hour"] / 1000000.0}

    i = 0
    while i <  dailyHarvestedFreshWeightPerHead.shape[0]:

        # if the date is not the harvest date, then skip.
        if dailyHarvestedFreshWeightPerHead[i] == 0.0:
            i += 1
            continue
        else:

            print ("averageDLIonEachCycle:{}".format(averageDLIonEachCycle[i]))
            print ("dailyHarvestedFreshWeightPerHead[i]:{}".format(dailyHarvestedFreshWeightPerHead[i]))
            print("getPenalizedUnitFreshWeight(averageDLIonEachCycle[i]):{}, i:{}".format(getPenalizedUnitFreshWeight(averageDLIonEachCycle[i]), i))

            if averageDLIonEachCycle[i] > optimumLightIntensityDLI["mol m-2 d-1"] and getPenalizedUnitFreshWeight(averageDLIonEachCycle[i]) > 0.0:
                # penalize the plant fresh weight
                print ("penaize the fresh weight, i:{}".format(i))
                penalizedUnitDailyHarvestedFreshWeight[i] = dailyHarvestedFreshWeightPerHead[i] - dailyHarvestedFreshWeightPerHead[i] / maximumYieldFW["g unit-1"] * (maximumYieldFW["g unit-1"] - getPenalizedUnitFreshWeight(averageDLIonEachCycle[i]))

                print("penalizedUnitDailyHarvestedFreshWeight[i]:{}".format(penalizedUnitDailyHarvestedFreshWeight[i]))
                print("unitDailyHarvestedFreshWeight[i]:{}".format(dailyHarvestedFreshWeightPerHead[i]))

            # if the penalty is too strong and the weight becomes zero
            elif averageDLIonEachCycle[i] > optimumLightIntensityDLI["mol m-2 d-1"] and getPenalizedUnitFreshWeight(averageDLIonEachCycle[i]) <= 0.0:
                print ("the light intensity may be too strong. The yield was penalized to zero")
                penalizedUnitDailyHarvestedFreshWeight[i] = 0.0

            # if no penalization occured
            else:
                penalizedUnitDailyHarvestedFreshWeight[i] = dailyHarvestedFreshWeightPerHead[i]
        i += 1

    return penalizedUnitDailyHarvestedFreshWeight


def getPenalizedUnitFreshWeight(lightIntensityDLI):
  '''
  The following parameters were derived from the soruce mentioned at penalizeDailyHarvestedFreshWeightPerHead

  '''
  a = -0.1563
  b = 8.3199
  c = 54.26
  return a * lightIntensityDLI**2 + b * lightIntensityDLI + c


def getWholeElectricityYieldEachOPVRatio(OPVAreaCoverageRatio, dailyJopvout, cropElectricityYieldSimulator1, greenhouseRoofArea = None):
    '''
    return the total electricity yield for a given period by the given OPV area(OPVAreaCoverageRatio * constant.greenhouseRoofArea)
    :param OPVAreaCoverageRatio: [-] proportionOPVAreaCoverageRatio
    :param dailyJopvout: [J/m^2] per day
    :return: dailyJopvout [J/m^2] by whole OPV area
    '''

    # get the OPV coverage ratio changing during the fallow period
    unfixedOPVCoverageRatio = OPVFilm.getDifferentOPVCoverageRatioInSummerPeriod(OPVAreaCoverageRatio, cropElectricityYieldSimulator1)
    # change the num of list from hourly data (365 * 24) to daily data (365)
    unfixedOPVCoverageRatio = unfixedOPVCoverageRatio[::24]

    if greenhouseRoofArea is None:
      return sum(dailyJopvout * unfixedOPVCoverageRatio * constant.greenhouseRoofArea)
    else:
      return sum(dailyJopvout * unfixedOPVCoverageRatio * greenhouseRoofArea)
    # # print "dailyJopvout:{}".format(dailyJopvout)
    # totalJopvout = sum(dailyJopvout)
    # if greenhouseRoofArea is None:
    #     return totalJopvout * unfixedOPVCoverageRatio * constant.greenhouseRoofArea
    # else:
    #     return totalJopvout * unfixedOPVCoverageRatio * greenhouseRoofArea


def getMonthlyElectricitySalesperArea(dailyJopvoutperArea, yearOfeachDay, monthOfeachDay, simulatorClass):
    '''

    :param dailyJopvoutperArea:
    :param yearOfeachDay:
    :param monthOfeachDay:
    :return:
    '''
    # unit: J/m^2/month
    monthlyElectricityYieldperArea = OPVFilm.getMonthlyElectricityProductionFromDailyData(dailyJopvoutperArea, yearOfeachDay, monthOfeachDay)
    # print("monthlyElectricityYieldperArea:{}".format(monthlyElectricityYieldperArea))

    # import the electricity sales price file: source (download the CSV file): https://www.eia.gov/electricity/data/browser/#/topic/7?agg=0,1&geo=0000000001&endsec=vg&freq=M&start=200101&end=201802&ctype=linechart&ltype=pin&rtype=s&maptype=0&rse=0&pin=
    fileName = constant.averageRetailPriceOfElectricityMonthly
    # import the file removing the header
    fileData = Util.readData(fileName, relativePath="", skip_header=1, d='\t')
    # print ("fileData:{}".format(fileData))

    simulatorClass.monthlyElectricityRetailPrice = fileData

    # print "monthlyElectricityYieldperArea.shape[0]:{}".format(monthlyElectricityYieldperArea.shape[0])
    # year = np.zeros(monthlyElectricityYieldperArea.shape[0])
    # month = np.zeros(monthlyElectricityYieldperArea.shape[0])
    monthlyResidentialElectricityPrice = np.zeros(monthlyElectricityYieldperArea.shape[0])

    index = 0
    for monthlyData in fileData:
        # exclude the data out of the set start month and end month
        # print("monthlyData:{}".format(monthlyData))
        if datetime.date(int(monthlyData[1]), int(monthlyData[0]), 1) + relativedelta(months=1) <= Util.getStartDateDateType() or \
                datetime.date(int(monthlyData[1]), int(monthlyData[0]), 1) > Util.getEndDateDateType():
            continue

        # year[index] = monthlyData[1]
        # month[index] = monthlyData[0]
        # take the residential electricity retail price
        monthlyResidentialElectricityPrice[index] = monthlyData[2]

        # print "monthlyData:{}".format(monthlyData)
        index += 1

    # print("monthlyResidentialElectricityPrice[Cents/kwh]:{}".format(monthlyResidentialElectricityPrice))

    # unit exchange: [J/m^2] -> [wh/m^2]
    monthlyWhopvoutperArea =Util.convertFromJouleToWattHour(monthlyElectricityYieldperArea)
    # unit exchange: [wh/m^2] -> [kwh/m^2]
    monthlyKWhopvoutperArea =Util.convertWhTokWh(monthlyWhopvoutperArea)
    # print("monthlyKWhopvoutperArea[kwh/m^2]:{}".format(monthlyKWhopvoutperArea))
    # [USD/month/m^2]
    monthlyElectricitySalesperArea = OPVFilm.getMonthlyElectricitySalesperArea(monthlyKWhopvoutperArea, monthlyResidentialElectricityPrice)
    # print "monthlyElectricitySalesperArea:{}".format(monthlyElectricitySalesperArea)

    return monthlyElectricitySalesperArea


def getMonthlyElectricitySales(OPVCoverage, monthlyElectricitySalesperArea, greenhouseRoofArea = None):
    '''
    return the monthly electricity saled given a cetain OPV coverage ratio

    :param OPVCoverageList:
    :param monthlyElectricitySalesperArea:
    :return:
    '''
    if greenhouseRoofArea is None:
        return monthlyElectricitySalesperArea * OPVCoverage * constant.greenhouseRoofArea
    else:
        return monthlyElectricitySalesperArea * OPVCoverage * greenhouseRoofArea

def getElectricitySalesperAreaEachOPVRatio():
    return 0

def getElectricityCostperArea():
    return 0


def getPlantSalesperSquareMeter(simulatorClass):
    """
    calculate the sales price of lettuce per square meter.
    The referred price is Lettuce, romaine, per lb. (453.6 gm) in U.S ( Northeast region: Connecticut, Maine, Massachusetts, New Hampshire, New Jersey, New York, Pennsylvania, Rhode Island, and Vermont.), city average, average price, not seasonally adjusted
    reference URL: https://data.bls.gov/timeseries/APU0000FL2101?amp%253bdata_tool=XGtable&output_view=data&include_graphs=true
    """

    # get the following data from the object
    totalDLIToPlants = simulatorClass.totalDLItoPlants

    #################### this conversion is not used any more ####################
    # # the price of lettuce per hundredweight [cwt]
    # priceperCwtEachHour = Lettuce.getLettucePricepercwt(year)
    # # unit conversion: cwt -> kg
    # priceperKgEachHour = priceperCwtEachHour / constant.kgpercwt * constant.plantPriceDiscountRatio_justForSimulation
    # # print "harvestedFreshWeightListperAreaKg:{}".format(harvestedFreshWeightListperAreaKg)
    # # print "dailyHarvestedFreshWeightListperAreaKg.shape:{}".format(dailyHarvestedFreshWeightListperAreaKg.shape)
    # # print "priceperKg:{}".format(priceperKg)
    # # convert the price each hour to the price each day
    # priceperKgEachDay = priceperKgEachHour[::24]
    # # print "priceperKgEachDay:{}".format(priceperKgEachDay)
    # # print "priceperKgEachDay.shape:{}".format(priceperKgEachDay.shape)
    #################################################################################0
    # get the retail price of lettuce harvested at each cycle
    # unit: USD/m^2/day
    plantSalesPerSquareMeter = Lettuce.getRetailPricePerArea(simulatorClass)
    # print ("plantSalesPerSquareMeter:{}".format(plantSalesPerSquareMeter))

    if constant.IfConsiderDiscountByTipburn  == True:
      # Tipburn discount
      # TODO: need to refine more
      plantSalesperSquareMeter = Lettuce.discountPlantSalesperSquareMeterByTipburn(plantSalesPerSquareMeter, totalDLIToPlants)

    return plantSalesPerSquareMeter

def getGreenhouseOperationCostForGrowingPlants(simulatorClass):
  '''
  This function estimates the economic cost for cooling and heating a greenhouse by simulating the energy balance model of a greenhouse.
  '''

  # get environment data to calculate the energy for cooling and heating

  # the energy for cooling and heating
  energyBalance.getGHEnergyConsumptionByCoolingHeating(simulatorClass)
  # unit: W
  Q_vW = simulatorClass.Q_vW
  # ############command to print out all array data
  # np.set_printoptions(threshold=np.inf)
  # print("Q_vW:{}".format(Q_vW))
  # np.set_printoptions(threshold=1000)
  # ############

  # if the energy balance is minus, we need to heat to maintain the internal temperature. [W m]
  requiredHeatingEnergyForPlants = np.array([-Q_vW["coolingOrHeatingEnergy W"][i] if Q_vW["coolingOrHeatingEnergy W"][i] < 0.0 else 0.0 for i in range (Q_vW["coolingOrHeatingEnergy W"].shape[0])])
  # if the energy balance is plus, we need to cool to maintain the internal temperature. [W m]
  requiredCoolingEnergyForPlants = np.array([Q_vW["coolingOrHeatingEnergy W"][i] if Q_vW["coolingOrHeatingEnergy W"][i] > 0.0 else 0.0 for i in range (Q_vW["coolingOrHeatingEnergy W"].shape[0])])


  # ############command to print out all array data
  # np.set_printoptions(threshold=np.inf)
  # print("requiredCoolingEnergyForPlants:{}".format(requiredCoolingEnergyForPlants))
  # np.set_printoptions(threshold=1000)
  # ############

  # unit: USD
  totalHeatingCostForPlants = energyBalance.getGHHeatingEnergyCostForPlants(requiredHeatingEnergyForPlants, simulatorClass)
  totalCoolingCostForPlants = energyBalance.getGHCoolingEnergyCostForPlants(requiredCoolingEnergyForPlants, simulatorClass)
  simulatorClass.totalHeatingCostForPlants = totalHeatingCostForPlants
  simulatorClass.totalCoolingCostForPlants = totalCoolingCostForPlants

  # unit: USD m-2
  totalHeatingCostForPlantsPerGHFloorArea = totalHeatingCostForPlants / constant.greenhouseFloorArea
  totalCoolingCostForPlantsPerGHFloorArea = totalCoolingCostForPlants / constant.greenhouseFloorArea
  simulatorClass.totalHeatingCostForPlantsPerGHFloorArea = totalHeatingCostForPlantsPerGHFloorArea
  simulatorClass.totalCoolingCostForPlantsPerGHFloorArea = totalCoolingCostForPlantsPerGHFloorArea


  return totalHeatingCostForPlants, totalCoolingCostForPlants


def getLaborCost(simulatorClass):
    """
    get the total labor cost during the simulation period
    :return:
    """

    harvestedShootFreshMassPerAreaKgPerDay = simulatorClass.harvestedShootFreshMassPerAreaKgPerDay
    # unit:kg
    totalHarvestedShootFreshMass = sum(harvestedShootFreshMassPerAreaKgPerDay) * constant.greenhouseCultivationFloorArea
    # print("totalHarvestedShootFreshMass:{}".format(totalHarvestedShootFreshMass))

    # source: https://onlinelibrary.wiley.com/doi/abs/10.1111/cjag.12161
    # unit: [labors/10000 kg yield]
    necessaryLaborPer10000kgYield = constant.necessaryLaborPer10000kgYield

    # source:https://www.bls.gov/regions/west/news-release/occupationalemploymentandwages_tucson.htm
    # unit:USD/labor/hour
    hourlyWagePerPerson = constant.hourlyWagePerPerson

    # unit:hour/day
    workingHourPerDay = constant.workingHourPerDay

    totalLaborCost = (totalHarvestedShootFreshMass / 10000.0) * necessaryLaborPer10000kgYield * workingHourPerDay * hourlyWagePerPerson * Util.getSimulationDaysInt()
    # print("totalLaborCost:{}".format(totalLaborCost))

    return totalLaborCost


def getPlantCostperSquareMeter(simulationDays):
    '''
    calculate the cost for plant cultivation for given period
    :param year:
    :return:
    '''
    # [USD/m^2]
    return constant.plantcostperSquaremeterperYear * simulationDays / constant.dayperYear






################################################# old code below################################

def calcOptimizedOPVAreaMaximizingtotalEconomicProfit(OPVAreaVector, totalEconomicProfitperYearVector):
    '''
    determine the best OPVArea maximizing the economic profit
    param:
        OPVAreaVector
        totalEconomicProfitperYearVector
    return:
        none
    '''
    maxtotalEconomicProfitperYear = np.max(totalEconomicProfitperYearVector)
    bestOPVArea = OPVAreaVector[np.argmax(totalEconomicProfitperYearVector)]
    print ("The OPV area maximizing the economic profit is {}m^2 the max economic profit is {}USD/year ".format(bestOPVArea, maxtotalEconomicProfitperYear))


def trainWeightsRLShadingCurtainDayStep(hasShadingCurtain, qLearningAgentsShadingCurtain=None, cropElectricityYieldSimulator1 = None):
  '''

  :param hasShadingCurtain:
  :param cropElectricityYieldSimulator1:
  :return:
  '''

  if hasShadingCurtain:

    # # set values necessary for RL training/testing
    # # for dLIEachdayThroughInnerStructure on a certain day
    # hourlyInnerLightIntensityPPFDThroughInnerStructure = cropElectricityYieldSimulator1.getHourlyInnerLightIntensityPPFDThroughInnerStructure()
    # # set dLIThroughInnerStructure to the object
    # dLIThroughInnerStructure = Util.convertFromHourlyPPFDWholeDayToDLI(hourlyInnerLightIntensityPPFDThroughInnerStructure)
    # qLearningAgentsShadingCurtain.setDLIThroughInnerStructure(dLIThroughInnerStructure)

    print ("training parameters: epsilon={}, gamma={}, alpha={}, period:{}".format(\
      qLearningAgentsShadingCurtain.epsilon, qLearningAgentsShadingCurtain.gamma, qLearningAgentsShadingCurtain.alpha,  constant.SimulationStartDate + "-" + constant.SimulationEndDate))
    for trainingIteration in range (0, qLearningAgentsShadingCurtain.numTraining):

      if trainingIteration % 100 == 0:
        # print("Iteration checkpoint: datetime.datetime.now():{}. trainingIteration:{}".format(datetime.datetime.now(), trainingIteration ))
        print("trainingIteration: {}, qLearningAgentsShadingCurtain.weights:{}, datetime.datetime.now():{}".format(\
          trainingIteration, qLearningAgentsShadingCurtain.weights, datetime.datetime.now()))

      # training the q value function
      for day in range (0, Util.getSimulationDaysInt()):

        state = day

        #########################################################################
        ############# set values necessary for RL training features##############
        #########################################################################
        # set day to the instance
        qLearningAgentsShadingCurtain.setDay(day)
        # dLIEachdayThroughInnerStructure on a certain day, necessary to cal DLI to PLants
        # qLearningAgentsShadingCurtain.setDLIEachDayThroughInnerStructure(dLIThroughInnerStructure[state])

        #set num of days from Jan 1st.
        daysFromJan1st = Util.getNumOfDaysFromJan1st(Util.getStartDateDateType() + datetime.timedelta(days=day))
        # date on a certain day
        qLearningAgentsShadingCurtain.setDaysFromJan1st(daysFromJan1st)

        # action = "openCurtain" or "closeCurtain"
        # if the state is at the terminal state, action is None.
        action = qLearningAgentsShadingCurtain.getAction(state)

        # if the q value is not initialized, initialize the q value. if initialized, just get the q value given state and action
        # state = qlearningAgentsShadingCurtain.getQValue(day, action)
        approximateQvalue = qLearningAgentsShadingCurtain.getApproximateQValue(state, action)
        # print ("approximateQvalue:{}".format(approximateQvalue))
        # set approximateQvalue to Q
        qLearningAgentsShadingCurtain.setApproximateQValue(approximateQvalue, state, action)


        # approximatedQvalueNextState = []
        # for action in qLearningAgentsShadingCurtain.getLegalActions(day):
        #   approximatedQvalueNextState.append(qLearningAgentsShadingCurtain.getApproximateQValue(day + 1, action))
        # approximateMaxQvalueNextState = max[approximatedQvalueNextState]

        # get the maximum q value in the next state
        if (state+1) == Util.getSimulationDaysInt():
          approximateMaxQvalueNextState = 0.0
        else:
          approximateMaxQvalueNextState = qLearningAgentsShadingCurtain.getApproximateValue(state + 1)

        # calc the difference between the current q value and the maximum q value in the next state, which is used for updating weights

        difference = (qLearningAgentsShadingCurtain.getReward(day) + approximateMaxQvalueNextState) - approximateQvalue
        # print ("qLearningAgentsShadingCurtain.getReward(day):{}".format(qLearningAgentsShadingCurtain.getReward(day)))
        # print ("approximateMaxQvalueNextState:{}".format(approximateMaxQvalueNextState))
        # print ("approximateQvalue:{}".format(approximateQvalue))
        # print ("difference:{}".format(difference))

        # update weight of the q learning function
        qLearningAgentsShadingCurtain.updateApproximateWeight(difference)

    # print ("qLearningAgentsShadingCurtain.weights:{}".format(qLearningAgentsShadingCurtain.weights))
    # print ("check trainingIteration:{}".format(trainingIteration))

    # print ("qLearningAgentsShadingCurtain.weights:{}".format(qLearningAgentsShadingCurtain.weights))
    print ("qLearningAgentsShadingCurtain.approximateQ:{}".format(qLearningAgentsShadingCurtain.approximateQ))

    return qLearningAgentsShadingCurtain
    # ####################################################################################################
    # Stop execution here...
    # sys.exit()
    # Move the above line to different parts of the assignment as you implement more of the functionality.
    # ####################################################################################################

def testWeightsRLShadingCurtainDayStep(hasShadingCurtain, qLearningAgentsShadingCurtain = None, cropElectricityYieldSimulator1=None):

  numTesting = qLearningAgentsShadingCurtain.numTesting

  if hasShadingCurtain:

    # change the exploration rate into zero because in testing, RL does not explore
    qLearningAgentsShadingCurtain.epsilon = 0.0
    # array to store the sales price at each iteration
    plantSalesperSquareMeterList = np.zeros(numTesting)

    for testingIteration in range(0, numTesting):

    # get values necessary for RL training, which was done at
    # hourlyInnerLightIntensityPPFDThroughInnerStructure = cropElectricityYieldSimulator1.getHourlyInnerLightIntensityPPFDThroughInnerStructure()
    # dLIThroughInnerStructure = Util.convertFromHourlyPPFDWholeDayToDLI(hourlyInnerLightIntensityPPFDThroughInnerStructure)
    # set dLIThroughInnerStructure to the object
    # qLearningAgentsShadingCurtain.setDLIThroughInnerStructure(dLIThroughInnerStructure)

      print("testingIteration: {}, qLearningAgentsShadingCurtain.weights:{}, datetime.datetime.now():{}, period:{}".format( \
        testingIteration, qLearningAgentsShadingCurtain.weights, datetime.datetime.now(), constant.SimulationStartDate + "-" + constant.SimulationEndDate ))

      # training the q value function
      for day in range(0, Util.getSimulationDaysInt()):

        state = day
        #########################################################################
        ############# set values necessary for RL training features##############
        #########################################################################
        # set day to the instance
        qLearningAgentsShadingCurtain.setDay(day)
        # dLIEachdayThroughInnerStructure on a certain day, necessary to cal DLI to PLants
        # qLearningAgentsShadingCurtain.setDLIEachDayThroughInnerStructure(dLIThroughInnerStructure[state])

        # set num of days from Jan 1st.
        daysFromJan1st = Util.getNumOfDaysFromJan1st(Util.getStartDateDateType() + datetime.timedelta(days=day))
        # date on a certain day
        qLearningAgentsShadingCurtain.setDaysFromJan1st(daysFromJan1st)

        # action = "openCurtain" or "closeCurtain"
        # if the state is at the terminal state, action is None.
        action = qLearningAgentsShadingCurtain.getPolicy(state)
        # store the action at each state at tuples in list for a record.
        qLearningAgentsShadingCurtain.policies[state] = action

        ################## calculate the daily plant yield start#####################

        #### calc the DLI on a certain state
        dLIEachDayThroughInnerStructure = qLearningAgentsShadingCurtain.getDLIThroughInnerStructureElement(state)

        dLIEachDayToPlants = 0.0
        if action == constant.openCurtainString:
          dLIEachDayToPlants = dLIEachDayThroughInnerStructure
        elif action == constant.closeCurtainString:
          dLIEachDayToPlants = dLIEachDayThroughInnerStructure * constant.shadingTransmittanceRatio

        #store the DLI ateach state by list for a record. since the sequence is important, not use a dictionary.
        qLearningAgentsShadingCurtain.dLIEachDayToPlants[day] = dLIEachDayToPlants

        ###### calc plant weight increase with a certain DLI
        # num of days from the latest seeding
        daysFromSeeding = state % constant.cultivationDaysperHarvest

        # if the calc method is A.J Both 2003 model
        if qLearningAgentsShadingCurtain.cropElectricityYieldSimulator1.getPlantGrowthModel() == constant.TaylorExpantionWithFluctuatingDLI:
          # daily [g/unit]
          unitDailyFreshWeightIncreaseElement = \
            Lettuce.calcUnitDailyFreshWeightIncreaseBoth2003Taylor(dLIEachDayToPlants, constant.cultivationDaysperHarvest, daysFromSeeding)
          # update the values to the instance
          qLearningAgentsShadingCurtain.setUnitDailyFreshWeightIncreaseElementShadingCurtain(unitDailyFreshWeightIncreaseElement, state)
          # print ("1 unitDailyFreshWeightIncrease [g/unit]:{}, state:{}".format(unitDailyFreshWeightIncreaseElement, state))

        else:
          print ("[test] error: feasture w_2 not considered. choosing un-existing plant growth model")
        ################## calculate the daily plant yield end#####################

      ################## calculate the total plant sales start#####################
      print ("DLI to plants at each day [mol/m^2/m^2]".format(qLearningAgentsShadingCurtain.dLIEachDayToPlants))

      unitPlantWeight = qLearningAgentsShadingCurtain.getUnitDailyFreshWeightIncreaseListShadingCurtain()
      print ("unitPlantWeight [g/unit]:{}".format(unitPlantWeight))
      totalUnitPlantWeight = sum(unitPlantWeight)


      # unit conversion; get the daily plant yield per given period per area: [g/unit] -> [g/m^2]
      unitPlantWeightperArea = Util.convertUnitShootFreshMassToShootFreshMassperArea(unitPlantWeight)
      # unit conversion:  [g/m^2] -> [kg/m^2]1
      unitPlantWeightperAreaKg = Util.convertFromgramTokilogram(unitPlantWeightperArea)

      # get the sales price of plant [USD/m^2]
      # if the average DLI during each harvest term is more than 17 mol/m^2/day, discount the price
      dailyPlantSalesperSquareMeter = getPlantSalesperSquareMeter(\
        cropElectricityYieldSimulator1.getYear(), unitPlantWeightperAreaKg, qLearningAgentsShadingCurtain.dLIEachDayToPlants)
      plantSalesperSquareMeter = sum(dailyPlantSalesperSquareMeter)
      plantSalesperSquareMeterList[testingIteration] = plantSalesperSquareMeter
      # print "dailyPlantSalesperSquareMeter.shape:{}".format(dailyPlantSalesperSquareMeter.shape)

      print ("plantSalesperSquareMeterList[{}]:{}".format(testingIteration, plantSalesperSquareMeterList))

      ################## calculate the total plant sakes end#####################

  else:
    print ("shading curtain assumed not to be given. the function without shading curtain will be made in the future")

  # return the average of testing results
  return plantSalesperSquareMeter






