#= ブランク関数計算プログラム
#
#= Blackbody emissive power spectrum calculation program
#
# Authors::   Hiroaki Saito
# Version::   $Id: gallery_BlackbodyEmissivePowerSpectrum.m,v 1.00 saito$
#  Tag Name::  $Name: dcrtm-20161126 $
#  Copyright:: Copyright (C) GFD Dennou Club, 2008-2016. All rights reserved.
#  License::   See COPYRIGHT[link:../../../COPYRIGHT]
hold off;
format long;
set(0, 'DefaultAxesLineWidth', 2, 'DefaultLineLineWidth', 3)
set(0, 'defaultAxesFontSize', 18);

# パラメタの設定
pi = 3.14159265359e0;
c =   299792458e0;
h = 6.626070040e-34;
kb = 1.38064852e-23;

# 波長域の設定
WaveLength(:,1) = [5e1:1e5]'.*1e-9;

# 温度のテーブル作成
Temp = [500 1000 2000 3000 5000 5800];

for i = 1:6
# プランク関数
  PlankEquation(:,i) = ...
  ((2.0e0*pi*h*c.^2)./(WaveLength(:,1)).^5) ...
  .*(1.0e0./(exp((h*c)./(WaveLength(:,1).*kb.*Temp(:,i))).-1.0e0));
end

# BlackbodyEmissivePowerSpectrum.nc データ出力
nccreate("BlackbodyEmissivePowerSpectrum.nc","WaveLength", ...
	 "Dimensions",{"WaveLength", 99951});
ncwrite("BlackbodyEmissivePowerSpectrum.nc","WaveLength",WaveLength(:,1).*1e6);
ncwriteatt("BlackbodyEmissivePowerSpectrum.nc","WaveLength", ...
	   "long_name", "WaveLength");
ncwriteatt("BlackbodyEmissivePowerSpectrum.nc","WaveLength", ...
	   "units", "um");

nccreate("BlackbodyEmissivePowerSpectrum.nc","Temp", ...
	 "Dimensions",{"Temp", 6});
ncwrite("BlackbodyEmissivePowerSpectrum.nc","Temp",Temp);
ncwriteatt("BlackbodyEmissivePowerSpectrum.nc","Temp", ...
	   "long_name", "Temperature");
ncwriteatt("BlackbodyEmissivePowerSpectrum.nc","Temp", ...
	   "units", "K");

nccreate("BlackbodyEmissivePowerSpectrum.nc","PlankEquation1", ...
	 "Dimensions",{"WaveLength", 99951,"Temp",6});
ncwrite("BlackbodyEmissivePowerSpectrum.nc","PlankEquation1", PlankEquation.*1e-6);
ncwriteatt("BlackbodyEmissivePowerSpectrum.nc","PlankEquation1", ...
	   "long_name", "Blackbody Emissive Power");
ncwriteatt("BlackbodyEmissivePowerSpectrum.nc","PlankEquation1", ...
	   "units", "W m-2 um-1");

# 各温度でのピーク波長のテーブル作成
MaxVal = ...
[1e6.*WaveLength((PlankEquation(:,1)==max(max(PlankEquation(:,1)))),1) ...
		max(PlankEquation(:,1)).*1e-6; ...
 1e6.*WaveLength((PlankEquation(:,2)==max(max(PlankEquation(:,2)))),1) ...
		max(PlankEquation(:,2)).*1e-6; ...
 1e6.*WaveLength((PlankEquation(:,3)==max(max(PlankEquation(:,3)))),1) ...
		max(PlankEquation(:,3)).*1e-6; ...
 1e6.*WaveLength((PlankEquation(:,4)==max(max(PlankEquation(:,4)))),1) ...
		max(PlankEquation(:,4)).*1e-6; ...
 1e6.*WaveLength((PlankEquation(:,5)==max(max(PlankEquation(:,5)))),1) ...
		max(PlankEquation(:,5)).*1e-6; ...
 1e6.*WaveLength((PlankEquation(:,6)==max(max(PlankEquation(:,6)))),1) ...
		max(PlankEquation(:,6)).*1e-6];

# PeakWaveLength.nc データ出力
nccreate("PeakWavelength.nc","MaxWaveLength", ...
	 "Dimensions",{"WaveLength", 6});
ncwrite("PeakWavelength.nc","MaxWaveLength",MaxVal(:,1));
ncwriteatt("PeakWavelength.nc","MaxWaveLength", ...
	   "long_name", "Peak WaveLength");
ncwriteatt("PeakWavelength.nc","MaxWaveLength", ...
	   "units", "um");

nccreate("PeakWavelength.nc","PeakBlackbodyEmission", ...
	 "Dimensions",{"WaveLength", 6});
ncwrite("PeakWavelength.nc","PeakBlackbodyEmission",MaxVal(:,2));
ncwriteatt("PeakWavelength.nc","PeakBlackbodyEmission", ...
	   "long_name", "Peak Blackbody Emissive Power");
ncwriteatt("PeakWavelength.nc","PeakBlackbodyEmission", ...
	   "units", "W m-2 um-1)");


ColorNum = 10000;                        
A = colormap(hot(ColorNum));

#描画;
loglog(WaveLength(:,1).*1e6,PlankEquation(:,1).*1e-6, "Color", A(500,:), ...
     WaveLength(:,1).*1e6,PlankEquation(:,2).*1e-6,  "Color", A(1000,:), ...
     WaveLength(:,1).*1e6,PlankEquation(:,3).*1e-6,  "Color", A(2000,:), ...
     WaveLength(:,1).*1e6,PlankEquation(:,4).*1e-6, "Color", A(3000,:), ...
     WaveLength(:,1).*1e6,PlankEquation(:,5).*1e-6, "Color", A(5000,:), ...
     WaveLength(:,1).*1e6,PlankEquation(:,6).*1e-6, "Color", A(5777,:))
hold on;

loglog(1e6.*WaveLength((PlankEquation(:,1)==max(max(PlankEquation(:,1)))),1), ...
       max(PlankEquation(:,1)).*1e-6, "ok", ...
       1e6.*WaveLength((PlankEquation(:,2)==max(max(PlankEquation(:,2)))),1), ...
       max(PlankEquation(:,2)).*1e-6, "ok", ...
       1e6.*WaveLength((PlankEquation(:,3)==max(max(PlankEquation(:,3)))),1), ...
       max(PlankEquation(:,3)).*1e-6, "ok", ...
       1e6.*WaveLength((PlankEquation(:,4)==max(max(PlankEquation(:,4)))),1), ...
       max(PlankEquation(:,4)).*1e-6, "ok", ...
       1e6.*WaveLength((PlankEquation(:,5)==max(max(PlankEquation(:,5)))),1), ...
       max(PlankEquation(:,5)).*1e-6, "ok", ...
       1e6.*WaveLength((PlankEquation(:,6)==max(max(PlankEquation(:,6)))),1), ...
       max(PlankEquation(:,6)).*1e-6, "ok")
loglog(MaxVal(:,1), MaxVal(:,2), "--k")
axis([5e-2 20 1e2 1e8])

title("Blackbody emissive power spectrum")
xlabel('Wavelength [\mum]','FontSize',20')
ylabel('Blackbody emissive power [W/m^{2}{\mu}m]','FontSize',20)

# Figure 保存;
print -dpng "BlackbodyEmissivePowerSpectrum.png" -color
#print -dpdf "BlackbodyEmissivePowerSpectrum.pdf" -color
