在 C 中实现具有非常小的数字的方程 - 普朗克定律生成黑体

Question

我有一个问题，经过多次头疼后，我认为与长双打中非常小的数字有关。

我正在尝试实现普朗克定律方程，以在给定波长范围和给定温度之间以 1nm 的间隔生成归一化黑体曲线。最终这将是一个接受输入的函数，现在它是 main()，变量固定并由 printf() 输出。

我在matlab 和python 中看到了示例，它们在类似的循环中执行与我相同的等式，完全没有问题。

这是equation：

我的代码生成了不正确的黑体曲线：

我已经独立测试了代码的关键部分。在尝试通过在 excel 中将方程分解成块来测试方程后，我注意到它确实会产生非常小的数字，我想知道我对大数字的实现是否会导致问题？有没有人对使用 C 来实现方程有任何见解？这对我来说是一个新领域，我发现数学比普通代码更难实现和调试。

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

/*********************更新 2017 年 3 月 24 日星期五 23:42**************** ********/

到目前为止，感谢您的建议，许多有用的指针特别是了解数字在 C 中的存储方式 (IEEE 754) 但我不认为这是这里的问题，因为它只适用于有效数字。我实施了大部分建议，但在这个问题上仍然没有进展。我怀疑 cmets 中的 Alexander 可能是对的，改变单位和操作顺序可能是我需要做的才能使方程像 matlab 或 python 示例一样工作，但我的数学知识不足以做到这一点。我把方程分解成几块，仔细看看它在做什么。

//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;
}

xcode 中运行时的方程变量值：

Answer 1

我注意到您的程序当前状态存在一些错误和/或可疑之处：

• 您已将nm_to_m 定义为 10^-9,，但您除以它。如果您的波长以纳米为单位，则应将其乘以 10^-9 以得到米。也就是说，如果hh 应该是您的波长（以米为单位），它大约是几个光小时。
• dd 显然也是如此。
• mm，作为负 1 的指数表达式，为零，这为您提供了无穷大的结果。这显然是因为双精度数中没有足够的数字来表示指数的重要部分。不要在此处使用 exp(...) - 1，而是尝试使用 expm1() 函数，该函数实现了一个定义明确的算法，用于计算指数负 1 而不会出现取消错误。
• 由于 interval 为 1，因此目前无关紧要，但如果将 interval 设置为其他值，您可能会看到结果与代码的含义不匹配。
• 除非您计划在未来对此进行更改，否则此程序不需要“保存”所有计算的值。您可以在运行它们时将它们打印出来。

另一方面，您似乎没有任何下溢或上溢的危险。您使用的最大和最小数字似乎与 10^±60 相差不远，这完全在普通doubles 可以处理的范围内，更不用说long doubles。话虽这么说，使用更多的标准化单位可能伤害，但在你目前显示的数量级上，我不会担心。

Answer 2

感谢 cmets 中的所有指针。对于在 C 中实现方程式时遇到类似问题的其他人，我在代码中遇到了一些愚蠢的错误：

• 写 6 而不是 9
• 当我应该乘法时除法
• 我的数组大小与 for() 循环的迭代次数相差一个错误
• 200 我指的是温度变量中的 2000

由于最后一个结果，特别是我没有得到我预期的结果（我的波长范围不适合绘制我正在计算的温度），这导致我假设在实施过程中出现了问题方程，特别是我在考虑 C 中的大/小数字，因为我不理解它们。事实并非如此。

总之，在代码中实现它之前，我应该确保我确切地知道在给定的测试条件下我的方程应该输出什么。我会努力让数学变得更舒服，尤其是algebra 和dimensional analysis。

以下是作为函数实现的工作代码，您可以随意将其用于任何用途，但显然不提供任何形式的保证等。

黑体.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;
}

黑体.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;
}

输出图：