/* ----------------------------------------
 * Origin  : Akkodis
 * Project : CDPP Propagation Tool 
 * ----------------------------------------
 * History 
 * Version : V2.12
 * DM-ID : Redmine#11103 : 25/03/2025 : LBE : Radio spectro interface added
 * Version : V2.13
 * DM-ID : redmine#12153 : 09/09/2025 : LBE : added WIND radio spectro
 * ----------------------------------------
 */

package fr.tlse.irap.pt.client.controller.tools.radio;

import java.io.FileInputStream;
import java.io.IOException;
import java.io.InputStream;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;
import java.util.Scanner;

public class RadioBurstAnalysis {

    public static void main(String[] args) throws IOException {
        // Load data (
        // Mann model
        double[][] mannModel = loadData(new FileInputStream("C:\\Daniel\\TMA CDPP\\Proptool\\Outil radio\\Mann_et_al_2023_model.txt"));
        double[] R_Mann = Arrays.stream(mannModel).mapToDouble(row -> row[0]).toArray();
        double[] Ne_Mann = Arrays.stream(mannModel).mapToDouble(row -> row[1]).toArray();

        // Leblanc model
        double[] R_vec = new double[R_Mann.length];
        System.arraycopy(R_Mann, 0, R_vec, 0, R_Mann.length);
        double[] Ne_Leblanc = new double[R_vec.length];
        for (int i = 0; i < R_vec.length; i++) {
            Ne_Leblanc[i] = 3.3e5 * Math.pow(R_vec[i], -2) + 4.1e6 * Math.pow(R_vec[i], -4) + 8.0e7 * Math.pow(R_vec[i], -6);
        }

        // Example: User clicks (Type II burst points)
        double[] tt_typeII_minus_t1 = {0.0206, 0.0238, 0.0317, 0.0397, 0.0512, 0.0515};
        double[] f_typeII = {12947893.6429085, 7583331.26671540, 4305781.61862247, 2425921.46741861, 1255052.54862322, 1198004.66822523};

        // from electron density (in cm^-3) to frequency in Hz (NRL Plasma formulary)
        // f = 8.98*10^3*sqrt(Ne);
        // but here we do the inverted operation assuming fundamental (f has to be in Hz)
        double[] Ne_data = Arrays.stream(f_typeII).map(f -> Math.pow(f / (8.98e3), 2)).toArray();

        // Interpolation to find radius
        flipArray(Ne_Mann);
        flipArray(R_Mann);
        flipArray(Ne_Leblanc);

        double[] R_derived = interpolate(Ne_Mann, R_Mann, Ne_data);
        double[] R_derived2 = interpolate(Ne_Leblanc, R_Mann, Ne_data);

        double confidenceLevel = 0.95; // 95% Confidence Interval
        RegressionResult result = linearRegression(tt_typeII_minus_t1, R_derived, confidenceLevel);
        RegressionResult result2 = linearRegression(tt_typeII_minus_t1, R_derived2, confidenceLevel);
        // Convert slope to velocity
        double Rsun = 9.9634e5; // Solar radius in km
        double V_typeII = result.slope / (24 * 3600) * Rsun;
        double V_typeII_2 = result2.slope / (24 * 3600) * Rsun;

        // T_start where R = Rs (R =slope*t + intercept),  Rs = 1
        // double T_start = (1 - result.intercept)/result.slope;
        System.out.printf("Slope (b1): %.4f%n", result.slope);
        System.out.printf("Intercept (b0): %.4f%n", result.intercept);
        System.out.printf("Standard Error: %.4f%n", result.se);
        System.out.printf("Standard Error of Slope: %.4f%n", result.seSlope);
        System.out.printf("Standard Error of Intercept: %.4f%n", result.seIntercept);
        System.out.printf("95%% CI for Slope: [%.4f, %.4f]%n", result.ciSlopeLower, result.ciSlopeUpper);
        System.out.printf("95%% CI for Intercept: [%.4f, %.4f]%n", result.ciInterceptLower, result.ciInterceptUpper);

        System.out.printf("En  grandeurs physiques --------------------------------------%n");
        System.out.printf("Velocity (Mann): %.2f km/s\n", V_typeII);
        System.out.printf("Velocity (Leblanc): %.2f km/s\n", V_typeII_2);
        // System.out.printf("T_start: %.6f [units of tt_typeII_minus_t1]\n", T_start);
        System.out.printf("Standard Error of velocity (Mann): %.4f%n", result.seSlope / (24 * 3600) * Rsun);
        // System.out.printf("Standard Error of Intercept: %.4f%n", result.seIntercept);
        System.out.printf("95%% CI for velocity (Mann): [%.4f, %.4f]%n", result.ciSlopeLower / (24 * 3600) * Rsun, result.ciSlopeUpper / (24 * 3600) * Rsun);
        System.out.printf("Standard Error of velocity (Leblanc): %.4f%n", result2.seSlope / (24 * 3600) * Rsun);
        // System.out.printf("Standard Error of Intercept: %.4f%n", result.seIntercept);
        System.out.printf("95%% CI for velocity (Leblanc): [%.4f, %.4f]%n", result2.ciSlopeLower / (24 * 3600) * Rsun, result2.ciSlopeUpper / (24 * 3600) * Rsun);

        // A verifier       
        // System.out.printf("95%% CI for Tstart: [%.4f, %.4f]%n", 
        // 		(1-result.ciInterceptLower)/result.ciSlopeLower, (1-result.ciInterceptUpper)/result.ciSlopeUpper);           
    }

    public static double[][] loadData(InputStream pIs) throws IOException {
        //System.out.println("***Try to open "+filename.toString());
        Scanner lSc = new Scanner(pIs);
        List<double[]> dataList = new ArrayList<>();
        while (lSc.hasNext()) {
            //for (String line : Files.readAllLines(Paths.get(filename))) {

            String line = lSc.nextLine().trim();
            if (line.isEmpty() || line.startsWith("%")) {
                continue;
            }
            String[] parts = line.split("\\s+");
            double[] values = Arrays.stream(parts).mapToDouble(Double::parseDouble).toArray();
            dataList.add(values);
        }
        lSc.close();
        return dataList.toArray(new double[0][]);
    }

    public static double[] interpolate(double[] x, double[] y, double[] newX) {
        double[] newY = new double[newX.length];

        // Handle cases where newX is outside the bounds of x
        for (int i = 0; i < newX.length; i++) {
            newY[i] = interpolate(x, y, newX[i]);
        }

        return newY;
    }

    public static double interpolate(double x[], double y[], double pNewX) {
        double xi = pNewX;
        double newY;

        
        //System.out.println("interp: Compare value "+pNewX+" to range ["+x[0]+"-"+x[x.length - 1]+"]");
        
        // If xi is less than the smallest value in x, return the first y value
        if (xi <= x[0]) {
            newY = y[0];
        } // If xi is greater than the largest value in x, return the last y value
        else if (xi >= x[x.length - 1]) {
            newY = y[y.length - 1];
        } // Otherwise, perform linear interpolation
        else {
            int index = 0;
            while (index < x.length - 1 && x[index + 1] < xi) {
                index++;
            }
            double x1 = x[index], x2 = x[index + 1];
            double y1 = y[index], y2 = y[index + 1];

            // Perform linear interpolation
            newY = y1 + (y2 - y1) * (xi - x1) / (x2 - x1);
        }
        return newY;
    }

    static void flipArray(double[] array) {
        for (int i = 0; i < array.length / 2; i++) {
            double temp = array[i];
            array[i] = array[array.length - 1 - i];
            array[array.length - 1 - i] = temp;
        }
    }

    public static RegressionResult linearRegression(double[] x, double[] y, double confidenceLevel) {
        int n = x.length;
        double sumX = 0, sumY = 0, sumXY = 0, sumX2 = 0;

        for (int i = 0; i < n; i++) {
            sumX += x[i];
            sumY += y[i];
            sumXY += x[i] * y[i];
            sumX2 += x[i] * x[i];
        }

        double xMean = sumX / n;
        double yMean = sumY / n;

        // Compute slope (b1) and intercept (b0)
        double numerator = (n * sumXY - sumX * sumY);
        double denominator = (n * sumX2 - sumX * sumX);
        if (denominator == 0) {
            throw new ArithmeticException("Denominator in slope calculation is zero, indicating perfect correlation.");
        }
        // double b1 = (n * sumXY - sumX * sumY) / (n * sumX2 - sumX * sumX);
        double b1 = numerator / denominator;
        double b0 = yMean - b1 * xMean;

        // Compute residual standard error (se)
        double sse = 0;
        for (int i = 0; i < n; i++) {
            double yPred = b0 + b1 * x[i];
            sse += Math.pow(y[i] - yPred, 2);
        }
        double se = Math.sqrt(sse / (n - 2)); // standard error of regression

        // Compute standard errors
        double seB1 = se / Math.sqrt(sumX2 - (sumX * sumX) / n);
        double seB0 = se * Math.sqrt((1.0 / n) + (xMean * xMean) / (sumX2 - (sumX * sumX) / n));

        // Get t-critical value for confidence interval 
        //use lookup table e.g. https://www.itl.nist.gov/div898/handbook/eda/section3/eda3672.htm
        double tValue = tCriticalValue.getTCriticalValue(n - 2, confidenceLevel, true);
//System.out.printf("--------- tValue (df = %d) = %f\n",n-2, tValue);
//System.out.printf("--------- seB1  = %f\n", seB1);
//System.out.printf("--------- tValue*seB1  = %f\n", tValue*seB1);
//System.out.printf("--------- slope b1  = %f\n", b1);
        // Compute confidence intervals
        double ciB1Lower = b1 - tValue * seB1; //intercept
        double ciB1Upper = b1 + tValue * seB1; //intercept
        double ciB0Lower = b0 - tValue * seB0; //slope
        double ciB0Upper = b0 + tValue * seB0; //slope
//System.out.printf("--------- slope b1  = %f\n", b1);
        return new RegressionResult(b1, b0, se, seB1, seB0, ciB1Lower, ciB1Upper, ciB0Lower, ciB0Upper);
    }

}

// class for regression results
class RegressionResult {

    double slope, intercept, se, seSlope, seIntercept;
    double ciSlopeLower, ciSlopeUpper;
    double ciInterceptLower, ciInterceptUpper;

    public RegressionResult(double slope, double intercept, double se, double seSlope, double seIntercept,
            double ciSlopeLower, double ciSlopeUpper,
            double ciInterceptLower, double ciInterceptUpper) {
        this.slope = slope;
        this.intercept = intercept;
        this.se = se; //residual standard error of fit
        this.seSlope = seSlope;
        this.seIntercept = seIntercept;
        this.ciSlopeLower = ciSlopeLower;
        this.ciSlopeUpper = ciSlopeUpper;
        this.ciInterceptLower = ciInterceptLower;
        this.ciInterceptUpper = ciInterceptUpper;
    }

}
