AE510_HW1_SI

#include <iostream>
#include <fstream>
#include <cmath>
#include <iomanip>
using namespace std;

// Define parameters
const double gamma_air = 1.4;   // Adiabatic coefficient for air
const double f = 0.02;      // Friction coefficient
const double Rspec = 286.9; // Air specific gas constant in J/kg*K
// Hydraulic diameter: Dh = 4 inches.

const double Dh = 4.0 * 0.0254;  // 4 inch in meters

// Total pipe length in feet
const double L = 70.0 * 0.3048;  // 70 ft in meters

// Function that returns dM/dx given M.
double dMdx(double M) {
    // Avoid division by zero (when 1 - gamma*M^2 is too close to zero).
    double denom = 1.0 - gamma_air * M * M;
    if (fabs(denom) < 1e-8) {
        cerr << "Denominator approaching zero. M = " << M << endl;
        return 0.0;
    }
    return (gamma_air * f / (2.0 * Dh)) * (M * M * M) / denom;
}

// Function that returns dP/dx given M and P.
double dPdx(double M, double P) {

    double denom = 1.0 - gamma_air * M * M;
    if (fabs(denom) < 1e-8) {
        cerr << "Denominator approaching zero. M = " << M << endl;
        return 0.0;
    }
    return (-1.0 * P * f / (2.0 * Dh) ) * ( ( gamma_air * M * M) / denom );
}

double dP0dx(double M, double P0) {

    double denom1 = gamma_air * M * M - 1.0;
    double denom2 = 1 + (gamma_air - 1.0) * M * M /2.0;
    double numer1 = gamma_air * M * M / 2.0;
    double numer2 = 1 - (gamma_air + 1.0) * M * M /2.0;

    if (fabs(denom1) < 1e-8) {
        cerr << "Denominator approaching zero. M = " << M << endl;
        return 0.0;
    }
    return (P0 * f / Dh) * numer1 * numer2 / (denom1 * denom2);
}

double dVdx(double M, double V) {

    double numer = gamma_air * M * M / 2.0;
    double denom = 1.0 - gamma_air * M * M;

    if (fabs(denom) < 1e-8) {
        cerr << "Denominator approaching zero. M = " << M << endl;
        return 0.0;
    }
    return (f/Dh) * (numer/denom) * V;
}

double dT0dx(double M, double T0) {

    double numer1 = gamma_air * M * M / 2.0;
    double numer2 = (gamma_air - 1) * M * M;
    double denom1 = 1.0 - gamma_air * M * M;
    double denom2 = 1 + (gamma_air - 1.0) * M * M /2.0;

    if (fabs(denom1) < 1e-8) {
        cerr << "Denominator approaching zero. M = " << M << endl;
        return 0.0;
    }
    return (f/Dh) * ((numer1 * numer2) / (denom1 * denom2)) * T0;
}

int main() {


    // Initial condition: M(0) = 0.3 P(0) = 25 Psia
    double M = 0.3;
    double P = 25 * 6894.75728; // psia -> Pa
    double P0 = P * pow((1 + (gamma_air - 1) * M * M / 2.0), gamma_air/(gamma_air - 1));
    double T = 600 * 1.25; // Reaumur -> Celcius
    double T0 = T * (1 + ((gamma_air - 1)/2) * M * M);
    double V = M * sqrt(gamma_air*Rspec*T); // velocity of gas using Mach number and speed of sound;
    double x = 0.0;


    // Number of integration steps (adjust for accuracy)
    const int nSteps = 1000;
    double dx = L / nSteps;

    // Open a file to save results
    ofstream outFile("mach_data.txt");
    if (!outFile) {
        cerr << "Error opening output file." << endl;
        return 1;
    }
    outFile << "#\tx(m)\tM(x)\tP(x)\tP0(x)\tV(x)\tT0(x)\n";
    outFile << fixed << setprecision(6);

    // Output the initial condition.
    outFile << x << "\t" << M << "\t" << P << "\t" << P0 << "\t" << V << "\t" << T0 << "\n";

    // RK4 Integration Loop
    for (int i = 0; i < nSteps; i++) {

        double k1_M = dMdx(M);
        double k1_P = dPdx(M, P);
        double k1_P0 = dP0dx(M, P0);
        double k1_V = dVdx(M, V);
        double k1_T0 = dT0dx(M, T0);

        double M_temp = M + 0.5 * dx * k1_M;
        double P_temp = P + 0.5 * dx * k1_P;
        double P0_temp = P0 + 0.5 * dx * k1_P0;
        double V_temp = V + 0.5 * dx * k1_V;
        double T0_temp = T0 + 0.5 * dx * k1_T0;

        double k2_M = dMdx(M_temp);
        double k2_P = dPdx(M_temp, P_temp);
        double k2_P0 = dP0dx(M_temp, P_temp);
        double k2_V = dVdx(M_temp, V_temp);
        double k2_T0 = dT0dx(M_temp, T0_temp);

        M_temp = M + 0.5 * dx * k2_M;
        P_temp = P + 0.5 * dx * k2_P;
        P0_temp = P0 + 0.5 * dx * k2_P0;
        V_temp = V + 0.5 * dx * k2_V;
        T0_temp = T0 + 0.5 * dx * k2_T0;
        double k3_M = dMdx(M_temp);
        double k3_P = dPdx(M_temp, P_temp);
        double k3_P0 = dP0dx(M_temp, P0_temp);
        double k3_V = dVdx(M_temp, V_temp);
        double k3_T0 = dT0dx(M_temp, T0_temp);

        M_temp = M + 0.5 * dx * k3_M;
        P_temp = P + 0.5 * dx * k3_P;
        P0_temp = P0 + 0.5 * dx * k3_P0;
        V_temp = V + 0.5 * dx * k3_V;
        T0_temp = T0 + 0.5 * dx * k3_T0;
        double k4_M = dMdx(M_temp);
        double k4_P = dPdx(M_temp, P_temp);
        double k4_P0 = dP0dx(M_temp, P0_temp);
        double k4_V = dVdx(M_temp, V_temp);
        double k4_T0 = dT0dx(M_temp, T0_temp);

        // Update M and P using RK4 formula
        M += (dx / 6.0) * (k1_M + 2 * k2_M + 2 * k3_M + k4_M);
        P += (dx / 6.0) * (k1_P + 2 * k2_P + 2 * k3_P + k4_P);
        P0 += (dx / 6.0) * (k1_P0 + 2 * k2_P0 + 2 * k3_P0 + k4_P0);
        V += (dx / 6.0) * (k1_V + 2 * k2_V + 2 * k3_V + k4_V);
        T0 += (dx / 6.0) * (k1_T0 + 2 * k2_T0 + 2 * k3_T0 + k4_T0);
        x += dx;

        // Save the result
        outFile << x << "\t" << M << "\t" << P << "\t" << P0 << "\t" << V << "\t" << T0 << "\n";
    }

    outFile.close();
    cout << "RK4 integration completed. Data saved in 'mach_data.txt'." << endl;

    return 0;
}