using MathNet.Numerics;
using MathNet.Numerics.IntegralTransforms;
using System;
using System.Linq;
using System.Numerics;

namespace 测试Scottplot
{
    internal class ClarkTaylorModel
    {
        public static void Construct()
        {
            // 1. 生成模拟数据
            double[] time = Enumerable.Range(0, 1000).Select(i => i * 0.001).ToArray();
            double[] laserPulse = GenerateGaussianPulse(time, 0.05, 0.02);
            double[] idealResponse = GenerateIdealResponse(time, 0.2); // 假设真实响应
            double[] measuredTemp = Convolve(idealResponse, laserPulse); // 模拟实测数据

            // 2. 反卷积
            double[] deconvolved = DeconvolveFFT(measuredTemp, laserPulse);

            // 3. 计算热扩散率
            double tHalf = FindHalfMaxTime(deconvolved, time);
            double alpha = CalculateThermalDiffusivity(1e-3, tHalf);
            Console.WriteLine($"Thermal Diffusivity: {alpha} m²/s");
        }
        //生成高斯脉冲
        static double[] GenerateGaussianPulse(double[] time, double peakTime, double width)
        {
            return time.Select(t => Math.Exp(-Math.Pow((t - peakTime) / width, 2))).ToArray();
        }

        // 生成理想响应（理论模型）
        static double[] GenerateIdealResponse(double[] time, double alpha)
        {
            double d = 1e-3; // 厚度1mm
            double t0 = 0.1388 * d * d / alpha; // 理论半峰时间
            return time.Select(t => t >= t0 ? Math.Exp(-0.5 * Math.Pow((t - t0) / (0.1 * t0), 2)) : 0).ToArray();
        }

        // 卷积运算（模拟实测信号）
        static double[] Convolve(double[] signal, double[] kernel)
        {
            int n = signal.Length;
            double[] result = new double[n];
            for (int i = 0; i < n; i++)
            {
                for (int j = 0; j <= i; j++)
                {
                    result[i] += signal[j] * (i - j < kernel.Length ? kernel[i - j] : 0);
                }
            }
            return result;
        }

        // FFT反卷积（含维纳滤波）
        static double[] DeconvolveFFT(double[] measured, double[] pulse)
        {
            int n = measured.Length;
            Complex[] fftMeasured = new Complex[n];
            Complex[] fftPulse = new Complex[n];

            for (int i = 0; i < n; i++)
            {
                fftMeasured[i] = new Complex(measured[i], 0);
                fftPulse[i] = new Complex(pulse[i], 0);
            }

            Fourier.Forward(fftMeasured);
            Fourier.Forward(fftPulse);

            double noiseLevel = 1e-6;
            Complex[] fftIdeal = new Complex[n];
            for (int i = 0; i < n; i++)
            {
                double denom = fftPulse[i].MagnitudeSquared() + noiseLevel;
                fftIdeal[i] = fftMeasured[i] * Complex.Conjugate(fftPulse[i]) / denom;
            }

            Fourier.Inverse(fftIdeal);
            return fftIdeal.Select(c => c.Real).ToArray();
        }

        // 找半峰时间
        static double FindHalfMaxTime(double[] response, double[] time)
        {
            double maxTemp = response.Max();
            int halfMaxIndex = Array.FindIndex(response, t => t >= maxTemp / 2);
            return time[halfMaxIndex];
        }

        // 计算热扩散率
        static double CalculateThermalDiffusivity(double thickness, double tHalf)
        {
            return 0.1388 * Math.Pow(thickness, 2) / tHalf;
        }
    }
    
  
    
}