[未测试]补充了仿真函数对应的单精度版本

BaseCommonLibrary/BaseTool/BaseConstVariable.h

@@ -148,6 +148,11 @@ struct Vector3 {
 };
 
 
+struct Vector3_single {
+	float x, y, z;
+};
+
+
 /*********************************************** FEKO仿真参数 ********************************************************************/
 struct SatellitePos {
 	double time;

Toolbox/SimulationSARTool/GPUBpSimulation.cu

@@ -21,10 +21,10 @@
 #include "GPUBPTool.cuh"
 #include "BPBasic0_CUDA.cuh"
 #include "GPUBpSimulation.cuh"
+#include "GPURFPC_single.cuh"
 #include "GPURFPC.cuh"
  
 
-
 double* getFreqPoints_mallocHost(double startFreq, double endFreq, long freqpoints)
 {
 	long double dfreq = (endFreq - startFreq) / (freqpoints - 1);
@@ -68,7 +68,7 @@ __global__ void kernel_RFPCProcess(
 	// ÈëÉä½Ç
 	Vector3 TS = vec_sub(S, T);//T-->S
 	double localIncAngle = angleBetweenVectors(TS, slp, false);// ÈëÉä½Ç
-	double sigma0 = GPU_getSigma0dB(p1, p2, p3, p4, p5, p6, localIncAngle);// 后向散射系数
+	double sigma0 = GPU_getSigma0dB_params(p1, p2, p3, p4, p5, p6, localIncAngle);// 后向散射系数
 	sigma0 = powf(10.0, sigma0 / 10.0);
 
 	// ¾àÀë
@@ -95,8 +95,6 @@ void RFPCProcess(double Tx, double Ty, double Tz,
 	double p1, double p2, double p3, double p4, double p5, double p6, 
 	GPUDATA& d_data)
 {
-
-
 	double* AntX = (double*)mallocCUDADevice(sizeof(double) * d_data.Np);
 	double* AntY = (double*)mallocCUDADevice(sizeof(double) * d_data.Np);
 	double* AntZ = (double*)mallocCUDADevice(sizeof(double) * d_data.Np);
@@ -128,3 +126,125 @@ void RFPCProcess(double Tx, double Ty, double Tz,
 
 
 
+/** 单精度成像测试 *****************************************************************************************************************************************/
+
+
+
+
+
+
+
+
+
+
+__global__ void kernel_RFPCProcess_single(
+	float* Sx, float* Sy, float* Sz,
+	float Tx, float Ty, float Tz,
+	float Tslx, float Tsly, float Tslz, // 目标的坡面向量
+	float p1, float p2, float p3, float p4, float p5, float p6,
+	long Np, long Nf,
+	float minF, float dFreq, float RefRange,
+	cuComplex* phdata
+) {
+	long prfid = blockIdx.x * blockDim.x + threadIdx.x; // 获取当前的线程编码
+	if (prfid >= Np || prfid < 0) { return; }
+	else {}
+
+	// 距离
+	Vector3_single S{ Sx[prfid],Sy[prfid],Sz[prfid] };
+	Vector3_single T{ Tx,Ty,Tz };
+	Vector3_single slp{ Tslx,Tsly,Tslz };
+	//printf("S=(%e,%e,%e)\n", S.x, S.y, S.z);
+	// 入射角
+	Vector3_single TS = vec_sub_single(S, T);//T-->S
+	float localIncAngle = angleBetweenVectors_single(TS, slp, false);// 入射角
+	float sigma0 = GPU_getSigma0dB_single(p1, p2, p3, p4, p5, p6, localIncAngle);// 后向散射系数
+	sigma0 = powf(10.0, sigma0 / 10.0);
+
+	// 距离
+	float R = sqrtf(vec_dot_single(TS, TS));
+
+	// 计算增益
+	float amp_echo = sigma0 / (powf(4 * LAMP_CUDA_PI, 2) * powf(R, 4)); // 反射强度
+	float phi = 0;
+	for (long fid = 0; fid < Nf; fid++) {
+		phi = -4.0 * PI / LIGHTSPEED * (minF + dFreq * fid) * (R - RefRange);
+		phdata[prfid * Nf + fid].x += amp_echo * cosf(phi);
+		phdata[prfid * Nf + fid].y += amp_echo * sinf(phi);
+		//printf("phdata(%d,%d)=complex(%e,%e)\n", prfid, fid, phdata[prfid * Nf + fid].x, phdata[prfid * Nf + fid].y);
+	}
+	//printf("Nf:%d\n", Nf);
+	//printf("amp_echo:%f\n", amp_echo);
+	//printf("sigma0:%f\n", sigma0);
+	//printf("R:%e\n", R);
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+float* getFreqPoints_mallocHost_single(float startFreq, float endFreq, long freqpoints)
+{
+	float dfreq = (endFreq - startFreq) / (freqpoints - 1);
+	float* freqlist = (float*)mallocCUDAHost(sizeof(float) * freqpoints);
+	for (long i = 0; i < freqpoints; i++) {
+		freqlist[i] = startFreq + dfreq * i;
+	}
+	return freqlist;
+}
+
+void RFPCProcess_single(
+	float Tx, float Ty, float Tz, 
+	float Tslx, float Tsly, float Tslz, 
+	float p1, float p2, float p3, float p4, float p5, float p6, 
+	GPUDATA_single& d_data)
+{
+	float* AntX = (float*)mallocCUDADevice(sizeof(float) * d_data.Np);
+	float* AntY = (float*)mallocCUDADevice(sizeof(float) * d_data.Np);
+	float* AntZ = (float*)mallocCUDADevice(sizeof(float) * d_data.Np);
+
+	HostToDevice(d_data.AntX, AntX, sizeof(float) * d_data.Np); printf("antX host to device finished!!\n");
+	HostToDevice(d_data.AntY, AntY, sizeof(float) * d_data.Np); printf("antY host to device finished!!\n");
+	HostToDevice(d_data.AntZ, AntZ, sizeof(float) * d_data.Np); printf("antZ host to device finished!!\n");
+	float minF = d_data.minF[0];
+	long grid_size = (d_data.Np + BLOCK_SIZE - 1) / BLOCK_SIZE;
+	float dfreq = d_data.deltaF;
+	float R0 = d_data.R0;
+	kernel_RFPCProcess_single << <grid_size, BLOCK_SIZE >> > (
+		AntX, AntY, AntZ,
+		Tx, Ty, Tz,
+		Tslx, Tsly, Tslz, // 目标的坡面向量
+		p1, p2, p3, p4, p5, p6,
+		d_data.Np, d_data.Nfft,
+		minF, dfreq, R0,
+		d_data.phdata
+		);
+
+	PrintLasterError("RFPCProcess");
+	FreeCUDADevice(AntX);
+	FreeCUDADevice(AntY);
+	FreeCUDADevice(AntZ);
+
+	return;
+
+}
+
+
+

Toolbox/SimulationSARTool/GPUBpSimulation.cuh

@@ -30,6 +30,17 @@ extern "C" void RFPCProcess(
 	GPUDATA& d_data
 	);
 
+extern "C" float* getFreqPoints_mallocHost_single(float startFreq, float endFreq, long freqpoints);
+ 
+extern "C" void RFPCProcess_single(
+	float Tx, float Ty, float Tz, // 目标点坐标
+	float Tslx, float Tsly, float Tslz, // 目标的坡面向量
+	float p1, float p2, float p3, float p4, float p5, float p6,// 地面目标后向散射系数与入射角关系 系数
+	GPUDATA_single& d_data
+);
+
+
+
 
 
 #endif

Toolbox/SimulationSARTool/PowerSimulationIncoherent/LookTableSimulationComputer.cu

@@ -13,19 +13,19 @@
 #include "LookTableSimulationComputer.cuh"
 
 
-extern  __device__ __host__ double getNumberDopplerCenterRate(double R, double r0, double r1, double r2, double r3, double r4, double reftime)
+extern  __device__  double getNumberDopplerCenterRate(double R, double r0, double r1, double r2, double r3, double r4, double reftime)
 {
 	double t=R / LIGHTSPEED - reftime;
 	double dopplerCenterRate = r0 + r1*t + r2*t*t + r3*t*t*t + r4*t*t*t*t;
 	return  dopplerCenterRate;
 }
 
-__device__ __host__ double getDopplerCenterRate(double Rx, double Ry, double Rz, double Vx, double Vy, double Vz, double fact_lamda)
+__device__  double getDopplerCenterRate(double Rx, double Ry, double Rz, double Vx, double Vy, double Vz, double fact_lamda)
 {
 	return -2*(Rx*Vx+Ry*Vy+Rz*Vz)* fact_lamda;
 }
 
-__device__ __host__ double getPolyfitNumber(double x, double a0, double a1, double a2, double a3, double a4, double a5)
+__device__  double getPolyfitNumber(double x, double a0, double a1, double a2, double a3, double a4, double a5)
 {
 	return a0 + a1 * x + a2 * x * x + a3 * x * x * x + a4 * x * x * x * x + a5 * x * x * x * x * x;
 }

Toolbox/SimulationSARTool/PowerSimulationIncoherent/LookTableSimulationComputer.cuh

@@ -8,11 +8,11 @@
 
 
 
-extern __device__ __host__ double getNumberDopplerCenterRate(double R, double r0, double r1, double r2, double r3, double r4,double reftime);
+extern __device__  double getNumberDopplerCenterRate(double R, double r0, double r1, double r2, double r3, double r4,double reftime);
 
-extern __device__ __host__ double getDopplerCenterRate(double Rx,double Ry,double Rz,double Vx,double Vy,double Vz,double fact_lamda); //fact_lamda lamda ľÄľšĘý
+extern __device__  double getDopplerCenterRate(double Rx,double Ry,double Rz,double Vx,double Vy,double Vz,double fact_lamda); //fact_lamda lamda ľÄľšĘý
 
-extern __device__ __host__ double getPolyfitNumber(double x, double a0, double a1, double a2, double a3, double a4,double a5);
+extern __device__  double getPolyfitNumber(double x, double a0, double a1, double a2, double a3, double a4,double a5);
 
 
 

Toolbox/SimulationSARTool/SimulationSAR/BPBasic0_CUDA.cu

@@ -276,6 +276,268 @@ void BPBasic0(GPUDATA& h_data)
     freeGPUData(d_data);
 }
 
+
+
+
+
+
+/** 单精度代码**********************************************************************************************/
+
+
+
+
+
+
+
+__global__ void phaseCompensationKernel_single(cufftComplex* phdata, const float* Freq, float r, int K, int Na) {
+    int freqIdx = blockIdx.x * blockDim.x + threadIdx.x;
+    int pulseIdx = blockIdx.y * blockDim.y + threadIdx.y;
+
+    if (freqIdx >= K || pulseIdx >= Na) return;
+
+    int idx = pulseIdx * K + freqIdx;
+    float phase = 4 * PI * Freq[freqIdx] * r / c;
+    float cos_phase = cosf(phase);
+    float sin_phase = sinf(phase);
+
+    cufftComplex ph = phdata[idx];
+    float new_real = ph.x * cos_phase - ph.y * sin_phase;
+    float new_imag = ph.x * sin_phase + ph.y * cos_phase;
+    phdata[idx] = make_cuComplex(new_real, new_imag);
+}
+
+__global__ void fftshiftKernel_single(cufftComplex* data, int Nfft, int Np) {
+    int pulse = blockIdx.x * blockDim.x + threadIdx.x;
+    if (pulse >= Np) return;
+
+    int half = Nfft / 2;
+    for (int i = 0; i < half; ++i) {
+        cufftComplex temp = data[pulse * Nfft + i];
+        data[pulse * Nfft + i] = data[pulse * Nfft + i + half];
+        data[pulse * Nfft + i + half] = temp;
+    }
+}
+
+__global__ void processPulseKernel_single(
+    long prfid,
+    int nx, int ny,
+    const float* x_mat, const float* y_mat, const float* z_mat,
+    float AntX, float AntY, float AntZ,
+    float R0, float minF,
+    const cufftComplex* rc_pulse,
+    const float r_start, const float dr, const int nR,
+    cufftComplex* im_final
+) {
+    //
+
+    long long idx = blockIdx.x * blockDim.x + threadIdx.x;
+    long long pixelcount = nx * ny;
+    if (idx >= pixelcount) return;
+
+    //printf("processPulseKernel start!!\n");
+
+    //if (x >= nx || y >= ny) return;
+    //int idx = x * ny + y;
+
+
+    float dx = AntX - x_mat[idx];
+    float dy = AntY - y_mat[idx];
+    float dz = AntZ - z_mat[idx];
+
+    //printf("processPulseKernel xmat !!\n");
+    float R = sqrtf(dx * dx + dy * dy + dz * dz);
+    float dR = R - R0;
+
+    if (dR < r_start || dR >= (r_start + dr * (nR - 1))) return;
+    // Linear interpolation
+    float pos = (dR - r_start) / dr;
+    int index = (int)floorf(pos);
+    float weight = pos - index;
+
+    if (index < 0 || index >= nR - 1) return;
+
+    cufftComplex rc_low = rc_pulse[prfid * nR + index];
+    cufftComplex rc_high = rc_pulse[prfid * nR + index + 1];
+    cufftComplex rc_interp;
+    rc_interp.x = rc_low.x * (1 - weight) + rc_high.x * weight;
+    rc_interp.y = rc_low.y * (1 - weight) + rc_high.y * weight;
+
+    // Phase correction
+    float phase = 4 * PI * minF / c * dR;
+    float cos_phase = cosf(phase);
+    float sin_phase = sinf(phase);
+
+    cufftComplex phCorr;
+    phCorr.x = rc_interp.x * cos_phase - rc_interp.y * sin_phase;
+    phCorr.y = rc_interp.x * sin_phase + rc_interp.y * cos_phase;
+
+    // Accumulate
+    im_final[idx].x += phCorr.x;
+    im_final[idx].y += phCorr.y;
+    //printf("r_start=%e;dr=%e;nR=%d\n", r_start, dr, nR);
+ 
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+void bpBasic0CUDA_single(GPUDATA_single& data, int flag, float* h_R)
+{
+    // Phase compensation
+    if (flag == 1) {
+        dim3 block(16, 16);
+        dim3 grid((data.K + 15) / 16, (data.Np + 15) / 16);
+        phaseCompensationKernel_single << <grid, block >> > (data.phdata, data.Freq, data.R0, data.K, data.Np);
+        PrintLasterError("bpBasic0CUDA Phase compensation");
+        //data.R0 = data.r;  // 假设data.r已正确设置
+    }
+
+    // FFT处理
+    cufftHandle plan;
+    cufftPlan1d(&plan, data.Nfft, CUFFT_C2C, data.Np);
+    cufftExecC2C(plan, data.phdata, data.phdata, CUFFT_INVERSE);
+    cufftDestroy(plan);
+
+    // FFT移位
+    dim3 blockShift(256);
+    dim3 gridShift((data.Np + 255) / 256);
+    fftshiftKernel_single << <gridShift, blockShift >> > (data.phdata, data.Nfft, data.Np);
+    PrintLasterError("bpBasic0CUDA Phase FFT Process");
+
+    printfinfo("fft finished!!\n");
+    // 图像重建
+
+
+
+    float r_start = data.r_vec[0];
+    float dr = (data.r_vec[data.Nfft - 1] - r_start) / (data.Nfft - 1);
+    printfinfo("dr = %f\n", dr);
+    long pixelcount = data.nx * data.ny;
+    long grid_size = (pixelcount + BLOCK_SIZE - 1) / BLOCK_SIZE;
+    printfinfo("grid finished!!\n");
+
+    //float* d_R = (float*)mallocCUDADevice(sizeof(float) * data.nx * data.ny);
+    printfinfo("r_start=%e;dr=%e;nR=%d\n", r_start, dr, data.Nfft);
+    printfinfo("BPimage .....\n");
+    for (long ii = 0; ii < data.Np; ++ii) {
+        processPulseKernel_single << <grid_size, BLOCK_SIZE >> > (
+            ii,
+            data.nx, data.ny,
+            data.x_mat, data.y_mat, data.z_mat,
+            data.AntX[ii], data.AntY[ii], data.AntZ[ii],
+            data.R0, data.minF[ii],
+            data.phdata,
+            r_start, dr, data.Nfft,
+            data.im_final
+            //,d_R
+            );
+        PrintLasterError("processPulseKernel");
+        if (ii % 1000 == 0) {
+            printfinfo("\rPRF(%f %)  %d / %d\t\t\t\t", (ii * 100.0 / data.Np), ii, data.Np);
+        }
+    }
+    //FreeCUDADevice(d_R);
+
+    PrintLasterError("bpBasic0CUDA Phase BPimage Process finished!!");
+}
+
+void initGPUData_single(GPUDATA_single& h_data, GPUDATA_single& d_data)
+{
+    d_data.AntX = h_data.AntX; //(double*)mallocCUDADevice(sizeof(double) * h_data.Np);
+    d_data.AntY = h_data.AntY;//(double*)mallocCUDADevice(sizeof(double) * h_data.Np);
+    d_data.AntZ = h_data.AntZ;// (double*)mallocCUDADevice(sizeof(double) * h_data.Np);
+    d_data.minF = h_data.minF;// (double*)mallocCUDADevice(sizeof(double) * h_data.Np);
+    d_data.x_mat = (float*)mallocCUDADevice(sizeof(float) * h_data.nx * h_data.ny);
+    d_data.y_mat = (float*)mallocCUDADevice(sizeof(float) * h_data.nx * h_data.ny);
+    d_data.z_mat = (float*)mallocCUDADevice(sizeof(float) * h_data.nx * h_data.ny);
+    d_data.r_vec = h_data.r_vec;// (double*)mallocCUDADevice(sizeof(double) * h_data.Nfft);
+    d_data.Freq = (float*)mallocCUDADevice(sizeof(float) * h_data.Nfft);
+    d_data.phdata = (cuComplex*)mallocCUDADevice(sizeof(cuComplex) * h_data.Nfft * h_data.Np);
+    d_data.im_final = (cuComplex*)mallocCUDADevice(sizeof(cuComplex) * h_data.nx * h_data.ny);
+
+    //HostToDevice(h_data.AntX, d_data.AntX,sizeof(double) * h_data.Np);
+    //HostToDevice(h_data.AntY, d_data.AntY,sizeof(double) * h_data.Np);
+    //HostToDevice(h_data.AntZ, d_data.AntZ,sizeof(double) * h_data.Np);
+    //HostToDevice(h_data.minF, d_data.minF,sizeof(double) * h_data.Np);
+    HostToDevice(h_data.x_mat, d_data.x_mat, sizeof(float) * h_data.nx * h_data.ny); printf("image X Copy finished!!!\n");
+    HostToDevice(h_data.y_mat, d_data.y_mat, sizeof(float) * h_data.nx * h_data.ny); printf("image Y Copy finished!!!\n");
+    HostToDevice(h_data.z_mat, d_data.z_mat, sizeof(float) * h_data.nx * h_data.ny); printf("image Z Copy finished!!!\n");
+    HostToDevice(h_data.Freq, d_data.Freq, sizeof(float) * h_data.Nfft);
+    //HostToDevice(h_data.r_vec, d_data.r_vec, sizeof(double) * h_data.Nfft);
+    HostToDevice(h_data.phdata, d_data.phdata, sizeof(cuComplex) * h_data.Nfft * h_data.Np); printf("image echo Copy finished!!!\n");
+    HostToDevice(h_data.im_final, d_data.im_final, sizeof(cuComplex) * h_data.nx * h_data.ny); printf("image data  Copy finished!!!\n");
+
+    // 拷贝标量参数
+    d_data.Nfft = h_data.Nfft;
+    d_data.K = h_data.K;
+    d_data.Np = h_data.Np;
+    d_data.nx = h_data.nx;
+    d_data.ny = h_data.ny;
+    d_data.R0 = h_data.R0;
+    d_data.deltaF = h_data.deltaF;
+}
+
+void freeGPUData_single(GPUDATA_single& d_data)
+{
+    //FreeCUDADevice((d_data.AntX));
+//FreeCUDADevice((d_data.AntY));
+//FreeCUDADevice((d_data.AntZ));
+//FreeCUDADevice((d_data.minF));
+    FreeCUDADevice((d_data.x_mat));
+    FreeCUDADevice((d_data.y_mat));
+    FreeCUDADevice((d_data.z_mat));
+    //FreeCUDADevice((d_data.r_vec));
+    FreeCUDADevice((d_data.Freq));
+    FreeCUDADevice((d_data.phdata));
+    FreeCUDADevice((d_data.im_final));
+}
+
+void freeHostData_single(GPUDATA_single& h_data)
+{
+    //FreeCUDAHost((h_data.AntX)); 
+//FreeCUDAHost((h_data.AntY));
+//FreeCUDAHost((h_data.AntZ));
+    FreeCUDAHost((h_data.minF));
+    //FreeCUDAHost((h_data.x_mat));
+    //FreeCUDAHost((h_data.y_mat));
+    //FreeCUDAHost((h_data.z_mat));
+    FreeCUDAHost((h_data.r_vec));
+    FreeCUDAHost((h_data.Freq));
+    FreeCUDAHost((h_data.phdata));
+    FreeCUDAHost((h_data.im_final));
+}
+
+void BPBasic0_single(GPUDATA_single& h_data)
+{
+    GPUDATA_single d_data;
+    initGPUData_single(h_data, d_data);
+    bpBasic0CUDA_single(d_data, 0);
+    DeviceToHost(h_data.im_final, d_data.im_final, sizeof(cuComplex) * h_data.nx * h_data.ny);
+    freeGPUData_single(d_data);
+}
+
 //int main() {
 //    GPUDATA h_data, d_data;
 //

Toolbox/SimulationSARTool/SimulationSAR/BPBasic0_CUDA.cuh

@@ -30,6 +30,21 @@ struct GPUDATA {
     double deltaF; // 频点范围
 };
 
+
+struct GPUDATA_single {
+    long Nfft, K, Np, nx, ny; // 傅里叶点数、频点数、脉冲数、图像列、图像行
+    float* AntX, * AntY, * AntZ, * minF; // 天线坐标、起始频率
+    float* x_mat, * y_mat, * z_mat;// 地面坐标
+    float* r_vec; // 坐标范围
+    float* Freq;// 频率
+    cuComplex* phdata;// 回波
+    cuComplex* im_final;// 图像
+    float R0; // 参考斜距
+    float deltaF; // 频点范围
+};
+
+
+
 extern "C" {
     void bpBasic0CUDA(GPUDATA& data, int flag,double* h_R=nullptr);
     void initGPUData(GPUDATA& h_data, GPUDATA& d_data);
@@ -38,4 +53,14 @@ extern "C" {
     void BPBasic0(GPUDATA& h_data);
 };
 
+
+extern "C" {
+    void bpBasic0CUDA_single(GPUDATA_single& data, int flag, float* h_R = nullptr);
+    void initGPUData_single(GPUDATA_single& h_data, GPUDATA_single& d_data);
+    void freeGPUData_single(GPUDATA_single& d_data);
+    void freeHostData_single(GPUDATA_single& d_data);
+    void BPBasic0_single(GPUDATA_single& h_data);
+};
+
+
 #endif

Toolbox/SimulationSARTool/SimulationSAR/GPUBPTool.cu

@@ -22,13 +22,13 @@
 #define M_PI 3.14159265358979323846
 #endif
 
-__device__ __host__  double angleBetweenVectors(Vector3 a, Vector3 b, bool returnDegrees ) {
+__device__   double angleBetweenVectors(Vector3 a, Vector3 b, bool returnDegrees ) {
 	// 计算点积
 	double dotProduct = a.x * b.x + a.y * b.y + a.z * b.z;
 
 	// 计算模长
-	double magA = std::sqrt(a.x * a.x + a.y * a.y + a.z * a.z);
+	double magA =  sqrt(a.x * a.x + a.y * a.y + a.z * a.z);
-	double magB = std::sqrt(b.x * b.x + b.y * b.y + b.z * b.z);
+	double magB =  sqrt(b.x * b.x + b.y * b.y + b.z * b.z);
 
 	// 处理零向量异常
 	if (magA == 0.0 || magB == 0.0) {
@@ -37,33 +37,33 @@ __device__ __host__  double angleBetweenVectors(Vector3 a, Vector3 b, bool retur
 
 	double cosTheta = dotProduct / (magA * magB);
 	cosTheta = cosTheta < -1 ? -1 : cosTheta>1 ? 1 : cosTheta;  // 截断到[-1, 1]
-	double angleRad = std::acos(cosTheta);
+	double angleRad = acos(cosTheta);
 	return returnDegrees ? angleRad * 180.0 / M_PI : angleRad;
 }
 
 // 向量运算
-__device__ __host__ Vector3 vec_sub(Vector3 a, Vector3 b) {
+__device__  Vector3 vec_sub(Vector3 a, Vector3 b) {
 	return { a.x - b.x, a.y - b.y, a.z - b.z };
 }
 
-__device__ __host__ double vec_dot(Vector3 a, Vector3 b) {
+__device__  double vec_dot(Vector3 a, Vector3 b) {
 	return a.x * b.x + a.y * b.y + a.z * b.z;
 }
 
-__device__ __host__ Vector3 vec_cross(Vector3 a, Vector3 b) {
+__device__  Vector3 vec_cross(Vector3 a, Vector3 b) {
 	return { a.y * b.z - a.z * b.y,
 			a.z * b.x - a.x * b.z,
 			a.x * b.y - a.y * b.x };
 }
 
-__device__ __host__ Vector3 vec_normalize(Vector3 v) {
+__device__  Vector3 vec_normalize(Vector3 v) {
 	double len = sqrt(vec_dot(v, v));
 	return (len > 1e-12) ? Vector3{ v.x / len, v.y / len, v.z / len } : v;
 }
 
 
 // 标准正交基底坐标计算
-__device__ __host__ Vector3 coordinates_orthonormal_basis(Vector3& A,
+__device__  Vector3 coordinates_orthonormal_basis(Vector3& A,
 	Vector3& e1,
 	Vector3& e2,
 	Vector3& e3) {
@@ -87,7 +87,7 @@ __device__ __host__ Vector3 coordinates_orthonormal_basis(Vector3& A,
 
 
 // 计算视线交点T
-extern __device__ __host__ Vector3 compute_T(Vector3 S, Vector3 ray, double H) {
+extern __device__  Vector3 compute_T(Vector3 S, Vector3 ray, double H) {
 	Vector3 dir = vec_normalize(ray);
 	double a_h = WGS84_A + H;
 
@@ -105,7 +105,7 @@ extern __device__ __host__ Vector3 compute_T(Vector3 S, Vector3 ray, double H) {
 }
 
 // 主计算函数A
-extern __device__ __host__ Vector3 compute_P(Vector3 S, Vector3 T, double R, double H) {
+extern __device__  Vector3 compute_P(Vector3 S, Vector3 T, double R, double H) {
 	Vector3 ex, ey, ez; // 平面基函数
 	Vector3 ST = vec_normalize(vec_sub(T, S));// S->T 
 	Vector3 SO = vec_normalize(vec_sub(Vector3{ 0, 0, 0 }, S)); // S->O 
@@ -203,6 +203,43 @@ extern __device__ __host__ Vector3 compute_P(Vector3 S, Vector3 T, double R, dou
 	return P;
 }
 
+__device__  double angleBetweenVectors_single(Vector3_single a, Vector3_single b, bool returnDegrees)
+{
+	// 计算点积
+	float dotProduct = a.x * b.x + a.y * b.y + a.z * b.z;
+
+	// 计算模长
+	float magA =  sqrtf(a.x * a.x + a.y * a.y + a.z * a.z);
+	float magB =  sqrtf(b.x * b.x + b.y * b.y + b.z * b.z);
+
+	// 处理零向量异常
+	if (magA == 0.0 || magB == 0.0) {
+		return NAN;
+	}
+
+	float cosTheta = dotProduct / (magA * magB);
+	cosTheta = cosTheta < -1 ? -1 : cosTheta>1 ? 1 : cosTheta;  // 截断到[-1, 1]
+	float angleRad =   acosf(cosTheta);
+	return returnDegrees ? angleRad * 180.0 / M_PI : angleRad;
+}
+
+__device__  Vector3_single vec_sub_single(Vector3_single a, Vector3_single b)
+{
+	return { a.x - b.x, a.y - b.y, a.z - b.z };
+}
+
+__device__  float vec_dot_single(Vector3_single a, Vector3_single b)
+{
+	return a.x * b.x + a.y * b.y + a.z * b.z;
+}
+
+__device__  Vector3_single vec_cross_single(Vector3_single a, Vector3_single b)
+{
+	return Vector3_single{ a.y * b.z - a.z * b.y,
+			a.z * b.x - a.x * b.z,
+			a.x * b.y - a.y * b.x };
+}
+
  
 
 //// 参数校验与主函数

Toolbox/SimulationSARTool/SimulationSAR/GPUBPTool.cuh

@@ -8,12 +8,16 @@
 
 
 
-extern __device__ __host__ double angleBetweenVectors(Vector3 a, Vector3 b, bool returnDegrees = false);
+extern __device__  double angleBetweenVectors(Vector3 a, Vector3 b, bool returnDegrees = false);
-extern __device__ __host__ Vector3 vec_sub(Vector3 a, Vector3 b);
+extern __device__  Vector3 vec_sub(Vector3 a, Vector3 b);
-extern __device__ __host__ double vec_dot(Vector3 a, Vector3 b);
+extern __device__  double vec_dot(Vector3 a, Vector3 b);
-extern __device__ __host__ Vector3 vec_cross(Vector3 a, Vector3 b);
+extern __device__  Vector3 vec_cross(Vector3 a, Vector3 b);
-extern __device__ __host__ Vector3 vec_normalize(Vector3 v);
+extern __device__  Vector3 vec_normalize(Vector3 v);
-extern __device__ __host__ Vector3 compute_T(Vector3 S, Vector3 ray_dir, double H);
+extern __device__  Vector3 compute_T(Vector3 S, Vector3 ray_dir, double H);
-// 
+extern __device__  Vector3 compute_P(Vector3 S, Vector3 T, double R, double H );
-extern __device__ __host__ Vector3 compute_P(Vector3 S, Vector3 T, double R, double H );
 
+
+extern __device__  double angleBetweenVectors_single(Vector3_single a, Vector3_single b, bool returnDegrees = false);
+extern __device__  Vector3_single vec_sub_single(Vector3_single a, Vector3_single b);
+extern __device__  float vec_dot_single(Vector3_single a, Vector3_single b);
+extern __device__  Vector3_single vec_cross_single(Vector3_single a, Vector3_single b);

Toolbox/SimulationSARTool/SimulationSAR/GPURFPC.cu

@@ -1,9 +1,4 @@
-
+#include <cuda.h>
-#include <time.h>
-#include <iostream>
-#include <memory>
-#include <cmath>
-#include <complex>
 #include <device_launch_parameters.h>
 #include <cuda_runtime.h>
 #include <cublas_v2.h>
@@ -19,12 +14,13 @@
 /* »úÆ÷º¯Êý  ****************************************************************************************************************************/
 
 
-extern __host__ __device__ double GPU_getSigma0dB(CUDASigmaParam param, double theta) {//ÏßÐÔÖµ
+extern  __device__ double GPU_getSigma0dB(CUDASigmaParam param, double theta) {//ÏßÐÔÖµ
 	double sigma = param.p1 + param.p2 * exp(-param.p3 * theta) + param.p4 * cos(param.p5 * theta + param.p6);
 	return  sigma;
 }
 
-extern __host__ __device__ double GPU_getSigma0dB(
+
+__device__ double GPU_getSigma0dB_params(
 	const double p1, const double p2, const  double p3, const  double p4, const  double p5, const  double p6, 
 	 double theta) {//ÏßÐÔÖµ
 	return  p1 + p2 * exp(-p3 * theta) + p4 * cos(p5 * theta + p6);
@@ -32,7 +28,7 @@ extern __host__ __device__ double GPU_getSigma0dB(
 
 
 
-extern __host__ __device__ CUDAVectorEllipsoidal GPU_SatelliteAntDirectNormal(
+extern  __device__ CUDAVectorEllipsoidal GPU_SatelliteAntDirectNormal(
 	double RstX, double RstY, double RstZ,
 	double  AntXaxisX, double  AntXaxisY, double  AntXaxisZ,
 	double  AntYaxisX, double  AntYaxisY, double  AntYaxisZ,
@@ -111,7 +107,7 @@ extern __host__ __device__ CUDAVectorEllipsoidal GPU_SatelliteAntDirectNormal(
 	return result;
 }
 
-extern __host__ __device__ double GPU_BillerInterpAntPattern(double* antpattern,
+extern  __device__ double GPU_BillerInterpAntPattern(double* antpattern,
 	double starttheta, double startphi, double dtheta, double dphi,
 	long thetapoints, long phipoints,
 	double searththeta, double searchphi) {

Toolbox/SimulationSARTool/SimulationSAR/GPURFPC.cuh

@@ -22,13 +22,13 @@ extern "C" struct  CUDASigmaParam {
 
  
 
-extern __host__ __device__ double GPU_getSigma0dB(
+extern  __device__ double GPU_getSigma0dB_params(
 	const double p1, const double p2, const  double p3, const  double p4, const  double p5, const  double p6,
 	double theta);
 
-extern __host__ __device__ double GPU_getSigma0dB(CUDASigmaParam param, double theta);
+extern  __device__ double GPU_getSigma0dB(CUDASigmaParam param, double theta);
 
-extern __host__ __device__ CUDAVectorEllipsoidal GPU_SatelliteAntDirectNormal(
+extern  __device__ CUDAVectorEllipsoidal GPU_SatelliteAntDirectNormal(
 	double RstX, double RstY, double RstZ,
 	double  AntXaxisX, double  AntXaxisY, double  AntXaxisZ,
 	double  AntYaxisX, double  AntYaxisY, double  AntYaxisZ,
@@ -36,7 +36,7 @@ extern __host__ __device__ CUDAVectorEllipsoidal GPU_SatelliteAntDirectNormal(
 	double  AntDirectX, double  AntDirectY, double  AntDirectZ
 );
 
-extern __host__ __device__ double GPU_BillerInterpAntPattern(double* antpattern,
+extern  __device__ double GPU_BillerInterpAntPattern(double* antpattern,
 	double starttheta, double startphi, double dtheta, double dphi,
 	long thetapoints, long phipoints,
 	double searththeta, double searchphi);

Toolbox/SimulationSARTool/SimulationSAR/GPURFPC_single.cu

@@ -0,0 +1,468 @@
+#include <cuda.h>
+#include <device_launch_parameters.h>
+#include <cuda_runtime.h>
+#include <cublas_v2.h>
+#include <cuComplex.h>
+
+#include "BaseConstVariable.h"
+#include "GPURFPC_single.cuh"
+
+
+#ifdef __CUDANVCC___
+
+
+/* 机器函数  ****************************************************************************************************************************/
+
+
+extern  __device__ float GPU_getSigma0dB_single(CUDASigmaParam_single param, float theta) {//线性值
+	float sigma = param.p1 + param.p2 * exp(-param.p3 * theta) + param.p4 * cos(param.p5 * theta + param.p6);
+	return  sigma;
+}
+
+
+extern  __device__ float GPU_getSigma0dB_single(const float p1, const float p2, const float p3, const float p4, const float p5, const float p6, float theta)
+{
+	return  p1 + p2 * expf(-p3 * theta) + p4 * cosf(p5 * theta + p6);
+}
+
+ 
+
+
+extern  __device__ CUDAVectorEllipsoidal GPU_SatelliteAntDirectNormal_single(
+	float RstX, float RstY, float RstZ,
+	float  AntXaxisX, float  AntXaxisY, float  AntXaxisZ,
+	float  AntYaxisX, float  AntYaxisY, float  AntYaxisZ,
+	float  AntZaxisX, float  AntZaxisY, float  AntZaxisZ,
+	float  AntDirectX, float  AntDirectY, float  AntDirectZ
+) {
+	CUDAVectorEllipsoidal result{ 0,0,-1 };
+
+	// 求解天线增益
+	float Xst = -1 * RstX; // 卫星 -->  地面
+	float Yst = -1 * RstY;
+	float Zst = -1 * RstZ;
+
+	// 归一化
+	float RstNorm = sqrtf(Xst * Xst + Yst * Yst + Zst * Zst);
+	float AntXaxisNorm = sqrtf(AntXaxisX * AntXaxisX + AntXaxisY * AntXaxisY + AntXaxisZ * AntXaxisZ);
+	float AntYaxisNorm = sqrtf(AntYaxisX * AntYaxisX + AntYaxisY * AntYaxisY + AntYaxisZ * AntYaxisZ);
+	float AntZaxisNorm = sqrtf(AntZaxisX * AntZaxisX + AntZaxisY * AntZaxisY + AntZaxisZ * AntZaxisZ);
+
+
+	float Rx = Xst / RstNorm;
+	float Ry = Yst / RstNorm;
+	float Rz = Zst / RstNorm;
+	float Xx = AntXaxisX / AntXaxisNorm;
+	float Xy = AntXaxisY / AntXaxisNorm;
+	float Xz = AntXaxisZ / AntXaxisNorm;
+	float Yx = AntYaxisX / AntYaxisNorm;
+	float Yy = AntYaxisY / AntYaxisNorm;
+	float Yz = AntYaxisZ / AntYaxisNorm;
+	float Zx = AntZaxisX / AntZaxisNorm;
+	float Zy = AntZaxisY / AntZaxisNorm;
+	float Zz = AntZaxisZ / AntZaxisNorm;
+
+	float Xant = (Rx * Yy * Zz - Rx * Yz * Zy - Ry * Yx * Zz + Ry * Yz * Zx + Rz * Yx * Zy - Rz * Yy * Zx) / (Xx * Yy * Zz - Xx * Yz * Zy - Xy * Yx * Zz + Xy * Yz * Zx + Xz * Yx * Zy - Xz * Yy * Zx);
+	float Yant = -(Rx * Xy * Zz - Rx * Xz * Zy - Ry * Xx * Zz + Ry * Xz * Zx + Rz * Xx * Zy - Rz * Xy * Zx) / (Xx * Yy * Zz - Xx * Yz * Zy - Xy * Yx * Zz + Xy * Yz * Zx + Xz * Yx * Zy - Xz * Yy * Zx);
+	float Zant = (Rx * Xy * Yz - Rx * Xz * Yy - Ry * Xx * Yz + Ry * Xz * Yx + Rz * Xx * Yy - Rz * Xy * Yx) / (Xx * Yy * Zz - Xx * Yz * Zy - Xy * Yx * Zz + Xy * Yz * Zx + Xz * Yx * Zy - Xz * Yy * Zx);
+
+
+
+	// 计算theta 与 phi
+	float Norm = sqrtf(Xant * Xant + Yant * Yant + Zant * Zant); // 计算 pho
+	float Zn = Zant / Norm;
+	float ThetaAnt = ( - 1 > Zn) ? PI : (Zn > 1 ? 0 : acos(Zn));//  acosf(Zant / Norm); // theta 与 Z轴的夹角
+	float PhiAnt = abs(Xant)<PRECISIONTOLERANCE ?0: atanf(Yant / Xant); // -pi/2 ~pi/2
+
+	if (abs(Yant) < PRECISIONTOLERANCE) { // X轴上
+		PhiAnt = 0;
+	}
+	else if (abs(Xant) < PRECISIONTOLERANCE) { // Y轴上，原点
+		if (Yant > 0) {
+			PhiAnt = PI / 2;
+		}
+		else {
+			PhiAnt = -PI / 2;
+		}
+	}
+	else if (Xant < 0) {
+		if (Yant > 0) {
+			PhiAnt = PI + PhiAnt;
+		}
+		else {
+			PhiAnt = -PI + PhiAnt;
+		}
+	}
+	else {  // Xant>0  X 正轴
+
+	}
+
+	if (isnan(PhiAnt)) {
+		printf("V=[%f,%f,%f];norm=%f;thetaAnt=%f;phiAnt=%f;\n", Xant, Yant, Zant, Norm, ThetaAnt, PhiAnt);
+	}
+
+	result.theta = ThetaAnt;
+	result.phi = PhiAnt;
+	result.Rho = Norm;
+	return result;
+}
+
+extern  __device__ float GPU_BillerInterpAntPattern_single(float* antpattern,
+	float starttheta, float startphi, float dtheta, float dphi,
+	long thetapoints, long phipoints,
+	float searththeta, float searchphi) {
+	float stheta = searththeta;
+	float sphi = searchphi;
+	if (stheta > 90) {
+		return 0;
+	}
+	else {}
+
+
+	float pthetaid = (stheta - starttheta) / dtheta;// 
+	float pphiid = (sphi - startphi) / dphi;
+
+	long lasttheta = floorf(pthetaid);
+	long nextTheta = lasttheta + 1;
+	long lastphi = floorf(pphiid);
+	long nextPhi = lastphi + 1;
+
+
+	if (lasttheta < 0 || nextTheta < 0 || lastphi < 0 || nextPhi < 0 ||
+		lasttheta >= thetapoints || nextTheta >= thetapoints || lastphi >= phipoints || nextPhi >= phipoints)
+	{
+		return 0;
+	}
+	else {
+		float x = stheta;
+		float y = sphi;
+
+		float x1 = lasttheta * dtheta + starttheta;
+		float x2 = nextTheta * dtheta + starttheta;
+		float y1 = lastphi * dphi + startphi;
+		float y2 = nextPhi * dphi + startphi;
+
+		float z11 = antpattern[lasttheta * phipoints + lastphi];
+		float z12 = antpattern[lasttheta * phipoints + nextPhi];
+		float z21 = antpattern[nextTheta * phipoints + lastphi];
+		float z22 = antpattern[nextTheta * phipoints + nextPhi];
+
+
+		//z11 = powf(10, z11 / 10); // dB-> 线性
+		//z12 = powf(10, z12 / 10);
+		//z21 = powf(10, z21 / 10);
+		//z22 = powf(10, z22 / 10);
+
+		float GainValue = (z11 * (x2 - x) * (y2 - y)
+			+ z21 * (x - x1) * (y2 - y)
+			+ z12 * (x2 - x) * (y - y1)
+			+ z22 * (x - x1) * (y - y1));
+		GainValue = GainValue / ((x2 - x1) * (y2 - y1));
+		return GainValue;
+	}
+}
+ 
+
+ 
+/* 核函数 ****************************************************************************************************************************/
+// 计算每块
+__global__ void CUDA_Kernel_Computer_R_amp_single(
+	float* antX, float* antY, float* antZ,
+	float* antXaxisX, float* antXaxisY, float* antXaxisZ,
+	float* antYaxisX, float* antYaxisY, float* antYaxisZ,
+	float* antZaxisX, float* antZaxisY, float* antZaxisZ,
+	float* antDirectX, float* antDirectY, float* antDirectZ,
+	long PRFCount, // 整体的脉冲数，
+	float* targetX, float* targetY, float* targetZ, long* demCls, 
+	float* demSlopeX, float* demSlopeY, float* demSlopeZ ,
+	long startPosId, long pixelcount,
+	CUDASigmaParam_single* sigma0Paramslist, long sigmaparamslistlen,
+	float Pt,
+	float refPhaseRange,
+	float* TransAntpattern,
+	float Transtarttheta, float Transstartphi, float Transdtheta, float Transdphi, int Transthetapoints, int Transphipoints,
+	float* ReceiveAntpattern,
+	float Receivestarttheta, float Receivestartphi, float Receivedtheta, float Receivedphi, int Receivethetapoints, int Receivephipoints,
+	float maxTransAntPatternValue, float maxReceiveAntPatternValue,
+	float NearR, float FarR,
+	float* d_temp_R, float* d_temp_amps// 计算输出
+) {
+	long idx = blockIdx.x * blockDim.x + threadIdx.x; // 获取当前的线程编码
+	long prfId = idx / SHAREMEMORY_FLOAT_HALF;
+	long posId = idx % SHAREMEMORY_FLOAT_HALF+ startPosId; // 当前线程对应的影像点
+	
+	if (prfId < PRFCount && posId < pixelcount) {
+		float RstX = antX[prfId] - targetX[posId]; // 计算坐标矢量
+		float RstY = antY[prfId] - targetY[posId];
+		float RstZ = antZ[prfId] - targetZ[posId];
+
+		float RstR = sqrt(RstX * RstX + RstY * RstY + RstZ * RstZ); // 矢量距离
+		if (RstR<NearR || RstR>FarR) {
+			d_temp_R[idx] = 0;
+			d_temp_amps[idx] = 0;
+			return;
+		}
+		else {
+			float slopeX = demSlopeX[posId];
+			float slopeY = demSlopeY[posId];
+			float slopeZ = demSlopeZ[posId];
+ 
+			float slopR = sqrtf(slopeX * slopeX + slopeY * slopeY + slopeZ * slopeZ); //  
+			if (abs(slopR - 0) > 1e-3) {
+				float dotAB = RstX * slopeX + RstY * slopeY + RstZ * slopeZ;
+				float localangle = acos(dotAB / (RstR * slopR));  
+
+				if (localangle < 0 || localangle >= LAMP_CUDA_PI / 2|| isnan(localangle)) {
+					d_temp_R[idx] = 0;
+					d_temp_amps[idx] = 0;
+					return;
+				}
+				else {}
+
+
+				float ampGain = 0;
+				// 求解天线方向图指向
+				CUDAVectorEllipsoidal antVector = GPU_SatelliteAntDirectNormal_single(
+					RstX, RstY, RstZ,
+					antXaxisX[prfId], antXaxisY[prfId], antXaxisZ[prfId],
+					antYaxisX[prfId], antYaxisY[prfId], antYaxisZ[prfId],
+					antZaxisX[prfId], antZaxisY[prfId], antZaxisZ[prfId],
+					antDirectX[prfId], antDirectY[prfId], antDirectZ[prfId]
+				);
+				antVector.theta = antVector.theta * r2d;
+				antVector.phi = antVector.phi * r2d;
+				//printf("theta: %f , phi: %f \n", antVector.theta, antVector.phi);
+				if (antVector.Rho > 0) {
+					//float TansantPatternGain = GPU_BillerInterpAntPattern(
+					//	TransAntpattern,
+					//	Transtarttheta, Transstartphi, Transdtheta, Transdphi, Transthetapoints, Transphipoints,
+					//	antVector.theta, antVector.phi);
+					//float antPatternGain = GPU_BillerInterpAntPattern(
+					//	ReceiveAntpattern,
+					//	Receivestarttheta, Receivestartphi, Receivedtheta, Receivedphi, Receivethetapoints, Receivephipoints,
+					//	antVector.theta, antVector.phi);
+
+					float sigma0 = 0;
+					{
+						long clsid = demCls[posId];
+						//printf("clsid=%d\n", clsid);
+						CUDASigmaParam_single tempsigma = sigma0Paramslist[clsid];
+
+						
+						if (abs(tempsigma.p1) < PRECISIONTOLERANCE &&
+							abs(tempsigma.p2) < PRECISIONTOLERANCE &&
+							abs(tempsigma.p3) < PRECISIONTOLERANCE &&
+							abs(tempsigma.p4) < PRECISIONTOLERANCE &&
+							abs(tempsigma.p5) < PRECISIONTOLERANCE &&
+							abs(tempsigma.p6) < PRECISIONTOLERANCE
+							) {
+							sigma0 = 0;
+						}
+						else {
+							float sigma = GPU_getSigma0dB_single(tempsigma, localangle);
+							sigma0 = powf(10.0, sigma / 10.0);
+						}
+					}
+					//ampGain = TansantPatternGain * antPatternGain;
+					ampGain = 1;
+					//if (10 * log10(ampGain / maxReceiveAntPatternValue / maxTransAntPatternValue) < -3) { // 小于-3dB
+					//	d_temp_R[idx] = 0;
+					//	d_temp_amps[idx] = 0;
+					//	return;
+					//}
+					//else {}
+
+
+
+
+					ampGain = ampGain / (powf(4 * LAMP_CUDA_PI, 2) * powf(RstR, 4)); // 反射强度
+					
+					float temp_amp = float(ampGain * Pt * sigma0);
+					float temp_R = float(RstR - refPhaseRange);
+
+					if (isnan(temp_amp) || isnan(temp_R)|| isinf(temp_amp) || isinf(temp_R)) {
+						printf("amp is nan or R is nan,amp=%f;R=%f; \n", temp_amp, temp_R);
+						d_temp_R[idx] = 0;
+						d_temp_amps[idx] = 0;
+						return;
+					}
+					else {}
+					
+					
+					d_temp_amps[idx] = temp_amp;
+					d_temp_R[idx] = temp_R;
+					return;
+				}
+				else {
+					d_temp_R[idx] = 0;
+					d_temp_amps[idx] = 0;
+					return;
+				}
+			}
+			else {
+				d_temp_R[idx] = 0;
+				d_temp_amps[idx] = 0;
+				return;
+			}
+		}
+	}
+}
+
+
+
+
+__global__ void CUDA_Kernel_Computer_echo_single(
+	float* d_temp_R, float* d_temp_amps, long posNum,
+	float f0, float dfreq, 
+	long FreqPoints, // 当前频率的分块
+	long maxfreqnum, // 最大脉冲值
+	float* d_temp_echo_real, float* d_temp_echo_imag,
+	long temp_PRF_Count
+) { 
+	__shared__ float s_R[SHAREMEMORY_FLOAT_HALF];  // 注意一个完整的block_size 共享相同内存
+	__shared__ float s_amp[SHAREMEMORY_FLOAT_HALF];
+
+	long tid = threadIdx.x;
+	long bid = blockIdx.x;
+	long idx = bid * blockDim.x + tid;
+	long prfId = idx / FreqPoints; // 脉冲ID
+	long fId = idx % FreqPoints;//频率ID
+
+	long psid = 0;
+	long pixelId = 0;
+	for (long ii = 0; ii < SHAREMEMORY_FLOAT_HALF_STEP; ii++) { // SHAREMEMORY_FLOAT_HALF_STEP * BLOCK_SIZE=SHAREMEMORY_FLOAT_HALF
+		psid = tid * SHAREMEMORY_FLOAT_HALF_STEP + ii;
+		pixelId = prfId * posNum + psid; // 
+		if (psid < posNum) {
+			s_R[psid] = d_temp_R[pixelId];
+			s_amp[psid] = d_temp_amps[pixelId];
+		}
+		else {
+			s_R[psid] = 0;
+			s_amp[psid] = 0;
+		}
+
+	}
+
+	__syncthreads(); // 确定所有待处理数据都已经进入程序中
+
+
+
+	if (fId < maxfreqnum && prfId < temp_PRF_Count) {
+		 
+		long echo_ID = prfId * maxfreqnum + fId; // 计算对应的回波位置
+		float factorjTemp = RFPCPIDIVLIGHT * (f0 + fId * dfreq);
+		float temp_real = 0;
+		float temp_imag = 0;
+		float temp_phi = 0;
+		float temp_amp = 0;
+		for (long dataid = 0; dataid < SHAREMEMORY_FLOAT_HALF; dataid++) {
+
+			temp_phi = s_R[dataid] * factorjTemp;
+			temp_amp = s_amp[dataid];
+			temp_real += (temp_amp * cosf(temp_phi));
+			temp_imag += (temp_amp * sinf(temp_phi));
+			//if (dataid > 5000) {
+			//	printf("echo_ID=%d; dataid=%d;ehodata=(%f,%f);R=%f;amp=%f;\n", echo_ID, dataid, temp_real, temp_imag, s_R[0], s_amp[0]);
+			//}
+			if (isnan(temp_phi) || isnan(temp_amp) || isnan(temp_real) || isnan(temp_imag)
+				|| isinf(temp_phi) || isinf(temp_amp) || isinf(temp_real) || isinf(temp_imag)
+				) {
+				printf("[amp,phi,real,imag]=[%f,%f,%f,%f];\n",temp_amp,temp_phi,temp_real,temp_imag);
+			}
+			
+		}
+		//printf("echo_ID=%d; ehodata=(%f,%f)\n", echo_ID, temp_real, temp_imag);
+		//printf("(%f %f %f) ", factorjTemp, s_amp[0], s_R[0]);
+		d_temp_echo_real[echo_ID] += /*d_temp_echo_real[echo_ID] + */temp_real;
+		d_temp_echo_imag[echo_ID] += /*d_temp_echo_imag[echo_ID] +*/ temp_imag;
+	}
+}
+
+
+
+/**
+* 分块计算主流程
+ */
+void CUDA_RFPC_MainProcess_single(
+	float* antX, float* antY, float* antZ,
+	float* antXaxisX, float* antXaxisY, float* antXaxisZ,
+	float* antYaxisX, float* antYaxisY, float* antYaxisZ,
+	float* antZaxisX, float* antZaxisY, float* antZaxisZ,
+	float* antDirectX, float* antDirectY, float* antDirectZ,
+	long PRFCount, long FreqNum,  
+	float f0, float dfreq,
+	float Pt,
+	float refPhaseRange,
+	float* TransAntpattern,
+	float Transtarttheta, float Transstartphi, float Transdtheta, float Transdphi, int Transthetapoints, int Transphipoints,
+	float* ReceiveAntpattern,
+	float Receivestarttheta, float Receivestartphi, float Receivedtheta, float Receivedphi, int Receivethetapoints, int Receivephipoints,
+	float maxTransAntPatternValue, float maxReceiveAntPatternValue,
+	float NearR, float FarR,
+	float* targetX, float* targetY, float* targetZ, long* demCls, long TargetNumber,
+	float* demSlopeX, float* demSlopeY, float* demSlopeZ, 
+	CUDASigmaParam_single* sigma0Paramslist, long sigmaparamslistlen,
+	float* out_echoReal, float* out_echoImag,
+	float* d_temp_R, float* d_temp_amp
+)
+{
+	long BLOCK_FREQNUM = NextBlockPad(FreqNum, BLOCK_SIZE); // 256*freqBlockID
+	long cudaBlocknum = 0;
+	long freqpoints = BLOCK_FREQNUM;
+	printf("freqpoints:%d\n", freqpoints);
+	long process = 0;
+	for (long sTi = 0; sTi < TargetNumber; sTi = sTi + SHAREMEMORY_FLOAT_HALF) {
+		cudaBlocknum = (PRFCount * SHAREMEMORY_FLOAT_HALF + BLOCK_SIZE - 1) / BLOCK_SIZE;
+		CUDA_Kernel_Computer_R_amp_single << <cudaBlocknum, BLOCK_SIZE >> > (
+			antX, antY, antZ,
+			antXaxisX, antXaxisY, antXaxisZ,
+			antYaxisX, antYaxisY, antYaxisZ,
+			antZaxisX, antZaxisY, antZaxisZ,
+			antDirectX, antDirectY, antDirectZ,
+			PRFCount, 
+			targetX, targetY, targetZ, demCls, 
+			demSlopeX, demSlopeY, demSlopeZ,  
+			sTi, TargetNumber,
+			sigma0Paramslist, sigmaparamslistlen,
+			Pt,
+			refPhaseRange,
+			TransAntpattern,
+			Transtarttheta, Transstartphi, Transdtheta, Transdphi, Transthetapoints, Transphipoints,
+			ReceiveAntpattern,
+			Receivestarttheta, Receivestartphi, Receivedtheta, Receivedphi, Receivethetapoints, Receivephipoints,
+			maxTransAntPatternValue, maxReceiveAntPatternValue,
+			NearR, FarR,
+			d_temp_R, d_temp_amp// 计算输出
+			);
+
+		PrintLasterError("CUDA_Kernel_Computer_R_amp");
+ 
+		
+		cudaBlocknum = (PRFCount * BLOCK_FREQNUM + BLOCK_SIZE - 1) / BLOCK_SIZE;
+		CUDA_Kernel_Computer_echo_single << <cudaBlocknum, BLOCK_SIZE >> > (
+			d_temp_R, d_temp_amp, SHAREMEMORY_FLOAT_HALF,
+			f0, dfreq, 
+			freqpoints, FreqNum,
+			out_echoReal, out_echoImag,
+			PRFCount
+			);
+		PrintLasterError("CUDA_Kernel_Computer_echo");
+
+		if ((sTi * 100.0 / TargetNumber ) - process >= 1) {
+			process = sTi * 100.0 / TargetNumber;
+			PRINT("TargetID [%f]: %d / %d finished\n", sTi*100.0/ TargetNumber,sTi, TargetNumber);
+		}
+		
+
+
+	}
+
+
+	cudaDeviceSynchronize();
+}
+
+
+#endif
+
+

Toolbox/SimulationSARTool/SimulationSAR/GPURFPC_single.cuh

@@ -0,0 +1,85 @@
+#ifndef _GPURFPC_SINGLE_H_
+#define _GPURFPC_SINGLE_H_
+
+#include "BaseConstVariable.h"
+#include "GPUTool.cuh"
+#include <cuda_runtime.h>
+#include <device_launch_parameters.h>
+#include <cublas_v2.h>
+#include <cuComplex.h>
+
+#define RFPCPIDIVLIGHT -4*PI/(LIGHTSPEED/1e9)
+
+ 
+extern "C" struct  CUDASigmaParam_single {
+	float p1;
+	float p2;
+	float p3;
+	float p4;
+	float p5;
+	float p6;
+};
+
+extern  __device__ float GPU_getSigma0dB_single(const float p1, const float p2, const float p3, const float p4, const float p5, const float p6, float theta)
+;
+
+extern  __device__ float GPU_getSigma0dB_single(CUDASigmaParam_single param, float theta);
+
+extern  __device__ CUDAVectorEllipsoidal GPU_SatelliteAntDirectNormal_single(
+	float RstX, float RstY, float RstZ,
+	float  AntXaxisX, float  AntXaxisY, float  AntXaxisZ,
+	float  AntYaxisX, float  AntYaxisY, float  AntYaxisZ,
+	float  AntZaxisX, float  AntZaxisY, float  AntZaxisZ,
+	float  AntDirectX, float  AntDirectY, float  AntDirectZ
+);
+
+extern  __device__ float GPU_BillerInterpAntPattern_single(float* antpattern,
+	float starttheta, float startphi, float dtheta, float dphi,
+	long thetapoints, long phipoints,
+	float searththeta, float searchphi);
+
+
+
+extern "C" void CUDA_RFPC_MainProcess_single(
+	// 天线
+	float* antX, float* antY, float* antZ, // 天线坐标
+	float* antXaxisX, float* antXaxisY, float* antXaxisZ, // 天线坐标系的X轴
+	float* antYaxisX, float* antYaxisY, float* antYaxisZ,// 天线坐标系的Y轴
+	float* antZaxisX, float* antZaxisY, float* antZaxisZ,// 天线坐标系的Z轴
+	float* antDirectX, float* antDirectY, float* antDirectZ,// 天线的指向
+	long PRFCount, long FreqNum,   // 脉冲数量，频率数量
+	float f0, float dfreq,// 起始频率，终止频率
+	float Pt,// 发射能量
+	float refPhaseRange,
+	// 天线方向图
+	float* TransAntpattern, float Transtarttheta, float Transstartphi, float Transdtheta, float Transdphi, int Transthetapoints, int Transphipoints, // 发射天线方向图
+	float* ReceiveAntpattern, float Receivestarttheta, float Receivestartphi, float Receivedtheta, float Receivedphi, int Receivethetapoints, int Receivephipoints,//接收天线方向图
+	float maxTransAntPatternValue,float maxReceiveAntPatternValue,
+	float NearR, float FarR, // 距离范围
+
+
+	// 地面
+	float* targetX, float* targetY, float* targetZ, long* demCls, long TargetPixelNumber, // 地面坐标、地表覆盖类型，像素数
+	float* demSlopeX, float* demSlopeY, float* demSlopeZ, // 地表坡度矢量
+	CUDASigmaParam_single* sigma0Paramslist, long sigmaparamslistlen,// 插值图像
+ 
+	float* out_echoReal, float* out_echoImag,// 输出回波
+
+	float* d_temp_R, float* d_temp_amp
+);
+
+
+
+
+
+
+
+
+
+
+
+#endif
+
+
+
+ 

Toolbox/SimulationSARTool/SimulationSARTool.vcxproj

@@ -41,7 +41,7 @@
     <CharacterSet>Unicode</CharacterSet>
   </PropertyGroup>
   <PropertyGroup Condition="'$(Configuration)|$(Platform)' == 'Release|x64'" Label="Configuration">
-    <ConfigurationType>DynamicLibrary</ConfigurationType>
+    <ConfigurationType>Application</ConfigurationType>
     <PlatformToolset>v143</PlatformToolset>
     <UseDebugLibraries>false</UseDebugLibraries>
     <WholeProgramOptimization>true</WholeProgramOptimization>
@@ -222,6 +222,7 @@
     <ClCompile Include="SimulationSAR\TBPImageAlgCls.cpp" />
     <ClCompile Include="UnitTestMain.cpp" />
     <CudaCompile Include="GPUBpSimulation.cu" />
+    <CudaCompile Include="SimulationSAR\GPURFPC_single.cu" />
     <CudaCompile Include="SimulationSAR\GPUTBPImage.cu" />
     <QtMoc Include="PowerSimulationIncoherent\QSimulationSARPolynomialOrbitModel.h" />
     <CudaCompile Include="PowerSimulationIncoherent\LookTableSimulationComputer.cuh" />
@@ -238,6 +239,7 @@
     <CudaCompile Include="SimulationSAR\BPBasic0_CUDA.cu" />
     <CudaCompile Include="SimulationSAR\GPUBPTool.cuh" />
     <ClInclude Include="SimulationSAR\BPBasic0_CUDA.cuh" />
+    <ClInclude Include="SimulationSAR\GPURFPC_single.cuh" />
     <ClInclude Include="SimulationSAR\GPUTBPImage.cuh" />
     <ClInclude Include="SimulationSAR\RFPCProcessCls.h" />
     <ClInclude Include="SimulationSAR\SARSatelliteSimulationAbstractCls.h" />

Toolbox/SimulationSARTool/SimulationSARTool.vcxproj.filters

@@ -68,6 +68,9 @@
     <ClInclude Include="GPUBpSimulation.cuh">
       <Filter>SimulationSAR</Filter>
     </ClInclude>
+    <ClInclude Include="SimulationSAR\GPURFPC_single.cuh">
+      <Filter>SimulationSAR</Filter>
+    </ClInclude>
   </ItemGroup>
   <ItemGroup>
     <ClCompile Include="SimulationSAR\QImageSARRFPC.cpp">
@@ -193,5 +196,8 @@
     <CudaCompile Include="GPUBpSimulation.cu">
       <Filter>SimulationSAR</Filter>
     </CudaCompile>
+    <CudaCompile Include="SimulationSAR\GPURFPC_single.cu">
+      <Filter>SimulationSAR</Filter>
+    </CudaCompile>
   </ItemGroup>
 </Project>

Toolbox/SimulationSARTool/UnitTestMain.cpp

@@ -44,7 +44,7 @@ void testSimualtionEchoPoint() {
 		.arg(h_data.AntZ[gpspoints - 1]);
 	/** 2. 成像网格 **************************************************************************************************/
 	qDebug() << u8"轨道文件读取结束\n2.成像网格读取。。。";
-	QString demxyzPath = u8"D:\\Programme\\vs2022\\RasterMergeTest\\simulationData\\demdataset\\demxyz.bin";
+	QString demxyzPath = u8"C:\\Users\\30453\\Desktop\\script\\已修改GF3场景\\data\\demxyz.bin";
 	gdalImage demgridimg(demxyzPath);
 	long dem_rowCount = demgridimg.height;
 	long dem_ColCount = demgridimg.width;
@@ -115,7 +115,7 @@ void testSimualtionEchoPoint() {
 	/** 7. 目标 **************************************************************************************************/
 	double Tx = -2029086.618142, Ty = 4139594.934504, Tz = 4392846.782027;
 	double Tslx = -2029086.618142, Tsly = 4139594.934504, Tslz = 4392846.782027;
-	double p1 = 33.3, p2 = 0, p3 = 0, p4 = 0, p5 = 0, p6 = 0;
+	double p1 = 1, p2 = 0, p3 = 0, p4 = 0, p5 = 0, p6 = 0;
 
 	/** 7. 构建回波 **************************************************************************************************/
 	GPUDATA d_data;
@@ -405,14 +405,405 @@ void testBpImage() {
 }
 
 
-//int main(int argc, char* argv[]) {
+
-//
+void testSimualtionEchoPoint_single() {
-//	QApplication a(argc, argv);
+	GPUDATA_single h_data{};
-//	
+	/** 1. 轨道 **************************************************************************************************/
-//	testBpImage();
+	qDebug() << u8"1.轨道文件读取中。。。";
-//
+	QString inGPSPath = u8"C:\\Users\\30453\\Desktop\\script\\data\\GF3_Simulation.gpspos.data";
-//	return 0;
+	long gpspoints = gdalImage(inGPSPath).height;
-//}
+	std::shared_ptr<SatelliteAntPos> antpos = SatelliteAntPosOperator::readAntPosFile(inGPSPath, gpspoints);
+	h_data.AntX = (float*)mallocCUDAHost(sizeof(float) * gpspoints);
+	h_data.AntY = (float*)mallocCUDAHost(sizeof(float) * gpspoints);
+	h_data.AntZ = (float*)mallocCUDAHost(sizeof(float) * gpspoints);
+
+	for (long i = 0; i < gpspoints; i = i + 1) {
+		h_data.AntX[i] = antpos.get()[i].Px;
+		h_data.AntY[i] = antpos.get()[i].Py;
+		h_data.AntZ[i] = antpos.get()[i].Pz;
+	}
+	//gpspoints = gpspoints / 2;
+	qDebug() << "GPS points Count:\t" << gpspoints;
+	qDebug() << "PRF 0:\t " << QString("%1,%2,%3").arg(h_data.AntX[0]).arg(h_data.AntY[0]).arg(h_data.AntZ[0]);
+	qDebug() << "PRF " << gpspoints - 1 << ":\t " << QString("%1,%2,%3")
+		.arg(h_data.AntX[gpspoints - 1])
+		.arg(h_data.AntY[gpspoints - 1])
+		.arg(h_data.AntZ[gpspoints - 1]);
+	/** 2. 成像网格 **************************************************************************************************/
+	qDebug() << u8"轨道文件读取结束\n2.成像网格读取。。。";
+	QString demxyzPath = u8"C:\\Users\\30453\\Desktop\\script\\已修改GF3场景\\data\\demxyz.bin";
+	gdalImage demgridimg(demxyzPath);
+	long dem_rowCount = demgridimg.height;
+	long dem_ColCount = demgridimg.width;
+	std::shared_ptr<float> demX = readDataArr<float>(demgridimg, 0, 0, dem_rowCount, dem_ColCount, 1, GDALREADARRCOPYMETHOD::VARIABLEMETHOD);
+	std::shared_ptr<float> demY = readDataArr<float>(demgridimg, 0, 0, dem_rowCount, dem_ColCount, 2, GDALREADARRCOPYMETHOD::VARIABLEMETHOD);
+	std::shared_ptr<float> demZ = readDataArr<float>(demgridimg, 0, 0, dem_rowCount, dem_ColCount, 3, GDALREADARRCOPYMETHOD::VARIABLEMETHOD);
+
+	h_data.x_mat = (float*)mallocCUDAHost(sizeof(float) * dem_rowCount * dem_ColCount);
+	h_data.y_mat = (float*)mallocCUDAHost(sizeof(float) * dem_rowCount * dem_ColCount);
+	h_data.z_mat = (float*)mallocCUDAHost(sizeof(float) * dem_rowCount * dem_ColCount);
+
+	for (long i = 0; i < dem_rowCount; i++) {
+		for (long j = 0; j < dem_ColCount; j++) {
+			h_data.x_mat[i * dem_ColCount + j] = demX.get()[i * dem_ColCount + j];
+			h_data.y_mat[i * dem_ColCount + j] = demY.get()[i * dem_ColCount + j];
+			h_data.z_mat[i * dem_ColCount + j] = demZ.get()[i * dem_ColCount + j];
+		}
+	}
+	qDebug() << "dem row Count:\t" << dem_rowCount;
+	qDebug() << "dem col Count:\t" << dem_ColCount;
+
+	qDebug() << u8"成像网格读取结束";
+	/** 3. 频率网格 **************************************************************************************************/
+	qDebug() << u8"3.处理频率";
+	float centerFreq = 5.3e9;
+	float bandwidth = 40e6;
+	long freqpoints = 2048;
+	std::shared_ptr<float> freqlist(getFreqPoints_mallocHost_single(centerFreq - bandwidth / 2, centerFreq + bandwidth / 2, freqpoints), FreeCUDAHost);
+	std::shared_ptr<float> minF(new float[gpspoints], delArrPtr);
+	for (long i = 0; i < gpspoints; i++) {
+		minF.get()[i] = freqlist.get()[0];
+	}
+	h_data.deltaF = bandwidth / (freqpoints - 1);
+	qDebug() << "start Freq:\t" << centerFreq - bandwidth / 2;
+	qDebug() << "end   Freq:\t" << centerFreq + bandwidth / 2;
+	qDebug() << "freq points:\t" << freqpoints;
+	qDebug() << "delta freq:\t" << freqlist.get()[1] - freqlist.get()[0];
+	qDebug() << u8"频率结束";
+	/** 4. 初始化回波 **************************************************************************************************/
+	qDebug() << u8"4.初始化回波";
+	std::shared_ptr<cuComplex> phdata(createEchoPhase_mallocHost(gpspoints, freqpoints), FreeCUDAHost);
+	qDebug() << "Azimuth Points:\t" << gpspoints;
+	qDebug() << "Range Points:\t" << freqpoints;
+	qDebug() << u8"初始化回波结束";
+	/** 5. 初始化图像 **************************************************************************************************/
+	qDebug() << u8"5.初始化图像";
+	h_data.im_final = (cuComplex*)mallocCUDAHost(sizeof(cuComplex) * dem_rowCount * dem_ColCount);
+	qDebug() << "Azimuth Points:\t" << gpspoints;
+	qDebug() << "Range Points:\t" << freqpoints;
+	qDebug() << u8"初始化图像结束";
+
+	/** 6. 模型参数初始化 **************************************************************************************************/
+	qDebug() << u8"6.模型参数初始化";
+
+	h_data.nx = dem_ColCount;
+	h_data.ny = dem_rowCount;
+
+	h_data.Np = gpspoints;
+	h_data.Freq = freqlist.get();
+	h_data.minF = minF.get();
+	h_data.Nfft = freqpoints;
+	h_data.K = h_data.Nfft;
+	h_data.phdata = phdata.get();
+	h_data.R0 = 900000;
+	qDebug() << u8"模型参数结束";
+
+
+	/** 7. 目标 **************************************************************************************************/
+	float Tx = -2029086.618142, Ty = 4139594.934504, Tz = 4392846.782027;
+	float Tslx = -2029086.618142, Tsly = 4139594.934504, Tslz = 4392846.782027;
+	float p1 = 1, p2 = 0, p3 = 0, p4 = 0, p5 = 0, p6 = 0;
+
+	/** 7. 构建回波 **************************************************************************************************/
+	GPUDATA_single d_data;
+	initGPUData_single(h_data, d_data);
+
+	RFPCProcess_single(Tx, Ty, Tz,
+		Tslx, Tsly, Tslz, // 目标的坡面向量
+		p1, p2, p3, p4, p5, p6,
+		d_data);
+
+
+	/** 8. 展示回波 **************************************************************************************************/
+	{
+		DeviceToHost(h_data.phdata, d_data.phdata, sizeof(cuComplex) * d_data.Np * d_data.Nfft);
+		cuComplex* h_ifftphdata = (cuComplex*)mallocCUDAHost(sizeof(cuComplex) * d_data.Np * d_data.Nfft);
+		cuComplex* d_ifftphdata = (cuComplex*)mallocCUDADevice(sizeof(cuComplex) * d_data.Np * d_data.Nfft);
+		CUDAIFFT(d_data.phdata, d_ifftphdata, d_data.Np, d_data.Nfft, d_data.Nfft);
+		FFTShift1D(d_ifftphdata, d_data.Np, d_data.Nfft);
+		DeviceToHost(h_ifftphdata, d_ifftphdata, sizeof(cuComplex) * d_data.Np * d_data.Nfft);
+		std::shared_ptr<std::complex<double>> ifftdata(new std::complex<double>[d_data.Np * d_data.Nfft], delArrPtr);
+		{
+			for (long i = 0; i < d_data.Np; i++) {
+				for (long j = 0; j < d_data.Nfft; j++) {
+					ifftdata.get()[i * d_data.Nfft + j] =
+						std::complex<double>(
+							h_ifftphdata[i * d_data.Nfft + j].x,
+							h_ifftphdata[i * d_data.Nfft + j].y);
+				}
+			}
+		}
+		testOutComplexDoubleArr(QString("echo_ifft_single.bin"), ifftdata.get(), d_data.Np, d_data.Nfft);
+
+		FreeCUDADevice(d_ifftphdata);
+		FreeCUDAHost(h_ifftphdata);
+
+	}
+	/** 9. 成像 **************************************************************************************************/
+
+	// 计算maxWr（需要先计算deltaF）
+	float deltaF = h_data.deltaF; // 从输入参数获取
+	float maxWr = 299792458.0f / (2.0f * deltaF);
+	qDebug() << "maxWr :\t" << maxWr;
+	float* r_vec_host = (float*)mallocCUDAHost(sizeof(float) * h_data.Nfft);//  new float[data.Nfft];
+	const float step = maxWr / h_data.Nfft;
+	const float start = -1 * h_data.Nfft / 2.0f * step;
+	printf("nfft=%d\n", h_data.Nfft);
+
+	for (int i = 0; i < h_data.Nfft; ++i) {
+		r_vec_host[i] = start + i * step;
+	}
+
+	h_data.r_vec = r_vec_host;
+	d_data.r_vec = h_data.r_vec;
+
+	qDebug() << "r_vec [0]:\t" << h_data.r_vec[0];
+	qDebug() << "r_vec step:\t" << step;
+
+	bpBasic0CUDA_single(d_data, 0);
+
+	DeviceToHost(h_data.im_final, d_data.im_final, sizeof(cuComplex) * d_data.nx * d_data.ny);
+
+	{
+		DeviceToHost(h_data.phdata, d_data.phdata, sizeof(cuComplex) * d_data.Np * d_data.Nfft);
+		std::shared_ptr<std::complex<double>> ifftdata(new std::complex<double>[d_data.Np * d_data.Nfft], delArrPtr);
+		{
+			for (long i = 0; i < d_data.Np; i++) {
+				for (long j = 0; j < d_data.Nfft; j++) {
+					ifftdata.get()[i * d_data.Nfft + j] =
+						std::complex<double>(
+							h_data.phdata[i * d_data.Nfft + j].x,
+							h_data.phdata[i * d_data.Nfft + j].y);
+				}
+			}
+		}
+		testOutComplexDoubleArr(QString("echo_ifft_BPBasic_single.bin"), ifftdata.get(), d_data.Np, d_data.Nfft);
+
+
+
+		std::shared_ptr<std::complex<double>> im_finals(new std::complex<double>[d_data.nx * d_data.ny], delArrPtr);
+		{
+			for (long i = 0; i < d_data.ny; i++) {
+				for (long j = 0; j < d_data.nx; j++) {
+					im_finals.get()[i * d_data.nx + j] = std::complex<double>(
+						h_data.im_final[i * d_data.nx + j].x,
+						h_data.im_final[i * d_data.nx + j].y);
+				}
+			}
+		}
+		testOutComplexDoubleArr(QString("im_finals_single.bin"), im_finals.get(), d_data.ny, d_data.nx);
+	}
+}
+
+
+void testBpImage_single() {
+
+	GPUDATA_single h_data{};
+	/** 0. 轨道 **************************************************************************************************/
+	QString echoPath = "D:\\Programme\\vs2022\\RasterMergeTest\\LAMPCAE_SCANE-all-scane\\GF3_Simulation.xml";
+	std::shared_ptr<EchoL0Dataset> echoL0ds(new EchoL0Dataset);
+	echoL0ds->Open(echoPath);
+	qDebug() << u8"加载回拨文件：\t" << echoPath;
+
+	/** 1. 轨道 **************************************************************************************************/
+	qDebug() << u8"1.轨道文件读取中。。。";
+	QString inGPSPath = echoL0ds->getGPSPointFilePath();
+	long gpspoints = gdalImage(inGPSPath).height;
+	std::shared_ptr<SatelliteAntPos> antpos = SatelliteAntPosOperator::readAntPosFile(inGPSPath, gpspoints);
+	h_data.AntX = (float*)mallocCUDAHost(sizeof(float) * gpspoints);
+	h_data.AntY = (float*)mallocCUDAHost(sizeof(float) * gpspoints);
+	h_data.AntZ = (float*)mallocCUDAHost(sizeof(float) * gpspoints);
+
+	for (long i = 0; i < gpspoints; i = i + 1) {
+		h_data.AntX[i] = antpos.get()[i].Px;
+		h_data.AntY[i] = antpos.get()[i].Py;
+		h_data.AntZ[i] = antpos.get()[i].Pz;
+	}
+	//gpspoints = gpspoints / 2;
+	qDebug() << "GPS points Count:\t" << gpspoints;
+	qDebug() << "PRF 0:\t " << QString("%1,%2,%3").arg(h_data.AntX[0]).arg(h_data.AntY[0]).arg(h_data.AntZ[0]);
+	qDebug() << "PRF " << gpspoints - 1 << ":\t " << QString("%1,%2,%3")
+		.arg(h_data.AntX[gpspoints - 1])
+		.arg(h_data.AntY[gpspoints - 1])
+		.arg(h_data.AntZ[gpspoints - 1]);
+	/** 2. 成像网格 **************************************************************************************************/
+	qDebug() << u8"轨道文件读取结束\n2.成像网格读取。。。";
+	QString demxyzPath = u8"D:\\Programme\\vs2022\\RasterMergeTest\\simulationData\\demdataset\\demxyz.bin";
+	gdalImage demgridimg(demxyzPath);
+	long dem_rowCount = demgridimg.height;
+	long dem_ColCount = demgridimg.width;
+	std::shared_ptr<double> demX = readDataArr<double>(demgridimg, 0, 0, dem_rowCount, dem_ColCount, 1, GDALREADARRCOPYMETHOD::VARIABLEMETHOD);
+	std::shared_ptr<double> demY = readDataArr<double>(demgridimg, 0, 0, dem_rowCount, dem_ColCount, 2, GDALREADARRCOPYMETHOD::VARIABLEMETHOD);
+	std::shared_ptr<double> demZ = readDataArr<double>(demgridimg, 0, 0, dem_rowCount, dem_ColCount, 3, GDALREADARRCOPYMETHOD::VARIABLEMETHOD);
+
+	h_data.x_mat = (float*)mallocCUDAHost(sizeof(float) * dem_rowCount * dem_ColCount);
+	h_data.y_mat = (float*)mallocCUDAHost(sizeof(float) * dem_rowCount * dem_ColCount);
+	h_data.z_mat = (float*)mallocCUDAHost(sizeof(float) * dem_rowCount * dem_ColCount);
+
+	for (long i = 0; i < dem_rowCount; i++) {
+		for (long j = 0; j < dem_ColCount; j++) {
+			h_data.x_mat[i * dem_ColCount + j] = demX.get()[i * dem_ColCount + j];
+			h_data.y_mat[i * dem_ColCount + j] = demY.get()[i * dem_ColCount + j];
+			h_data.z_mat[i * dem_ColCount + j] = demZ.get()[i * dem_ColCount + j];
+		}
+	}
+	qDebug() << "dem row Count:\t" << dem_rowCount;
+	qDebug() << "dem col Count:\t" << dem_ColCount;
+
+	qDebug() << u8"成像网格读取结束";
+	/** 3. 频率网格 **************************************************************************************************/
+	qDebug() << u8"3.处理频率";
+	float centerFreq = echoL0ds->getCenterFreq();
+	float bandwidth = echoL0ds->getBandwidth();
+	size_t freqpoints = echoL0ds->getPlusePoints();
+	std::shared_ptr<float> freqlist(getFreqPoints_mallocHost_single(centerFreq - bandwidth / 2, centerFreq + bandwidth / 2, freqpoints), FreeCUDAHost);
+	std::shared_ptr<float> minF(new float[gpspoints], delArrPtr);
+	for (long i = 0; i < gpspoints; i++) {
+		minF.get()[i] = freqlist.get()[0];
+	}
+	h_data.deltaF = bandwidth / (freqpoints - 1);
+	qDebug() << "start Freq:\t" << centerFreq - bandwidth / 2;
+	qDebug() << "end   Freq:\t" << centerFreq + bandwidth / 2;
+	qDebug() << "freq points:\t" << freqpoints;
+	qDebug() << "delta freq:\t" << freqlist.get()[1] - freqlist.get()[0];
+	qDebug() << u8"频率结束";
+	/** 4. 初始化回波 **************************************************************************************************/
+	qDebug() << u8"4.初始化回波";
+	std::shared_ptr<std::complex<double>> echoData = echoL0ds->getEchoArr();
+	size_t echosize = sizeof(cuComplex) * gpspoints * freqpoints;
+	qDebug() << "echo data size (byte) :\t" << echosize;
+	h_data.phdata = (cuComplex*)mallocCUDAHost(echosize);
+	shared_complexPtrToHostCuComplex(echoData.get(), h_data.phdata, gpspoints * freqpoints);
+
+	qDebug() << "Azimuth Points:\t" << gpspoints;
+	qDebug() << "Range Points:\t" << freqpoints;
+	qDebug() << u8"初始化回波结束";
+	/** 5. 初始化图像 **************************************************************************************************/
+	qDebug() << u8"5.初始化图像";
+	h_data.im_final = (cuComplex*)mallocCUDAHost(sizeof(cuComplex) * dem_rowCount * dem_ColCount);
+	qDebug() << "Azimuth Points:\t" << gpspoints;
+	qDebug() << "Range Points:\t" << freqpoints;
+	qDebug() << u8"初始化图像结束";
+
+	/** 6. 模型参数初始化 **************************************************************************************************/
+	qDebug() << u8"6.模型参数初始化";
+
+	h_data.nx = dem_ColCount;
+	h_data.ny = dem_rowCount;
+
+	h_data.Np = gpspoints;
+	h_data.Freq = freqlist.get();
+	h_data.minF = minF.get();
+	h_data.Nfft = freqpoints;
+	h_data.K = h_data.Nfft;
+
+	h_data.R0 = echoL0ds->getRefPhaseRange();
+	qDebug() << u8"模型参数结束";
+
+
+	/** 7. 目标 **************************************************************************************************/
+	double Tx = -2029086.618142, Ty = 4139594.934504, Tz = 4392846.782027;
+	double Tslx = -2029086.618142, Tsly = 4139594.934504, Tslz = 4392846.782027;
+	double p1 = 33.3, p2 = 0, p3 = 0, p4 = 0, p5 = 0, p6 = 0;
+
+	/** 7. 构建回波 **************************************************************************************************/
+	GPUDATA_single d_data;
+	initGPUData_single(h_data, d_data);
+
+	/** 8. 展示回波 **************************************************************************************************/
+	{
+		DeviceToHost(h_data.phdata, d_data.phdata, sizeof(cuComplex) * d_data.Np * d_data.Nfft);
+		cuComplex* h_ifftphdata = (cuComplex*)mallocCUDAHost(sizeof(cuComplex) * d_data.Np * d_data.Nfft);
+		cuComplex* d_ifftphdata = (cuComplex*)mallocCUDADevice(sizeof(cuComplex) * d_data.Np * d_data.Nfft);
+		CUDAIFFT(d_data.phdata, d_ifftphdata, d_data.Np, d_data.Nfft, d_data.Nfft);
+		FFTShift1D(d_ifftphdata, d_data.Np, d_data.Nfft);
+		DeviceToHost(h_ifftphdata, d_ifftphdata, sizeof(cuComplex) * d_data.Np * d_data.Nfft);
+		std::shared_ptr<std::complex<double>> ifftdata(new std::complex<double>[d_data.Np * d_data.Nfft], delArrPtr);
+		{
+			for (long i = 0; i < d_data.Np; i++) {
+				for (long j = 0; j < d_data.Nfft; j++) {
+					ifftdata.get()[i * d_data.Nfft + j] =
+						std::complex<double>(
+							h_ifftphdata[i * d_data.Nfft + j].x,
+							h_ifftphdata[i * d_data.Nfft + j].y);
+				}
+			}
+		}
+		testOutComplexDoubleArr(QString("echo_ifft.bin"), ifftdata.get(), d_data.Np, d_data.Nfft);
+
+		FreeCUDADevice(d_ifftphdata);
+		FreeCUDAHost(h_ifftphdata);
+
+	}
+	/** 9. 成像 **************************************************************************************************/
+
+	// 计算maxWr（需要先计算deltaF）
+	double deltaF = h_data.deltaF; // 从输入参数获取
+	double maxWr = 299792458.0f / (2.0f * deltaF);
+	qDebug() << "maxWr :\t" << maxWr;
+	float* r_vec_host = (float*)mallocCUDAHost(sizeof(float) * h_data.Nfft);//  new double[data.Nfft];
+	const double step = maxWr / h_data.Nfft;
+	const double start = -1 * h_data.Nfft / 2.0f * step;
+	printf("nfft=%d\n", h_data.Nfft);
+
+	for (int i = 0; i < h_data.Nfft; ++i) {
+		r_vec_host[i] = start + i * step;
+	}
+
+	h_data.r_vec = r_vec_host;
+	d_data.r_vec = h_data.r_vec;
+
+	qDebug() << "r_vec [0]:\t" << h_data.r_vec[0];
+	qDebug() << "r_vec step:\t" << step;
+
+	bpBasic0CUDA_single(d_data, 0);
+
+	DeviceToHost(h_data.im_final, d_data.im_final, sizeof(cuComplex) * d_data.nx * d_data.ny);
+	{
+		DeviceToHost(h_data.phdata, d_data.phdata, sizeof(cuComplex) * d_data.Np * d_data.Nfft);
+		std::shared_ptr<std::complex<double>> ifftdata(new std::complex<double>[d_data.Np * d_data.Nfft], delArrPtr);
+		{
+			for (long i = 0; i < d_data.Np; i++) {
+				for (long j = 0; j < d_data.Nfft; j++) {
+					ifftdata.get()[i * d_data.Nfft + j] =
+						std::complex<double>(
+							h_data.phdata[i * d_data.Nfft + j].x,
+							h_data.phdata[i * d_data.Nfft + j].y);
+				}
+			}
+		}
+		testOutComplexDoubleArr(QString("echo_ifft_BPBasic.bin"), ifftdata.get(), d_data.Np, d_data.Nfft);
+
+		std::shared_ptr<std::complex<double>> im_finals(new std::complex<double>[d_data.nx * d_data.ny], delArrPtr);
+		{
+			for (long i = 0; i < d_data.ny; i++) {
+				for (long j = 0; j < d_data.nx; j++) {
+					im_finals.get()[i * d_data.nx + j] = std::complex<double>(
+						h_data.im_final[i * d_data.nx + j].x,
+						h_data.im_final[i * d_data.nx + j].y);
+				}
+			}
+		}
+		testOutComplexDoubleArr(QString("im_finals.bin"), im_finals.get(), d_data.ny, d_data.nx);
+	}
+
+
+
+}
+
+
+
+
+
+int main(int argc, char* argv[]) {
+
+	QApplication a(argc, argv);
+	
+	//testSimualtionEchoPoint();
+	//testSimualtionEchoPoint_single();
+	void testBpImage_single();
+	return 0;
+}