[feature request] dewpoint calculation

[b3-4r]

Hi, I was wondering- would it be possible to add a dewpoint display (and perhaps logging) based on the drybulb and RH measurements? It's a pretty useful number to know to predict fogging and surface condensation.

Thanks a lot!

[NerdyProjects]

Logging is not neccessary, as you can easily calculate it later. For live display, I use MijiaTemp which can also calculate the dew point. But yep, display would be very nice!

There was a similar issue and apparently the ALU power of the CPU is really bad - you would need some fix point math lib, maybe it is a problem of flash space as well... e.g. would suffer the recording buffer.

[b3-4r]

Ah right, low power silicon. I'll figure out a way to do it offboard. Perhaps the bangle.js app could do the calculation and optionally beam it back to be displayed but that's a topic for another forum.

[cultur98]

Hi there, thanks for bringing this up! It actually refers to a question I already asked a while ago. Victor provided a off device download solution which works fine.

But I was interested in live values. So I wrote myself an ATC Thermometer to MQTT Gateway based on ESP32 and Arduino. I also used this code from Victor. If you are interested in the code I can share it. (But be aware that is terrible .. spaghetti .. hardly commented .. but it works ... somehow.)

However I still like the idea the same as you @b3-4r to show the dew point on the device. This would be a very cheap, usable and portable solution! A anti-mould-ometer in the pocket. So to say.

To calculate the dew point on the thingy internally I can think of three possible solutions:

1. Use a lookup table. This needs flash memory but does not need a lot of computational power. One probably would need to sacrifice a few values for the logger. (Which is also a great feature of this project!!)
2. Run the calculation on the device with fixed precision math. Does need some computational power and battery energy. But the dewpoint must not be calculated with every new value from the sensor.
3. Use the approximation which is fine for relative humidity values higher than 50%.

My favorite is option 2. although (and maybe probably) I don't know whether it works.

[cultur98]

Holiday and mediocre weather: I managed to test the dew point calculation on one of my DEVICE_LYWSD03MMC The implementation is quite a hack hijacking the show_clock() function. But it works!

I use functions from libfixmath and the first formula mentioned in the Wikipedia article for the dew point.

It shows temperature/humidity and the dew point (and "-1" below) alternating.

[cultur98]

Code for the dew point calculation the log() eats the time

// the fixed point number functions have been copied from
// the library https://github.com/PetteriAimonen/libfixmath
// Great thanks !!

typedef int32_t fix16_t;

static const fix16_t fix16_one = 0x00010000; /*!< fix16_t value of 1 */
static const fix16_t fix16_e   = 178145;     /*!< fix16_t value of e */

static const fix16_t fix16_maximum  = 0x7FFFFFFF; /*!< the maximum value of fix16_t */
static const fix16_t fix16_minimum  = 0x80000000; /*!< the minimum value of fix16_t */
static const fix16_t fix16_overflow = 0x80000000; /*!< the value used to indicate overflows when FIXMATH_NO_OVERFLOW is not specified */

static inline fix16_t fix16_from_int(int a)     { return a * fix16_one; }
static inline float   fix16_to_float(fix16_t a) { return (float)a / fix16_one; }

static inline int fix16_to_int(fix16_t a)
{
#ifdef FIXMATH_NO_ROUNDING
    return (a >> 16);
#else
    if (a >= 0)
        return (a + (fix16_one >> 1)) / fix16_one;
    return (a - (fix16_one >> 1)) / fix16_one;
#endif
}

static inline fix16_t fix16_from_float(float a)
{
    float temp = a * fix16_one;
#ifndef FIXMATH_NO_ROUNDING
    temp += (temp >= 0) ? 0.5f : -0.5f;
#endif
    return (fix16_t)temp;
}

static inline uint32_t fix_abs(fix16_t in)
{
    if(in == fix16_minimum)
    {
        // minimum negative number has same representation as
        // its absolute value in unsigned
        return 0x80000000;
    }
    else
    {
        return ((in >= 0)?(in):(-in));
    }
}

fix16_t fix16_add(fix16_t a, fix16_t b)
{
    // Use unsigned integers because overflow with signed integers is
    // an undefined operation (http://www.airs.com/blog/archives/120).
    uint32_t _a = a;
    uint32_t _b = b;
    uint32_t sum = _a + _b;

    // Overflow can only happen if sign of a == sign of b, and then
    // it causes sign of sum != sign of a.
    if (!((_a ^ _b) & 0x80000000) && ((_a ^ sum) & 0x80000000))
        return fix16_overflow;

    return sum;
}

fix16_t fix16_sub(fix16_t a, fix16_t b)
{
    uint32_t _a = a;
    uint32_t _b = b;
    uint32_t diff = _a - _b;

    // Overflow can only happen if sign of a != sign of b, and then
    // it causes sign of diff != sign of a.
    if (((_a ^ _b) & 0x80000000) && ((_a ^ diff) & 0x80000000))
        return fix16_overflow;

    return diff;
}

/* Saturating arithmetic */
fix16_t fix16_sadd(fix16_t a, fix16_t b)
{
    fix16_t result = fix16_add(a, b);

    if (result == fix16_overflow)
        return (a >= 0) ? fix16_maximum : fix16_minimum;

    return result;
}   

fix16_t fix16_ssub(fix16_t a, fix16_t b)
{
    fix16_t result = fix16_sub(a, b);

    if (result == fix16_overflow)
        return (a >= 0) ? fix16_maximum : fix16_minimum;

    return result;
}

fix16_t fix16_mul(fix16_t inArg0, fix16_t inArg1)
{
    // Each argument is divided to 16-bit parts.
    //                  AB
    //          *    CD
    // -----------
    //                  BD  16 * 16 -> 32 bit products
    //               CB
    //               AD
    //              AC
    //           |----| 64 bit product
    int32_t A = (inArg0 >> 16), C = (inArg1 >> 16);
    uint32_t B = (inArg0 & 0xFFFF), D = (inArg1 & 0xFFFF);

    int32_t AC = A*C;
    int32_t AD_CB = A*D + C*B;
    uint32_t BD = B*D;

    int32_t product_hi = AC + (AD_CB >> 16);

    // Handle carry from lower 32 bits to upper part of result.
    uint32_t ad_cb_temp = AD_CB << 16;
    uint32_t product_lo = BD + ad_cb_temp;
    if (product_lo < BD)
        product_hi++;

#ifndef FIXMATH_NO_OVERFLOW
    // The upper 17 bits should all be the same (the sign).
    if (product_hi >> 31 != product_hi >> 15)
        return fix16_overflow;
#endif

#ifdef FIXMATH_NO_ROUNDING
    return (product_hi << 16) | (product_lo >> 16);
#else
    // Subtracting 0x8000 (= 0.5) and then using signed right shift
    // achieves proper rounding to result-1, except in the corner
    // case of negative numbers and lowest word = 0x8000.
    // To handle that, we also have to subtract 1 for negative numbers.
    uint32_t product_lo_tmp = product_lo;
    product_lo -= 0x8000;
    product_lo -= (uint32_t)product_hi >> 31;
    if (product_lo > product_lo_tmp)
        product_hi--;

    // Discard the lowest 16 bits. Note that this is not exactly the same
    // as dividing by 0x10000. For example if product = -1, result will
    // also be -1 and not 0. This is compensated by adding +1 to the result
    // and compensating this in turn in the rounding above.
    fix16_t result = (product_hi << 16) | (product_lo >> 16);
    result += 1;
    return result;
#endif
}

fix16_t fix16_div(fix16_t a, fix16_t b)
{
    // This uses the basic binary restoring division algorithm.
    // It appears to be faster to do the whole division manually than
    // trying to compose a 64-bit divide out of 32-bit divisions on
    // platforms without hardware divide.

    if (b == 0)
        return fix16_minimum;

    uint32_t remainder = fix_abs(a);
    uint32_t divider = fix_abs(b);

    uint32_t quotient = 0;
    uint32_t bit = 0x10000;

    /* The algorithm requires D >= R */
    while (divider < remainder)
    {
        divider <<= 1;
        bit <<= 1;
    }

    #ifndef FIXMATH_NO_OVERFLOW
    if (!bit)
        return fix16_overflow;
    #endif

    if (divider & 0x80000000)
    {
        // Perform one step manually to avoid overflows later.
        // We know that divider's bottom bit is 0 here.
        if (remainder >= divider)
        {
                quotient |= bit;
                remainder -= divider;
        }
        divider >>= 1;
        bit >>= 1;
    }

    /* Main division loop */
    while (bit && remainder)
    {
        if (remainder >= divider)
        {
                quotient |= bit;
                remainder -= divider;
        }

        remainder <<= 1;
        bit >>= 1;
    }    

    #ifndef FIXMATH_NO_ROUNDING
    if (remainder >= divider)
    {
        quotient++;
    }
    #endif

    fix16_t result = quotient;

    /* Figure out the sign of result */
    if ((a ^ b) & 0x80000000)
    {
        #ifndef FIXMATH_NO_OVERFLOW
        if (result == fix16_minimum)
                return fix16_overflow;
        #endif

        result = -result;
    }

    return result;
}

fix16_t fix16_exp(fix16_t inValue) {
    if(inValue == 0        ) return fix16_one;
    if(inValue == fix16_one) return fix16_e;
    if(inValue >= 681391   ) return fix16_maximum;
    if(inValue <= -772243  ) return 0;

    /* The algorithm is based on the power series for exp(x):
     * http://en.wikipedia.org/wiki/Exponential_function#Formal_definition
     * 
     * From term n, we get term n+1 by multiplying with x/n.
     * When the sum term drops to zero, we can stop summing.
     */

    // The power-series converges much faster on positive values
    // and exp(-x) = 1/exp(x).
    bool neg = (inValue < 0);
    if (neg) inValue = -inValue;

    fix16_t result = inValue + fix16_one;
    fix16_t term = inValue;

    uint_fast8_t i;        
    for (i = 2; i < 30; i++)
    {
        term = fix16_mul(term, fix16_div(inValue, fix16_from_int(i)));
        result += term;

        if ((term < 500) && ((i > 15) || (term < 20)))
            break;
    }

    if (neg) result = fix16_div(fix16_one, result);
    return result;
}

fix16_t fix16_log(fix16_t inValue)
{
    fix16_t guess = fix16_from_int(2);
    fix16_t delta;
    int scaling = 0;
    int count = 0;

    if (inValue <= 0)
        return fix16_minimum;

    // Bring the value to the most accurate range (1 < x < 100)
    const fix16_t e_to_fourth = 3578144;
    while (inValue > fix16_from_int(100))
    {
        inValue = fix16_div(inValue, e_to_fourth);
        scaling += 4;
    }

    while (inValue < fix16_one)
    {
        inValue = fix16_mul(inValue, e_to_fourth);
        scaling -= 4;
    }

    do
    {
        // Solving e(x) = y using Newton's method
        // f(x) = e(x) - y
        // f'(x) = e(x)
        fix16_t e = fix16_exp(guess);
        delta = fix16_div(inValue - e, e);

        // It's unlikely that logarithm is very large, so avoid overshooting.
        if (delta > fix16_from_int(3))
            delta = fix16_from_int(3);

        guess += delta;
    } while ((count++ < 10)
        && ((delta > 1) || (delta < -1)));

    return guess + fix16_from_int(scaling);
}

int16_t calc_dew_f1(int16_t t_e2, int16_t rh_e2) {
    // the next four initializations have to be done only once 
    // could run in an init function to safe time
  fix16_t f_100 = fix16_from_int(100);
  fix16_t f_10000 = fix16_from_int(10000);
  fix16_t f_b = fix16_from_float(18.678);
  fix16_t f_c = fix16_from_float(257.14);

    // convert input values to fix16_t
  fix16_t T  = fix16_from_int(t_e2);
  T = fix16_div(T, f_100);
  fix16_t RH = fix16_from_int(rh_e2);
  fix16_t RH_100 = fix16_div(RH, f_10000);

    // the actual calculation of the dewpoint
  fix16_t gamma = fix16_add(fix16_log(RH_100), fix16_div(fix16_mul(f_b, T), fix16_add(f_c, T)));
  fix16_t Tdp = fix16_div(fix16_mul(f_c, gamma), fix16_sub(f_b, gamma));
  fix16_t i_Tdp = fix16_mul(Tdp, f_100);
  return (fix16_to_int(i_Tdp)); // in Celsius
}

[CraigHutchinson]

@cultur98 Great work.