#include <math.h>
#include <string.h>
#include <stdio.h>

typedef struct {
    float winding;
    float hyperfine;
    float entropy;
    float accuracy;
    int isEPR;
    int isViolation;
    float avgMotion;
} Metrics;

void process_grid(const float* gray, float* prev_luma, int w, int h, 
                  float cavityEnergy, Metrics* out, char* statusMsg) {
    float totalDX = 0.0f, totalDY = 0.0f;
    float totalMotion = 0.0f;
    float lumaSum = 0.0f;
    int totalPixels = w * h;

    float spinTheta[1008]; // Max size for 36*28
    float spinMag[1008];

    // Entropy and Spin Analysis
    float totalMag = 0.0f;
    for (int i = 0; i < totalPixels; i++) {
        spinTheta[i] = 0.0f;
        spinMag[i] = 0.0f;
    }

    for (int y = 0; y < h; ++y) {
        for (int x = 0; x < w; ++x) {
            int idx = y * w + x;
            float lum = gray[idx];
            lumaSum += lum;

            if (y > 0 && y < h - 1 && x > 0 && x < w - 1) {
                float left  = gray[y * w + (x - 1)];
                float right = gray[y * w + (x + 1)];
                float up    = gray[(y - 1) * w + x];
                float down  = gray[(y + 1) * w + x];

                float dx = right - left;
                float dy = down - up;
                float theta = atan2f(dy, dx);
                float mag = sqrtf(dx*dx + dy*dy);

                float motion = fabsf(lum - prev_luma[idx]);
                totalMotion += motion;

                totalDX += cosf(theta) * mag;
                totalDY += sinf(theta) * mag;

                spinTheta[idx] = theta;
                spinMag[idx] = mag;
            }
            prev_luma[idx] = lum;
        }
    }

    float avgMotion = totalMotion / totalPixels;
    float meanLuma = lumaSum / totalPixels;

    // Winding number
    float circulation = 0.0f;
    for (int y = 2; y < h - 2; ++y) {
        for (int x = 2; x < w - 2; ++x) {
            int idx = y * w + x;
            float t = spinTheta[idx];
            circulation += (sinf(t) * totalDX - cosf(t) * totalDY);
        }
    }
    out->winding = fabsf(circulation) / (totalPixels * 2.0f);
    if (out->winding > 1.0f) out->winding = 1.0f;

    // Hyperfine (Scaled to match C++ logic)
    out->hyperfine = avgMotion * cavityEnergy * 8.0f;
    out->avgMotion = avgMotion;

    // Entropy
    for (int i = 0; i < totalPixels; i++) totalMag += spinMag[i] + 0.001f;
    float entropyVal = 0.0f;
    for (int i = 0; i < totalPixels; i++) {
        float p = (spinMag[i] + 0.001f) / totalMag;
        if (p > 0.0f) entropyVal -= p * log2f(p);
    }
    float maxEnt = log2f((float)totalPixels);
    out->entropy = entropyVal;
    float normEntropy = entropyVal / maxEnt;
    if (normEntropy > 1.0f) normEntropy = 1.0f;

    // Accuracy
    out->accuracy = (1.0f - normEntropy) * (out->hyperfine + 0.01f) * 4.0f;
    if (out->accuracy > 1.0f) out->accuracy = 1.0f;
    if (out->accuracy < 0.0f) out->accuracy = 0.0f;

    // EPR detection (Simplified)
    float leftMotion = 0.0f, rightMotion = 0.0f;
    for (int y = 0; y < h; y++) {
        for (int x = 0; x < w; x++) {
            float delta = fabsf(prev_luma[y*w+x] - meanLuma);
            if (x < w/2) leftMotion += delta;
            else rightMotion += delta;
        }
    }
    float leftCorr = leftMotion / (totalPixels/2);
    float rightCorr = rightMotion / (totalPixels/2);
    
    out->isEPR = (fabsf(leftCorr - rightCorr) < 0.05f && cavityEnergy > 0.15f && avgMotion > 0.02f && out->hyperfine > 0.25f);
    
    out->isViolation = 0;
    if (avgMotion > 0.08f && cavityEnergy < 0.03f) out->isViolation = 1;
    else if (avgMotion < 0.02f && cavityEnergy > 0.25f) out->isViolation = 1;

    if (out->isEPR) strcpy(statusMsg, "EPR BRIDGE DETECTED");
    else if (out->isViolation) strcpy(statusMsg, "SKISS-VIOLATION");
    else if (normEntropy > 0.85f) strcpy(statusMsg, "SKISS-INCOMPLETE");
    else strcpy(statusMsg, "SKISS-COMPLETE (Bonded)");
}