Implementing equations with very small numbers in C - Plank's Law generating blackbody

Question

I have a problem that, after much head scratching, I think is to do with very small numbers in a long-double.

I am trying to implement Planck's law equation to generate a normalised blackbody curve at 1nm intervals between a given wavelength range and for a given temperature. Ultimately this will be a function accepting inputs, for now it is main() with the variables fixed and outputting by printf().

I see examples in matlab and python, and they are implementing the same equation as me in a similar loop with no trouble at all.

This is the equation:

My code generates an incorrect blackbody curve:

I have tested key parts of the code independently. After trying to test the equation by breaking it into blocks in excel I noticed that it does result in very small numbers and I wonder if my implementation of large numbers could be causing the issue? Does anyone have any insight into using C to implement equations? This a new area to me and I have found the maths much harder to implement and debug than normal code.

#include <stdio.h>
#include <math.h>
#include <stdlib.h>  

//global variables
const double H = 6.626070040e-34; //Planck's constant (Joule-seconds)
const double C = 299800000;       //Speed of light in vacume (meters per second)
const double K = 1.3806488e-23;   //Boltzmann's constant (Joules per Kelvin)
const double nm_to_m = 1e-6;      //conversion between nm and m
const int interval = 1;           //wavelength interval to caculate at (nm)

//typedef structure to hold results
typedef struct {
    int *wavelength;
    long double *radiance;
    long double *normalised;
} results;

int main() {
    int min = 100 , max = 3000;            //wavelength bounds to caculate between, later to be swaped to function inputs
    double temprature = 200;               //temprature in kelvin, later to be swaped to function input
    double new_valu, old_valu = 0;
    static results SPD_data, *SPD;         //setup a static results structure and a pointer to point to it
    SPD = &SPD_data;
    SPD->wavelength = malloc(sizeof(int) * (max - min));          //allocate memory based on wavelength bounds
    SPD->radiance = malloc(sizeof(long double) * (max - min));
    SPD->normalised = malloc(sizeof(long double) * (max - min));

    for (int i = 0; i <= (max - min); i++) {
        //Fill wavelength vector
        SPD->wavelength[i] = min + (interval * i);

        //Computes radiance for every wavelength of blackbody of given temprature
        SPD->radiance[i] = ((2 * H * pow(C, 2)) / (pow((SPD->wavelength[i] / nm_to_m), 5))) * (1 / (exp((H * C) / ((SPD->wavelength[i] / nm_to_m) * K * temprature))-1));

        //Copy SPD->radiance to SPD->normalised
        SPD->normalised[i] = SPD->radiance[i];

        //Find largest value
        if (i <= 0) {
            old_valu = SPD->normalised[0];
        } else if (i > 0){
            new_valu = SPD->normalised[i];
            if (new_valu > old_valu) {
                old_valu = new_valu;
            }
        }
    }
    //for debug perposes
    printf("wavelength(nm)  radiance(Watts per steradian per meter squared)  normalised radiance\n");

    for (int i = 0; i <= (max - min); i++) {
        //Normalise SPD
        SPD->normalised[i] = SPD->normalised[i] / old_valu;

        //for debug perposes
        printf("%d %Le %Lf\n", SPD->wavelength[i], SPD->radiance[i], SPD->normalised[i]);
    }
    return 0; //later to be swaped to 'return SPD';
}

/*********************UPDATE Friday 24th Mar 2017 23:42*************************/

Thank you for the suggestions so far, lots of useful pointers especially understanding the way numbers are stored in C (IEEE 754) but I don't think that is the issue here as it only applies to significant digits. I implemented most of the suggestions but still no progress on the problem. I suspect Alexander in the comments is probably right, changing the units and order of operations is likely what I need to do to make the equation work like the matlab or python examples, but my knowledge of maths is not good enough to do this. I broke the equation down into chunks to take a closer look at what it was doing.

//global variables
const double H = 6.6260700e-34;   //Planck's constant (Joule-seconds) 6.626070040e-34
const double C = 299792458;         //Speed of light in vacume (meters per second)
const double K = 1.3806488e-23;     //Boltzmann's constant (Joules per Kelvin) 1.3806488e-23
const double nm_to_m = 1e-9;        //conversion between nm and m
const int interval = 1;             //wavelength interval to caculate at (nm)
const int min = 100, max = 3000;    //max and min wavelengths to caculate between (nm)
const double temprature = 200;      //temprature (K)

//typedef structure to hold results
typedef struct {
    int *wavelength;
    long double *radiance;
    long double *normalised;
} results;

//main program
int main()
{
    //setup a static results structure and a pointer to point to it
    static results SPD_data, *SPD;
    SPD = &SPD_data;
    //allocate memory based on wavelength bounds
    SPD->wavelength = malloc(sizeof(int) * (max - min));
    SPD->radiance = malloc(sizeof(long double) * (max - min));
    SPD->normalised = malloc(sizeof(long double) * (max - min));

    //break equasion into visible parts for debuging
    long double aa, bb, cc, dd, ee, ff, gg, hh, ii, jj, kk, ll, mm, nn, oo;

    for (int i = 0; i < (max - min); i++) {
        //Computes radiance at every wavelength interval for blackbody of given temprature
        SPD->wavelength[i] = min + (interval * i);

        aa = 2 * H;
        bb = pow(C, 2);
        cc = aa * bb;
        dd = pow((SPD->wavelength[i] / nm_to_m), 5);
        ee = cc / dd;

        ff = 1;
        gg = H * C;
        hh = SPD->wavelength[i] / nm_to_m;
        ii = K * temprature;
        jj = hh * ii;
        kk = gg / jj;
        ll = exp(kk);
        mm = ll - 1;
        nn = ff / mm;

        oo = ee * nn;

        SPD->radiance[i] = oo;
    }

    //for debug perposes
    printf("wavelength(nm) | radiance(Watts per steradian per meter squared)\n");
    for (int i = 0; i < (max - min); i++) {
        printf("%d %Le\n", SPD->wavelength[i], SPD->radiance[i]);
    }
    return 0;
}

Equation variable values during runtime in xcode:

Answer 1

I notice a couple of things that are wrong and/or suspicious about the current state of your program:

• You have defined nm_to_m as 10^-9,, yet you divide by it. If your wavelength is measured in nanometers, you should multiply it by 10^-9 to get it in meters. To wit, if hh is supposed to be your wavelength in meters, it is on the order of several light-hours.
• The same is obviously true for dd as well.
• mm, being the exponential expression minus 1, is zero, which gives you infinity in the results deriving from it. This is apparently because you don't have enough digits in a double to represent the significant part of the exponential. Instead of using exp(...) - 1 here, try using the expm1() function instead, which implements a well-defined algorithm for calculating exponentials minus 1 without cancellation errors.
• Since interval is 1, it doesn't currently matter, but you can probably see that your results wouldn't match the meaning of the code if you set interval to something else.
• Unless you plan to change something about this in the future, there shouldn't be a need for this program to "save" the values of all calculations. You could just print them out as you run them.

On the other hand, you don't seem to be in any danger of underflow or overflow. The largest and smallest numbers you use don't seem to be a far way from 10^±60, which is well within what ordinary doubles can deal with, let alone long doubles. The being said, it might not hurt to use more normalized units, but at the magnitudes you currently display, I wouldn't worry about it.

Answer 2

Thanks for all the pointers in the comments. For anyone else running into a similar problem with implementing equations in C, I had a few silly errors in the code:

• writing a 6 not a 9
• dividing when I should be multiplying
• an off by one error with the size of my array vs the iterations of for() loop
• 200 when I meant 2000 in the temperature variable

As a result of the last one particularly I was not getting the results I expected (my wavelength range was not right for plotting the temperature I was calculating) and this was leading me to the assumption that something was wrong in the implementation of the equation, specifically I was thinking about big/small numbers in C because I did not understand them. This was not the case.

In summary, I should have made sure I knew exactly what my equation should be outputting for given test conditions before implementing it in code. I will work on getting more comfortable with maths, particularly algebra and dimensional analysis.

Below is the working code, implemented as a function, feel free to use it for anything but obviously no warranty of any kind etc.

blackbody.c

//
//  Computes radiance for every wavelength of blackbody of given temprature
//
//  INPUTS: int min wavelength to begin calculation from (nm), int max wavelength to end calculation at (nm), int temperature (kelvin)
//  OUTPUTS: pointer to structure containing:
//              - spectral radiance (Watts per steradian per meter squared per wavelength at 1nm intervals)
//              - normalised radiance
//

//include & define
#include "blackbody.h"

//global variables
const double H = 6.626070040e-34;   //Planck's constant (Joule-seconds) 6.626070040e-34
const double C = 299792458;         //Speed of light in vacuum (meters per second)
const double K = 1.3806488e-23;     //Boltzmann's constant (Joules per Kelvin) 1.3806488e-23
const double nm_to_m = 1e-9;        //conversion between nm and m
const int interval = 1;             //wavelength interval to calculate at (nm), to change this line 45 also need to be changed

bbresults* blackbody(int min, int max, double temperature) {
    double new_valu, old_valu = 0;      //variables for normalising result
    bbresults *SPD;
    SPD = malloc(sizeof(bbresults));
    //allocate memory based on wavelength bounds
    SPD->wavelength = malloc(sizeof(int) * (max - min));
    SPD->radiance = malloc(sizeof(long double) * (max - min));
    SPD->normalised = malloc(sizeof(long double) * (max - min));

    for (int i = 0; i < (max - min); i++) {
        //Computes radiance for every wavelength of blackbody of given temperature
        SPD->wavelength[i] = min + (interval * i);
        SPD->radiance[i] = ((2 * H * pow(C, 2)) / (pow((SPD->wavelength[i] * nm_to_m), 5))) * (1 / (expm1((H * C) / ((SPD->wavelength[i] * nm_to_m) * K * temperature))));

        //Copy SPD->radiance to SPD->normalised
        SPD->normalised[i] = SPD->radiance[i];

        //Find largest value
        if (i <= 0) {
            old_valu = SPD->normalised[0];
        } else if (i > 0){
            new_valu = SPD->normalised[i];
            if (new_valu > old_valu) {
                old_valu = new_valu;
            }
        }
    }

    for (int i = 0; i < (max - min); i++) {
        //Normalise SPD
        SPD->normalised[i] = SPD->normalised[i] / old_valu;
    }

    return SPD;
}

blackbody.h

#ifndef blackbody_h
#define blackbody_h

#include <stdio.h>
#include <math.h>
#include <stdlib.h>

//typedef structure to hold results
typedef struct {
    int *wavelength;
    long double *radiance;
    long double *normalised;
} bbresults;

//function declarations
bbresults* blackbody(int, int, double);

#endif /* blackbody_h */

main.c

#include <stdio.h>
#include "blackbody.h"

int main() {
    bbresults *TEST;
    int min = 100, max = 3000, temp = 5000;

    TEST = blackbody(min, max, temp);

    printf("wavelength | normalised radiance |           radiance               |\n");
    printf("   (nm)    |          -          | (W per meter squr per steradian) |\n");
    for (int i = 0; i < (max - min); i++) {
        printf("%4d %Lf %Le\n", TEST->wavelength[i], TEST->normalised[i], TEST->radiance[i]);
    }

    free(TEST);
    free(TEST->wavelength);
    free(TEST->radiance);
    free(TEST->normalised);
    return 0;
}

Plot of output: