Search code examples
cmathlong-double

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

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:

Plancks Law

My code generates an incorrect blackbody curve: Blackbody code incorrect output

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:

variable values during runtime in xcode

Solution

  • 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:

    plot of output