Code: Select all
Status: dzVents: !Info: SOLAR 3: Function math.pow(x, y) has been deprecated in Lua 5.3. Please consider changing code to x^y
Moderator: leecollings
Code: Select all
Status: dzVents: !Info: SOLAR 3: Function math.pow(x, y) has been deprecated in Lua 5.3. Please consider changing code to x^y
Code: Select all
function math.pow(x, y) -- Function math.pow(x, y) has been deprecated in Lua 5.3.
return x^y
end
Code: Select all
--[[
W: Weather script
Prerequisites
==================================
Requires dz v2020 or later
SolarData: real-time solar data
from:
https://easydomoticz.com/forum/viewtopic.php?f=17&t=1340
https://www.domoticz.com/forum/viewtopic.php?p=148904#p148904
https://domoticz.com/forum/viewtopic.php?f=72&t=35541
https://www.domoticz.com/wiki/index.php?title=Lua_dzVents_-_Solar_Data:_Azimuth_Altitude_Lux
-- Authors ----------------------------------------------------------------
V1.0 - Sébastien Joly - Great original work
V1.1 - Neutrino - Adaptation to dz
V1.2 - Jmleglise - An acceptable approximation of the lux below 1° altitude for Dawn and dusk + translation + several changes to be more userfriendly.
V1.3 - Jmleglise - No update of the Lux data when <=0 to get the sunset and sunrise with lastUpdate
V1.4 - use the API instead of updateDevice to update the data of the virtual sensor to be able of using devicechanged['Lux'] in our scripts. (Due to a bug in dz that doesn't catch the devicechanged event of the virtual sensor)
V1.5 - xces - UTC time calculation.
V2.0 - BakSeeDaa - Converted to dzVents and changed quite many things.
V2.4.1-DarkSky - oredin - Use Dark Sky API instead of WU API
19/05/2019 - hestia - sources independent: use sensor with input values for CloudCover and Pressure
06/11/2019 - hestia - function math.pow(x, y) has been deprecated in Lua 5.3. => changing code to x^y
10/01/2002 - hestia - latitude and Longitude from settings - option to stop the script if lux and radiation not needed
07/02/2021 - hestia - modify the logging
22/02/2021 - hestia - from a cloud cover device instead of okta, remove the stop of the script (see 10/01/2002)
05/11/2021 - hestia - add a radiation device for the radiation no cloud, on the floor
25/02/2021 - hestia - checked with V3.1 (14/02/2021) - Jmleglise = Merge the different fork. Some clean Up, comment and wiki
19/08/2021 - hestia - add several radiation devices to have direct and scattered, with (weighted) or w/o cloud impact radiation separate, in order to calcultate the radiation on a surface with any angle from the ground.
The result given previously for the radiation was on the ground only, so not usable for a window (vertical) or a roof (w an angle).
Change the formula of the coefficient of mitigation M, to avoid the requirement of the pressure (no longer needed this weather info: more simple)
https://www.ftexploring.com/solar-energy/air-mass-and-insolation1.htm#fn1
+ refactoring of some parts
02/03/2021 - hestia - some cleaning et comments
--]]
-- Variables to customize ------------------------------------------------
-- Triggers of the script to choose
local idxCloudCover = 1913 -- device holding cloud coverage in percentage. Ex : 60% (for e.g. a device from Openweathermap...)
local intervalMins = 5 -- the interval of running this script. No need to be faster than the data source. (For example it is 10 min)
--local TEST = 204 -- a dummy switch for testing w/o waiting minutes / remove comment to use / comment to ignore
------
------ Dummy devices to save result (idx or 'name') (to create w dz GUI) ------
-- Custom Sensor SubType
local idxSolarAzimuth = 337 -- Azimuth of the sun
local idxSolarAltitude = 338 -- Altitude of the sun
-- total radiation with cloud impact, on the floor => the values saved with the V3.1
local idxRadiation = 304 -- Solar Radiation SubType (W/m2)
local idxLux = 305 -- Lux Type (Lux)
-- additional values from V3.1 - Solar Radiation SubType (W/m2)
local idxRadiationDirect = 1914 -- direct radiation with cloud impact, if no incidence
local idxRadiationScattered = 1915 -- scattered radiation with cloud impact
local idxRadiation0 = 1916 -- total radiation with cloud impact, if no incidence
-- radiation without cloud impact, on the floor (with incidence = 90° - sun altitude)
local idxRadiationNoCloud = 1970 -- solar Radiation (in Watt/m2)
--local idxRadiationNoCloudCustom = 1764 -- custom sensor for Solar Radiation (in Watt/m2) => to choose one of both
------
-- Parameters to give
-- Latitude and Longitude: custon values for this script ; if empty (nul) values from domoticz settings are taken
local script_latitude -- Latitude (Decimal number) Decimal Degrees ; something like 51.748485
local script_longitude -- Longitude (Decimal number) Decimal Degrees ; something like 5.629728
local LocalAltitude = 87 -- Altitude. (Integer) Meters above sea level
-- How to convert Lux to W/m2
-- "There is no simple conversion...it depends on the wavelength or color of the light
-- However, for the SUN, there is an approximate conversion of 0.0079 W/m2 per Lux"
-- see InverseSquareLawPresentation.pdf https://studyres.com/doc/16229818/light-vs.-distance p32
local LuxToWm2 = 0.0079
-- Logging to choose
local LOG_DEBUG = 0 -- 0 =>ERROR / 1 => FORCE / 2 => DEBUG
local LOG_LEVEL
local LOGGING
if LOG_DEBUG == 2 then
LOGGING = domoticz.LOG_DEBUG
LOG_LEVEL = domoticz.LOG_DEBUG
elseif LOG_DEBUG == 1 then
LOGGING = domoticz.LOG_FORCE
LOG_LEVEL = domoticz.LOG_FORCE
else
LOGGING = domoticz.LOG_ERROR
LOG_LEVEL = domoticz.LOG_INFO
end
return {
logging = {
level = LOGGING
},
on = {
--devices = {idxCloudCover},
devices = {TEST},
timer = {'every '..tostring(intervalMins)..' minutes'}
},
execute = function(dz, item, triggerInfo)
_G.logMarker = dz.moduleLabel -- set logmarker to scriptname
local _u = dz.utils
local function leapYear(year)
return year%4==0 and (year%100~=0 or year%400==0)
end
local latitude
if script_latitude ~= nil then
latitude = script_latitude
else
latitude = dz.settings.location.latitude
end
local longitude
if longitude ~= nil then
longitude = script_longitude
else
longitude = dz.settings.location.longitude
end
local arbitraryTwilightLux = 6.32 -- W/m² egal 800 Lux (the theoritical value is 4.74 but I have more accurate result with 6.32...) /!\
local constantSolarRadiation = 1361 -- Solar Constant W/m²
local year = os.date('%Y')
local numOfDay = os.date('%j')
local nbDaysInYear = (leapYear(year) and 366 or 365)
local angularSpeed = 360/365.25
local declination = math.deg(math.asin(0.3978 * math.sin(math.rad(angularSpeed) *(numOfDay - (81 - 2 * math.sin((math.rad(angularSpeed) * (numOfDay - 2))))))))
local timeDecimal = (os.date('!%H') + os.date('!%M') / 60) -- Coordinated Universal Time (UTC)
local solarHour = timeDecimal + (4 * longitude / 60 ) -- The solar Hour
local hourlyAngle = 15 * ( 12 - solarHour ) -- hourly Angle of the sun
local sunAltitude = math.deg(math.asin(math.sin(math.rad(latitude))* math.sin(math.rad(declination)) + math.cos(math.rad(latitude)) * math.cos(math.rad(declination)) * math.cos(math.rad(hourlyAngle))))-- the height of the sun in degree, compared with the horizon
local azimuth = math.acos((math.sin(math.rad(declination)) - math.sin(math.rad(latitude)) * math.sin(math.rad(sunAltitude))) / (math.cos(math.rad(latitude)) * math.cos(math.rad(sunAltitude) ))) * 180 / math.pi -- deviation of the sun from the North, in degree
local sinAzimuth = (math.cos(math.rad(declination)) * math.sin(math.rad(hourlyAngle))) / math.cos(math.rad(sunAltitude))
if sinAzimuth<0 then azimuth=360-azimuth end
local sunstrokeDuration = math.deg(2/15 * math.acos(- math.tan(math.rad(latitude)) * math.tan(math.rad(declination)))) -- duration of sunstroke in the day . Not used in this calculation.
local RadiationAtm = constantSolarRadiation * (1 +0.034 * math.cos( math.rad( 360 * numOfDay / nbDaysInYear ))) -- Sun radiation (in W/m²) in the entrance of atmosphere.
-- Coefficient of mitigation M
MM0 = math.exp(-0.0001184 * LocalAltitude)
--dz.log('MM0: ' .. MM0, LOG_LEVEL)
-- https://www.ftexploring.com/solar-energy/air-mass-and-insolation1.htm#fn1
local sinusSunAltitude = math.sin(math.rad(sunAltitude))
local M0 = math.sqrt(1229 + (614 * sinusSunAltitude)^2) - 614 * sinusSunAltitude
local M = M0 / MM0
dz.log('triggerInfo = ' .. triggerInfo.type, LOG_LEVEL)
dz.log('Latitude = ' .. latitude .. ', Longitude: ' .. longitude, LOG_LEVEL)
dz.log('Local Altitude = ' .. tostring(LocalAltitude) .. ' m', LOG_LEVEL)
dz.log('Angular Speed = ' .. angularSpeed .. ' per day', LOG_LEVEL)
dz.log('Declination = ' .. declination .. '°', LOG_LEVEL)
dz.log('Universal Coordinated Time (UTC) = '.. timeDecimal ..' H.dd', LOG_LEVEL)
dz.log('Solar Hour = '.. solarHour ..' H.dd', LOG_LEVEL)
dz.log('Altitude of the sun = ' .. _u.round(sunAltitude,0) .. '°', dz.LOG_FORCE)
dz.log('Angular hourly = '.. hourlyAngle .. '°', LOG_LEVEL)
dz.log('Azimuth of the sun = ' .. _u.round(azimuth,0) .. '°', dz.LOG_FORCE)
dz.log('Duration of the sun stroke of the day = ' .. _u.round(sunstrokeDuration,2) ..' H.dd', LOG_LEVEL)
dz.log('Radiation max in atmosphere = ' .. _u.round(RadiationAtm,2) .. ' W/m²', LOG_LEVEL)
dz.log('Coefficient of mitigation M = ' .. M ..' M0 = '..M0, LOG_LEVEL)
-- Save the results : sun azimuth and altitude
dz.devices(idxSolarAzimuth).updateCustomSensor(_u.round(azimuth,0))
dz.devices(idxSolarAltitude).updateCustomSensor(_u.round(sunAltitude,0))
local cloudPercentage = dz.devices(idxCloudCover).percentage
local okta = cloudPercentage / 12.5
local Kc = 1 - 0.75 * ((okta / 8)^3.4) -- Factor of mitigation for the cloud layer
dz.log('Okta = '.. okta, LOG_LEVEL)
dz.log('Kc = ' .. Kc, LOG_LEVEL)
local directRadiation, directRadiationFloor, scatteredRadiation
if sunAltitude > 1 then -- Below 1° of Altitude , the formulae reach their limit of precision.
directRadiation = RadiationAtm * (0.6^M) -- no take into account the slope of the target
directRadiationFloor = directRadiation * sinusSunAltitude -- sinusSunAltitude = the cosinus of incidence (the floor is horizontal)
scatteredRadiation = RadiationAtm * (0.271 - 0.294 * (0.6^M)) * sinusSunAltitude -- this sinus is not the incidence, no incidence effect: "scattered"
elseif sunAltitude <= 1 and sunAltitude >= -7 then -- apply theoretical Lux of twilight
directRadiation = 0
directRadiationFloor = 0
scatteredRadiation = arbitraryTwilightLux-(1-sunAltitude)/8*arbitraryTwilightLux -- in scatteredRadiation to simplify the code
elseif sunAltitude < -7 then -- no management of nautical and astronomical twilight...
directRadiation = 0
directRadiationFloor = 0
scatteredRadiation = 0
end
local totalRadiationFloor, weightedTotalRadiationFloor, weightedDirectRadiation, weightedScatteredRadiation, weightedTotalRadiation
totalRadiationFloor = directRadiationFloor + scatteredRadiation -- radiation on the floor if no cloud
weightedTotalRadiationFloor = totalRadiationFloor * Kc -- with cloud impact: what could give a weather station (sensor horizontal)
local weightedLux = _u.round(weightedTotalRadiationFloor / LuxToWm2, 0) -- Radiation in Lux
weightedTotalRadiationFloor = _u.round(weightedTotalRadiationFloor,0)
totalRadiationFloor = _u.round(totalRadiationFloor, 0) -- radiation on the floor if no cloud
weightedDirectRadiation = _u.round(directRadiation * Kc, 0) -- cloud, no incidence
weightedScatteredRadiation = _u.round(scatteredRadiation * Kc,0) -- cloud, incidence no matters
weightedTotalRadiation = _u.round((directRadiation + scatteredRadiation) * Kc, 0) -- cloud, no incidence
-- with no incidence with the target
dz.log('Scattered Radiation = '.. _u.round(scatteredRadiation,0) ..' W/m²', LOG_LEVEL)
dz.log('Scattered Radiation with Cloud impact = '.. weightedScatteredRadiation..' W/m²', LOG_LEVEL) -- (Scattered Radiation) * Kc
dz.devices(idxRadiationScattered).updateRadiation(weightedScatteredRadiation)
dz.log('Direct Radiation = '.. _u.round(directRadiation,0) ..' W/m²', LOG_LEVEL)
dz.log('Direct Radiation with Cloud impact = '.. weightedDirectRadiation..' W/m²', LOG_LEVEL) -- (Direct Radiation Floor) * Kc
dz.devices(idxRadiationDirect).updateRadiation(weightedDirectRadiation)
dz.log('Total Radiation with Cloud impact = '.. weightedTotalRadiation..' W/m²', dz.LOG_FORCE) -- (Direct Radiation + Scattered Radiation) * Kc
dz.devices(idxRadiation0).updateRadiation(weightedTotalRadiation)
-- projection on the floor
dz.log('Direct Radiation Floor = '.. _u.round(directRadiationFloor,0) ..' W/m²', LOG_LEVEL) -- Direct Radiation * sinus(SunAltitude)
dz.log('Total Radiation Floor = ' .. totalRadiationFloor ..' W/m²', LOG_LEVEL) -- Direct Radiation Floor + Scattered Radiation
if idxRadiationNoCloud ~= nil then
dz.devices(idxRadiationNoCloud).updateRadiation(totalRadiationFloor)
end
if idxRadiationNoCloudCustom ~= nil then
dz.devices(idxRadiationNoCloudCustom).updateCustomSensor(totalRadiationFloor)
end
-- projection on the floor with cloud impact (weighted by Kc)
dz.log('Total Radiation Floor with Cloud impact = '.. weightedTotalRadiationFloor..' W/m²', LOG_LEVEL) -- (Total Radiation Floor) * Kc
dz.devices(idxRadiation).updateRadiation(weightedTotalRadiationFloor)
dz.log('Total Radiation Floor with Cloud impact = '.. weightedLux ..' Lux', LOG_LEVEL) -- (Total Radiation Floor) * Kc converted in Lux
-- No update if Lux is already 0. So lastUpdate of the Lux sensor will keep the time when Lux has reached 0.
-- (Kind of timeofday['SunsetInMinutes'])
if dz.devices(idxLux).lux + weightedLux > 0 then
dz.devices(idxLux).updateLux(weightedLux)
end
end
}
If your info is from a weather station, from what I saw on the doc of some of them, the sensor in horizontal, so it is not the direct normal radiance, but the full radiance on a horizontal area.I also have a device with "Sunpower" value in Watt/m2 from a nearby weatherstation, which I assume is the Direct Normal Irradiance DNI, measured on a surface perpendicular to the sun beam, and therefore already incorporating cloud cover effect. Using your sun altitude and azimuth calculation and the angle and azimuth of my PV installation, I can easily calculate the beam irradiance.
you could have a look at this scripttheoretical performance of my PVpanel installation and to then compare this to actual performance
Code: Select all
math.max(_u.round(F_solarRadiation * F_cosIncidence + F_solarRadiationScattered, 0), 0)
Agreed, so far I'm using info from WU weather stations near me that give the full radiation (on an horizontal area) and make a calculation of cloud cover from the solar radiation script that give the radiation w/o cloud(I am so far not using the openweathermap cloud coverage since it is very inaccurate)
I'm not at the inception of the script, it's a long story... but I'm almost sure it is from an observer on the surface of the earth.willemd wrote: ↑Thursday 11 August 2022 17:46 ...
Does the Azimuth and Altitude calculation give the position of the Sun vs. a point on the earth surface as it is at that moment in the universe (i.e. in terms of positions and angles), or as it is seen at that moment from earth (in terms of observed altitude above the horizon)? If it is the position in the universe then we will only see the sun at that position 8 minutes and 20 seconds later, right? And we will only see and feel the sun power at the max when the actual position is already past the highest altitude angle.
...
I don't know and I don't have such a sensor. It could be relative to the position of the sensor and the environment.I see a definitive time shift between the theoretical solar radiation calculated using this altitude and the actual solar radiation from a local weather station on a clear day (we had quite a few recently). The max in sunpower (actualGHI) occurs later than max theoretical GHI
Found it, so simple: there is a timestamp in the json file received from the buienradar weatherstation. The timestamp is not shown/used in the standard domoticz buienradar interface, but when monitoring the update I can see there is a new timestamp every 10 minutes and when it is updated the timestamp is approx. 20 minutes old, so there is a 20 minute delay and it remain that way for 10 minutes, so average delay is 25 minutes, which is exactly what I see when comparing the sinus curve of the theoretical GHI and the actual sunpower measured by the weatherstation. (maybe average does not make sense, it just depends on when exactly the buienradar interface runs)willemd wrote: ↑Thursday 11 August 2022 17:46
Another explanation could also be that it is because the weather station provides an average sunpower measurement of the past 10 minutes instead of an actual value at that moment and is retrieved only every 10 minutes.
Or maybe both effects? The time shift seems to be approx 25 to 30 minutes. The max sunpower measured is approx 25 minutes after solar noon.
Code: Select all
import socketserver
import urllib.request
import math
from datetime import datetime, timedelta
from suncalc import get_position, get_times
date = datetime.utcnow() + timedelta(hours=2)
lon=4.xxxxxxxxxxxxxxx
lat=52.xxxxxxxxxxxxxx
calc = get_position(date, lon, lat)
azi = round(calc.get("azimuth") * (180/math.pi) + 180,1)
alt = round(calc.get("altitude") * (180/math.pi),1)
print(azi)
print(alt)
try:
url_domo = "http://[DOMO-IP:PORT]/json.htm?type=command¶m=udevice&idx=[SENSOR-IDX]&nvalue=0&svalue="
urllib.request.urlopen(url_domo + str(azi),None,0.1)
url_domo2 = "http://[DOMO-IP:PORT]/json.htm?type=command¶m=udevice&idx=[SENSOR-IDX]&nvalue=0&svalue="
urllib.request.urlopen(url_domo2 + str(alt),None,0.1)
except:
pass
exit()
Hi,ouiouiblog wrote: ↑Sunday 07 April 2024 14:44 I was wondering if someone on the forum understood the formulas that are in that script or had some sources to explain them.
Users browsing this forum: No registered users and 0 guests