#include "WindInverserToolGPU.cuh"
#include "GPUTool.cuh"


// 使用 __device__ 限定符，使函数在GPU上执行且只能被其他GPU函数调用
__device__ float windspeedcmd5_device(float angle, float incidenceangle, float sigma0) {
    // CMOD5模型系数数组 - 常量内存通常更适合此类数据
    const float cmod5_c[30] = {
        1.0f, -0.688f, -0.793f, 0.338f, -0.173f, 0.0f, 0.004f,
        0.111f, 0.0162f, 6.34f, 2.57f, -2.18f, 0.4f,
        -0.6f, 0.045f, 0.007f, 0.33f, 0.012f, 22.0f,
        1.95f, 3.0f, 8.39f, -3.44f, 1.36f, 5.35f, 1.99f, 0.29f, 3.80f, 1.53f
    };

    // 确保入射角不小于16度
    if (incidenceangle + 0.5f <= 16.0f) {
        incidenceangle = 16.0f;
    }

    // 在0-26m/s范围内搜索风速
    for (float U = 0.0f; U <= 26.0f; U += 0.01f) {
        // 计算归一化参数
        float y0 = cmod5_c[19];  // 索引调整：C数组从0开始
        float n = cmod5_c[20];
        float a = y0 - (y0 - 1.0f) / n;
        float b = 1.0f / (n * powf(y0 - 1.0f, n - 1.0f));  // 使用powf而非^

        // 归一化入射角
        float x = (incidenceangle - 40.0f) / 25.0f;

        // 计算基础后向散射系数
        float a0 = cmod5_c[1] + cmod5_c[2] * x + cmod5_c[3] * x * x + cmod5_c[4] * x * x * x;
        float a1 = cmod5_c[5] + cmod5_c[6] * x;
        float a2 = cmod5_c[7] + cmod5_c[8] * x;
        float gamma = cmod5_c[9] + cmod5_c[10] * x + cmod5_c[11] * x * x;
        float s0 = cmod5_c[12] + cmod5_c[13] * x;
        float s = a2 * U;

        // 计算饱和函数
        float f;
        if (s < s0) {
            float g = 1.0f / (1.0f + expf(-s0));  // 使用expf
            float alpha = s0 * (1.0f - g);
            f = powf(s / s0, alpha) * g;  // 使用powf
        }
        else {
            f = 1.0f / (1.0f + expf(-s));
        }

        // 计算各向同性部分
        float b0 = powf(f, gamma) * powf(10.0f, a0 + a1 * U);  // 使用powf

        // 计算一阶调和系数
        float b1_term = cmod5_c[15] * U * (0.5f + x - tanhf(4.0f * (x + cmod5_c[16] + cmod5_c[17] * U)));
        float b1 = (cmod5_c[14] * (1.0f + x) - b1_term) / (expf(0.34f * (U - cmod5_c[18])) + 1.0f);

        // 计算二阶调和系数
        float v0 = cmod5_c[21] + cmod5_c[22] * x + cmod5_c[23] * x * x;
        float d1 = cmod5_c[24] + cmod5_c[25] * x + cmod5_c[26] * x * x;
        float d2 = cmod5_c[27] + cmod5_c[28] * x;
        float y = (U / v0 + 1.0f);

        float v2;
        if (y < y0) {
            v2 = a + b * powf(y - 1.0f, n);  // 使用powf
        }
        else {
            v2 = y;
        }
        float b2 = (-d1 + d2 * v2) * expf(-v2);

        // 计算模拟后向散射系数并与实测值比较
        float angle_rad = angle * 3.14159265f / 180.0f;
        float sigma0_model = 10.0f * log10f(b0 * powf(1.0f + b1 * cosf(angle_rad) + b2 * cosf(2.0f * angle_rad), 1.6f));

        if (sigma0_model - sigma0 >= 0.0f) {
            return U;
        }
    }

    return 26.0f; // 如果未找到，返回最大风速
}








__global__ void CUDA_WindSpeedInversered_Kernel(float* in_wind_dir_U_arr_cuda,
												float* in_wind_dir_V_arr_cuda, 
												float* in_inc_angle_arr_cuda, 
												float* in_sar_sigma0_arr_cuda, 
												float* out_wind_speed_arr_cuda, int64_t pixelCount) {
	int64_t  idx = blockIdx.x * blockDim.x + threadIdx.x;
	if (idx < pixelCount) {
		float sar_sigma0 = in_sar_sigma0_arr_cuda[idx];
        if (isinf(sar_sigma0) || isnan(sar_sigma0)) {
            out_wind_speed_arr_cuda[idx] = -1;
        }
        else {
            float U10 = in_wind_dir_U_arr_cuda[idx];
            float V10 = in_wind_dir_V_arr_cuda[idx];
            float inc_angle = in_inc_angle_arr_cuda[idx];

            float wind_dir_angle = 270 - atan2(V10, U10) * 180.0 / PI;
            wind_dir_angle = fmodf(wind_dir_angle, 360.0f);
            float out_wind_speed = windspeedcmd5_device(wind_dir_angle, inc_angle, sar_sigma0);

            //printf("windspeedcmd5(%f, %f, %f)-%f \n", wind_dir_angle, inc_angle, sar_sigma0, out_wind_speed);

            out_wind_speed_arr_cuda[idx] = out_wind_speed;
        }
	}
}





void WindSpeedInverseredFunc(float* in_wind_dir_U_arr_cuda, float* in_wind_dir_V_arr_cuda, float* in_inc_angle_arr_cuda, float* in_sar_sigma0_arr_cuda, float* out_wind_speed_arr_cuda, int64_t pixelCount)
{
	// 设置 CUDA 核函数的网格和块的尺寸
	int64_t blockSize = 256; // 每个块的线程数
	int64_t numBlocks = (pixelCount + blockSize - 1) / blockSize; // 根据 pixelcount 计算网格大小
	// 调用 CUDA 核函数
	CUDA_WindSpeedInversered_Kernel << <numBlocks, blockSize >> > (in_wind_dir_U_arr_cuda, in_wind_dir_V_arr_cuda, in_inc_angle_arr_cuda, in_sar_sigma0_arr_cuda, out_wind_speed_arr_cuda, pixelCount);

#ifdef __CUDADEBUG__
	cudaError_t err = cudaGetLastError();
	if (err != cudaSuccess) {
		printf("CUDAdistanceAB CUDA Error: %s\n", cudaGetErrorString(err));
		exit(2);
	}
#endif // __CUDADEBUG__
	cudaDeviceSynchronize();
}