#include <cstdio>
#include <cmath>
int main() {
    const double p = 4.0;                 // 極対数
    const double R = 0.05;                // 巻線抵抗 [Ohm]
    const double Ld = 1.0e-3;             // d軸インダクタンス [H]
    const double Lq = 2.0e-3;             // q軸インダクタンス [H]
    const double psia = 0.0601835;        // 永久磁石磁束鎖交数 [Wb]
    const double I_rms = 50.0;            // インバータ最大相電流実効値 [Arms]
    const double V_dc = 200.0;            // 直流リンク電圧 [V]
    const double PI = 3.14159265358979323846;
    // 絶対変換を考慮したdq電流の最大値
    const double Imax = I_rms * std::sqrt(3.0);
    // 今回は正弦波PWMとして、相電圧最大値を直流電圧の1/2とする
    const double Vmax = V_dc / 2.0;
    std::printf("Imax = %.4f A\n", Imax);
    std::printf("Vmax = %.4f V\n", Vmax);
    std::printf("N[rpm], MaxTorque[Nm]\n");

    for (int N = 0; N <= 20000; N += 200) {
        double omega_e = 2.0 * PI * ((double)N / 60.0) * p;
        double max_T = 0.0;
        for (double Ia = 0.0; Ia <= Imax; Ia += 1.0) {
            for (double beta_deg = 0.0; beta_deg <= 90.0; beta_deg += 0.5) {
                double beta = beta_deg * PI / 180.0;
                double id = -Ia * std::sin(beta);
                double iq =  Ia * std::cos(beta);
                double vd = R * id - omega_e * Lq * iq;
                double vq = R * iq + omega_e * (Ld * id + psia);
                // 絶対変換ベースのdq電圧を相電圧振幅へ換算
                double Va = std::sqrt(2.0 / 3.0) * std::sqrt(vd * vd + vq * vq);
                // 電圧制限を満たす場合のみトルクを計算
                if (Va <= Vmax) {
                    double T = p * (psia * iq + (Ld - Lq) * id * iq);
                    if (T > max_T) {
                        max_T = T;
                    }
                }
            }
        }
        std::printf("%d, %.4f\n", N, max_T);
    }
    return 0;
}