在 C 中用非常小的数字实现方程 - 普朗克定律生成黑体
Implementing equations with very small numbers in C - Plank's Law generating blackbody
我有一个问题，经过多次摸索之后，我认为与 long-double 中非常小的数字有关。
我正在尝试实施普朗克定律方程，以在给定波长范围和给定温度之间以 1 纳米的间隔生成归一化黑体曲线。最终这将是一个接受输入的函数，现在它是 main()，变量固定并由 printf() 输出。
我在 matlab and 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)，但我认为这不是这里的问题，因为它只适用于重要的数字。我实施了大部分建议，但问题仍然没有进展。我怀疑评论中的 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 中运行时的方程变量值：
我注意到一些错误的地方 and/or 怀疑你程序的当前状态：
- 您已将
nm_to_m 定义为 10-9,，但您除以它。如果您的波长是以纳米为单位测量的，则应将其乘以 10-9 以获得以米为单位的波长。也就是说，如果 hh 应该是以米为单位的波长，那么它大约是几个光小时。
dd显然也是如此。mm，作为指数表达式负 1，为零，由此得出的结果为无穷大。这显然是因为双精度数中没有足够的数字来表示指数的重要部分。不要在这里使用 exp(...) - 1，而是尝试使用 expm1() 函数，它实现了一个定义明确的算法来计算指数减 1 而没有取消错误。- 由于
interval 为 1，目前它并不重要，但如果将 interval 设置为其他内容，您可能会发现结果与代码的含义不符。
- 除非您打算将来对此进行更改，否则本程序不需要 "save" 所有计算的值。你可以直接把它们打印出来 运行。
另一方面，您似乎没有任何下溢或上溢的危险。你用的最大和最小的数字好像都和10±60相差不远，在普通的double可以应付的范围内，更不用说long double秒。话虽这么说，使用更多的标准化单位可能不会 伤害 ，但在您当前显示的幅度下，我不会担心。
感谢评论中的所有指点。对于其他 运行 在 C 中实现方程式的类似问题，我在代码中有一些愚蠢的错误：
- 写 6 而不是 9
- 在应该乘法的时候除法
- 我的数组大小与 for() 循环的迭代相差一个错误
- 200，我指的是温度变量中的 2000
由于最后一个结果，特别是我没有得到我预期的结果（我的波长范围不适合绘制我正在计算的温度），这让我假设方程式的实现，特别是我在考虑 C 中的 big/small 数字，因为我不理解它们。事实并非如此。
总而言之，在用代码实现之前，我应该确保我确切地知道在给定的测试条件下我的方程式应该输出什么。我会努力学习数学，尤其是 algebra and dimensional analysis.
下面是作为函数实现的工作代码，可以随意将其用于任何用途，但显然不提供任何类型的保证等。
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;
}
输出图：
我有一个问题，经过多次摸索之后，我认为与 long-double 中非常小的数字有关。
我正在尝试实施普朗克定律方程，以在给定波长范围和给定温度之间以 1 纳米的间隔生成归一化黑体曲线。最终这将是一个接受输入的函数，现在它是 main()，变量固定并由 printf() 输出。
我在 matlab and 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)，但我认为这不是这里的问题，因为它只适用于重要的数字。我实施了大部分建议，但问题仍然没有进展。我怀疑评论中的 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 中运行时的方程变量值：
我注意到一些错误的地方 and/or 怀疑你程序的当前状态：
- 您已将
nm_to_m定义为 10-9,，但您除以它。如果您的波长是以纳米为单位测量的，则应将其乘以 10-9 以获得以米为单位的波长。也就是说，如果hh应该是以米为单位的波长，那么它大约是几个光小时。 dd显然也是如此。mm，作为指数表达式负 1，为零，由此得出的结果为无穷大。这显然是因为双精度数中没有足够的数字来表示指数的重要部分。不要在这里使用exp(...) - 1，而是尝试使用expm1()函数，它实现了一个定义明确的算法来计算指数减 1 而没有取消错误。- 由于
interval为 1，目前它并不重要，但如果将interval设置为其他内容，您可能会发现结果与代码的含义不符。 - 除非您打算将来对此进行更改，否则本程序不需要 "save" 所有计算的值。你可以直接把它们打印出来 运行。
另一方面，您似乎没有任何下溢或上溢的危险。你用的最大和最小的数字好像都和10±60相差不远，在普通的double可以应付的范围内，更不用说long double秒。话虽这么说，使用更多的标准化单位可能不会 伤害 ，但在您当前显示的幅度下，我不会担心。
感谢评论中的所有指点。对于其他 运行 在 C 中实现方程式的类似问题，我在代码中有一些愚蠢的错误：
- 写 6 而不是 9
- 在应该乘法的时候除法
- 我的数组大小与 for() 循环的迭代相差一个错误
- 200，我指的是温度变量中的 2000
由于最后一个结果，特别是我没有得到我预期的结果（我的波长范围不适合绘制我正在计算的温度），这让我假设方程式的实现，特别是我在考虑 C 中的 big/small 数字，因为我不理解它们。事实并非如此。
总而言之，在用代码实现之前，我应该确保我确切地知道在给定的测试条件下我的方程式应该输出什么。我会努力学习数学，尤其是 algebra and dimensional analysis.
下面是作为函数实现的工作代码，可以随意将其用于任何用途，但显然不提供任何类型的保证等。
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;
}
输出图：