mirror of
https://github.com/fhem/fhem-mirror.git
synced 2025-01-31 06:39:11 +00:00
1380 lines
48 KiB
Perl
1380 lines
48 KiB
Perl
# $Id$
|
||
##############################################################################################################
|
||
#
|
||
# 98_feels_like.pm
|
||
#
|
||
# Module that listens to events generated by weather sensors in order to compute
|
||
# - radiation
|
||
# - clouds coverage
|
||
# - equivalent temperature according to UTCI and Apparent Temperature (AT)
|
||
#
|
||
# using the solar radiation sensor of Fine Offset / Ambient Weather / WH2601 ... weather stations.
|
||
#
|
||
# See documentation at end.
|
||
#
|
||
# written by holger.a.teutsch@gmail.com
|
||
#
|
||
|
||
use strict;
|
||
use warnings;
|
||
|
||
# so "perl -w 98_feels_like.pm" should produce no syntax errors
|
||
use vars qw($readingFnAttributes %defs);
|
||
|
||
sub feels_like_Initialize($$)
|
||
{
|
||
my ($hash) = @_;
|
||
$hash->{DefFn} = "feels_like_Define";
|
||
$hash->{NotifyFn} = "feels_like_Notify";
|
||
$hash->{AttrFn} = "feels_like_Attr";
|
||
$hash->{NotifyOrderPrefix} = "48-"; # Want to be called earlier than default
|
||
$hash->{AttrList} = "disable:0,1 maxTimediff sensorSkyViewFactor skyViewFactor " .
|
||
"bowenRatio linkeTurbity directRadiationRatio coverageCallback utciWindCutoff " .
|
||
"sunVisibility alphaMin alphaMax sensorType:wh2601,generic $readingFnAttributes";
|
||
}
|
||
|
||
##########################
|
||
sub
|
||
feels_like_Define($$)
|
||
{
|
||
my ($hash, $def) = @_;
|
||
my @a = split(/\s+/, $def);
|
||
|
||
return "wrong syntax: " .
|
||
"define <name> feels_like devicename temperature humidity " .
|
||
"[wind_speed [solar_radiation [pressure]]]" if (@a < 4 || @a> 8);
|
||
|
||
my $name = $a[0];
|
||
my $devname = $a[2];
|
||
|
||
my $match = 0;
|
||
my $i = 3;
|
||
foreach my $k ('T', 'H', 'W', 'S', 'P') {
|
||
my $rk = $k . '_READING';
|
||
my $dk = $k . '_DEV';
|
||
|
||
# in case it's a redefine
|
||
delete($hash->{$dk});
|
||
delete($hash->{$rk});
|
||
|
||
if($#a >= $i) {
|
||
my @x = split(/:/, $a[$i]);
|
||
if (@x == 2) {
|
||
$hash->{$dk} = $x[0];
|
||
$hash->{$rk} = $x[1];
|
||
} else {
|
||
$hash->{$dk} = $devname;
|
||
$hash->{$rk} = $x[0];
|
||
$match++;
|
||
}
|
||
}
|
||
|
||
$i++;
|
||
}
|
||
|
||
return "feels_like: at least one reading device must match $devname" if ($match == 0);
|
||
|
||
$hash->{DEVNAME} = $devname;
|
||
|
||
# set NOTIFYDEV
|
||
$hash->{NOTIFYDEV} = $devname;
|
||
$hash->{STATE} = "active";
|
||
return undef;
|
||
}
|
||
|
||
##########################
|
||
sub
|
||
feels_like_Attr(@)
|
||
{
|
||
my ($cmd, $name, $a_name, $a_val) = @_;
|
||
my $hash = $defs{$name};
|
||
|
||
if ($cmd eq "set" && $a_name eq "sunVisibility") {
|
||
my @x = split(/, */, $a_val);
|
||
if ((@x & 1) == 0 || $x[0] != 0 || $x[$#x] != 360) {
|
||
return "Value of sunVisibility must must be a list odd # of elements "
|
||
. "starting with 0 and ending with 360";
|
||
}
|
||
}
|
||
|
||
if ($a_name eq 'disable') {
|
||
if ($cmd eq 'set') {
|
||
$hash->{DISABLED} = $a_val;
|
||
$hash->{STATE} = $a_val == 1 ? 'disabled' : 'active';
|
||
}
|
||
elsif ($cmd eq 'del') {
|
||
$hash->{DISABLED} = 0;
|
||
$hash->{STATE} = 'active';
|
||
}
|
||
}
|
||
|
||
return undef;
|
||
}
|
||
|
||
##########################
|
||
sub
|
||
feels_like_Notify($$)
|
||
{
|
||
my ($hash, $dev) = @_;
|
||
|
||
my $hash_name = $hash->{NAME};
|
||
my $dev_name = $dev->{NAME};
|
||
|
||
return undef if ($hash->{DISABLED});
|
||
|
||
Log3($hash_name, 4, "feels_like_Notify: $hash_name for $dev_name");
|
||
|
||
my $max_timediff = AttrVal($hash_name, "maxTimediff", 10 * 60);
|
||
|
||
# extract readings from event queue
|
||
my @R = ('', '', '', '', '');
|
||
|
||
my $events = deviceEvents($dev, undef);
|
||
return undef if (!$events);
|
||
|
||
my $have_one = 0;
|
||
|
||
foreach my $ev (@{$events}) {
|
||
next if (!defined($ev));
|
||
|
||
my ($ev_name, $val, $rest) = split(/: +/, $ev, 3);
|
||
next if (!defined($ev_name) || !defined($val));
|
||
Log3($hash_name, 5, "feels_like_Notify ev_name='$ev_name' val='$val'");
|
||
|
||
my $i = 0;
|
||
foreach my $k ('T', 'H', 'W', 'S', 'P') {
|
||
$i++;
|
||
my $dk = $k . '_DEV';
|
||
next if (!exists($hash->{$dk}) || ($hash->{$dk} ne $dev_name));
|
||
|
||
my $rk = $k . '_READING';
|
||
next if (!exists($hash->{$rk}) || ($hash->{$rk} ne $ev_name));
|
||
$R[$i-1] = $val;
|
||
$have_one = 1;
|
||
}
|
||
}
|
||
|
||
# at least one match required
|
||
return undef unless $have_one;
|
||
|
||
my ($T, $H, $W, $S, $P) = @R;
|
||
Log3($hash_name, 4, "feels_like_Notify $dev_name events T: $T, H: $H, W: $W, S: $S, P: $P");
|
||
|
||
# now get readings that were not filled by events
|
||
my $i = 0;
|
||
foreach my $k ('T', 'H', 'W', 'S', 'P') {
|
||
$i++;
|
||
my $rk = $k . '_READING';
|
||
next if (!exists($hash->{$rk}) || $R[$i-1]);
|
||
|
||
my $dk = $k . '_DEV';
|
||
my $val = ReadingsNum($hash->{$dk}, $hash->{$rk}, undef);
|
||
if (!defined($val)) {
|
||
Log(1, "feels_like: reading $hash->{$dk}:$hash->{$rk} does not exist, sorry!");
|
||
return undef;
|
||
}
|
||
|
||
my $ra = ReadingsAge($hash->{$dk}, $hash->{$rk}, undef);
|
||
if (!defined($ra) || $ra > $max_timediff) {
|
||
Log(1, "feels_like: reading $hash->{$dk}:$hash->{$rk} is too old, sorry!");
|
||
return undef;
|
||
}
|
||
|
||
$R[$i-1] = $val;
|
||
}
|
||
|
||
($T, $H, $W, $S, $P) = @R;
|
||
Log3($hash_name, 4, "feels_like_Notify final T: $T, H: $H, W: $W, S: $S, P: $P");
|
||
|
||
|
||
# running as device, save name for log level determination
|
||
# in library calls
|
||
$feels_like::dev_instance = $hash_name;
|
||
|
||
readingsBeginUpdate($dev);
|
||
my ($T_tmrt, $ER, $coverage, $oktas, $clouds, $alt, $azimuth);
|
||
|
||
$W = 0 if ($W eq '');
|
||
$P = 1013 if ($P eq '');
|
||
|
||
if ($S ne '') {
|
||
readingsBeginUpdate($hash);
|
||
($T_tmrt, $ER) = mean_radiant_temperature($hash_name, $T, $H, $W, $P, $S);
|
||
|
||
readingsBulkUpdate($hash, "sun_visible", $hash->{helper}{sun_visible});
|
||
readingsBulkUpdate($hash, "sr_invalid", $hash->{helper}{sr_invalid});
|
||
|
||
my $x = $hash->{helper}{alpha};
|
||
readingsBulkUpdate($hash, "alpha", defined($x) ? round($x, 2) : '');
|
||
$x = $hash->{helper}{nr_stddev};
|
||
readingsBulkUpdate($hash, "nr_stddev", defined($x) ? round($x, 3) : '');
|
||
$x = $hash->{helper}{nr_current};
|
||
readingsBulkUpdate($hash, "nr_current", defined($x) ? round($x, 3) : '');
|
||
|
||
my ($oktas, $clouds);
|
||
my $coverage = $hash->{helper}{nr_coverage};
|
||
if (defined($coverage)) {
|
||
$coverage = round($coverage, 2);
|
||
$oktas = round($coverage * 8, 0);
|
||
$clouds = ('SKC', 'FEW', 'FEW', 'SCT', 'SCT', 'BKN', 'BKN', 'BKN', 'OVC')[$oktas];
|
||
} else {
|
||
$coverage = $oktas = $clouds = '';
|
||
}
|
||
|
||
readingsBulkUpdate($hash, "nr_coverage", $coverage);
|
||
readingsBulkUpdate($hash, "nr_oktas", $oktas);
|
||
readingsBulkUpdate($hash, "nr_clouds", $clouds);
|
||
|
||
$coverage = $hash->{helper}{coverage};
|
||
if (defined($coverage)) {
|
||
$coverage = round($coverage, 2);
|
||
$oktas = round($coverage * 8, 0);
|
||
$clouds = ('SKC', 'FEW', 'FEW', 'SCT', 'SCT', 'BKN', 'BKN', 'BKN', 'OVC')[$oktas];
|
||
} else {
|
||
$coverage = $oktas = $clouds = '';
|
||
}
|
||
|
||
readingsBulkUpdate($dev, "extra_radiation", round($ER, 1));
|
||
readingsBulkUpdate($dev, "coverage", $coverage);
|
||
readingsBulkUpdate($dev, "oktas", $oktas);
|
||
readingsBulkUpdate($dev, "clouds", $clouds);
|
||
|
||
readingsEndUpdate($hash, 1);
|
||
} else {
|
||
$T_tmrt = $T;
|
||
$ER = 0;
|
||
}
|
||
|
||
# for higher winds UTCI seems a bit odd
|
||
my $max_W = AttrVal($hash_name, "utciWindCutoff", 3.0);
|
||
my $T_utci = T_UTCI($T, $H, min($W, $max_W), $P, $T_tmrt);
|
||
my $AT = T_apparent_temperature($T, $H, $W, $ER);
|
||
|
||
readingsBulkUpdate($dev, "temperature_utci", round($T_utci, 1));
|
||
readingsBulkUpdate($dev, "temperature_mrt", round($T_tmrt, 1));
|
||
readingsBulkUpdate($dev, "temperature_at", round($AT, 1));
|
||
readingsEndUpdate($dev, 1);
|
||
|
||
$feels_like::dev_instance = undef;
|
||
return undef;
|
||
}
|
||
|
||
# don't pollute main namespace
|
||
package feels_like;
|
||
|
||
use POSIX;
|
||
use Math::Trig qw (deg2rad rad2deg);
|
||
|
||
# this stuff is needed so the core functions can be called (and debugged) outside of FHEM
|
||
# without too much pain
|
||
our $verbose = 1;
|
||
our $dev_instance = undef;
|
||
|
||
sub Log($$)
|
||
{
|
||
my ($level, $str) = @_;
|
||
if (defined($dev_instance)) {
|
||
main::Log3($dev_instance, $level, "feels_like_core " . $str);
|
||
} else {
|
||
return unless($level <= $verbose);
|
||
main::Log(1, "feels_like_lib " . $str);
|
||
}
|
||
}
|
||
|
||
|
||
# Derived by Holger Teutsch from -- UTCI_a002.f90 -- http://www.utci.org
|
||
|
||
#~ value is the UTCI in degree Celsius
|
||
#~ computed by a 6th order approximating polynomial from the 4 Input paramters
|
||
#~
|
||
#~ Input parameters
|
||
#~ - Ta : air temperature, degree Celsius
|
||
#~ - ehPa : water $vapour presure, hPa=hecto Pascal
|
||
#~ - Tmrt : mean radiant temperature, degree Celsius
|
||
#~ - va : wind speed 10 m above ground level in m/s
|
||
#~
|
||
#~ UTCI_approx, Version a 0.002, October 2009
|
||
#~ Copyright (C) 2009 Peter Broede
|
||
|
||
sub UTCI_approx($$$$)
|
||
{
|
||
my ($Ta, $ehPa, $Tmrt, $va) = @_;
|
||
|
||
|
||
my $D_Tmrt = $Tmrt - $Ta;
|
||
|
||
# approx is only valid for va >= 0.5 m/s
|
||
$va = 0.5 if ($va < 0.5);
|
||
|
||
my $Pa = $ehPa/10.0; #~ use $vapour pressure in kPa
|
||
|
||
#~ calculate 6th order polynomial as approximation
|
||
|
||
my $Ta_2 = $Ta*$Ta;
|
||
my $Ta_3 = $Ta_2*$Ta;
|
||
my $Ta_4 = $Ta_3*$Ta;
|
||
my $Ta_5 = $Ta_4*$Ta;
|
||
|
||
my $va_2 = $va*$va;
|
||
my $va_3 = $va_2*$va;
|
||
my $va_4 = $va_3*$va;
|
||
my $va_5 = $va_4*$va;
|
||
|
||
my $D_Tmrt_2 = $D_Tmrt*$D_Tmrt;
|
||
my $D_Tmrt_3 = $D_Tmrt_2*$D_Tmrt;
|
||
my $D_Tmrt_4 = $D_Tmrt_3*$D_Tmrt;
|
||
my $D_Tmrt_5 = $D_Tmrt_4*$D_Tmrt;
|
||
|
||
my $Pa_2 = $Pa*$Pa;
|
||
my $Pa_3 = $Pa_2*$Pa;
|
||
my $Pa_4 = $Pa_3*$Pa;
|
||
my $Pa_5 = $Pa_4*$Pa;
|
||
|
||
my $UTCI_approx = $Ta +
|
||
6.07562052E-01 +
|
||
-2.27712343E-02 * $Ta +
|
||
8.06470249E-04 * $Ta_2 +
|
||
-1.54271372E-04 * $Ta_3 +
|
||
-3.24651735E-06 * $Ta_4 +
|
||
7.32602852E-08 * $Ta_5 +
|
||
1.35959073E-09 * $Ta_5*$Ta +
|
||
-2.25836520E+00 * $va +
|
||
8.80326035E-02 * $Ta*$va +
|
||
2.16844454E-03 * $Ta_2*$va +
|
||
-1.53347087E-05 * $Ta_3*$va +
|
||
-5.72983704E-07 * $Ta_4*$va +
|
||
-2.55090145E-09 * $Ta_5*$va +
|
||
-7.51269505E-01 * $va_2 +
|
||
-4.08350271E-03 * $Ta*$va_2 +
|
||
-5.21670675E-05 * $Ta_2*$va_2 +
|
||
1.94544667E-06 * $Ta_3*$va_2 +
|
||
1.14099531E-08 * $Ta_4*$va_2 +
|
||
1.58137256E-01 * $va_3 +
|
||
-6.57263143E-05 * $Ta*$va_3 +
|
||
2.22697524E-07 * $Ta_2*$va_3 +
|
||
-4.16117031E-08 * $Ta_3*$va_3 +
|
||
-1.27762753E-02 * $va_4 +
|
||
9.66891875E-06 * $Ta*$va_4 +
|
||
2.52785852E-09 * $Ta_2*$va_4 +
|
||
4.56306672E-04 * $va_5 +
|
||
-1.74202546E-07 * $Ta*$va_5 +
|
||
-5.91491269E-06 * $va_5*$va +
|
||
3.98374029E-01 * $D_Tmrt +
|
||
1.83945314E-04 * $Ta*$D_Tmrt +
|
||
-1.73754510E-04 * $Ta_2*$D_Tmrt +
|
||
-7.60781159E-07 * $Ta_3*$D_Tmrt +
|
||
3.77830287E-08 * $Ta_4*$D_Tmrt +
|
||
5.43079673E-10 * $Ta_5*$D_Tmrt +
|
||
-2.00518269E-02 * $va*$D_Tmrt +
|
||
8.92859837E-04 * $Ta*$va*$D_Tmrt +
|
||
3.45433048E-06 * $Ta_2*$va*$D_Tmrt +
|
||
-3.77925774E-07 * $Ta_3*$va*$D_Tmrt +
|
||
-1.69699377E-09 * $Ta_4*$va*$D_Tmrt +
|
||
1.69992415E-04 * $va_2*$D_Tmrt +
|
||
-4.99204314E-05 * $Ta*$va_2*$D_Tmrt +
|
||
2.47417178E-07 * $Ta_2*$va_2*$D_Tmrt +
|
||
1.07596466E-08 * $Ta_3*$va_2*$D_Tmrt +
|
||
8.49242932E-05 * $va_3*$D_Tmrt +
|
||
1.35191328E-06 * $Ta*$va_3*$D_Tmrt +
|
||
-6.21531254E-09 * $Ta_2*$va_3*$D_Tmrt +
|
||
-4.99410301E-06 * $va_4*$D_Tmrt +
|
||
-1.89489258E-08 * $Ta*$va_4*$D_Tmrt +
|
||
8.15300114E-08 * $va_5*$D_Tmrt +
|
||
7.55043090E-04 * $D_Tmrt_2 +
|
||
-5.65095215E-05 * $Ta*$D_Tmrt_2 +
|
||
-4.52166564E-07 * $Ta_2*$D_Tmrt_2 +
|
||
2.46688878E-08 * $Ta_3*$D_Tmrt_2 +
|
||
2.42674348E-10 * $Ta_4*$D_Tmrt_2 +
|
||
1.54547250E-04 * $va*$D_Tmrt_2 +
|
||
5.24110970E-06 * $Ta*$va*$D_Tmrt_2 +
|
||
-8.75874982E-08 * $Ta_2*$va*$D_Tmrt_2 +
|
||
-1.50743064E-09 * $Ta_3*$va*$D_Tmrt_2 +
|
||
-1.56236307E-05 * $va_2*$D_Tmrt_2 +
|
||
-1.33895614E-07 * $Ta*$va_2*$D_Tmrt_2 +
|
||
2.49709824E-09 * $Ta_2*$va_2*$D_Tmrt_2 +
|
||
6.51711721E-07 * $va_3*$D_Tmrt_2 +
|
||
1.94960053E-09 * $Ta*$va_3*$D_Tmrt_2 +
|
||
-1.00361113E-08 * $va_4*$D_Tmrt_2 +
|
||
-1.21206673E-05 * $D_Tmrt_3 +
|
||
-2.18203660E-07 * $Ta*$D_Tmrt_3 +
|
||
7.51269482E-09 * $Ta_2*$D_Tmrt_3 +
|
||
9.79063848E-11 * $Ta_3*$D_Tmrt_3 +
|
||
1.25006734E-06 * $va*$D_Tmrt_3 +
|
||
-1.81584736E-09 * $Ta*$va*$D_Tmrt_3 +
|
||
-3.52197671E-10 * $Ta_2*$va*$D_Tmrt_3 +
|
||
-3.36514630E-08 * $va_2*$D_Tmrt_3 +
|
||
1.35908359E-10 * $Ta*$va_2*$D_Tmrt_3 +
|
||
4.17032620E-10 * $va_3*$D_Tmrt_3 +
|
||
-1.30369025E-09 * $D_Tmrt_4 +
|
||
4.13908461E-10 * $Ta*$D_Tmrt_4 +
|
||
9.22652254E-12 * $Ta_2*$D_Tmrt_4 +
|
||
-5.08220384E-09 * $va*$D_Tmrt_4 +
|
||
-2.24730961E-11 * $Ta*$va*$D_Tmrt_4 +
|
||
1.17139133E-10 * $va_2*$D_Tmrt_4 +
|
||
6.62154879E-10 * $D_Tmrt_5 +
|
||
4.03863260E-13 * $Ta*$D_Tmrt_5 +
|
||
1.95087203E-12 * $va*$D_Tmrt_5 +
|
||
-4.73602469E-12 * $D_Tmrt_5*$D_Tmrt +
|
||
5.12733497E+00 * $Pa +
|
||
-3.12788561E-01 * $Ta*$Pa +
|
||
-1.96701861E-02 * $Ta_2*$Pa +
|
||
9.99690870E-04 * $Ta_3*$Pa +
|
||
9.51738512E-06 * $Ta_4*$Pa +
|
||
-4.66426341E-07 * $Ta_5*$Pa +
|
||
5.48050612E-01 * $va*$Pa +
|
||
-3.30552823E-03 * $Ta*$va*$Pa +
|
||
-1.64119440E-03 * $Ta_2*$va*$Pa +
|
||
-5.16670694E-06 * $Ta_3*$va*$Pa +
|
||
9.52692432E-07 * $Ta_4*$va*$Pa +
|
||
-4.29223622E-02 * $va_2*$Pa +
|
||
5.00845667E-03 * $Ta*$va_2*$Pa +
|
||
1.00601257E-06 * $Ta_2*$va_2*$Pa +
|
||
-1.81748644E-06 * $Ta_3*$va_2*$Pa +
|
||
-1.25813502E-03 * $va_3*$Pa +
|
||
-1.79330391E-04 * $Ta*$va_3*$Pa +
|
||
2.34994441E-06 * $Ta_2*$va_3*$Pa +
|
||
1.29735808E-04 * $va_4*$Pa +
|
||
1.29064870E-06 * $Ta*$va_4*$Pa +
|
||
-2.28558686E-06 * $va_5*$Pa +
|
||
-3.69476348E-02 * $D_Tmrt*$Pa +
|
||
1.62325322E-03 * $Ta*$D_Tmrt*$Pa +
|
||
-3.14279680E-05 * $Ta_2*$D_Tmrt*$Pa +
|
||
2.59835559E-06 * $Ta_3*$D_Tmrt*$Pa +
|
||
-4.77136523E-08 * $Ta_4*$D_Tmrt*$Pa +
|
||
8.64203390E-03 * $va*$D_Tmrt*$Pa +
|
||
-6.87405181E-04 * $Ta*$va*$D_Tmrt*$Pa +
|
||
-9.13863872E-06 * $Ta_2*$va*$D_Tmrt*$Pa +
|
||
5.15916806E-07 * $Ta_3*$va*$D_Tmrt*$Pa +
|
||
-3.59217476E-05 * $va_2*$D_Tmrt*$Pa +
|
||
3.28696511E-05 * $Ta*$va_2*$D_Tmrt*$Pa +
|
||
-7.10542454E-07 * $Ta_2*$va_2*$D_Tmrt*$Pa +
|
||
-1.24382300E-05 * $va_3*$D_Tmrt*$Pa +
|
||
-7.38584400E-09 * $Ta*$va_3*$D_Tmrt*$Pa +
|
||
2.20609296E-07 * $va_4*$D_Tmrt*$Pa +
|
||
-7.32469180E-04 * $D_Tmrt_2*$Pa +
|
||
-1.87381964E-05 * $Ta*$D_Tmrt_2*$Pa +
|
||
4.80925239E-06 * $Ta_2*$D_Tmrt_2*$Pa +
|
||
-8.75492040E-08 * $Ta_3*$D_Tmrt_2*$Pa +
|
||
2.77862930E-05 * $va*$D_Tmrt_2*$Pa +
|
||
-5.06004592E-06 * $Ta*$va*$D_Tmrt_2*$Pa +
|
||
1.14325367E-07 * $Ta_2*$va*$D_Tmrt_2*$Pa +
|
||
2.53016723E-06 * $va_2*$D_Tmrt_2*$Pa +
|
||
-1.72857035E-08 * $Ta*$va_2*$D_Tmrt_2*$Pa +
|
||
-3.95079398E-08 * $va_3*$D_Tmrt_2*$Pa +
|
||
-3.59413173E-07 * $D_Tmrt_3*$Pa +
|
||
7.04388046E-07 * $Ta*$D_Tmrt_3*$Pa +
|
||
-1.89309167E-08 * $Ta_2*$D_Tmrt_3*$Pa +
|
||
-4.79768731E-07 * $va*$D_Tmrt_3*$Pa +
|
||
7.96079978E-09 * $Ta*$va*$D_Tmrt_3*$Pa +
|
||
1.62897058E-09 * $va_2*$D_Tmrt_3*$Pa +
|
||
3.94367674E-08 * $D_Tmrt_4*$Pa +
|
||
-1.18566247E-09 * $Ta*$D_Tmrt_4*$Pa +
|
||
3.34678041E-10 * $va*$D_Tmrt_4*$Pa +
|
||
-1.15606447E-10 * $D_Tmrt_5*$Pa +
|
||
-2.80626406E+00 * $Pa_2 +
|
||
5.48712484E-01 * $Ta*$Pa_2 +
|
||
-3.99428410E-03 * $Ta_2*$Pa_2 +
|
||
-9.54009191E-04 * $Ta_3*$Pa_2 +
|
||
1.93090978E-05 * $Ta_4*$Pa_2 +
|
||
-3.08806365E-01 * $va*$Pa_2 +
|
||
1.16952364E-02 * $Ta*$va*$Pa_2 +
|
||
4.95271903E-04 * $Ta_2*$va*$Pa_2 +
|
||
-1.90710882E-05 * $Ta_3*$va*$Pa_2 +
|
||
2.10787756E-03 * $va_2*$Pa_2 +
|
||
-6.98445738E-04 * $Ta*$va_2*$Pa_2 +
|
||
2.30109073E-05 * $Ta_2*$va_2*$Pa_2 +
|
||
4.17856590E-04 * $va_3*$Pa_2 +
|
||
-1.27043871E-05 * $Ta*$va_3*$Pa_2 +
|
||
-3.04620472E-06 * $va_4*$Pa_2 +
|
||
5.14507424E-02 * $D_Tmrt*$Pa_2 +
|
||
-4.32510997E-03 * $Ta*$D_Tmrt*$Pa_2 +
|
||
8.99281156E-05 * $Ta_2*$D_Tmrt*$Pa_2 +
|
||
-7.14663943E-07 * $Ta_3*$D_Tmrt*$Pa_2 +
|
||
-2.66016305E-04 * $va*$D_Tmrt*$Pa_2 +
|
||
2.63789586E-04 * $Ta*$va*$D_Tmrt*$Pa_2 +
|
||
-7.01199003E-06 * $Ta_2*$va*$D_Tmrt*$Pa_2 +
|
||
-1.06823306E-04 * $va_2*$D_Tmrt*$Pa_2 +
|
||
3.61341136E-06 * $Ta*$va_2*$D_Tmrt*$Pa_2 +
|
||
2.29748967E-07 * $va_3*$D_Tmrt*$Pa_2 +
|
||
3.04788893E-04 * $D_Tmrt_2*$Pa_2 +
|
||
-6.42070836E-05 * $Ta*$D_Tmrt_2*$Pa_2 +
|
||
1.16257971E-06 * $Ta_2*$D_Tmrt_2*$Pa_2 +
|
||
7.68023384E-06 * $va*$D_Tmrt_2*$Pa_2 +
|
||
-5.47446896E-07 * $Ta*$va*$D_Tmrt_2*$Pa_2 +
|
||
-3.59937910E-08 * $va_2*$D_Tmrt_2*$Pa_2 +
|
||
-4.36497725E-06 * $D_Tmrt_3*$Pa_2 +
|
||
1.68737969E-07 * $Ta*$D_Tmrt_3*$Pa_2 +
|
||
2.67489271E-08 * $va*$D_Tmrt_3*$Pa_2 +
|
||
3.23926897E-09 * $D_Tmrt_4*$Pa_2 +
|
||
-3.53874123E-02 * $Pa_3 +
|
||
-2.21201190E-01 * $Ta*$Pa_3 +
|
||
1.55126038E-02 * $Ta_2*$Pa_3 +
|
||
-2.63917279E-04 * $Ta_3*$Pa_3 +
|
||
4.53433455E-02 * $va*$Pa_3 +
|
||
-4.32943862E-03 * $Ta*$va*$Pa_3 +
|
||
1.45389826E-04 * $Ta_2*$va*$Pa_3 +
|
||
2.17508610E-04 * $va_2*$Pa_3 +
|
||
-6.66724702E-05 * $Ta*$va_2*$Pa_3 +
|
||
3.33217140E-05 * $va_3*$Pa_3 +
|
||
-2.26921615E-03 * $D_Tmrt*$Pa_3 +
|
||
3.80261982E-04 * $Ta*$D_Tmrt*$Pa_3 +
|
||
-5.45314314E-09 * $Ta_2*$D_Tmrt*$Pa_3 +
|
||
-7.96355448E-04 * $va*$D_Tmrt*$Pa_3 +
|
||
2.53458034E-05 * $Ta*$va*$D_Tmrt*$Pa_3 +
|
||
-6.31223658E-06 * $va_2*$D_Tmrt*$Pa_3 +
|
||
3.02122035E-04 * $D_Tmrt_2*$Pa_3 +
|
||
-4.77403547E-06 * $Ta*$D_Tmrt_2*$Pa_3 +
|
||
1.73825715E-06 * $va*$D_Tmrt_2*$Pa_3 +
|
||
-4.09087898E-07 * $D_Tmrt_3*$Pa_3 +
|
||
6.14155345E-01 * $Pa_4 +
|
||
-6.16755931E-02 * $Ta*$Pa_4 +
|
||
1.33374846E-03 * $Ta_2*$Pa_4 +
|
||
3.55375387E-03 * $va*$Pa_4 +
|
||
-5.13027851E-04 * $Ta*$va*$Pa_4 +
|
||
1.02449757E-04 * $va_2*$Pa_4 +
|
||
-1.48526421E-03 * $D_Tmrt*$Pa_4 +
|
||
-4.11469183E-05 * $Ta*$D_Tmrt*$Pa_4 +
|
||
-6.80434415E-06 * $va*$D_Tmrt*$Pa_4 +
|
||
-9.77675906E-06 * $D_Tmrt_2*$Pa_4 +
|
||
8.82773108E-02 * $Pa_5 +
|
||
-3.01859306E-03 * $Ta*$Pa_5 +
|
||
1.04452989E-03 * $va*$Pa_5 +
|
||
2.47090539E-04 * $D_Tmrt*$Pa_5 +
|
||
1.48348065E-03 * $Pa_5*$Pa;
|
||
return $UTCI_approx;
|
||
}
|
||
|
||
#~ calculates saturation vapour pressure over water in hPa for input air temperature (ta) in celsius according to:
|
||
#~ Hardy, R.; ITS-90 Formulations for Vapor Pressure, Frostpoint Temperature,
|
||
# Dewpoint Temperature and Enhancement Factors in the Range -100 to 100 °C;
|
||
#~ Proceedings of Third International Symposium on Humidity and Moisture;
|
||
# edited by National Physical Laboratory (NPL), London, 1998, pp. 214-221
|
||
#~ http://www.thunderscientific.com/tech_info/reflibrary/its90formulas.pdf (retrieved 2008-10-01)
|
||
sub es($)
|
||
{
|
||
my ($ta) = @_;
|
||
|
||
my @g = (-2.8365744E3, -6.028076559E3, 1.954263612E1, -2.737830188E-2, 1.6261698E-5, 7.0229056E-10,
|
||
-1.8680009E-13, 2.7150305);
|
||
|
||
my $tk = $ta + 273.15; # air temp in K
|
||
my $res = $g[7] * log($tk);
|
||
for (my $i = 0; $i < 7; $i++) {
|
||
$res += $g[$i] *$tk**($i-2);
|
||
}
|
||
$res=exp($res)*0.01; # *0.01: convert Pa to hPa
|
||
return $res;
|
||
}
|
||
|
||
# http://www.bom.gov.au/info/thermal_stress/#atapproximation
|
||
#
|
||
# [STDM 1994] : Steadman: Norms of apparent temperature in Australia
|
||
# http://reg.bom.gov.au/amm/docs/1994/steadman.pdf
|
||
#
|
||
sub apparent_temperature($$$$)
|
||
{
|
||
# ambient Temperature (C), humidity (%), windspeed in 10 m (m/s)
|
||
# effective radiation (W/m^2)
|
||
my ($Ta, $rh, $ws, $ER) = @_;
|
||
|
||
# vapour pressure
|
||
my $e = $rh / 100.0 * es($Ta);
|
||
|
||
# apparent temperature
|
||
my $AT = $Ta + 0.348 * $e -0.7 * $ws + 0.7 * $ER / ($ws + 10) - 4.25;
|
||
|
||
return $AT;
|
||
}
|
||
|
||
sub cosD($)
|
||
{
|
||
my ($phi) = @_;
|
||
return cos($phi * 0.0174532925);
|
||
}
|
||
|
||
sub sinD($)
|
||
{
|
||
my ($phi) = @_;
|
||
return sin($phi * 0.0174532925);
|
||
}
|
||
|
||
sub asinD($)
|
||
{
|
||
my ($x) = @_;
|
||
return POSIX::asin($x) / 0.0174532925;
|
||
}
|
||
|
||
|
||
sub acosD($)
|
||
{
|
||
my ($x) = @_;
|
||
return POSIX::acos($x) / 0.0174532925;
|
||
}
|
||
|
||
|
||
# tweakable constants
|
||
our @T_L = (3.4, 3.6, 4.3, 4.2, 4.7, 4.6, 4.7, 4.7, 4.3, 3.8, 3.4, 3.4); # Linke turbity for each month
|
||
our $f_svf_fl = 0.7; # sky view factor for computation of T_mrt
|
||
our $f_svf_sensor = 0.6; # sky view factor of sensor
|
||
our $bowen = 0.6; # bowen ratio (default = suburb)
|
||
|
||
# given constants
|
||
my $I_0 = 1350; # global insolation W/m²
|
||
my $sigma = 5.6E-8; # Stefan-Boltzmann
|
||
my $sw_abs = 0.72; # short wave absoption of body
|
||
my $lw_abs = 0.95; # long wave
|
||
|
||
##################################
|
||
# In the morning hours and or low altitudes of the sun the radiation sensor is obstructed by the rain gauge.
|
||
# Therefore the radiation curve is heavily skewed to the afternoon.
|
||
# Using several samplings af clear day radiation over the year a correction factor for the sensor was developed.
|
||
# The assumption is that the sensor always measures the diffuse radiation and that the direct or normal
|
||
# part is to be corrected.
|
||
#
|
||
# A polynom regression created with Python's sciphi.curve_fit is used here.
|
||
#
|
||
# I_h = poly(azimuth, sday) * (I_measured - D)
|
||
# sday is the distance to June 21st in days
|
||
#
|
||
my $azimuth_min = 80;
|
||
my $azimuth_max = 280;
|
||
my $alpha_min = 0.28;
|
||
my $alpha_max = 3.20;
|
||
my $nx = 7;
|
||
my $ny = 11;
|
||
my @poly = (-106.92078815088152,1.1559639672742617,-0.10641701193735834,
|
||
-0.003289004759141109,0.00024341610319517833,-7.393885292177471e-06,
|
||
1.1036469480601365e-07,-8.433653719722458e-10,4.6816105400180725e-12,
|
||
-2.2915699724588908e-14,5.1558807208283076e-17,-1.079309274936794e-21,
|
||
4.475568377940087,-0.02979654850540337,0.004700175671662923,
|
||
-2.0556721233089508e-05,-6.038914362044792e-08,2.006751148263329e-08,
|
||
-3.111327815305771e-10,-1.374995798484272e-12,1.9735294352791308e-14,
|
||
2.1301007272254039e-19,-1.1889565579254942e-22,-6.300423390960505e-22,
|
||
-0.07727338274859613,0.00020814642590810098,-6.394225596258378e-05,
|
||
7.868630252758761e-09,-7.464104467557937e-09,1.1747261979615144e-10,
|
||
6.942150907013742e-13,3.4379341050323286e-15,6.450994446627346e-18,
|
||
-4.444113964025086e-19,6.797271157288228e-22,-2.4187064817431376e-26,
|
||
0.0007233278231093484,-2.738439691372625e-07,5.245059294961144e-07,
|
||
3.0960830052172135e-09,-8.666235563108327e-12,-8.543981121106956e-13,
|
||
-3.0921133374358243e-16,-3.403768199365476e-18,2.0279200191426042e-21,
|
||
-4.248840692118507e-21,1.7318873922658786e-23,8.574899408718153e-28,
|
||
-3.97231902915211e-06,-5.352290412779365e-10,-3.192650063883535e-09,
|
||
-8.74015720283199e-12,2.283911280963708e-13,7.341513792108321e-16,
|
||
-1.646709199019733e-19,1.3337349967893863e-19,-1.768592696157379e-24,
|
||
1.8081165441415707e-23,-9.255858370952993e-28,-6.417530710544851e-29,
|
||
1.2816665080935244e-08,-1.1922083845200477e-11,1.3421147775930614e-11,
|
||
-8.568582018392473e-14,1.9981576527458447e-16,-3.2157555342893155e-22,
|
||
6.697464022398249e-20,-1.7492784034332204e-21,2.9982339842976173e-26,
|
||
3.1608055013268483e-26,-4.358123093981532e-28,-4.283467678330636e-31,
|
||
-2.2523634293967502e-11,7.160123195453301e-14,-3.16054616229181e-14,
|
||
4.4391029509107566e-16,-3.1617351024523765e-18,-1.106231164369411e-20,
|
||
1.3620891422547596e-26,4.542693150881559e-24,-1.1189308112788874e-26,
|
||
-2.7544251659486844e-29,1.8059972620664227e-36,6.576920263313157e-33,
|
||
1.6644569253648212e-14,-1.0182652809916555e-16,2.986157294663078e-17,
|
||
-5.245136086661749e-19,1.8660756935566013e-21,1.0544047345884695e-22,
|
||
-1.6410689680364017e-24,1.0506963847508829e-26,-5.70549250231417e-29,
|
||
3.5859454340809775e-32,2.027506749740463e-33,-1.4181289166467323e-35);
|
||
# numdp: 3936, mean error: 0.06565
|
||
|
||
# azimuth, sday
|
||
sub sensor_factor($$) {
|
||
my ($x, $y) = @_;
|
||
|
||
# precompute powers
|
||
my @yp;
|
||
|
||
my $t = 1.0;
|
||
for my $i (0 .. $ny) {
|
||
push(@yp, $t);
|
||
$t *= $y;
|
||
}
|
||
|
||
my $f = 0.0;
|
||
my $xp = 1.0;
|
||
for my $i (0 .. $nx) {
|
||
for my $j (0 .. $ny) {
|
||
$f += $poly[$i * ($ny + 1) +$j] * $xp * $yp[$j];
|
||
}
|
||
$xp *= $x;
|
||
}
|
||
|
||
return $f;
|
||
}
|
||
|
||
|
||
##########################
|
||
sub Julian_Date($$$$)
|
||
{
|
||
my ($yr, $mn, $mday, $hr) = @_;
|
||
|
||
# http://aa.usno.navy.mil/faq/docs/JD_Formula.php
|
||
my $JD = 367 * $yr - int(7 * ($yr + int(($mn+9)/12)) / 4) + int(275 * $mn/9) + $mday + 1721013.5 + $hr / 24;
|
||
}
|
||
|
||
# get (azimuth, elevation) of sun
|
||
# https://de.wikipedia.org/wiki/Sonnenstand#Astronomische_Zusammenh%C3%A4nge
|
||
sub sun_position($$$)
|
||
{
|
||
my ($LAT, $LON, $TS) = @_;
|
||
my ($sec, $min, $hr, $mday, $mn, $yr) = gmtime($TS);
|
||
|
||
$yr += 1900;
|
||
$mn += 1;
|
||
$hr += $min / 60;
|
||
|
||
my $JD = Julian_Date($yr, $mn, $mday, $hr);
|
||
my $n = $JD - 2451545.0;
|
||
|
||
my $L = fmod(280.460 + 0.9856474 * $n, 360);
|
||
my $g = fmod(357.528 + 0.9856003 * $n, 360);
|
||
my $Lambda = fmod($L + 1.915 * sinD($g) + 0.01997 * sinD(2 * $g), 360);
|
||
|
||
my $eps = 23.439 - 0.0000004 * $n;
|
||
|
||
my $sL = sinD($Lambda);
|
||
my $cL = cosD($Lambda);
|
||
my $alpha = rad2deg(atan(cosD($eps) * $sL / $cL));
|
||
|
||
if ($cL < 0) {
|
||
$alpha += 180;
|
||
}
|
||
|
||
my $delta = asinD(sinD($eps) * sinD($Lambda));
|
||
|
||
my $T0 = (Julian_Date($yr, $mn, $mday, 0) - 2451545.0) / 36525;
|
||
|
||
my $Theta_hG = fmod(6.697376 + 2400.05134 * $T0 + 1.002738 * $hr, 24);
|
||
my $Theta = fmod($Theta_hG * 15 + $LON, 360);
|
||
my $tau = $Theta - $alpha;
|
||
|
||
my $sd = sinD($delta);
|
||
my $cd = cosD($delta);
|
||
my $X = cosD($tau) * sinD($LAT) - $sd / $cd * cosD($LAT);
|
||
my $a = rad2deg(atan(sinD($tau) / $X));
|
||
$a += 180 if ($X < 0);
|
||
$a = fmod($a + 180, 360);
|
||
my $h = asinD($cd * cosD($tau) * cosD($LAT) + $sd * sinD($LAT));
|
||
|
||
my $R = 1.02 / tan(deg2rad($h + 10.3 / ($h + 5.11)));
|
||
my $hR = $h + $R / 60;
|
||
|
||
Log(5, "JD: $JD, n: $n, L: $L, g: $g, Lambda: $Lambda, Eps: $eps, alpha: $alpha, delta: $delta\n");
|
||
Log(5, "T0: $T0, Theta_hG: $Theta_hG, Theta: $Theta, a: $a, h: $h, h: $hR\n");
|
||
return ($a, $h);
|
||
}
|
||
|
||
# compute radiation components for sun altitude, pressure, and sky view factor
|
||
sub radiation($$$$)
|
||
{
|
||
my ($alt, $pabs, $f_svf, $f_direct) = @_;
|
||
|
||
my ($month, $day_of_year) = (gmtime(time()))[4, 7];
|
||
my $t_l = $T_L[$month];
|
||
|
||
# Steadman 1994, eq. (18)
|
||
my $I_0_d = ($I_0 - 46 * sinD(($day_of_year - 94) / 365 * 360));
|
||
|
||
# Matzerakis 2010, eq. (3) .. (10)
|
||
my $z = 90 - $alt;
|
||
my $f = exp(-0.027 * $pabs / 1013.2 * $t_l / cosD($z));
|
||
my $G_0 = 0.84 * $I_0_d * cosD($z) * $f;
|
||
|
||
my $m_R0 = 1.0 / (sinD($alt) + 0.50572 * ($alt + 6.07995)**-1.6364);
|
||
my $delta_R0 = 1.0 / (0.9 * $m_R0 + 9.4);
|
||
my $tau = exp(-$t_l * $delta_R0 * $m_R0 * $pabs / 1013.2);
|
||
|
||
my $I_d = $I_0 * $tau;
|
||
my $I_h = $I_d * cosD($z);
|
||
|
||
my $x = ($G_0 - $I_h);
|
||
|
||
# anisotropic radiation is =! 0 only if the sun is visible
|
||
# so multiply with $f_direct
|
||
my $D_aniso = $x * $tau * $f_direct;
|
||
my $D_iso = $x * (1.0 - $tau) * $f_svf;
|
||
|
||
my $D_0 = $D_aniso + $D_iso;
|
||
my $D_8 = 0.28 * $G_0 * $f_svf;
|
||
|
||
Log(4, sprintf("D_aniso: %4.1f, D_iso: %4.1f, D_0: %4.1f, D_8: %4.1f", $D_aniso, $D_iso, $D_0, $D_8));
|
||
return ($G_0, $I_h, $I_d, $D_0, $D_8);
|
||
}
|
||
|
||
# long wave surface radiation
|
||
sub L_surface($$$$)
|
||
{
|
||
my ($Ta, $va, $G, $A) = @_;
|
||
my $eps = 0.8;
|
||
|
||
# Matzerakis 2010, eq. (12), (13)
|
||
|
||
# iterate for surface temperature, (converges rapidly)
|
||
my $Ts = $Ta;
|
||
foreach my $i (1..10) {
|
||
my $Ts_o = $Ts;
|
||
my $E = $eps * $sigma * ($Ts + 273)**4 + (1 - $eps) * $A;
|
||
my $Q = $sw_abs * $G + $A - $E;
|
||
my $B = ($Q >= 0) ? -0.19 * $Q : -0.32 * $Q;
|
||
$Ts = $Ta + ($Q + $B) / ((6.2 + 4.2 * $va) * (1 + 1 / $bowen));
|
||
Log(5, sprintf("E = %4.1f, Q = %4.1f, B = %4.1f, Ts = %5.2f", $E, $Q, $B, $Ts));
|
||
|
||
last if abs($Ts_o - $Ts) < 0.05;
|
||
}
|
||
|
||
# Radiation from surface
|
||
my $E = $eps * $sigma * ($Ts + 273)**4 + (1 - $eps) * $A;
|
||
return ($E, $Ts);
|
||
}
|
||
|
||
# compute T_tmrt, mean radiant temperature
|
||
sub T_RMT($$$$$$)
|
||
{
|
||
my ($I_d_n, $D, $LW, $E, $alt, $f_direct) = @_;
|
||
|
||
# Soil reflectance
|
||
my $Rsoil = 0.2;
|
||
|
||
# projection factor for cylindrical body
|
||
my $fp = 0.308 * cosD($alt * (1 - $alt * $alt / 48402));
|
||
|
||
# upper hemisphere
|
||
my $R = 0.5 * ($LW + $sw_abs / $lw_abs * $D);
|
||
|
||
# lower hemisphere
|
||
$R += 0.5 * $sw_abs / $lw_abs * ($I_d_n * sinD($alt) + $D) * $Rsoil;
|
||
$R += 0.5 * $E;
|
||
|
||
my $T_mrt_ambient = ($R / $sigma)**0.25 - 273;
|
||
|
||
# direct
|
||
my $R_d = $sw_abs / $lw_abs * $fp * $I_d_n;
|
||
|
||
my $T_mrt_full = (($R + $R_d)/ $sigma)**0.25 - 273;
|
||
my $T_mrt = (($f_direct * $R_d + $R) / $sigma)**0.25 - 273;
|
||
Log(4, sprintf("T_mrt_ambient: %4.1f, T_mrt_full: %4.1f, T_mrt: %4.1f", $T_mrt_ambient, $T_mrt_full, $T_mrt));
|
||
|
||
return $T_mrt;
|
||
}
|
||
|
||
# some defaults for stand alone testing
|
||
our $LAT = 50.138804;
|
||
our $LON = 8.501993;
|
||
our $altitude = 146;
|
||
|
||
# should be callable outside of FHEM for development purposes
|
||
sub _mean_radiant_temperature($$$$$$$$$$)
|
||
{
|
||
my ($azimuth, $alt, $Ta, $Hr, $va, $pabs, $f_direct, $Isensor, $coverage, $T_sky) = @_;
|
||
|
||
my $ehPa = es($Ta) * $Hr * 0.01;
|
||
Log(4, sprintf("Isensor: %3.1f, Pa %3.1f", $Isensor, $ehPa));
|
||
Log(4, sprintf("Elevation: %3.1f, Azimuth: %3.1f", $alt, $azimuth));
|
||
|
||
my ($G_0, $I_d, $I_h, $D_0, $D_8);
|
||
|
||
if ($alt >= 5) {
|
||
# get radiation values with f_svf_fl for 'feels like' temperature
|
||
($G_0, $I_h, $I_d, $D_0, $D_8) = radiation($alt, $pabs, $f_svf_fl, $f_direct);
|
||
|
||
} elsif (0 < $alt && $alt < 5) {
|
||
# twilight mode
|
||
$I_d = $I_h = 0.0;
|
||
$D_0 = $D_8 = $Isensor;
|
||
$coverage = 0.5;
|
||
} else {
|
||
# night mode
|
||
$I_d = $I_h = $D_0 = $D_8 = 0.0;
|
||
$coverage = 0.75;
|
||
}
|
||
|
||
Log(4, sprintf("Coverage used for computation: %3.2f", $coverage));
|
||
|
||
# compute coverage dependent values
|
||
|
||
# direct radiation, normal
|
||
my $I_d_n = $I_d * (1 - $coverage);
|
||
|
||
# direct radiation, horizontal
|
||
my $I_h_n = $I_h * (1 - $coverage);
|
||
|
||
# diffuse radiation, Matzerakis 2010, eq. (7)
|
||
my $D = $D_0 + ($D_8 - $D_0) * $coverage;
|
||
|
||
# for very thick clouds everything measured is diffuse
|
||
$D = $Isensor if ($D > $Isensor);
|
||
|
||
my $I_lw;
|
||
if (defined($T_sky)) {
|
||
$I_lw = $sigma * ($T_sky + 273)**4;
|
||
Log(4, sprintf("Direct_n: %4.0f, Diffuse: %4.0f Lw: %4.0f T_sky (IR): %4.1f", $I_d_n, $D, $I_lw, $T_sky));
|
||
} else {
|
||
# air long wave radiation, Matzerakis 2010, eq. (11)
|
||
$I_lw = $sigma * (273 + $Ta)**4 * (0.82 -0.25 * 10**(-0.0945 * $ehPa)) * (1 + 0.21 * $coverage**2.5);
|
||
|
||
# sky temperature just for reference
|
||
my $T_s = ($I_lw / $sigma)**0.25 - 273;
|
||
|
||
Log(4, sprintf("Direct_n: %4.0f, Diffuse: %4.0f Lw: %4.0f T_sky (comp): %4.1f", $I_d_n, $D, $I_lw, $T_s));
|
||
}
|
||
|
||
# Shouldn't that be multiplied by sky_view ?
|
||
|
||
# surface long wave radiation
|
||
my ($E, $Tsoil) = L_surface($Ta, $va, $I_h_n + $D, $I_lw);
|
||
Log(4, sprintf("Tsoil = %4.1f, E = %4.1f", $Tsoil, $E));
|
||
|
||
my $T_mrt = T_RMT($I_d_n, $D, $I_lw, $E, $alt, $f_direct);
|
||
|
||
# ER = 'absorbed extra radiation' for apparent temperature
|
||
# 0.725 = effective radiating part of body, 6 = radiating heat xfer coefficient
|
||
my $ER = 0.725 * 6 * ($T_mrt - $Ta);
|
||
return ($T_mrt, $ER);
|
||
}
|
||
|
||
# compute coverage by analysis of Isensor's time series
|
||
# General method:
|
||
# - compute (correction factor adjusted) ratio of actual radiation to full radiation
|
||
# called 'normalized radiation, nr_..' ideally in (0,1)
|
||
# - store nr_.. values in array
|
||
# - count points in for bands of interval (0,1)
|
||
# - derive coverage by from relative ratios of bands
|
||
# - allow adjustment using e.g. an IR thermometer by calling an optional callback
|
||
#
|
||
sub compute_coverage($$$$$$)
|
||
{
|
||
my ($hash, $az, $alt, $pabs, $Isensor, $cb) = @_;
|
||
|
||
my $max_nr_history = 225; # 60 minutes
|
||
|
||
# create the device's nr_history array
|
||
if (!exists($hash->{helper}{nr_history})) {
|
||
$hash->{helper}{nr_history} = [];
|
||
}
|
||
|
||
my $nr_history = $hash->{helper}{nr_history};
|
||
|
||
my ($alpha, $nr_coverage, $nr_stddev, $nr_band_0, $nr_band_1, $nr_band_2, $nr_band_3);
|
||
|
||
my $name = $hash->{NAME};
|
||
|
||
my $sr_invalid = 0;
|
||
my $sun_visible = 0;
|
||
$hash->{helper}{nr_current} = undef;
|
||
|
||
# check whether sun is visible
|
||
my @sv = split(/, */, main::AttrVal($name, "sunVisibility", "0,5,360"));
|
||
my $az_low = shift(@sv);
|
||
|
||
while (@sv > 1) {
|
||
my $alt_min = shift(@sv);
|
||
my $az_high = shift(@sv);
|
||
if ($az_low <= $az && $az < $az_high) {
|
||
$sun_visible = $alt >= $alt_min ? 1 : 0;
|
||
last;
|
||
}
|
||
$az_low = $az_high;
|
||
}
|
||
|
||
Log(4, "sun_visible: $sun_visible");
|
||
|
||
if ($alt >= 5 && $sun_visible) {
|
||
Log(4, sprintf("compute_coverage: Isensor: %3.1f, pabs: %3.1f", $Isensor, $pabs));
|
||
Log(4, sprintf("az: %3.1f, alt: %3.1f", $az, $alt));
|
||
|
||
# get radiation with f_svf for sensor
|
||
my ($G_0, $I_h, $I_d, $D_0, $D_8) = radiation($alt, $pabs, $f_svf_sensor, 1.0);
|
||
Log(4, sprintf("G_0: %4.0f, I_d: %5.1f, I_h: %5.1f, D_0: %4.0f, D_8: %4.0f",
|
||
$G_0, $I_d, $I_h, $D_0, $D_8));
|
||
|
||
my $sday; # distance from June 21 in days
|
||
my $day_of_year = (gmtime)[7];
|
||
if ($day_of_year > 344) {
|
||
$sday = (365 - $day_of_year) - 172;
|
||
} else {
|
||
$sday = $day_of_year - 172;
|
||
}
|
||
|
||
$sday = -$sday if ($sday < 0);
|
||
|
||
my $sensor_type = main::AttrVal($name, "sensorType", "wh2601");
|
||
|
||
if ($sensor_type eq "wh2601") {
|
||
# apply sensor specific corrections
|
||
$alpha = sensor_factor($az, $sday);
|
||
Log(4, sprintf("sday: $sday, alpha: %3.2f", $alpha));
|
||
|
||
# likely to be 100% wrong if out of range
|
||
my $amin = main::AttrVal($name, "alphaMin", $alpha_min);
|
||
my $amax = main::AttrVal($name, "alphaMax", $alpha_max);
|
||
$amin = $amin > $alpha_min ? $amin : $alpha_min;
|
||
$amax = $amax < $alpha_max ? $amax : $alpha_max;
|
||
|
||
if ($alpha > $amax || $alpha < $amin
|
||
|| $az < $azimuth_min || $az > $azimuth_max) {
|
||
$sr_invalid = 1;
|
||
Log(4, "Alpha or sun's position is out of range");
|
||
}
|
||
} else {
|
||
$alpha = 1.0;
|
||
$sr_invalid = 0;
|
||
}
|
||
|
||
if (! $sr_invalid) {
|
||
# use last coverage computed as estimation
|
||
my $last_coverage = $hash->{helper}{coverage};
|
||
$last_coverage = defined($last_coverage) ? $last_coverage : 0.5;
|
||
my $D = (1 - $last_coverage) * $D_0 + $last_coverage * $D_8;
|
||
my $I_h_m = $alpha * ($Isensor - $D);
|
||
my $x = $G_0 - $D;
|
||
my $nr = $I_h_m / $x;
|
||
$hash->{helper}{nr_current} = $nr;
|
||
Log(4, sprintf("I_h_m: %3.1f, Limit: %3.1f, nr: %3.2f", $I_h_m, $x, $nr));
|
||
|
||
my $n_nr_history = push(@$nr_history, $nr);
|
||
|
||
while ($n_nr_history > $max_nr_history) {
|
||
shift(@$nr_history);
|
||
$n_nr_history--;
|
||
}
|
||
|
||
# at least ~ 10 minutes of data
|
||
if ($n_nr_history > 40) {
|
||
$nr_band_0 = $nr_band_1 = $nr_band_2 = $nr_band_3 = 0;
|
||
foreach my $nr (@$nr_history) {
|
||
if ($nr > 0.8) {
|
||
$nr_band_0++;
|
||
} elsif (0.6 < $nr) {
|
||
$nr_band_1++;
|
||
} elsif (0.2 < $nr) {
|
||
$nr_band_2++;
|
||
} else {
|
||
$nr_band_3++;
|
||
}
|
||
}
|
||
|
||
$nr_band_0 /= $n_nr_history;
|
||
$nr_band_1 /= $n_nr_history;
|
||
$nr_band_2 /= $n_nr_history;
|
||
$nr_band_3 /= $n_nr_history;
|
||
|
||
# if 100% clear force coverage to be rounded down for oktas
|
||
if ($nr_band_0 == 1) {
|
||
$nr_coverage = 0.05;
|
||
} else {
|
||
# meaning:
|
||
# for alternate high and med it's statistical
|
||
# for 100% low it will be rounded up to OVC
|
||
$nr_coverage = $nr_band_0 * 0.5/8 + $nr_band_1 * 3.5/8
|
||
+ $nr_band_2 * 6.5/8 + $nr_band_3 * 7.5/8;
|
||
}
|
||
|
||
$hash->{helper}{nr_coverage_raw} = $nr_coverage;
|
||
|
||
# experimental stuff for low coverage
|
||
# around 30 minutes
|
||
my $n_sample = 120;
|
||
if ($n_nr_history > $n_sample) {
|
||
my $sum = 0;
|
||
my $sum2 = 0;
|
||
foreach my $i (($n_nr_history-$n_sample)..($n_nr_history-1)) {
|
||
my $x = $nr_history->[$i];
|
||
$sum += $x;
|
||
$sum2 += $x * $x;
|
||
}
|
||
|
||
$sum /= $n_sample;
|
||
$sum2 /= $n_sample;
|
||
|
||
$nr_stddev = sqrt($sum2 - $sum * $sum);
|
||
|
||
if ($nr_coverage <= 4/8 && $nr_stddev > 0.2) {
|
||
$nr_coverage = 4/8;
|
||
} elsif ($nr_coverage < 1/8) {
|
||
if ($nr_stddev >= 0.1) {
|
||
$nr_coverage = 2/8;
|
||
} elsif ($nr_stddev > 0.015) {
|
||
$nr_coverage = 1/8;
|
||
}
|
||
}
|
||
}
|
||
}
|
||
}
|
||
}
|
||
|
||
# reset history on obstruction
|
||
if (!$sun_visible || $sr_invalid) {
|
||
@{$nr_history} = ();
|
||
$nr_coverage = undef;
|
||
}
|
||
|
||
my $T_sky;
|
||
|
||
# call the coverage callback
|
||
my $coverage = $nr_coverage;
|
||
if (defined($cb)) {
|
||
Log(4, "cb: calling with coverage: " . (defined($coverage) ? $coverage : 'undef'));
|
||
no strict "refs";
|
||
eval {
|
||
($coverage, $T_sky) = &{"main::$cb"}($hash, $az, $alt, $coverage);
|
||
};
|
||
|
||
if ($@) {
|
||
Log(0, "Eval of $cb failed: $@");
|
||
}
|
||
|
||
use strict "refs";
|
||
|
||
Log(4, "cb: return coverage: " . (defined($coverage) ? $coverage : 'undef')
|
||
. ", T_Sky: " . ($T_sky || 'undef'));
|
||
}
|
||
|
||
# store in helper so we can inspect it with list or create readings
|
||
$hash->{helper}{nr_band_0} = $nr_band_0;
|
||
$hash->{helper}{nr_band_1} = $nr_band_1;
|
||
$hash->{helper}{nr_band_2} = $nr_band_2;
|
||
$hash->{helper}{nr_band_3} = $nr_band_3;
|
||
|
||
$hash->{helper}{sr_invalid} = $sr_invalid;
|
||
$hash->{helper}{sun_visible} = $sun_visible;
|
||
$hash->{helper}{alpha} = $alpha;
|
||
$hash->{helper}{nr_stddev} = $nr_stddev;
|
||
$hash->{helper}{nr_coverage} = $nr_coverage;
|
||
|
||
$hash->{helper}{coverage} = $coverage;
|
||
|
||
$coverage = defined($coverage) ? $coverage : 0.5;
|
||
return ($coverage, $T_sky);
|
||
}
|
||
|
||
##########################################################################################################
|
||
# functions below are "public"
|
||
package main;
|
||
|
||
sub mean_radiant_temperature($$$$$$)
|
||
{
|
||
my ($devname, $Ta, $Hr, $va, $p, $Isensor) = @_;
|
||
my $hash = $main::defs{$devname};
|
||
if (!defined($hash)) {
|
||
Log(1, 'feels_like: $devname is invalid');
|
||
return undef;
|
||
}
|
||
|
||
$feels_like::LAT = AttrVal("global", "latitude", undef);
|
||
$feels_like::LON = AttrVal("global", "longitude", undef);
|
||
$feels_like::altitude = AttrVal("global", "altitude", undef);
|
||
|
||
Log(1, 'feels_like needs global attributes "latitude", "longitude", and "altitude"')
|
||
unless (defined($LAT) && defined($LON) && defined($altitude));
|
||
|
||
# set tweakable constants
|
||
my $x = AttrVal($devname, "sensorSkyViewFactor", undef);
|
||
$feels_like::f_svf_sensor = $x if (defined($x));
|
||
|
||
$x = AttrVal($devname, "skyViewFactor", undef);
|
||
$feels_like::f_svf_fl = $x if (defined($x));
|
||
|
||
$x = AttrVal($devname, "bowenRatio", undef);
|
||
$feels_like::bowen = $x if (defined($x));
|
||
|
||
$x = AttrVal($devname, "linkeTurbity", undef);
|
||
if (defined($x)) {
|
||
my @y;
|
||
$x = "\@y = $x";
|
||
eval($x);
|
||
if ($@) {
|
||
Log(1, "Invalid syntax for linkeTurbity -> $@");
|
||
} else {
|
||
@feels_like::T_L = @y;
|
||
}
|
||
}
|
||
|
||
# part of direct radiation that should be used for T_mrt
|
||
my $f_direct = AttrVal($devname, "directRadiationRatio", 0.4);
|
||
|
||
my ($az, $alt) = feels_like::sun_position($LAT, $LON, time());
|
||
$alt = 0 if ($alt < 0);
|
||
|
||
readingsBulkUpdate($hash, "azimuth", round($az, 1));
|
||
readingsBulkUpdate($hash, "altitude", round($alt, 1));
|
||
|
||
my $pabs = $p - $altitude / 8;
|
||
my $cb = AttrVal($devname, "coverageCallback", undef);
|
||
|
||
my ($coverage, $T_sky) = feels_like::compute_coverage($hash, $az, $alt, $pabs, $Isensor, $cb);
|
||
|
||
my @res = feels_like::_mean_radiant_temperature($az, $alt, $Ta, $Hr, $va, $p,
|
||
$f_direct, $Isensor, $coverage, $T_sky);
|
||
return (@res);
|
||
}
|
||
|
||
##########################
|
||
sub T_UTCI($$$$$)
|
||
{
|
||
my ($Ta, $Hr, $va, $p, $T_mrt) = @_;
|
||
|
||
my $ehPa = feels_like::es($Ta) * $Hr * 0.01;
|
||
|
||
# UTCI seems to heavily overestimate for very low wind speeds ($va >= 0.5 is already implied in UTCI_approx)
|
||
if ($va < 2.0) {
|
||
$va = 1.0 + 0.5 * $va;
|
||
}
|
||
|
||
return feels_like::UTCI_approx($Ta, $ehPa, $T_mrt, $va);
|
||
}
|
||
|
||
sub T_apparent_temperature($$$$)
|
||
{
|
||
my ($Ta, $Hr, $va, $ER) = @_;
|
||
|
||
return feels_like::apparent_temperature($Ta, $Hr, $va, $ER);
|
||
}
|
||
|
||
1;
|
||
|
||
=pod
|
||
=item helper
|
||
=item summary compute a 'feels like' temperature according to UTCI or AT and cloud coverage
|
||
=begin html
|
||
|
||
<a name="feels_like"></a>
|
||
<h3>feels_like</h3>
|
||
|
||
<ul>
|
||
Computation of a 'feels like' temperature based on ambient temperature, humidity, wind speed and solar radiation.<br>
|
||
While for a simple computation temperature, humidity, wind speed are sufficient this module is specifically tailored
|
||
to support stations with the Fine Offset sensor array marketed under several names like<br>
|
||
<ul>
|
||
<br>
|
||
<li>Ambient Weather WS-1200 / IP-1400 / ...</li>
|
||
<li>Froggit WS2601</li>
|
||
<li>Renkforce WH2600</li>
|
||
<br>
|
||
</ul>
|
||
Most of them are supported by the module <a href="commandref.html#HP1000">HP1000</a>.<br>
|
||
<br>
|
||
Algorithms providing a 'feels like' temperature work the following way:
|
||
<ul>
|
||
<li>Compute the <i>net extra radiation</i> on the human body using environmental data, the position of the sun and the amount of cloud cover<br>
|
||
resulting in the <a href = "https://en.wikipedia.org/wiki/Mean_radiant_temperature">Mean Radiant Temperature</a></li>
|
||
<li>Make a complete heat balance of the human body including this extra radiation.</li>
|
||
<li>Compare this to the heat balance of a walking person in the shadow with no wind and average humidity.</li>
|
||
<li>The temperature where these two match defines the <i>feels like</i>, <i>equivalent</i>, or <i>apparent</i> temperature.</li>
|
||
</ul>
|
||
</ul>
|
||
|
||
<ul>
|
||
This implementation computes the short and long wave radiation parts based on the position of the sun.<br>
|
||
Taking into account the <i>solar radiation reading</i> the amount of cloud cover is computed in <a href="https://en.wikipedia.org/wiki/Okta">Oktas</a>.<br>
|
||
Using this a <i>Mean Radiant Temperature</i> is computed that is consistent with the provided <i>solar radiation reading</i>.<br><br>
|
||
|
||
Currently this module provides two different <i>equivalent temperatures</i>:<br>
|
||
<ul>
|
||
<li>
|
||
<b><a href="http://www.utci.org">Universal Temperature Climate Index</a></b><br>
|
||
Currently considered as the gold standard. It was created 2007 by a joint effort of world wide experts.
|
||
See the pdf <a href="http://www.utci.org/utci_poster.pdf">Poster</a> for a quick introduction.<br>
|
||
</li>
|
||
|
||
<li>
|
||
<b><a href="http://www.bom.gov.au/info/thermal_stress/#apparent">Apparent Temperature</a></b><br>
|
||
An older standard developed by R. B. Steadman in the 1970s to 1990s and still in use in Australia and USA.<br>
|
||
</li>
|
||
|
||
</ul>
|
||
</ul>
|
||
|
||
<a name="feels_likedefine"></a>
|
||
<b>Define</b>
|
||
<ul>
|
||
<code>define <name> feels_like <out-device> <temperature> <humidity> [<wind_speed> [<solarradiation> [<pressure>]]]</code><br>
|
||
<br>
|
||
<ul>
|
||
Calculates a 'feels like' temperatures for specified readings and writes them into new readings for device <out-device>.<br>
|
||
Input readings can be specified as <device>:<reading>. At least one input reading must belong to <out-device>.<br><br>
|
||
|
||
<i>Input Readings</i>
|
||
<table>
|
||
<tr><td>temperature</td> <td>ambient temperature [°C]</td></tr>
|
||
<tr><td>humidity</td> <td>ambient relative humidity [%]</td></tr>
|
||
<tr><td>wind_speed</td> <td>speed of wind 10 m above ground (so be sure to calibrate accordingly) [m/s]</td></tr>
|
||
<tr><td>solarradiation </td> <td>solar radiation [W/m^2]</td></tr>
|
||
<tr><td>pressure</td> <td>pressure [hPa]</td></tr>
|
||
</table><br><br>
|
||
|
||
<i>Generated Readings</i>
|
||
<table>
|
||
<tr><td>temperature_utci</td> <td>temperature according to UTCI [°C]</td></tr>
|
||
<tr><td>temperature_at</td> <td>temperature according to Steadman [°C]</td></tr>
|
||
<tr><td>temperature_mrt</td> <td>mean radiant temperature [°C]</td></tr>
|
||
<tr><td>extra_radiation</td> <td>extra radiation on human body [W/m^2]</td></tr>
|
||
</table><br><br>
|
||
|
||
If <solarradiation> is specified as input parameter these readings are generated in addition:<br><br>
|
||
<table>
|
||
<tr><td>coverage</td> <td>cloud cover as value 0.0 -> 1.0</td></tr>
|
||
<tr><td>oktas</td> <td>cloud cover in oktas 1 -> 8</td></tr>
|
||
<tr><td>clouds</td> <td>letter code for coverage SKC -> OVC</td></tr>
|
||
</table><br>
|
||
The value for <i>clouds</i> is compatible with Wunderground and can be uploaded with field code "clouds".<br><br>
|
||
<b>Note:</b>For cloud detection it is essential that an event for <solarradiation> is generated for each transmission
|
||
of the weather station, i.e. each 16 seconds.
|
||
</ul>
|
||
<br>
|
||
|
||
Example:<PRE>
|
||
# Use full functionality for a wh2601 weather station
|
||
define wh2601_fl feels_like wh2601 temperature humidity wind_speed solarradiation pressure
|
||
|
||
# Use different sensors
|
||
define Thermo_1_fl feels_like Thermo_1 temperature Thermo_2:humidity WS_Sensor:wind_speed
|
||
</PRE>
|
||
</ul>
|
||
<ul>
|
||
<br>
|
||
<b>Attributes</b><br>
|
||
<table>
|
||
<tr><td>latitude, longitude, altitude </td> <td>required</td> <td>must be defined in <code>global</code></td></tr>
|
||
<tr><td>maxTimediff</td> <td>optional, should match installation </td> <td>Maximum allowed time difference in seconds for readings of different sensors to be considered coherent (default: 600)</td></tr>
|
||
<tr><td>sensorSkyViewFactor</td> <td>optional, should match installation </td> <td>Fraction of sky that is visible to the sensor (default: 0.6)</td></tr>
|
||
<tr><td>skyViewFactor</td> <td>optional, your personal choice</td> <td>Fraction of visible sky to be used for computation of T_mrt (default: 0.7)</td></tr>
|
||
<tr><td>directRadiationRatio</td> <td>optional, your personal choice</td> <td>Fraction of direct insolation to be used for computation of T_mrt (default 0.4)</td></tr>
|
||
<tr><td>utciWindCutoff</td> <td>optional, your personal choice</td> <td>Max wind speed for UTCI computation (default: 3 m/s)</td></tr>
|
||
<tr><td>bowenRatio</td> <td>optional, technical</td> <td><a href="https://en.wikipedia.org/wiki/Bowen_ratio">Bowen ratio</a> of environment (default 0.6 = grassland)</td></tr>
|
||
<tr><td>linkeTurbity</td> <td>optional, technical</td> <td><a href="http://www.soda-pro.com/help/general-knowledge/linke-turbidity-factor">Linke Turbity</a> of atmosphere.<br>
|
||
Value should be an array of 12 values ( = one value / month) (default: values for Frankfurt)<br>
|
||
You can find values for other locations <a href="http://www.soda-pro.com/web-services/atmosphere/turbidity-linke-2003">here</a></td></tr>
|
||
<tr><td>sunVisibility</td> <td>optional, should match installation</td> <td>Comma separated list of azimuth/altitude values to describe whether the sun is visible (e.g. due to buildings, trees, ...).<br>Must start with 0 and end with 360.<br>
|
||
E.g. 0,5,75,20.5,130,10,360 -> sun must be > 5° between 0 and 75° azimuth, then > 20.5° between 75° and 130° etc.</td></tr>
|
||
<tr><td>coverageCallback</td> <td>optional, technical</td> <td>Callback function for use with an additional infrared thermometer</td></tr>
|
||
<tr><td>sensorType</td> <td>optional, technical</td> <td>Type of sensor for applying specific corrections, values:wh2601,generic (default: wh2601)</td></tr>
|
||
</table>
|
||
|
||
<br><br>
|
||
If your environment is a<br>
|
||
large paved parking lot use <code>skyViewFactor = 1.0, directRadiationRatio = 1.0, bowenRatio = 5</code>,<br>
|
||
shadowy wetland use <code>skyViewFactor = 0.2, directRadiationRatio = 0.2, bowenRatio = 0.2</code>,<br>
|
||
suburb, golf course... use the defaults.
|
||
</ul>
|
||
|
||
<br>
|
||
<ul>
|
||
<b>Bibliography</b><br>
|
||
P. Bröde et al. (2010): The Universal Thermal Climate Index UTCI in Operational Use<br>
|
||
Proceedings of Conference: Adapting to Change: New Thinking on Comfort Cumberland
|
||
Lodge, Windsor, UK, 9-11 April 2010
|
||
<a href="http://www.utci.org/isb/documents/windsor_vers05.pdf">download</a>
|
||
<br><br>
|
||
Andreas Matzarakis, Frank Rutz, Helmut Mayer (2010): Modelling radiation fluxes in simple and complex
|
||
environments: basics of the RayMan model<br>
|
||
Int J Biometeorol (2010) 54:131–139 <a href="http://www.urbanclimate.net/rayman/rayman12.zip">download</a>
|
||
<br><br>
|
||
Robert G. Steadman. (1994): Norms of apparent temperature in Australia.<br>
|
||
Aust. Met. Mag., Vol 43, 1-16. <a href="http://reg.bom.gov.au/amm/docs/1994/steadman.pdf">download</a>
|
||
<br><br>
|
||
John L. Monteith† and Mike H. Unsworth: Principles of Environmental Physics, Fourth Edition<br>
|
||
Elsevier, Kidlington, UK, 2013
|
||
</ul>
|
||
|
||
=end html
|
||
=cut
|