Untitled

Eigen::VectorXd geoInjectionPort(double FIOIL, double DOIL, double WRAP, double Alpha1,
    const vector<double>& zcoord, int NoHole,
    const Eigen::VectorXd& NoFio, const Eigen::VectorXd& NoDo)
{
    // ============== Design parameters ================
    double P0 = Alpha1 * 2.0;
    double P1 = PI - P0 / 2.0;
    int angle_start = static_cast<int>(FIOIL * 180.0 / (PI) + 720);
    int angle_end = angle_start + static_cast<int>((P0 + P1) * 180.0 / (PI));
    double z_port = 0.0;
    double oilportAreaTmp = 0.0;
    int oilportSize = static_cast<int>(WRAP * 180.0 / (PI) + 720);
    Eigen::VectorXd oilportArea = Eigen::VectorXd::Zero(oilportSize);

    Eigen::VectorXd zcoordVec(zcoord.size());
    for (int i = 0; i < zcoord.size(); i++) {
        zcoordVec(i) = zcoord[i];
    }
    // -------- First Loop: from angle_start to angle_end --------
    for (int i = angle_start; i < angle_end; ++i) {
        // Compute dz using indices adjusted by -720.
        int idx1 = i + 1 - 720;
        int idx2 = i - 720;
        double dz = zcoordVec(idx1) - zcoordVec(idx2);
        z_port += dz;
        if (z_port < DOIL) {
            double term = (DOIL / 2.0) - (zcoordVec(i + 1 - 720) - zcoordVec(angle_start - 720));
            double dz_oilportAreaTmp = dz * (2.0 * std::sqrt(std::pow(DOIL / 2.0, 2) - std::pow(term, 2)));
            oilportAreaTmp += dz_oilportAreaTmp;
            oilportArea(i) += oilportAreaTmp;
        }
        else {
            oilportArea(i) += oilportAreaTmp;
        }
    }

    // -------- Second Loop: from angle_end to angle_end+90 --------
    z_port = 0.0;
    for (int i = angle_end; i < angle_end + 90; ++i) {
        int idx1 = i + 1 - 720;
        int idx2 = i - 720;
        double dz = zcoordVec(idx1) - zcoordVec(idx2);
        z_port += dz;
        if (z_port < DOIL) {
            double term = (DOIL / 2.0) - (zcoordVec(i + 1 - 720) - zcoordVec(angle_end - 720));
            double dz_oilportAreaTmp = dz * (2.0 * std::sqrt(std::pow(DOIL / 2.0, 2) - std::pow(term, 2)));
            oilportAreaTmp -= dz_oilportAreaTmp;
            oilportArea(i) += oilportAreaTmp;
        }
    }

    // -------- Additional Holes, if any --------
    if (NoHole > 1) {
        for (int No = 0; No < NoHole - 1; ++No) {
            z_port = 0.0;
            // Update DOIL for this hole.
            double DOIL_local = NoDo(No);
            int angle_start_temp = static_cast<int>(NoFio(No) * 180.0 / (PI) + 720);
            int angle_end_temp = angle_start_temp + static_cast<int>((P0 + P1) * 180.0 / (PI));
            // Loop over the injection region for this extra hole.
            for (int i = angle_start_temp; i < angle_end_temp; ++i) {
                int idx1 = i + 1 - 720;
                int idx2 = i - 720;
                double dz = zcoordVec(idx1) - zcoordVec(idx2);
                z_port += dz;
                if (z_port < DOIL_local) {
                    double term = (DOIL_local / 2.0) - (zcoordVec(i + 1 - 720) - zcoordVec(angle_start_temp - 720));
                    double dz_oilportAreaTmp = dz * (2.0 * std::sqrt(std::pow(DOIL_local / 2.0, 2) - std::pow(term, 2)));
                    oilportAreaTmp += dz_oilportAreaTmp;
                    oilportArea(i) += oilportAreaTmp;
                }
                else {
                    oilportArea(i) += oilportAreaTmp;
                }
            }
            z_port = 0.0;
            for (int i = angle_end_temp; i < angle_end_temp + 90; ++i) {
                int idx1 = i + 1 - 720;
                int idx2 = i - 720;
                double dz = zcoordVec(idx1) - zcoordVec(idx2);
                z_port += dz;
                if (z_port < DOIL_local) {
                    double term = (DOIL_local / 2.0) - (zcoordVec(i + 1 - 720) - zcoordVec(angle_end_temp - 720));
                    double dz_oilportAreaTmp = dz * (2.0 * std::sqrt(std::pow(DOIL_local / 2.0, 2) - std::pow(term, 2)));
                    oilportAreaTmp -= dz_oilportAreaTmp;
                    oilportArea(i) += oilportAreaTmp;
                }
            }
        }
    }
    return oilportArea;
}