プランクの放射法則
４．任意の温度に対する完全導体中の電磁波シミュレーション
量子力学の創成に重要な影響を与えたプランクの放射法則のシミュレーションを行います。
本稿では、完全導体で囲まれた真空中における電磁波の基本モードを用いた各種シミュレーションを次のステップで行います。
１．プランクの法則の解析解
２．完全導体で囲まれた空洞における電磁波の基本モード
３．完全導体で囲まれた空洞における各電磁波モード振幅の温度依存性
４．任意の温度に対する完全導体中の電磁波シミュレーション
本稿の目的である「プランクの放射法則」の理解の最後として、任意の温度に対する完全導体中の電場の時間発展を計算します。
完全導体中の電磁波復習ついでに、計算までの道筋を示します。
（１）完全導体中の電磁波のモード
は互いに独立ではなく、
電磁波の波数ベクトル
によって決まります。ただし、電磁波の振幅は任意です。
と
の関係は、
「２．完全導体で囲まれた空洞における電磁波の基本モード」
で示したとおりです。
（２）エネルギーの量子化による電場強度の制限
「３．完全導体で囲まれた空洞における各電磁波モード振幅の温度依存性」で示したとおり、 エネルギーの量子仮説を課すことで、電場強度は
となります。
ここで重要なのは、平均励起数
は温度
と角振動数
に依存することです。つまり、エネルギーの量子仮説は電磁波の振幅に制限を加えることを意味します。
ただし、
は電磁波の波数ベクトル
で決まります。
（３）電場強度の決定と各電場モードの重ね合わせ
完全導体中の任意の電場
は、
電磁波の波数ベクトル
と垂直で、互いにも垂直な
、
に分解することができます。
さらに、熱によって励起された電磁波の各電場振幅
は（等重率の原理的に？）等しいと仮定します。つまり、
と表せます。くどいですが、
,
は、波数ベクトル
に依存します。式(1)に対して、
「２．完全導体で囲まれた空洞における電磁波の基本モード」
で示した
の関係を踏まえて、
各モード
の電磁波は
となります。
も同様です。
任意の温度
に対する最終的な電場は、
各モード
に対する電場の和で得られます。また、
も加えます。つまり、
となります。
任意の温度に対する完全導体中の電磁波シミュレーション
一辺がL=5000[nm]の完全導体中の電磁波Ezの時間発展（温度 T=300[K]）。
横軸 x[nm]
縦軸 y[nm]
カラーバー：電場強度 [a.u.]
時間間隔：Δt=10^{-16}[s]
※z=L/2 の断面における電場を描画
では、
「２．完全導体で囲まれた空洞における各電磁波モード振幅の温度依存性」で示したとおり、ほとんど最低次のモードしか励起していないため、時間的に周期的な振る舞いとなっています。
では、
「２．完全導体で囲まれた空洞における各電磁波モード振幅の温度依存性」で示したとおり、最低次のモードとあと２つのモードしか励起していません。ただし、
と比べて、電場の振幅は約800倍です。
では、
「２．完全導体で囲まれた空洞における各電磁波モード振幅の温度依存性」で示したとおり、様々なモードが励起されています。それにともなって、
電場の振幅は約10000倍です。
C言語プログラムソース
/*
完全導体中の電場の時間依存性
（2012.01.06公開）
*/
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include <fstream>
#include <sstream>
#include <iostream>
#include <string>
#include <cstdio>
#include <iomanip>
#include <stdio.h>
#include <complex>
#if defined(_OPENMP)
#include <omp.h>
#endif
#if defined(_MSC_VER)
#include <direct.h> // Windowsフォルダ作成用
#elif defined(__GNUC__)
#include <sys/stat.h> // UNIX系ディレクトリ作成用
#endif
using namespace std;
const double PI = acos(-1.0);
const double e = 2.7182818284590452354;
const complex<double> I = complex<double>(0.0,1.0);
const double c = 2.99792458E+8;
const double lambda0 = 500.0E-9;
const double mu0 = 4.0*PI*1.0E-7;
const double epsilon0 = 1.0/(4.0*PI*c*c)*1.0E+7;
const double h = 6.6260896 * 1.0E-34;
const double hbar = h/(2.0*PI);
const double me = 9.10938215 * 1.0E-31;
const double eV = 1.60217733 * 1.0E-19;
const double kB = 1.3806503 * 1.0E-23;
const int Lx=5000, Ly=5000, Ld=50;
const double dz = 1.0E-9; //XY平面を描画する際には dz = 1.0E-8
const double dt = 1.0E-16;
const double L = Lx * dz;
const double V = pow(L,3);
const int ts = 0, te = 1000; //<------------ステップ数を指定
string folder = "BlackBody3";//作成するフォルダ名
char str[200];
int main(){
#if defined(_MSC_VER)
_mkdir(folder.c_str()); // Windowsフォルダ作成
#elif defined(__GNUC__)
mkdir(folder.c_str(), S_IRWXU | S_IRWXG | S_IRWXO); // UNIX系のディレクトリ作成
#endif
#if defined(_OPENMP)
omp_set_num_threads(7);
cout << "OpenMPを利用します（最大スレッド数：" << omp_get_max_threads() << "）" << endl;
#endif
const int n_max = 10;
double z = L/2.0;
for(double T=1000.0; T<=1000.0; T+=200.0){
cout << T << endl;
char str[200];
sprintf(str, "/%d", int(T));
string folder_n = folder + str;
#if defined(_MSC_VER)
_mkdir(folder_n.c_str()); // Windowsフォルダ作成
#elif defined(__GNUC__)
mkdir(folder_n.c_str(), S_IRWXU | S_IRWXG | S_IRWXO); // UNIX系のディレクトリ作成
#endif
#pragma omp parallel
{
#pragma omp for //schedule(dynamic, 100)
for(int tn=ts; tn<=te; tn++)
{
double t_real = dt * double(tn);
//cout << tn << endl;
char str[200];
string str1;
ofstream fout_s, fout_p;
sprintf(str, "/%d/s%d.data", int(T) , tn); str1 = folder + str;
fout_s.open(str1.c_str());
for(int ix=0; ix<=Lx ;ix=ix+Ld)
{
for(int iy=0 ; iy<=Lx ;iy=iy+Ld)
{
double x = ix *dz, y = iy*dz;
complex<double> EE1[3] ={0.0}, EE2[3] = {0.0};
for(int nux=0; nux<=n_max; nux++ )
{
for(int nuy=0; nuy<=n_max; nuy++ )
{
for(int nuz=0; nuz<=n_max; nuz++ )
{
if((nux==0 && nuy==0) || (nuy==0 && nuz==0) || (nuz==0 && nux==0) ) continue;
double kx = PI/L * double(nux);
double ky = PI/L * double(nuy);
double kz = PI/L * double(nuz);
double k = sqrt(pow(kx,2)+pow(ky,2)+pow(kz,2));
double omega = c * k;
double nbar = 1.0/(exp((hbar*omega)/(kB*T))-1.0);
double EE = sqrt(1.0/V * 16.0/epsilon0 * nbar * hbar * omega)/sqrt(2.0);
double E1x,E1y,E1z,E2x,E2y,E2z;
if(nux==0){
E1x = EE;
E1y = 0;
E1z = 0;
E2x = 0;
E2y = EE /sqrt(pow(1.0, 2)+pow(-ky/kz, 2));
E2z = EE*(-ky/kz)/sqrt(pow(1.0, 2)+pow(-ky/kz, 2));
}else if(nuy==0){
E1x = 0;
E1y = EE;
E1z = 0;
E2x = EE /sqrt(pow(1.0, 2)+pow(-kx/kz, 2));
E2y = 0;
E2z = EE*(-kx/kz)/sqrt(pow(1.0, 2)+pow(-kx/kz, 2));
}else if(nuz==0){
E1x = 0;
E1y = 0;
E1z = EE;
E2x = EE /sqrt(pow(1.0, 2)+pow(-kx/ky, 2));
E2y = EE*(-kx/ky)/sqrt(pow(1.0, 2)+pow(-kx/ky, 2));
E2z = 0;
}else{
E1x = EE /sqrt(pow(1.0, 2)+pow(-kx/ky, 2));
E1y = EE*(-kx/ky)/sqrt(pow(1.0, 2)+pow(-kx/ky, 2));
E1z = 0;
E2x = EE /sqrt(pow(1.0, 2)+pow(ky/kx, 2)+pow(-kx/kz - ky*ky/(kz*kx), 2));
E2y = EE*(ky/kx) /sqrt(pow(1.0, 2)+pow(ky/kx, 2)+pow(-kx/kz - ky*ky/(kz*kx), 2));
E2z = EE*(-kx/kz-ky*ky/(kz*kx))/sqrt(pow(1.0, 2)+pow(ky/kx, 2)+pow(-kx/kz - ky*ky/(kz*kx), 2));
}
EE1[0] += E1x*cos(kx*x)*sin(ky*y)*sin(kz*z) * exp(-I*omega*t_real);
EE1[1] += E1y*sin(kx*x)*cos(ky*y)*sin(kz*z) * exp(-I*omega*t_real);
EE1[2] += E1z*sin(kx*x)*sin(ky*y)*cos(kz*z) * exp(-I*omega*t_real);
EE2[0] += E2x*cos(kx*x)*sin(ky*y)*sin(kz*z) * exp(-I*omega*t_real);
EE2[1] += E2y*sin(kx*x)*cos(ky*y)*sin(kz*z) * exp(-I*omega*t_real);
EE2[2] += E2z*sin(kx*x)*sin(ky*y)*cos(kz*z) * exp(-I*omega*t_real);
}
}
}
fout_s << x/dz << " " << y/dz << " " << EE1[0].real()+EE2[0].real() << " " << EE1[1].real() +EE2[1].real()<< " " << EE1[2].real() + EE2[2].real() << " "<< sqrt(pow(abs(EE1[0]+EE2[0]),2)+pow(abs(EE1[1]+EE2[1]),2)+pow(abs(EE1[2]+EE2[2]),2)) << endl;
}
fout_s <<""<< endl;
}
fout_s.close();
}
}
/*
for(int i=0; i <= n_max*n_max*n_max; i++ ){
if(count[i]!=0){
double omega = c * PI/L * sqrt(double(i));
double lambda = (2.0*PI*c)/omega;
double EE = E[i] * count[i];
fout_s << lambda << " " << EE << endl;
}
}
*/
}
}
