#include <iostream>
#include <cmath>
#include <vector>
#include <iomanip>
#include <fstream>

// 物理定数（SI単位）
const double h = 6.62607015e-34;      // プランク定数 [J·s]
const double k = 1.380649e-23;        // ボルツマン定数 [J/K]
const double c = 2.99792458e8;        // 光速 [m/s]
const double eV = 1.602176634e-19;    // 電子ボルト [J]

// ν0 = 13.6 eV に対応する周波数
const double E0 = 13.6 * eV;
const double nu0 = E0 / h;

// 温度ごとに設定
double T = 1e4; // デフォルト値

// Planck関数 J_ν
double J_nu(double nu) {
    double numerator = 2.0 * h * std::pow(nu, 3) / (c * c);
    double exponent = h * nu / (k * T);
    double denominator = std::exp(exponent) - 1.0;
    return numerator / denominator;
}

// 被積分関数 f1 = 分子（ν - ν0）付き
double integrand1(double nu) {
    if (nu <= nu0) return 0.0;
    return 4.0 * M_PI * J_nu(nu) * (nu - nu0) / (h * std::pow(nu, 4));
}

// 被積分関数 f2 = 分母（ν - ν0 なし）
double integrand2(double nu) {
    if (nu <= nu0) return 0.0;
    return 4.0 * M_PI * J_nu(nu) / (h * std::pow(nu, 4));
}

// 台形則による数値積分
double trapezoidal_integrate(double (*f)(double), double a, double b, int n) {
    double h = (b - a) / n;
    double sum = 0.5 * (f(a) + f(b));
    for (int i = 1; i < n; ++i) {
        sum += f(a + i * h);
    }
    return sum * h;
}

int main() {
    // 出力ファイル
    std::ofstream file("output.csv");


    // 調べたい温度リスト
    std::vector<double> temperatures = {1e4,1.5e4,2e4,2.5e4, 3e4,3.5e4, 4e4,4.5e4, 5e4};

    // 積分設定
    int N = 100000;                   // 分割数
    double nu_max = nu0 * 1000.0;     // 上限（十分高い周波数）

    std::cout << std::fixed << std::setprecision(6);
    std::cout << " T [K]       I1/I2 (平均ν - ν0)\n";
    std::cout << "-------------------------------\n";

    file << "T,Ti\n";  // ヘッダ行

    for (double temp : temperatures) {
        T = temp; // 温度を更新

        double I1 = trapezoidal_integrate(integrand1, nu0, nu_max, N);
        double I2 = trapezoidal_integrate(integrand2, nu0, nu_max, N);
        double ratio = I1 / I2;

        std::cout << std::setw(8) << T << "    " << ratio * h / k << " Hz\n";

        file << std::setw(8) << T << "," << ratio * h/ k << "\n";
    }

    file.close();
    std::cout << "output.csv に出力しました。" << std::endl;

    return 0;
}