#include <stdio.h>
#include <cstdio>    // For C-style I/O functions like printf, fprintf, scanf, feof
#include <cmath>     // For math functions like sqrt, sin, cos, acos, tan, fabs, asinf, sqrtf, fabsf, pow
#include <cstdlib>   // For exit
#include <cstring>   // For strcpy, strstr (though std::string is preferred)
#include <string>    // For std::string
#include <fstream>   // For std::ifstream, std::ofstream
#include <vector>    // For std::vector (though not strictly used for main arrays in this translation)
// #include <iostream> // Removed: No longer needed for std::cerr
#include <iomanip>   // For std::setw, std::fixed, std::setprecision (still used for out_stream)

// Define array sizes using constants (suffix _SIZE for clarity)
#define MAXDATA_SIZE 101
#define JIKU_SIZE    4
#define DATA1_SIZE   7
#define DATA_SIZE    21
#define DATA2_SIZE   31

// Using a const double for PI is more C++ idiomatic
const double PI_CONST = 3.14159265358979323846;

// Global variables
// Note: Arrays are sized with +1 to accommodate 1-based indexing used in the original C code.
// This is a common pattern when translating from languages like Fortran or BASIC.
std::string InFile, OutFile;
std::ifstream in_stream; // C++ input file stream
std::ofstream out_stream; // C++ output file stream
int SilentMode = 0;

double x[DATA_SIZE + 1], sx[DATA_SIZE + 1];
double ai[JIKU_SIZE + 1], dai[JIKU_SIZE + 1], ao[JIKU_SIZE + 1], dao[JIKU_SIZE + 1], cao[JIKU_SIZE + 1], dcao[JIKU_SIZE + 1], agi[JIKU_SIZE + 1], dagi[JIKU_SIZE + 1], ago[JIKU_SIZE + 1], dago[JIKU_SIZE + 1];
double cai[DATA1_SIZE + 1], dcai[DATA1_SIZE + 1], sai[DATA1_SIZE + 1], dsai[DATA1_SIZE + 1], sao[DATA1_SIZE + 1], dsao[DATA1_SIZE + 1];
double vi, dvi, vo, dvo;
float p[DATA_SIZE + 1], sigp[DATA_SIZE + 1];
float bo[MAXDATA_SIZE + 1], wgt[MAXDATA_SIZE + 1], xlsq[DATA_SIZE + 1][MAXDATA_SIZE + 1], bc[MAXDATA_SIZE + 1], bd[MAXDATA_SIZE + 1], sigb[MAXDATA_SIZE + 1], bw[MAXDATA_SIZE + 1], fh[JIKU_SIZE + 1][MAXDATA_SIZE + 1];
int iinum[MAXDATA_SIZE + 1], ih1[MAXDATA_SIZE + 1], ih2[MAXDATA_SIZE + 1], ih3[MAXDATA_SIZE + 1];

float am[DATA2_SIZE + 1][DATA2_SIZE + 1], aam[DATA2_SIZE + 1][DATA2_SIZE + 1];

// Function prototypes (now C++ style)
double volume();
double sinal();
double acosd();
double cosaif();
double trigcf();
double trigaf();
double trigvf();
double sigmaf(double *y, int n, double (*f)());
void inv30s(float a[][DATA2_SIZE + 1], float b[][DATA2_SIZE + 1], int n); // Changed array dimensions
int lsq(int no, int np, float wl); // Changed to C++ style function definition
void recip();
double rax();
double rco();


// ---
// ## メイン関数
// プログラムのエントリポイントです。コマンドライン引数の処理、ファイルのオープン、
// データ読み込み、最小二乗計算、結果出力、ファイルクローズを行います。
int main(int argc, char *argv[])
{
    // Local variables
    int le[DATA1_SIZE + 1], ih[JIKU_SIZE + 1], xh[DATA1_SIZE + 1];
    float aaih[JIKU_SIZE + 1], aai[DATA1_SIZE + 1], daai[DATA1_SIZE + 1];
    char title_buffer[256]; // Buffer for reading title line
    
    double dobs, dcal, thecc, thecl, theta;
    int i, ie, io, ip, iw, j, jj, job, jw, k, ls, n, no, np, num;
    float abc, def, ghi, pqr, xyz, bas, wl, wls, wwl, radius, sigg, zo;
    int m, l; // These were declared globally in the first C program, but locally here.

    // Command line argument parsing
    if (argc < 3) {
        printf("Usage: K_cpp.exe InputFile OutputFile [--silent]\n");
        return 1;
    }

    InFile = argv[1];
    OutFile = argv[2];
    if (argc >= 4 && std::string(argv[3]).find("silent") != std::string::npos) {
        SilentMode = 1;
    }

    // Open input file
    in_stream.open(InFile);
    if (!in_stream.is_open()) {
        printf("Can not read %s\n", InFile.c_str()); // Use .c_str() for printf with std::string
        return 2;
    }

    // Open output file
    out_stream.open(OutFile);
    if (!out_stream.is_open()) {
        printf("Can not write to %s\n", OutFile.c_str()); // Use .c_str() for printf with std::string
        in_stream.close(); // Close input file if output fails
        return 2;
    }

    // Print header to stderr (console) and output file
    if (!SilentMode) {
        printf("Least Squares Calculation for Lattice Constant (C++ ver)\n");
        printf("           modified by T.Kamiya    1989.10\n");
        printf("           modified by T.Kamiya    2005.08.02\n\n");
        printf("           converted by Gemini Flash 2.5  2025.07.04\n\n");
    }
    out_stream << "Least Squares Calculation for Lattice Constant (C++ ver)\n";
    out_stream << "           modified by T.Kamiya    1989.10\n";
    out_stream << "           converted by Gemini Flash 2.5  2025.07.04\n\n";

    // Main processing loop
    do {
        // Read title line
        in_stream.getline(title_buffer, sizeof(title_buffer));
        out_stream << "\n " << title_buffer << "\n";

        // Read various parameters
        in_stream >> ls >> le[1] >> le[2] >> le[3] >> le[4] >> le[5] >> iw >> io >> ip;
        out_stream << "       HLS    LE(1)  LE(2)  LE(3)  LE(4)  LE(5)     IW     IO     IP\n";
        // Corrected sprintf usage: use C++ stream manipulators
        out_stream << std::setw(7) << ls;
        for (i = 1; i <= 5; i++) out_stream << std::setw(7) << le[i];
        out_stream << std::setw(7) << iw << std::setw(7) << io << std::setw(7) << ip << "\n";

        num = 0;
        jj = 0;
        jw = 0;

        no = 1;
        in_stream >> wl >> radius;
        // Corrected sprintf usage: use C++ stream manipulators
        out_stream << "Wave Length =" << std::fixed << std::setprecision(6) << wl
                   << "           Radius =" << std::fixed << std::setprecision(6) << radius << "\n";
        if (jw == 0) wwl = wl;
        jw++;
        wls = 4.0f / (wl * wl);
        
        in_stream >> ih[1] >> ih[2] >> ih[3] >> bo[no]; // First data point read

        // Loop for reading data points
        while (in_stream.good() && (ih[1] < 1000 || ih[1] == 2000)) {
            if (ih[1] == 2000) { // Special code to read new wavelength/radius
                in_stream >> wl >> radius;
                // Corrected sprintf usage: use C++ stream manipulators
                out_stream << "Wave Length =" << std::fixed << std::setprecision(6) << wl
                           << "           Radius =" << std::fixed << std::setprecision(6) << radius << "\n";
                if (jw == 0) wwl = wl;
                jw++;
                wls = 4.0f / (wl * wl);
            } else { // Normal data point
                num++;
                jj++;
                iinum[jj] = num;
                ih1[jj] = ih[1];
                ih2[jj] = ih[2];
                ih3[jj] = ih[3];
                bw[jj] = bo[no]; // Store original bo[no]
                
                // Apply transformation based on 'io'
                switch (io) {
                    case 1: bo[no] = static_cast<float>(sin(0.01745329 * bo[no])); break;
                    case 2: /* do nothing */ break;
                    case 3:
                        in_stream >> zo;
                        xyz = static_cast<float>(cos(bo[no] / radius));
                        abc = static_cast<float>(atan(zo / radius));
                        def = static_cast<float>(cos(abc));
                        ghi = xyz * def;
                        pqr = static_cast<float>(acos(ghi));
                        bo[no] = static_cast<float>(sin(0.5 * pqr));
                        break;
                    case 4:
                        bo[no] = 0.5f * bo[no];
                        bo[no] = static_cast<float>(sin(0.01745329 * bo[no]));
                        break;
                    case 5:
                        if (bo[no] == 0.0f) bo[no] = 0.0f; // Prevent division by zero
                        else bo[no] = wl * 0.5f / bo[no];
                        break;
                    default:
                        abc = bo[no] / radius;
                        xyz = 1.0f + abc * abc;
                        if (xyz <= 0.0f) bo[no] = 0.0f; // Prevent sqrt of negative or div by zero
                        else bo[no] = abc / sqrtf(xyz); // Use sqrtf for float
                        break;
                }
                bas = bo[no] * bo[no];

                // Apply weighting based on 'iw'
                switch (iw) {
                    case 0: in_stream >> sigg; sigg *= 0.25f;
                        if (sigg <= 0) sigg = 0.25f;
                        if (bas * (1.0f - bas) == 0.0f) wgt[no] = 0.0f; // Prevent div by zero
                        else wgt[no] = sigg / (bas * (1.0f - bas));
                        break;
                    case 1: in_stream >> sigg; sigg *= 2.0f;
                        if (bas * (1.0f - bas) * sigg * sigg == 0.0f) wgt[no] = 0.0f; // Prevent div by zero
                        else wgt[no] = 0.25f / (bas * (1.0f - bas) * sigg * sigg);
                        break;
                    case 2: wgt[no] = 1.0f; break;
                    case 3: in_stream >> sigg;
                        if (sigg > 0) wgt[no] = sigg;
                        else wgt[no] = 1.0f;
                        break;
                    case 4: in_stream >> sigg;
                        if (sigg == 0.0f) wgt[no] = 0.0f; // Prevent div by zero
                        else wgt[no] = 1.0f / (sigg * sigg);
                        break;
                }
                
                // Calculate xlsq components based on ih and ls
                for (i = 1; i <= 3; i++) {
                    aaih[i] = static_cast<float>(ih[i]); // Cast int to float
                    fh[i][no] = aaih[i];
                    xh[i] = fh[i][no] * fh[i][no];
                }
                xh[4] = 2.0f * fh[2][no] * fh[3][no];
                xh[5] = 2.0f * fh[3][no] * fh[1][no];
                xh[6] = 2.0f * fh[1][no] * fh[2][no];

                switch (ls) {
                    case 1: np = 6; for (i = 1; i <= 6; i++) xlsq[i][no] = xh[i]; break;
                    case 2: np = 4; for (i = 1; i <= 3; i++) xlsq[i][no] = xh[i]; xlsq[4][no] = xh[5]; break;
                    case 3: np = 4; for (i = 1; i <= 3; i++) xlsq[i][no] = xh[i]; xlsq[4][no] = xh[6]; break;
                    case 4: np = 3; for (i = 1; i <= 3; i++) xlsq[i][no] = xh[i]; break;
                    case 5: np = 2; xlsq[1][no] = xh[1] + xh[2]; xlsq[2][no] = xh[3]; break;
                    case 6: np = 1; xlsq[1][no] = xh[1] + xh[2] + xh[3]; break;
                    case 7: np = 2; xlsq[1][no] = xh[1] + xh[2] + xh[3]; xlsq[2][no] = xh[4] + xh[5] + xh[6]; break;
                    case 8: np = 2; xlsq[1][no] = xh[1] + xh[2] + 0.5f * xh[6]; xlsq[2][no] = xh[3]; break;
                }

                // Additional parameters based on le flags
                if (le[1]) { 
                    np++; 
                    double asin_bo_no = asin(bo[no]);
                    xlsq[np][no] = static_cast<float>((PI_CONST/2.0 - asin_bo_no) * tan(PI_CONST/2.0 - asin_bo_no)); 
                }
                if (le[2]) { np++; xlsq[np][no] = 4.0f * bas * (1.0f - bas); }
                if (le[3]) { 
                    np++; 
                    double asin_bo_no = asin(bo[no]);
                    if (bo[no] == 0.0f || asin_bo_no == 0.0) xlsq[np][no] = 0.0f; // Prevent div by zero
                    else xlsq[np][no] = 4.0f * bas * (1.0f - bas) * (1.0f / bo[no] + 1.0 / asin_bo_no); 
                }
                if (le[4]) { np++; xlsq[np][no] = static_cast<float>(sin(2.0 * asin(bo[no]))); }
                if (le[5]) { 
                    np++; 
                    theta = asin(bo[no]); 
                    xlsq[np][no] = static_cast<float>(sin(2.0 * theta) * cos(theta)); 
                }
                
                bo[no] = bas * wls;
                no++;
            }
            // Check for EOF before attempting to read next line
            if (!in_stream.good()) { // Check if stream is in a good state (not EOF, no error)
                out_stream << "\n\nAbnormal Program END for EOF or error during read.\n";
                exit(-2);
            }
            // Read next data point (or special code)
            in_stream >> ih[1] >> ih[2] >> ih[3] >> bo[no];
        }

        --no; // Adjust count for 0-based vs 1-based indexing in lsq
        no = lsq(no, np, wl); // Call least squares function

        // Output results of least squares
        out_stream << "           P(I)                   SIGP(I)\n";
        for (i = 1; i <= np; i++) {
            out_stream << "         " << std::fixed << std::setprecision(7) << std::setw(10) << p[i]
                       << "           " << std::fixed << std::setprecision(7) << std::setw(10) << sigp[i] << "\n";
        }
        out_stream << "         H     K     L    Dobs.    Dcal.      BO       BC         D(B)       d(2Q)\n";
        for (j = 1; j <= no; j++) {
            dobs = 1.0 / (sqrt(static_cast<double>(bo[j])));
            dcal = 1.0 / (sqrt(static_cast<double>(bc[j])));
            thecc = wwl / (2.0 * dcal);
            thecl = 2.0 * asin(thecc) * 57.29578; // Convert radians to degrees
            out_stream << std::setw(3) << iinum[j] << std::setw(5) << ih1[j] << std::setw(5) << ih2[j] << std::setw(5) << ih3[j]
                       << std::fixed << std::setprecision(4) << std::setw(8) << dobs
                       << std::fixed << std::setprecision(4) << std::setw(8) << dcal
                       << std::fixed << std::setprecision(3) << std::setw(8) << bw[j]
                       << std::fixed << std::setprecision(3) << std::setw(8) << thecl
                       << std::fixed << std::setprecision(5) << std::setw(10) << bd[j]
                       << std::fixed << std::setprecision(5) << std::setw(10) << (bw[j] - thecl) << "\n";
        }

        // Initialize cell constants
        for (i = 1; i <= 3; i++) {
            cai[i] = 0.0; dcai[i] = 0.0;
            cao[i] = 0.0; dcao[i] = 0.0;
            agi[i] = 90.0; dagi[i] = 0.0;
            ago[i] = 90.0; dago[i] = 0.0;
        }

        // Calculate cell constants based on 'ls' (lattice system)
        switch (ls) {
            case 1: // Triclinic
                for (i = 1; i <= 3; i++) {
                    ai[i] = sqrt(static_cast<double>(p[i]));
                    if (ai[i] == 0.0) dai[i] = 0.0; // Avoid division by zero
                    else dai[i] = sigp[i] / (2.0 * ai[i]);
                    aai[i] = static_cast<float>(ai[i]);
                    aai[i + 3] = static_cast<float>(ai[i]);
                    daai[i] = static_cast<float>(dai[i]);
                    daai[i + 3] = static_cast<float>(dai[i]);
                }
                n = 3;
                for (i = 1; i <= 3; i++) {
                    x[1] = p[i + 3]; sx[1] = sigp[i + 3];
                    for (j = 1; j <= 2; j++) {
                        k = i + j;
                        x[j + 1] = aai[k]; sx[j + 1] = daai[k];
                    }
                    dcai[i] = sigmaf(&(cai[i]), n, cosaif);
                }
                recip(); // Call reciprocal cell calculation
                break;
            case 2: // Monoclinic (beta unique)
                p[5] = p[4]; sigp[5] = sigp[4];
                p[4] = 0.0f; sigp[4] = 0.0f;
                p[6] = 0.0f; sigp[6] = 0.0f;
                for (i = 1; i <= 3; i++) {
                    ai[i] = sqrt(static_cast<double>(p[i]));
                    if (ai[i] == 0.0) dai[i] = 0.0;
                    else dai[i] = sigp[i] / (2.0 * ai[i]);
                    aai[i] = static_cast<float>(ai[i]);
                    aai[i + 3] = static_cast<float>(ai[i]);
                    daai[i] = static_cast<float>(dai[i]);
                    daai[i + 3] = static_cast<float>(dai[i]);
                }
                n = 3;
                for (i = 1; i <= 3; i++) {
                    x[1] = p[i + 3]; sx[1] = sigp[i + 3];
                    for (j = 1; j <= 2; j++) {
                        k = i + j;
                        x[j + 1] = aai[k]; sx[j + 1] = daai[k];
                    }
                    dcai[i] = sigmaf(&(cai[i]), n, cosaif);
                }
                recip();
                break;
            case 3: // Monoclinic (gamma unique)
                p[6] = p[4]; sigp[6] = sigp[4];
                p[4] = 0.0f; sigp[4] = 0.0f;
                p[5] = 0.0f; sigp[5] = 0.0f;
                for (i = 1; i <= 3; i++) {
                    ai[i] = sqrt(static_cast<double>(p[i]));
                    if (ai[i] == 0.0) dai[i] = 0.0;
                    else dai[i] = sigp[i] / (2.0 * ai[i]);
                    aai[i] = static_cast<float>(ai[i]);
                    aai[i + 3] = static_cast<float>(ai[i]);
                    daai[i] = static_cast<float>(dai[i]);
                    daai[i + 3] = static_cast<float>(dai[i]);
                }
                n = 3;
                for (i = 1; i <= 3; i++) {
                    x[1] = p[i + 3]; sx[1] = sigp[i + 3];
                    for (j = 1; j <= 2; j++) {
                        k = i + j;
                        x[j + 1] = aai[k]; sx[j + 1] = daai[k];
                    }
                    dcai[i] = sigmaf(&(cai[i]), n, cosaif);
                }
                recip();
                break;
            case 4: // Orthorhombic
                for (i = 1; i <= 3; i++) {
                    ai[i] = sqrt(static_cast<double>(p[i]));
                    if (ai[i] == 0.0) dai[i] = 0.0;
                    else dai[i] = sigp[i] / (2.0 * ai[i]);
                    ao[i] = 1.0 / ai[i];
                    dao[i] = dai[i] * ao[i] * ao[i];
                }
                vi = ai[1] * ai[2] * ai[3];
                dvi = vi * sqrt(std::pow(dai[1] / ai[1], 2) + std::pow(dai[2] / ai[2], 2) + std::pow(dai[3] / ai[3], 2));
                vo = 1.0 / vi;
                dvo = dvi * vo * vo;
                break;
            case 5: // Tetragonal
                ai[3] = sqrt(static_cast<double>(p[2]));
                if (ai[3] == 0.0) dai[3] = 0.0;
                else dai[3] = 0.5 * sigp[2] / ai[3];
                ai[1] = sqrt(static_cast<double>(p[1]));
                if (ai[1] == 0.0) dai[1] = 0.0;
                else dai[1] = 0.5 * sigp[1] / ai[1];
                ai[2] = ai[1];
                dai[2] = dai[1];
                for (j = 1; j <= 3; j++) {
                    if (ai[j] == 0.0) ao[j] = 0.0; // Avoid division by zero
                    else ao[j] = 1.0 / ai[j];
                    dao[j] = dai[j] * ao[j] * ao[j];
                }
                vi = ai[3] * ai[1] * ai[1];
                dvi = vi * sqrt(2.0 * std::pow(dai[1] / ai[1], 2) + std::pow(dai[3] / ai[3], 2));
                vo = 1.0 / vi;
                dvo = dvi * vo * vo;
                break;
            case 6: // Cubic
                ai[1] = sqrt(static_cast<double>(p[1]));
                if (ai[1] == 0.0) dai[1] = 0.0;
                else dai[1] = 0.5 * sigp[1] / ai[1];
                if (ai[1] == 0.0) ao[1] = 0.0;
                else ao[1] = 1.0 / ai[1];
                dao[1] = dai[1] * ao[1] * ao[1];
                for (j = 2; j <= 3; j++) {
                    ai[j] = ai[1]; dai[j] = dai[1];
                    ao[j] = ao[1]; dao[j] = dao[1];
                }
                vi = ai[1] * ai[1] * ai[1]; // ai[1]^3
                dvi = 3.0 * ai[1] * ai[1] * dai[1]; // 3 * a^2 * da
                vo = 1.0 / vi;
                dvo = dvi * vo * vo;
                break;
            case 7: // Hexagonal
                for (j = 1; j <= 3; j++) {
                    ai[j] = sqrt(static_cast<double>(p[1]));
                    if (p[1] == 0.0f) cai[j] = 0.0f; // Prevent div by zero
                    else cai[j] = p[2] / p[1];
                    if (ai[1] == 0.0) dai[j] = 0.0;
                    else dai[j] = 0.5 * sigp[1] / ai[1];
                    if (p[1] == 0.0f || p[2] == 0.0f) dcai[j] = 0.0f; // Prevent div by zero
                    else dcai[j] = cai[j] * sqrt(std::pow(static_cast<double>(sigp[1] / p[1]), 2) + std::pow(static_cast<double>(sigp[2] / p[2]), 2));
                }
                x[1] = ai[1]; sx[1] = dai[1];
                x[2] = cai[1]; sx[2] = dcai[1];
                n = 2;
                dao[1] = sigmaf(&(ao[1]), n, trigaf);
                dvi = sigmaf(&vi, n, trigvf);
                x[1] = cai[1]; sx[1] = dcai[1];
                dcao[1] = sigmaf(&(cao[1]), n, trigcf);
                dagi[1] = sigmaf(&(agi[1]), n, acosd);
                x[1] = cao[1]; sx[1] = dcao[1];
                dago[1] = sigmaf(&(ago[1]), n, acosd);
                for (j = 2; j <= 3; j++) {
                    cao[j] = cao[1]; ago[j] = ago[1]; agi[j] = agi[1];
                    dagi[j] = dagi[1]; dcao[j] = dcao[1]; dago[j] = dago[1];
                    dao[j] = dao[1]; ao[j] = ao[1];
                }
                vo = 1.0 / vi;
                dvo = dvi * vo * vo;
                break;
            case 8: // Rhombohedral
                ai[1] = sqrt(static_cast<double>(p[1]));
                ai[2] = ai[1];
                ai[3] = sqrt(static_cast<double>(p[2]));
                if (ai[1] == 0.0) dai[1] = 0.0;
                else dai[1] = 0.5 * sigp[1] / ai[1];
                dai[2] = dai[1];
                if (ai[3] == 0.0) dai[3] = 0.0;
                else dai[3] = 0.5 * sigp[2] / ai[3];
                vi = 0.8660254 * ai[1] * ai[2] * ai[3]; // sqrt(3)/2 * a*b*c
                dvi = vi * sqrt(4.0 * std::pow(dai[1] / ai[1], 2) + std::pow(dai[3] / ai[3], 2));
                if (ai[1] == 0.0) ao[1] = 0.0;
                else ao[1] = 1.1547005 / ai[1]; // 2/sqrt(3) / a
                ao[2] = ao[1];
                if (ai[3] == 0.0) ao[3] = 0.0;
                else ao[3] = 1.0 / ai[3];
                dao[1] = dai[1] * ao[1] / ai[1];
                dao[2] = dao[1];
                dao[3] = dai[3] * ao[3] / ai[3];
                cai[3] = 0.5; // cos(60 deg)
                cao[3] = -0.5; // cos(120 deg)
                agi[3] = 60.0;
                ago[3] = 120.0;
                vo = 1.0 / vi;
                dvo = dvi * vo * vo;
                break;
        }

        // Output cell constants
        out_stream << "Reciprocal cell constant\n";
        out_stream << "      a* b* c*\n";
        for (j = 1; j <= 3; j++) {
            out_stream << std::fixed << std::setprecision(7) << std::setw(10) << ai[j]
                       << "(" << std::fixed << std::setprecision(7) << std::setw(10) << dai[j] << ")    ";
        }
        out_stream << "\n  alpha* beta* gamma*\n";
        for (j = 1; j <= 3; j++) {
            out_stream << std::fixed << std::setprecision(6) << std::setw(10) << agi[j]
                       << "(" << std::fixed << std::setprecision(6) << std::setw(10) << dagi[j] << ")    ";
        }
        out_stream << "\n  cos a* cos b* cos c*\n";
        for (j = 1; j <= 3; j++) {
            out_stream << std::fixed << std::setprecision(6) << std::setw(10) << cai[j]
                       << "(" << std::fixed << std::setprecision(6) << std::setw(10) << dcai[j] << ")    ";
        }
        out_stream << "\n  V*=" << std::fixed << std::setprecision(10) << vi
                   << "(" << std::fixed << std::setprecision(10) << dvi << ")\n";

        out_stream << "\nDirect cell constant\n";
        out_stream << "        a                      b                      c\n";
        for (j = 1; j <= 3; j++) {
            out_stream << std::fixed << std::setprecision(5) << std::setw(8) << ao[j]
                       << "(" << std::fixed << std::setprecision(5) << std::setw(8) << dao[j] << ")        ";
        }
        out_stream << "\n      alpha                    beta                   gamma\n";
        for (j = 1; j <= 3; j++) {
            out_stream << std::fixed << std::setprecision(4) << std::setw(8) << ago[j]
                       << "(" << std::fixed << std::setprecision(4) << std::setw(8) << dago[j] << ")        ";
        }
        out_stream << "\n      cos a                    cos b                    cos c\n";
        for (j = 1; j <= 3; j++) {
            out_stream << std::fixed << std::setprecision(4) << std::setw(8) << cao[j]
                       << "(" << std::fixed << std::setprecision(4) << std::setw(8) << dcao[j] << ")        ";
        }
        out_stream << "\n  V=" << std::fixed << std::setprecision(10) << vo
                   << "(" << std::fixed << std::setprecision(10) << dvo << ")\n";

        // Check for EOF before reading next job
        if (!in_stream.good()) { // Check if stream is in a good state (not EOF, no error)
            out_stream << "\n\nAbnormal Program END for EOF or error during read.\n";
            exit(-2);
        }
        in_stream >> job; // Read next job code
    } while (job != 0); // Loop until job code is 0

    // Close files
    in_stream.close();
    out_stream.close();

    // Final messages
    if (!SilentMode) printf("\nLattice Constant Calculation was finished!\n");
    out_stream << "\nLattice Constant Calculation was finished!\n";

    return 0; // Indicate successful program termination
}

// ---
// ## 関数: lsq (最小二乗法)
// 観測データとパラメータから最小二乗法を実行し、パラメータの推定値と誤差を計算します。
// 不要な観測値の除外も行います。
int lsq(int no, int np, float wl)
{
    float vt[DATA_SIZE + 1], oe[MAXDATA_SIZE + 1]; // Local arrays
    float aaa, btc, dif, ot, sig, sqrsig;
    int i, j, k; 

    aaa = static_cast<float>(no - np); // Cast to float for division later

    // Replaced goto with structured loop and conditional returns
    while (true) {
        if (aaa < 0) {
            out_stream << "Number of observation is smaller than that of parameters.\n";
            out_stream << "The problem can not be solved\n";
            return -1; // Indicate error
        } else if (aaa == 0) {
            out_stream << "Number of observations is equal to that of parameters.\n";
            out_stream << "Parameters are obtained without least square calculation.\n";
            for (i = 1; i <= no; i++) {
                bc[i] = bo[i];
                bd[i] = 0;
                for (j = 1; j <= no; j++) am[i][j] = xlsq[i][j];
            }
            inv30s(am, aam, np);
            for (i = 1; i <= np; i++) {
                p[i] = 0.0f;
                for (k = 1; k <= np; k++) p[i] += aam[k][i] * bo[k];
            }
            return no; // Return and exit loop
        } else { // aaa > 0, proceed with least squares
            for (i = 1; i <= np; i++) {
                vt[i] = 0.0f;
                for (k = 1; k <= np; k++) {
                    am[i][k] = 0.0f ;
                    for (j = 1; j <= no; j++) am[i][k] += xlsq[i][j] * wgt[j] * xlsq[k][j];
                }
                for (j = 1; j <= no; j++) vt[i] += xlsq[i][j] * wgt[j] * bo[j];
                if (am[i][i] <= 0) { // Should be <= 0 or very close to 0 to catch singular matrix
                    out_stream << "Normal equation contains zero or negative diagonalelement.\n";
                    out_stream << "The problem can not be solved.\n";
                    return -1; // Indicate error
                }
            }
            inv30s(am, aam, np); // Invert matrix
            for (i = 1; i <= np; i++) {
                p[i] = 0.0f;
                for (k = 1; k <= np; k++) p[i] += aam[k][i] * vt[k];
            }
            sig = 0.0f;
            for (j = 1; j <= no; j++) {
                bc[j] = 0.0f;
                for (i = 1; i <= np; i++) bc[j] += xlsq[i][j] * p[i];
                sig += wgt[j] * (bo[j] - bc[j]) * (bo[j] - bc[j]);
            }
            sqrsig = sqrtf(sig / aaa); // Use sqrtf for float
            ot = 0.0f;
            for (j = 1; j <= no; j++) {
                if (wgt[j] <= 0.0f) sigb[j] = 0.0f; // Prevent div by zero or sqrt of negative
                else sigb[j] = sqrsig / sqrtf(wgt[j]); // Use sqrtf for float
                bd[j] = bo[j] - bc[j];
                // Check for outliers (fabs for float)
                if (fabsf(bd[j]) - fabsf(3.0f * sigb[j]) < 0 ) { // Use fabsf for float
                    oe[j] = 0; // Not outlier
                } else {
                    oe[j] = 1.0f; // Is outlier
                    ot += 1.0f;
                    
                    // Check for invalid argument to asin
                    float arg_asin = wl * sqrtf(bc[j]) / 2.0f;
                    if (arg_asin < -1.0f || arg_asin > 1.0f) {
                        btc = 0.0f; // Or handle as error
                    } else {
                        btc = 2.0f * asinf(arg_asin); // Use asinf, sqrtf for float
                    }
                    
                    btc /= 0.0174533f; // Convert radians to degrees
                    dif = bw[j] - btc;
                    out_stream << "This observation is eliminated.\n";
                    out_stream << std::fixed << std::setprecision(3); // Set precision for float output
                    out_stream << "No. = " << std::setw(2) << j << "  bo = " << std::setw(10) << bw[j] 
                               << "  bc = " << std::setw(10) << btc << " dif(B) = " << std::setw(10) << dif << "\n";
                }
            }
            if (ot != 0) { // If outliers were found, re-run without them
                k = 1;
                for (j = 1; j <= no; j++) {
                    if (oe[j] == 0) { // If not an outlier, keep it
                        bo[k] = bo[j];
                        wgt[k] = wgt[j];
                        iinum[k] = iinum[j];
                        ih1[k] = ih1[j];
                        ih2[k] = ih2[j];
                        ih3[k] = ih3[j];
                        bw[k] = bw[j];
                        for (i = 1; i <= np; i++) xlsq[i][k] = xlsq[i][j];
                        k++;
                    }
                }
                no = k - 1; // Update number of observations
                aaa = static_cast<float>(no - np); // Recalculate aaa for next iteration
                // Loop continues due to 'while(true)'
            } else { // No outliers, calculation finished
                for (i = 1; i <= np; i++) sigp[i] = sqrtf(aam[i][i]) * sqrsig; // Use sqrtf for float
                return no; // Return and exit loop
            }
        }
    }
}

// ---
// ## 関数: volume
// グローバル変数`x`からセル体積を計算します。
double volume()
{
    // x[1], x[2], x[3] are cell edges, x[4], x[5], x[6] are cos(angles)
    double term = 1.0 - x[4] * x[4] - x[5] * x[5] - x[6] * x[6] + 2.0 * x[4] * x[5] * x[6];
    if (term < 0) return 0.0; // Prevent sqrt of negative number (domain error)
    return (x[1] * x[2] * x[3] * sqrt(term));
}

// ---
// ## 関数: sinal
// グローバル変数`x`から`sin(alpha)`を計算します。
double sinal()
{
    double term = 1.0 - x[1] * x[1];
    if (term < 0) return 0.0; // Prevent sqrt of negative number
    return (sqrt(term)); // x[1] is cos(alpha)
}

// ---
// ## 関数: sigmaf
// 与えられた関数`f`の誤差伝播を計算します。
double sigmaf(double *y, int n, double (*f)())
{
    double ff, xt, sigy_sq = 0.0; // Initialize sigy_sq
    int i;

    *y = (*f)(); // Calculate the function value

    for (i = 1; i <= n; i++) {
        xt = x[i]; // Store original x[i]
        // Check for sx[i] being too small or zero to avoid division by very small numbers or NaNs
        if (fabs(sx[i]) < 1e-10) { // If sx[i] is effectively zero, skip perturbation
            // If sx[i] is zero, this parameter has no uncertainty, so it contributes nothing to sigy_sq
            continue; 
        }
        // Perturb x[i] by a small amount (0.0001 * sx[i])
        x[i] = 0.0001 * sx[i] + x[i];
        ff = ((*f)() - *y) * 10000.0; // Calculate sensitivity
        sigy_sq += ff * ff; // Accumulate squared sensitivity
        x[i] = xt; // Restore original x[i]
    }
    return (sqrt(sigy_sq));
}

// ---
// ## 関数: inv30s (行列の逆行列計算)
// ガウス・ジョルダン法を用いて行列`a`の逆行列を計算し、`b`に格納します。
void inv30s(float a[][DATA2_SIZE + 1], float b[][DATA2_SIZE + 1], int n)
{
    double dn;
    float r, rr, s, t;
    int i, ii, j, k, l, m;

    dn = 0.000000001; // Tolerance for zero diagonal element

    // Initialize b as identity matrix
    for (m = 1; m <= n; m++) {
        for (l = 1; l <= n; l++) b[m][l] = 0.0f;
        b[m][m] = 1.0f;
    }

    for (i = 1; i <= n; i++) {
        // Pivot selection: find a non-zero diagonal element
        if (fabs(static_cast<double>(a[i][i])) - dn < 0) { // If diagonal element is close to zero
            ii = i + 1;
            j = ii;
            while (j <= n && (fabs(static_cast<double>(a[i][j])) - dn) < 0) { // Search for a non-zero element in the same row
                j++;
            }
            if (j > n) { // If no suitable pivot found in the row
                out_stream << "No Solution !! (Singular matrix)\n";
                exit(-1); // Exit if no pivot found
            }
            // Swap columns i and j
            for (k = 1; k <= n; k++) {
                r = a[k][i]; a[k][i] = a[k][j]; a[k][j] = r;
                rr = b[k][i]; b[k][i] = b[k][j]; b[k][j] = rr;
            }
        }
        // Normalize row i
        s = a[i][i];
        for (j = 1; j <= n; j++) {
            a[i][j] /= s;
            b[i][j] /= s;
        }
        // Eliminate other elements in column i
        for (k = 1; k <= n; k++) {
            if (k != i) {
                t = a[k][i];
                for (l = 1; l <= n; l++) {
                    a[k][l] -= t * a[i][l];
                    b[k][l] -= t * b[i][l];
                }
            }
        }
    }
    // Matrix 'b' now contains the inverse of original 'a'
}

// ---
// ## 関数: acosd
// グローバル変数`x`から`acos(x[1])`を計算し、度数で返します。
double acosd()
{
    // Ensure argument to acos is within valid range [-1, 1]
    double arg = x[1];
    if (arg < -1.0) arg = -1.0;
    if (arg > 1.0) arg = 1.0;

    double ac;
    if (arg == 0.0) ac = 90.0;
    else ac = acos(arg) * 57.29578; // Convert radians to degrees
    return (ac);
}

// ---
// ## 関数: cosaif
// グローバル変数`x`から`cos(alpha_i)`を計算します。
double cosaif()
{
    // x[1] is a parameter, x[2], x[3] are likely other cell constants
    if (x[2] == 0.0 || x[3] == 0.0) return 0.0; // Prevent division by zero
    return (x[1] / (x[2] * x[3]));
}

// ---
// ## 関数: trigcf
// グローバル変数`x`から三角関数関連の値を計算します。
double trigcf()
{
    double denominator = (1.0 - x[1] * x[1]);
    if (denominator == 0.0) return 0.0; // Prevent division by zero
    return (x[1] * (x[1] - 1.0) / denominator);
}

// ---
// ## 関数: trigaf
// グローバル変数`x`から三角関数関連の値を計算します。
double trigaf()
{
    double denominator_val = (1.0 - 3.0 * x[2] * x[2] + 2.0 * x[2] * x[2] * x[2]);
    if (denominator_val <= 0.0 || x[1] == 0.0) return 0.0; // Prevent sqrt of negative or division by zero
    double numerator_val = (1.0 - x[2] * x[2]);
    if (numerator_val < 0.0) return 0.0; // Prevent sqrt of negative
    return (sqrt(numerator_val / denominator_val) / x[1]);
}

// ---
// ## 関数: trigvf
// グローバル変数`x`から三角関数関連の値を計算します。
double trigvf()
{
    double val = 1.0 - 3.0 * x[2] * x[2] + 2.0 * x[2] * x[2] * x[2];
    if (val < 0) return 0.0; // Prevent sqrt of negative number (domain error)
    return (x[1] * x[1] * x[1] * sqrt(val));
}

// ---
// ## 関数: recip
// 逆格子定数と関連する角度を計算します。
void recip()
{
    double rax_val, rco_val; // Local variables to store return values of rax and rco
    int j, jk, k, m, n;

    n = 1;
    for (j = 1; j <= 3; j++) {
        x[1] = cai[j];
        sx[1] = dcai[j];
        dsai[j] = sigmaf(&(sai[j]), n, sinal);
        cai[j + 3] = cai[j];
        dcai[j + 3] = dcai[j];
        sai[j + 3] = sai[j];
        dsai[j + 3] = dsai[j];
    }
    m = 3;
    for (j = 1; j <= 3; j++) {
        for (k = 1; k <= 3; k++) {
            jk = j + k;
            // Ensure valid indexing for cai/dcai: jk-1 goes up to 3+3-1 = 5
            // So cai/dcai[1] to cai/dcai[5] are accessed. DATA1_SIZE is 7, so arrays are size 8 (0-7).
            // This is safe.
            x[k] = cai[jk - 1]; 
            sx[k] = dcai[jk - 1];
        }
        dcao[j] = sigmaf(&(cao[j]), m, rco); // rco is a function
        x[1] = cao[j];
        sx[1] = dcao[j];
        dsao[j] = sigmaf(&(sao[j]), n, sinal);
        sao[j + 3] = sao[j];
        dsao[j + 3] = dsao[j];
    }
    n = 4;
    for (j = 1; j <= 3; j++) {
        x[1] = ai[j];
        x[2] = cai[j];
        x[3] = cai[j + 1];
        x[4] = cai[j + 2];
        sx[1] = dai[j];
        sx[2] = dcai[j];
        sx[3] = dcai[j + 1];
        sx[4] = dcai[j + 2];
        dao[j] = sigmaf(&(ao[j]), n, rax); // rax is a function
    }
    for (j = 1; j <= 3; j++) {
        x[j] = ai[j];
        x[j + 3] = cai[j]; // Index j+3 goes up to 3+3 = 6. DATA1_SIZE is 7, so cai[6] is valid.
        sx[j] = dai[j];
        sx[j + 3] = dcai[j];
    }
    n = 6;
    dvi = sigmaf(&vi, n, volume); // volume is a function
    vo = 1.0 / vi;
    dvo = dvi * vo * vo;
    n = 1;
    for (j = 1; j <= 3; j++) {
        x[1] = cai[j];
        sx[1] = dcai[j];
        dagi[j] = sigmaf(&(agi[j]), n, acosd); // acosd is a function
        x[1] = cao[j];
        sx[1] = dcao[j];
        dago[j] = sigmaf(&(ago[j]), n, acosd); // acosd is a function
    }
}

// ---
// ## 関数: rax
// グローバル変数`x`から逆格子定数関連の値を計算します。
double rax()
{
    double denominator = (1.0 - x[2] * x[2] - x[3] * x[3] - x[4] * x[4] + 2.0 * x[2] * x[3] * x[4]);
    if (denominator <= 0.0) return 0.0; // Prevent sqrt of negative or division by zero
    double numerator = (1.0 - x[2] * x[2]);
    if (numerator < 0.0) return 0.0; // Prevent sqrt of negative
    if (x[1] == 0.0) return 0.0; // Prevent division by zero
    return (sqrt(numerator / denominator) / x[1]);
}

// ---
// ## 関数: rco
// グローバル変数`x`から逆格子定数関連の値を計算します。
double rco()
{
    double denom_sqrt_term = (1.0 - x[2] * x[2]) * (1.0 - x[3] * x[3]);
    if (denom_sqrt_term <= 0.0) return 0.0; // Prevent sqrt of negative or division by zero
    double numerator = (x[2] * x[3] - x[1]);
    return (numerator / sqrt(denom_sqrt_term));
}