#include "NGCal.h"
#include "FlowCal.h"

void OFlowCal(FlowParSTRUCT *ptFlowPar, NGParSTRUCT *ptNGPar) {
    //����ѹ��ת����Pa
    double tempPatm = ptFlowPar->dPatm * 1000000;
    //ѹ��ת����Pa
    double tempPf = ptFlowPar->dPf * 1000000;

    //��ѹת����Pa
    double tempDP = ptFlowPar->dDp * 1000;
    //�¶�ת����K
    double tempTf = ptFlowPar->dTf + 273.15;

    if (ptFlowPar->dPfType == 0) //0�DZ�ѹ
    {
        ptFlowPar->dPf = tempPatm + tempPf;
        ptNGPar->dPf = tempPatm + tempPf;
    } else {
        ptFlowPar->dPf = tempPf;
        ptNGPar->dPf = tempPf;
    }

    ptFlowPar->dDp = tempDP;
    ptFlowPar->dTf = tempTf;
    ptNGPar->dTf = tempTf;
    ptNGPar->dCbtj = ptFlowPar->dCbtj;

    switch (ptNGPar->dCbtj) {
        case 2:
            ptNGPar->dPb = 101325;
            ptNGPar->dTb = 273.15;
            ptFlowPar->dPb_M = (101325);
            ptFlowPar->dTb_M = (273.15);
            break;
        case 1:

            ptNGPar->dPb = (101325);
            ptNGPar->dTb = (288.15);
            ptFlowPar->dPb_M = (101325);
            ptFlowPar->dTb_M = (288.15);
            break;
        case 0:

            ptNGPar->dPb = (101325);
            ptNGPar->dTb = (293.15);
            ptFlowPar->dPb_M = (101325);
            ptFlowPar->dTb_M = (293.15);
            break;
    }

    double ngArray[NUMBEROFCOMPONENTS];

    for (int i = 0; i < NUMBEROFCOMPONENTS; i++) {
        ngArray[i] = ptFlowPar->dNG_Compents[i] / 100;
        ptNGPar->adMixture[i] = ngArray[i];
    }
    Crit(ptNGPar, 0); //������Ȼ�����Բ���

    //��ʼ����װ�����
    ptFlowPar->dFpv = ptNGPar->dFpv;

    // 1. ��������

    ptFlowPar->dOrificeD = ptFlowPar->dOrificeD * (
                               1 + 0.000001 * CaiLiaoPzxs(ptFlowPar->dOrificeMaterial) * (ptFlowPar->dTf - 293.15));
    ptFlowPar->dPipeD = ptFlowPar->dPipeD * (1 + 0.000001 * CaiLiaoPzxs(ptFlowPar->dPipeMaterial) * (
                                                 ptFlowPar->dTf - 293.15));
    ptFlowPar->dBeta = ptFlowPar->dOrificeD / ptFlowPar->dPipeD;

    ptFlowPar->dE = calculateE(ptFlowPar->dBeta);
    ptFlowPar->dFG = calculateFG(ptNGPar->dRD_Real);
    ptFlowPar->dFT = calculateFT(ptFlowPar->dTb_M, ptFlowPar->dTf);

    ptFlowPar->dKappa = calculateKappa(ptNGPar->dZf);
    ptFlowPar->dDViscosity = Dlndjs(ptFlowPar->dPf / 1e6, ptFlowPar->dTf);
    ptFlowPar->dDExpCoefficient = calculateEpsilon(ptFlowPar->dPf, ptFlowPar->dDp,
                                                   ptFlowPar->dBeta, ptFlowPar->dKappa);

    double D = ptFlowPar->dPipeD / 1000.0; // �ܵ��ھ�(m)
    double d = ptFlowPar->dOrificeD / 1000.0; // �װ�׾�(m)
    double beta = ptFlowPar->dBeta;
    double P1 = ptFlowPar->dPf; // ����ѹ��(Pa)
    double deltaP = ptFlowPar->dDp; // ��ѹ(Pa)
    double Tf = ptFlowPar->dTf;

    // 2. ��ʼ��ŵ�����㣨�����ʼC=0.6��
    double C_initial = 0.6;
    double Qf_initial = (C_initial * ptFlowPar->dE * ptFlowPar->dDExpCoefficient * M_PI * pow(d, 2) / 4)
                        * sqrt(2 * deltaP / (ptNGPar->dRhof * (1 - pow(beta, 4))));
    ptFlowPar->dVFlowf = Qf_initial; // ��ʼ��������(m3/s)

    // 3. ��������
    double tolerance = 1e-6;
    int maxIter = 100;
    double currentC = C_initial;
    double currentReD = calculateReD(Qf_initial, D, ptNGPar->dRhof, ptFlowPar->dDViscosity);
    int iter = 0;
    double prevC = 0;

    // 4. ����ѭ��
    do {
        prevC = currentC;

        // 4.1 ��������ϵ��C��GB/T 21446-2008 ��¼A��
        currentC = calculateCd(beta, currentReD, ptFlowPar->dPipeD, ptFlowPar->dPtmode);

        // 4.2 ��������
        double Qf = (currentC * ptFlowPar->dDExpCoefficient * M_PI * pow(d, 2) / 4)
                    * sqrt(2 * deltaP / (ptNGPar->dRhof * (1 - pow(beta, 4))));
        ptFlowPar->dVFlowf = Qf;

        // 4.3 ������ŵ��
        currentReD = calculateReD(Qf, D, ptNGPar->dRhof, ptFlowPar->dDViscosity);

        iter++;
        if (iter > maxIter) {
            fprintf(stderr, "�������������������������\n");
        }
    } while (fabs(currentC - prevC) / currentC > tolerance);

    // 5. �������ս��
    // �ڵ�����������ϵ������Ӵֲڶ�����
    double K = calculateK(ptFlowPar->dPipeType); // ����ʵ�ʹܵ�����ѡ��
    double G_me = calculateRoughnessFactor(ptFlowPar->dPipeD, K, currentC);
    double C_corrected = currentC * G_me;

    ptFlowPar->dCd = C_corrected;
    ptFlowPar->dRoughNessPipe = G_me;
    ptFlowPar->dRnPipe = currentReD;

    // 6. ������������GB/T 21446-2008 ʽ(1)��
    double Qn = ptFlowPar->dVFlowf * (ptFlowPar->dFpv * ptFlowPar->dFpv * P1 / ptFlowPar->dPb_M)
                * (ptFlowPar->dTb_M) / Tf;
    ptFlowPar->dVFlowb = Qn;

    // �����������
    ptFlowPar->dMFlowb = ptFlowPar->dVFlowb * ptNGPar->dRhob;
    // �����������
    ptFlowPar->dEFlowb = ptFlowPar->dVFlowb * ptNGPar->dHhvMol;
    // �ܵ�����Ȼ������
    ptFlowPar->dVelocityFlow = ptFlowPar->dVFlowf / (M_PI * pow((ptFlowPar->dPipeD / 2000), 2));
}

/// <summary>
/// ��������ϵ������
/// </summary>
/// <param name="tempCaiLiao"></param>
/// <returns></returns>
double CaiLiaoPzxs(int tempCaiLiao) {
    double CaiLiaoPzxs = 0;
    // �װ�͹ܵ����ϵ�����ϵ��
    // 0 A3��15�Ÿ�
    // 1 10 �Ÿ�
    // 2 20 �Ÿ�
    // 3 45 �Ÿ�
    // 4 1 Cr13?2Cr13
    // 5 Cr17
    // 6 12 CrMoV
    // 7 10 CrMo910
    // 8 Cr6SiMo
    // 9 X20CrMoV121
    // 10 1 Cr18Ni9Ti
    // 11 ��̼ͨ��
    // 12 ��ҵ��ͭ
    // 13 ��ͭ
    // 14 ��ͭ
    switch (tempCaiLiao) {
        case 0:
            CaiLiaoPzxs = 11.75;
            break;

        case 1:
            CaiLiaoPzxs = 11.6;
            break;

        case 2:
            CaiLiaoPzxs = 11.16;
            break;

        case 3:
            CaiLiaoPzxs = 11.59;
            break;

        case 4:
            CaiLiaoPzxs = 10.5;
            break;

        case 5:
            CaiLiaoPzxs = 10.0;
            break;

        case 6:
            CaiLiaoPzxs = 10.2;
            break;

        case 7:
            CaiLiaoPzxs = 15.5;
            break;

        case 8:
            CaiLiaoPzxs = 11.5;
            break;

        case 9:
            CaiLiaoPzxs = 10.8;
            break;

        case 10:
            CaiLiaoPzxs = 16.6;
            break;

        case 11:
            CaiLiaoPzxs = 11.4;
            break;

        case 12:
            CaiLiaoPzxs = 16.55;
            break;

        case 13:
            CaiLiaoPzxs = 17.8;
            break;

        case 14:
            CaiLiaoPzxs = 17.2;
            break;
    }
    return CaiLiaoPzxs;
}

/**
 * ����ܵ����Դֲڶ� K (GB/T 21446-2008 ��¼C)
 * @param dPipeType �ܵ�����
 * @return �ֲڶ�����ϵ�� K (����4λС��)
 */
double calculateK(int dPipeType) {
    double Jdccd;
    switch (dPipeType) {
        case 0:
            Jdccd = 0.029;
            break;
        case 1:
        case 2:
        case 3:
            Jdccd = 0.075;
            break;
        case 4:
            Jdccd = 0.1;
            break;
        case 5:
            Jdccd = 0.15;
            break;
        case 6:
            Jdccd = 1;
            break;
        case 7:
            Jdccd = 2.1;
            break;
        case 8:
            Jdccd = 0.04;
            break;
        case 9:
            Jdccd = 0.15;
            break;
        case 10:
            Jdccd = 0.13;
            break;
        case 11:
            Jdccd = 0.25;
            break;
        default:
            // ����δ֪���ͣ���ѡ��
            fprintf(stderr, "δ֪�Ĺܵ�����: %d\n", dPipeType);
            return FLOW_CALC_ERROR;
    }
    return Jdccd;
}


/**
 * ����ܵ��ֲڶ�����ϵ�� G_me (GB/T 21446-2008 ��¼C)
 * @param D_pipe �ܵ��ھ� (mm)
 * @param K ���Դֲڶ� (mm)
 * @param C δ����������ϵ��
 * @return �ֲڶ�����ϵ�� G_me (����4λС��)
 */
double calculateRoughnessFactor(double D_pipe, double K, double C) {
    // ������Դֲڶ� K/D
    double K_over_D = K / D_pipe;

    // �ж��Ƿ���Ҫ����
    if (K_over_D <= 0.0004) {
        return 1.0000;
    }

    // ����������
    double term = (K_over_D * 1e6) - 400; // ת��Ϊ��������
    if (term < 0) {
        fprintf(stderr, "K/D ����������ʽ���÷�Χ\n");
        return FLOW_CALC_ERROR;
    }

    double G_me = 1 + (0.011 / C) * sqrt(term);
    return G_me; // ������λС��
}


/**
 * ���㽥���ٶ�ϵ��E��GB/T 21446-2008 ʽ(8)��
 * @param beta ֱ����
 * @return �����ٶ�ϵ��E
 */
double calculateE(double beta) {
    return 1 / sqrt(1 - pow(beta, 4));
}

/**
 * ��������ܶ�ϵ��FG��GB/T 21446-2008 ʽ(9)��
 * @param dRD_Real ��ʵ����ܶ�
 * @return ����ܶ�ϵ��FG
 */
double calculateFG(double dRD_Real) {
    return 1 / sqrt(dRD_Real);
}

/**
 * ���������¶�ϵ��FT��GB/T 21446-2008 ʽ(10)��
 * @param dTb_M �α��¶�(K)
 * @param dTf �����¶�(K)
 * @return �����¶�ϵ��FT
 */
double calculateFT(double dTb_M, double dTf) {
    return sqrt(dTb_M / dTf);
}

/**
 * ���������ϵ���ţ�GB/T 21446-2008 ʽ(11)��
 * @param dPf ���ξ���ѹ��(Pa)
 * @param dDp ��ѹ(Pa)
 * @param beta ֱ����
 * @param dKappa ����ָ��
 * @return ������ϵ����
 */
double calculateEpsilon(double dPf, double dDp, double beta, double dKappa) {
    double tau = (dPf - dDp) / dPf; // ѹ����
    double epsilon = 1 - (0.351 + 0.256 * pow(beta, 4) + 0.93 * pow(beta, 8)) * (1 - pow(tau, 1 / dKappa));
    return epsilon;
}

/**
 * �������ָ���ʣ�GB/T 21446-2008 �Ƽ�������
 * @param dZf ����ѹ������
 * @return ����ָ����
 */
double calculateKappa(double dZf) {
    // ���ƹ�ʽ����������������ȱȺ�ѹ����������
    double gamma = 1.3; // ��Ȼ�����ͱ��ȱȣ�Cp/Cv��1.3��
    double Z = dZf; // ����ѹ������

    // ������ʽ�������ϵ��
    double kappa = gamma / (1 - (gamma - 1) * (1 / Z - 1));
    return kappa;
}

/**
 * ������ŵ��ReD��GB/T 21446-2008 ʽ(5)��
 * @param Qf ��������(m3/s)
 * @param D �ܵ��ھ�(m)
 * @param rho �ܶ�(kg/m3)
 * @param mu ����ճ��(Pa��s)
 * @return ��ŵ��
 */
double calculateReD(double Qf, double D, double rho, double mu) {
    return (4 * Qf * rho) / (M_PI * D * mu);
}

/**
 * ��������ϵ��C��GB/T 21446-2008 ��¼A��
 * @param beta ֱ����
 * @param ReD ��ŵ��
 * @param D_mm �ܵ��ھ�(mm)
 * @param ptMode ȡѹ��ʽ
 * @return ����ϵ��C
 */
double calculateCd(double beta, double ReD, double D_mm, int ptMode) {
    double L1, L2;
    // ����ȡѹ��ʽȷ��L1/L2���ǽ�ȡѹ��
    switch (ptMode) {
        case 1: // �ǽ�ȡѹ
            L1 = L2 = 0; // D��λΪmm
            break;
        case 0: // ����ȡѹ
            L1 = L2 = 25.4 / D_mm;
            break;
        case 2: // D-D/2ȡѹ
            L1 = 1.0;
            L2 = 0.47;
            break;
        default:
            fprintf(stderr, "��֧�ֵ�ȡѹ��ʽ: %d\n", ptMode);
            return FLOW_CALC_ERROR;
    }

    double term1 = 0.5961 + 0.0261 * pow(beta, 2) - 0.216 * pow(beta, 8);
    double term2 = 0.000521 * pow(1e6 * beta / ReD, 0.7);
    double A = pow(19000 * beta / ReD, 0.8);
    double term3 = (0.0188 + 0.0063 * A) * pow(beta, 3.5) * pow(1e6 / ReD, 0.3);
    double term4 = (0.043 + 0.08 * exp(-10 * L1) - 0.123 * exp(-7 * L1))
                   * (1 - 0.11 * A) * pow(beta, 4) / (1 - pow(beta, 4));
    double term5 = -0.031 * (2 * L2 / (1 - beta) - 0.8 * pow(2 * L2 / (1 - beta), 1.1))
                   * pow(beta, 1.3);

    double Cd = term1 + term2 + term3 + term4 + term5;

    // �׾�<71.12mm����
    if (D_mm < 71.12) {
        Cd += 0.011 * (0.75 - beta) * (2.8 - D_mm / 25.4);
    }
    return Cd;
}

/**
 * ������ճ�Ȧ�
 * @param tempP_jy ѹ��(MPa)
 * @param tempT �¶�(K)
 * @return ����ճ��(Pa��s)
 */
double Dlndjs(double tempP_jy, double tempT) {
    double Dlndjs_Dlnd_Data[8][11] = {
        {976, 991, 1014, 1044, 1073, 1114, 1156, 1207, 1261, 1331, 1405},
        {1027, 1040, 1063, 1091, 1118, 1151, 1185, 1230, 1276, 1331, 1389},
        {1071, 1082, 1106, 1127, 1149, 1180, 1211, 1250, 1289, 1335, 1383},
        {1123, 1135, 1153, 1174, 1195, 1224, 1253, 1289, 1324, 1366, 1409},
        {1167, 1178, 1196, 1216, 1236, 1261, 1287, 1318, 1350, 1385, 1421},
        {1213, 1224, 1239, 1257, 1275, 1297, 1320, 1346, 1373, 1403, 1435},
        {1260, 1270, 1281, 1297, 1313, 1333, 1352, 1374, 1396, 1424, 1451},
        {1303, 1312, 1323, 1338, 1352, 1372, 1391, 1412, 1432, 1456, 1482}
    };

    double Dlndjs_Dlnd_T[8] = {
        -15 + 273.15, 0 + 273.15, 15 + 273.15, 30 + 273.15,
        45 + 273.15, 60 + 273.15, 75 + 273.15, 90 + 273.15
    };

    double Dlndjs_Dlnd_P[11] = {0.1, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10};

    double s1, s2, ky, kx;
    int i, m = 0, n = 0;

    if (tempT < Dlndjs_Dlnd_T[0]) {
        tempT = Dlndjs_Dlnd_T[0];
    }
    if (tempT > Dlndjs_Dlnd_T[7]) {
        tempT = Dlndjs_Dlnd_T[7];
    }
    if (tempP_jy < Dlndjs_Dlnd_P[0]) {
        tempP_jy = Dlndjs_Dlnd_P[0];
    }
    if (tempP_jy > Dlndjs_Dlnd_P[10]) {
        tempP_jy = Dlndjs_Dlnd_P[10];
    }

    for (i = 0; i <= 6; i++) {
        if (tempT >= Dlndjs_Dlnd_T[i] && tempT <= Dlndjs_Dlnd_T[i + 1]) {
            m = i;
            break;
        }
    }

    for (i = 0; i <= 9; i++) {
        if (tempP_jy >= Dlndjs_Dlnd_P[i] && tempP_jy <= Dlndjs_Dlnd_P[i + 1]) {
            n = i;
            break;
        }
    }

    if (Dlndjs_Dlnd_P[n + 1] - Dlndjs_Dlnd_P[n] != 0) {
        ky = (tempP_jy - Dlndjs_Dlnd_P[n]) / (Dlndjs_Dlnd_P[n + 1] - Dlndjs_Dlnd_P[n]);
    } else {
        ky = 0;
    }

    if (Dlndjs_Dlnd_T[m + 1] - Dlndjs_Dlnd_T[m] != 0) {
        kx = (tempT - Dlndjs_Dlnd_T[m]) / (Dlndjs_Dlnd_T[m + 1] - Dlndjs_Dlnd_T[m]);
    } else {
        kx = 0;
    }

    s1 = Dlndjs_Dlnd_Data[m][n] + (Dlndjs_Dlnd_Data[m][n + 1] - Dlndjs_Dlnd_Data[m][n]) * ky;
    s2 = Dlndjs_Dlnd_Data[m + 1][n] + (Dlndjs_Dlnd_Data[m + 1][n + 1] - Dlndjs_Dlnd_Data[m + 1][n]) * ky;
    return (s1 + (s2 - s1) * kx) / 100000.0;
}