生成黑体的C-平面定律中极小数方程的实现

Question

我有一个问题，在经历了大量的挠头之后，我认为这是一个很小的数字在一个长-双。

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

我在matlab和python中看到了一些例子，它们在一个类似的循环中实现了与我相同的等式，一点也不麻烦。

这是方程式

​

​

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

​

我已经独立测试了代码的关键部分。在尝试通过在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';
}

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

到目前为止，感谢您的建议，许多有用的指针，特别是理解数字存储在C (IEEE 754)中的方式，但我不认为这是问题所在，因为它只适用于重要的数字。我执行了大部分建议，但在这个问题上仍未取得进展。我怀疑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

感谢评论中的所有指点。对于其他在用C实现等式时遇到类似问题的人，我在代码中出现了一些愚蠢的错误：

• 写6而不是9
• 当我应该乘积的时候除法
• 与数组的大小与for()循环的迭代量对应的off错误
• 当我的意思是温度变量为2000时

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

总之，在用代码实现它之前，我应该确保确切地知道在给定的测试条件下应该输出什么公式。我将努力使数学变得更舒服，特别是代数和量纲分析。

下面是作为一个函数实现的工作代码，您可以随意使用它来做任何事情，但是显然没有任何类型的保证等等。

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

产出图：

​

​