# $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 $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';

        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 $hash->{$dk} ne $dev_name;

            my $rk = $k . '_READING';
            next if $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, 2) : '');

        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 = 10;
my @poly = (-80.23280867846618,6.145065782244877,-0.7313678902154087,
         0.03250849135687465,-0.0008978201231421577,1.2347982929985227e-05,
         -6.695047898526558e-08,-6.925834867866924e-13,1.118146677702321e-12,
         -3.761823314379515e-15,8.31708341485985e-18,2.889991801034238,
         -0.021047044910027124,0.0036703670873613507,3.1069176573053456e-07,
         2.6470054821477945e-06,-5.508774210575208e-08,-1.2439092637735542e-11,
         3.7634863724249864e-12,-1.3255609900004629e-14,4.339293538942833e-20,
         -4.68151013627075e-20,-0.04712919605086813,6.893886246929368e-06,
         -2.4691399433555226e-05,-2.7337830810319805e-06,3.726252525682523e-08,
         -1.0515675063167236e-13,1.591454248835739e-13,3.053932162293519e-18,
         -9.360829215949115e-17,4.002296426843836e-19,-1.4894216811964654e-26,
         0.0004197684396043355,-1.8732722222493324e-06,6.033160136276518e-07,
         1.2176711956423159e-08,-2.062820321383044e-10,-1.1291543376052022e-13,
         -2.439658121332727e-16,1.7915626144639213e-19,-8.842129866028896e-20,
         1.6742219332109106e-21,4.263995859323078e-27,-2.138344965559571e-06,
         3.214635151772535e-10,-4.334635876244454e-09,-2.5406627767533693e-13,
         -4.663374881131278e-16,5.985383017939705e-15,5.032395312757073e-18,
         2.8018146883583936e-22,-1.6611650939055894e-22,5.358372694525906e-24,
         -8.462210368880217e-26,6.152870110997017e-09,1.347360353950362e-10,
         9.431561152562534e-12,2.0079746039341004e-14,1.5719724367552726e-16,
         -5.2446655819462686e-18,-2.0705070511073696e-19,-2.480922327826887e-22,
         1.6382316916191638e-23,-1.4051484459006766e-25,6.577693048301765e-28,
         -9.252604312148363e-12,-6.014074995390795e-13,6.689614000017945e-17,
         -2.602920051152143e-16,-9.332333141173395e-22,1.4543725660072884e-20,
         4.2294696816900514e-22,-7.190623685255967e-25,-1.8751070930782814e-26,
         2.3268376697378167e-28,-1.330640955851776e-30,5.664842198596399e-15,
         7.143889573805286e-16,-9.369499427606346e-18,8.11159678374021e-20,
         9.545333699396679e-21,-1.717759863496281e-22,5.874702834998392e-25,
         6.326001764152421e-27,-1.009832457124965e-28,4.329510333618922e-31,
         -3.330400033066543e-35);

# 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;

    # 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);

        $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 {
            # 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;
            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 * 2.5/8
                                    + $nr_band_2 * 6.5/8 + $nr_band_3 * 7.5/8;
                }

                # 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 (0 && $nr_coverage < 1/8) {
                        if ($nr_stddev > 0.005) {
                            $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 &lt;name&gt; feels_like &lt;out-device&gt; &lt;temperature&gt; &lt;humidity&gt; [&lt;wind_speed&gt; [&lt;solar_radiation&gt; [&lt;pressure&gt;]]]</code><br>
    <br>
    <ul>
        Calculates a 'feels like' temperatures for specified readings and writes them into new readings for device &lt;out-device&gt;.<br>
        Input readings can be specified as &lt;device&gt;:&lt;reading&gt. At least one input reading must belong to &lt;out-device&gt;.<br><br>

<i>Input Readings</i>
<table>
<tr><td>temperature</td> <td>ambient temperature [&deg;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>solar_radiation&nbsp;</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 [&deg;C]</td></tr>
<tr><td>temperature_at</td> <td>temperature according to Steadman [&deg;C]</td></tr>
<tr><td>temperature_mrt</td> <td>mean radiant temperature [&deg;C]</td></tr>
<tr><td>extra_radiation</td> <td>extra radiation on human body [W/m^2]</td></tr>
</table><br><br>

If &lt;solar_radiation&gt; 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 &lt;solar_radiation&gt; 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 solar_radiation 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 &nbsp;</td> <td>required</td> <td>must be defined in <code>global</code></td></tr>
<tr><td>sensorSkyViewFactor</td> <td>optional, should match installation &nbsp;</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>
</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&ouml;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