float变换为double,float 精度不够

BaseTool/BaseConstVariable.h

@@ -14,8 +14,8 @@
 
 #define __CUDANVCC___ // 定义CUDA函数
 
-//#define __PRFDEBUG__
+#define __PRFDEBUG__
-//#define __PRFDEBUG_PRFINF__
+#define __PRFDEBUG_PRFINF__
 //#define __ECHOTIMEDEBUG__
 
 #define __TBPIMAGEDEBUG__

BaseTool/EchoDataFormat.cpp

@@ -52,9 +52,9 @@ std::shared_ptr<PluseData> CreatePluseDataArr(long pluseCount)
     return std::shared_ptr<PluseData>(new PluseData[pluseCount],delArrPtr);
 }
 
-std::shared_ptr<std::complex<float>> CreateEchoData(long plusePoints)
+std::shared_ptr<std::complex<double>> CreateEchoData(long plusePoints)
 {
-    return std::shared_ptr<std::complex<float>>(new std::complex<float>[plusePoints],delArrPtr);
+    return std::shared_ptr<std::complex<double>>(new std::complex<double>[plusePoints],delArrPtr);
 }
 
 EchoL0Dataset::EchoL0Dataset()
@@ -164,7 +164,7 @@ ErrorCode EchoL0Dataset::OpenOrNew(QString folder, QString filename, long PluseC
         CPLSetConfigOption("GDAL_FILENAME_IS_UTF8", "YES");
  
         GDALDriver* poDriver = GetGDALDriverManager()->GetDriverByName("ENVI");
-        GDALDataset* poDstDS = (poDriver->Create(this->echoDataFilePath.toUtf8().constData(), PlusePoints, PluseCount, 1, GDT_CFloat32, NULL));
+        GDALDataset* poDstDS = (poDriver->Create(this->echoDataFilePath.toUtf8().constData(), PlusePoints, PluseCount, 1, GDT_CFloat64, NULL));
         GDALFlushCache((GDALDatasetH)poDstDS);
         GDALClose(poDstDS);
         //poDstDS.reset();
@@ -251,14 +251,14 @@ QString EchoL0Dataset::getEchoDataFilename()
     return GPSPointFilePath;
 }
 
-void EchoL0Dataset::initEchoArr(std::complex<float> init0)
+void EchoL0Dataset::initEchoArr(std::complex<double> init0)
 {
     long blockline = Memory1MB * 2000 / 8 / 2 / this->PlusePoints;
     
     long start = 0;
     for (start = 0; start < this->PluseCount; start = start + blockline) {
         long templine = start + blockline < this->PluseCount ? blockline : this->PluseCount - start;
-        std::shared_ptr<std::complex<float>> echotemp = this->getEchoArr(start, templine);
+        std::shared_ptr<std::complex<double>> echotemp = this->getEchoArr(start, templine);
         for (long i = 0; i < templine; i++) {
             for (long j = 0; j < this->PlusePoints; j++) {
                 echotemp.get()[i * this->PlusePoints + j] = init0;
@@ -533,7 +533,7 @@ std::shared_ptr<double> EchoL0Dataset::getAntPos()
     return temp;
 }
 
-std::shared_ptr<std::complex<float>> EchoL0Dataset::getEchoArr(long startPRF, long PRFLen)
+std::shared_ptr<std::complex<double>> EchoL0Dataset::getEchoArr(long startPRF, long PRFLen)
 {
     if (!(startPRF < this->PluseCount)) {
         qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::ECHO_L0DATA_PRFIDXOUTRANGE))<<startPRF<<"  "<<this->PluseCount;
@@ -560,14 +560,14 @@ std::shared_ptr<std::complex<float>> EchoL0Dataset::getEchoArr(long startPRF, lo
     long width = rasterDataset->GetRasterXSize();
     long height = rasterDataset->GetRasterYSize();
     long band_num = rasterDataset->GetRasterCount();
-    std::shared_ptr<std::complex<float>>  temp = nullptr;
+    std::shared_ptr<std::complex<double>>  temp = nullptr;
 	if (height != this->PluseCount || width != this->PlusePoints) {
 		qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::ECHO_L0DATA_ECHOFILEFORMATERROR));
 	}
 	else {
-		if (gdal_datatype == GDT_CFloat32) {
+		if (gdal_datatype == GDT_CFloat64) {
-            temp= std::shared_ptr<std::complex<float>>(new std::complex<float>[PRFLen * width ], delArrPtr);
+            temp= std::shared_ptr<std::complex<double>>(new std::complex<double>[PRFLen * width ], delArrPtr);
-            poBand->RasterIO(GF_Read, 0, startPRF, width, PRFLen, temp.get(), width, PRFLen, GDT_CFloat32, 0, 0);
+            poBand->RasterIO(GF_Read, 0, startPRF, width, PRFLen, temp.get(), width, PRFLen, GDT_CFloat64, 0, 0);
             GDALFlushCache((GDALDatasetH)rasterDataset.get());
 		}
 		else {
@@ -581,7 +581,7 @@ std::shared_ptr<std::complex<float>> EchoL0Dataset::getEchoArr(long startPRF, lo
     return temp;
 }
 
-std::shared_ptr<std::complex<float>> EchoL0Dataset::getEchoArr()
+std::shared_ptr<std::complex<double>> EchoL0Dataset::getEchoArr()
 {
     return  this->getEchoArr(0,this->PluseCount);
 }
@@ -616,7 +616,7 @@ ErrorCode EchoL0Dataset::saveAntPos(std::shared_ptr<double> ptr)
     return ErrorCode::SUCCESS;
 }
 
-ErrorCode EchoL0Dataset::saveEchoArr(std::shared_ptr<std::complex<float>> echoPtr, long startPRF, long PRFLen)
+ErrorCode EchoL0Dataset::saveEchoArr(std::shared_ptr<std::complex<double>> echoPtr, long startPRF, long PRFLen)
 {
     if (!(startPRF < this->PluseCount)) {
         qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::ECHO_L0DATA_PRFIDXOUTRANGE));
@@ -655,8 +655,8 @@ ErrorCode EchoL0Dataset::saveEchoArr(std::shared_ptr<std::complex<float>> echoPt
         return ErrorCode::ECHO_L0DATA_ECHOFILEFORMATERROR;
     }
     else {
-        if (gdal_datatype == GDT_CFloat32) {
+        if (gdal_datatype == GDT_CFloat64) {
-            poBand->RasterIO(GF_Write, 0, startPRF, width, PRFLen,  echoPtr.get(), width, PRFLen, GDT_CFloat32, 0, 0);
+            poBand->RasterIO(GF_Write, 0, startPRF, width, PRFLen,  echoPtr.get(), width, PRFLen, GDT_CFloat64, 0, 0);
             GDALFlushCache((GDALDatasetH)rasterDataset);
         }
         else {

BaseTool/EchoDataFormat.h

@@ -72,7 +72,7 @@ struct PluseData {
 long getPluseDataSize(PluseData& pluseData);
 ErrorCode getPluseDataFromBuffer(char* buffer, PluseData& data);
 std::shared_ptr<PluseData> CreatePluseDataArr(long pluseCount);
-std::shared_ptr<std::complex<float>> CreateEchoData(long plusePoints);
+std::shared_ptr<std::complex<double>> CreateEchoData(long plusePoints);
 
 
 
@@ -108,7 +108,7 @@ public:
 	QString getGPSPointFilename();
 	QString getEchoDataFilename();
 
-	void initEchoArr(std::complex<float> init0);
+	void initEchoArr(std::complex<double> init0);
 
 
 private: // 产品名称设置
@@ -184,11 +184,11 @@ public: 	//
 public:
 	// 读取文件
 	std::shared_ptr<double> getAntPos();
-	std::shared_ptr<std::complex<float>> getEchoArr(long startPRF, long PRFLen);
+	std::shared_ptr<std::complex<double>> getEchoArr(long startPRF, long PRFLen);
-	std::shared_ptr<std::complex<float>> getEchoArr();
+	std::shared_ptr<std::complex<double>> getEchoArr();
 	//保存文件
 	ErrorCode saveAntPos(std::shared_ptr<double> ptr); // 注意这个方法很危险，请写入前检查数据是否正确
-	ErrorCode saveEchoArr(std::shared_ptr<std::complex<float>> echoPtr, long startPRF, long PRFLen);
+	ErrorCode saveEchoArr(std::shared_ptr<std::complex<double>> echoPtr, long startPRF, long PRFLen);
 
 };
 

BaseTool/ImageOperatorBase.cpp

@@ -290,7 +290,7 @@ ReadComplexMatrixData(int start_line, int width, int line_num,
 		data_mat.col(1) = (complex_short_mat.imag().array().cast<double>()).array();
 		_flag			= true;
 	} else if(gdal_datatype == GDT_CFloat32) {
-		Eigen::MatrixX<std::complex<float>> complex_short_mat(line_num * width, 1);
+		Eigen::MatrixX<std::complex<double>> complex_short_mat(line_num * width, 1);
 		rasterDataset->GetRasterBand(1)->RasterIO(GF_Read, 0, start_line, width, line_num,
 												  complex_short_mat.data(), width, line_num,
 												  gdal_datatype, 0, 0); // real

BaseTool/ImageOperatorBase.h

@@ -429,7 +429,7 @@ std::shared_ptr<T> readDataArr(gdalImage& imgds, int start_row, int start_col, i
 		delete[] temp;
 	}
 	//else if (gdal_datatype == GDT_CFloat32) {
-	//	if (std::is_same<T, std::complex<double>>::value || std::is_same<T, std::complex<float>>::value) {
+	//	if (std::is_same<T, std::complex<double>>::value || std::is_same<T, std::complex<double>>::value) {
 	//		float* temp = new float[rows_count * cols_count * 2];
 	//		demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, rows_count, gdal_datatype, 0, 0);
 
@@ -451,7 +451,7 @@ std::shared_ptr<T> readDataArr(gdalImage& imgds, int start_row, int start_col, i
 	//	}
 	//}
 	//else if (gdal_datatype == GDT_CFloat64 ) {
-	//	if (std::is_same<T, std::complex<double>>::value || std::is_same<T, std::complex<float>>::value) {
+	//	if (std::is_same<T, std::complex<double>>::value || std::is_same<T, std::complex<double>>::value) {
 	//		double* temp = new double[rows_count * cols_count * 2];
 	//		demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, rows_count, gdal_datatype, 0, 0);
 
@@ -502,7 +502,7 @@ std::shared_ptr<T> readDataArrComplex(gdalImageComplex& imgds, int start_row, in
 	//Eigen::MatrixXd datamatrix(rows_count, cols_count);
 
 	if (gdal_datatype == GDT_CFloat32) {
-		if (std::is_same<T, std::complex<double>>::value || std::is_same<T, std::complex<float>>::value) {
+		if (std::is_same<T, std::complex<double>>::value || std::is_same<T, std::complex<double>>::value) {
 			float* temp = new float[rows_count * cols_count * 2];
 			demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, rows_count, gdal_datatype, 0, 0);
 
@@ -524,7 +524,7 @@ std::shared_ptr<T> readDataArrComplex(gdalImageComplex& imgds, int start_row, in
 		}
 	}
 	else if (gdal_datatype == GDT_CFloat64) {
-		if (std::is_same<T, std::complex<double>>::value || std::is_same<T, std::complex<float>>::value) {
+		if (std::is_same<T, std::complex<double>>::value || std::is_same<T, std::complex<double>>::value) {
 			double* temp = new double[rows_count * cols_count * 2];
 			demBand->RasterIO(GF_Read, start_col, start_row, cols_count, rows_count, temp, cols_count, rows_count, gdal_datatype, 0, 0);
 

BaseTool/SARSimulationImageL1.cpp

@@ -605,7 +605,7 @@ ErrorCode SARSimulationImageL1Dataset::saveAntPos(std::shared_ptr<double> ptr)
     return ErrorCode::SUCCESS;
 }
 
-std::shared_ptr<std::complex<float>> SARSimulationImageL1Dataset::getImageRaster()
+std::shared_ptr<std::complex<double>> SARSimulationImageL1Dataset::getImageRaster()
 {
     if (this->Rasterlevel != RasterLevel::RasterSLC) {
         return nullptr;
@@ -614,7 +614,7 @@ std::shared_ptr<std::complex<float>> SARSimulationImageL1Dataset::getImageRaster
     return  this->getImageRaster(0, this->rowCount);
 }
 
-ErrorCode SARSimulationImageL1Dataset::saveImageRaster(std::shared_ptr<std::complex<float>> echoPtr, long startPRF, long PRFLen)
+ErrorCode SARSimulationImageL1Dataset::saveImageRaster(std::shared_ptr<std::complex<double>> echoPtr, long startPRF, long PRFLen)
 {
     if (!(startPRF < this->rowCount)) {
         qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::ECHO_L0DATA_PRFIDXOUTRANGE));
@@ -663,7 +663,7 @@ ErrorCode SARSimulationImageL1Dataset::saveImageRaster(std::shared_ptr<std::comp
     return ErrorCode::SUCCESS;
 }
 
-std::shared_ptr<std::complex<float>> SARSimulationImageL1Dataset::getImageRaster(long startPRF, long PRFLen)
+std::shared_ptr<std::complex<double>> SARSimulationImageL1Dataset::getImageRaster(long startPRF, long PRFLen)
 {
     if (this->Rasterlevel != RasterLevel::RasterSLC) {
         return nullptr;
@@ -695,13 +695,13 @@ std::shared_ptr<std::complex<float>> SARSimulationImageL1Dataset::getImageRaster
     long width = rasterDataset->GetRasterXSize();
     long height = rasterDataset->GetRasterYSize();
     long band_num = rasterDataset->GetRasterCount();
-    std::shared_ptr<std::complex<float>>  temp = nullptr;
+    std::shared_ptr<std::complex<double>>  temp = nullptr;
     if (height != this->rowCount || width != this->colCount) {
         qDebug() << QString::fromStdString(errorCode2errInfo(ErrorCode::ECHO_L0DATA_ECHOFILEFORMATERROR));
     }
     else {
         if (gdal_datatype == GDT_CFloat32) {
-            temp = std::shared_ptr<std::complex<float>>(new std::complex<float>[PRFLen * width], delArrPtr);
+            temp = std::shared_ptr<std::complex<double>>(new std::complex<double>[PRFLen * width], delArrPtr);
             poBand->RasterIO(GF_Read, 0, startPRF, width, PRFLen, temp.get(), width, PRFLen, GDT_CFloat32, 0, 0);
             GDALFlushCache((GDALDatasetH)rasterDataset.get());
         }
@@ -723,7 +723,7 @@ Eigen::MatrixXcd SARSimulationImageL1Dataset::getImageRasterMatrix()
     }
 
 
-    std::shared_ptr<std::complex<float>> data = this->getImageRaster();
+    std::shared_ptr<std::complex<double>> data = this->getImageRaster();
     Eigen::MatrixXcd dataimg = Eigen::MatrixXcd::Zero(this->rowCount, this->colCount);
     for (long i = 0; i < this->rowCount; i++) {
         for (long j = 0; j < this->colCount; j++) {
@@ -735,7 +735,7 @@ Eigen::MatrixXcd SARSimulationImageL1Dataset::getImageRasterMatrix()
 
 ErrorCode SARSimulationImageL1Dataset::saveImageRaster(Eigen::MatrixXcd& dataimg, long startPRF)
 {
-    std::shared_ptr<std::complex<float>> data(new std::complex<float>[dataimg.rows()* dataimg.cols()]);
+    std::shared_ptr<std::complex<double>> data(new std::complex<double>[dataimg.rows()* dataimg.cols()]);
     long dataimgrows = dataimg.rows();
     long dataimgcols = dataimg.cols();
     for (long i = 0; i < dataimgrows; i++) {

BaseTool/SARSimulationImageL1.h

@@ -57,9 +57,9 @@ public:
 	std::shared_ptr<double> getAntPos();
 	ErrorCode saveAntPos(std::shared_ptr<double> ptr); // 注意这个方法很危险，请写入前检查数据是否正确
 	
-	std::shared_ptr<std::complex<float>> getImageRaster();
+	std::shared_ptr<std::complex<double>> getImageRaster();
-	ErrorCode saveImageRaster(std::shared_ptr<std::complex<float>> echoPtr, long startPRF, long PRFLen);
+	ErrorCode saveImageRaster(std::shared_ptr<std::complex<double>> echoPtr, long startPRF, long PRFLen);
-	std::shared_ptr<std::complex<float>> getImageRaster(long startPRF, long PRFLen);
+	std::shared_ptr<std::complex<double>> getImageRaster(long startPRF, long PRFLen);
 	
 	Eigen::MatrixXcd getImageRasterMatrix();
 	ErrorCode saveImageRaster(Eigen::MatrixXcd& data, long startPRF);

GPURFPCKernel.cu

@@ -0,0 +1,599 @@
+#include <iostream>
+#include <memory>
+#include <cmath>
+#include <complex>
+#include <device_launch_parameters.h>
+#include <cuda_runtime.h>
+#include <cublas_v2.h>
+#include <cuComplex.h>
+
+#include "BaseConstVariable.h"
+#include "GPURFPC.cuh"
+
+
+#ifdef __CUDANVCC___
+
+
+
+
+
+__device__ float GPU_getSigma0dB(CUDASigmaParam param, float theta) {//线性值
+	float sigma = param.p1 + param.p2 * exp(-param.p3 * theta) + param.p4 * cos(param.p5 * theta + param.p6);
+	return  sigma;
+}
+
+
+__device__ CUDAVectorEllipsoidal GPU_SatelliteAntDirectNormal(
+	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 AntXaxisX = antXaxisX;
+	float AntXaxisY = antXaxisY;
+	float AntXaxisZ = antXaxisZ;
+	float AntYaxisX = antYaxisX;
+	float AntYaxisY = antYaxisY;
+	float AntYaxisZ = antYaxisZ;
+	float AntZaxisX = antZaxisX;
+	float AntZaxisY = antZaxisY;
+	float AntZaxisZ = antZaxisZ;
+
+	// 归一化
+	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 ThetaAnt = acosf(Zant / Norm); // theta 与 Z轴的夹角
+	float PhiAnt = 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;
+}
+
+__device__ float GPU_BillerInterpAntPattern(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;
+	}
+}
+
+__device__ cuComplex  GPU_calculationEcho(float sigma0, float TransAnt, float ReciveAnt,
+	float localangle, float R, float slopeangle, float Pt, float lamda) {
+	float amp = Pt * TransAnt * ReciveAnt;
+	amp = amp * sigma0;
+	amp = amp / (powf(4 * LAMP_CUDA_PI, 2) * powf(R, 4)); // 反射强度
+	float phi = (-4 * LAMP_CUDA_PI / lamda) * R;
+	cuComplex echophi = make_cuComplex(0, phi);
+	cuComplex echophiexp = cuCexpf(echophi);
+	cuComplex echo = make_cuComplex(echophiexp.x * amp, echophiexp.y * amp);
+	return echo;
+}
+
+
+__global__ void CUDA_SatelliteAntDirectNormal(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,
+	float* thetaAnt, float* phiAnt
+	, long len) {
+	long  idx = blockIdx.x * blockDim.x + threadIdx.x;
+	if (idx < len) {
+		float Xst = -1 * RstX[idx]; // 卫星 -->  地面
+		float Yst = -1 * RstY[idx];
+		float Zst = -1 * RstZ[idx];
+		float AntXaxisX = antXaxisX;
+		float AntXaxisY = antXaxisY;
+		float AntXaxisZ = antXaxisZ;
+		float AntYaxisX = antYaxisX;
+		float AntYaxisY = antYaxisY;
+		float AntYaxisZ = antYaxisZ;
+		float AntZaxisX = antZaxisX;
+		float AntZaxisY = antZaxisY;
+		float AntZaxisZ = antZaxisZ;
+
+		// 归一化
+		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 ThetaAnt = acosf(Zant / Norm); // theta 与 Z轴的夹角
+		float PhiAnt = 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);
+		}
+
+		//if (abs(ThetaAnt - 0) < PRECISIONTOLERANCE) {
+		//	PhiAnt = 0;
+		//}
+		//else {}
+
+
+		thetaAnt[idx] = ThetaAnt * r2d;
+		phiAnt[idx] = PhiAnt * r2d;
+		//printf("Rst=[%f,%f,%f];AntXaxis = [%f, %f, %f];AntYaxis=[%f,%f,%f];AntZaxis=[%f,%f,%f];phiAnt=%f;thetaAnt=%f;\n", Xst, Yst, Zst
+		//	, AntXaxisX, AntXaxisY, AntXaxisZ
+		//	, AntYaxisX, AntYaxisY, AntYaxisZ
+		//	, AntZaxisX, AntZaxisY, AntZaxisZ
+		//	, phiAnt[idx]
+		//	, thetaAnt[idx]
+		//);
+	}
+}
+
+__global__ void CUDA_BillerInterpAntPattern(float* antpattern,
+	float starttheta, float startphi, float dtheta, float dphi,
+	long thetapoints, long phipoints,
+	float* searththeta, float* searchphi, float* searchantpattern,
+	long len) {
+	long  idx = blockIdx.x * blockDim.x + threadIdx.x;
+	if (idx < len) {
+		float stheta = searththeta[idx];
+		float sphi = searchphi[idx];
+		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)
+		{
+			searchantpattern[idx] = 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);
+			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));
+			searchantpattern[idx] = GainValue;
+		}
+	}
+}
+
+__global__ void CUDA_calculationEcho(float* sigma0, float* TransAnt, float* ReciveAnt,
+	float* localangle, float* R, float* slopeangle,
+	float nearRange, float Fs, float Pt, float lamda, long FreqIDmax,
+	cuComplex* echoArr, long* FreqID,
+	long len) {
+	long  idx = blockIdx.x * blockDim.x + threadIdx.x;
+	if (idx < len) {
+		float r = R[idx];
+		float amp = Pt * TransAnt[idx] * ReciveAnt[idx];
+		amp = amp * sigma0[idx];
+		amp = amp / (powf(4 * LAMP_CUDA_PI, 2) * powf(r, 4)); // 反射强度
+
+		// 处理相位
+		float phi = (-4 * LAMP_CUDA_PI / lamda) * r;
+		cuComplex echophi = make_cuComplex(0, phi);
+		cuComplex echophiexp = cuCexpf(echophi);
+
+		float timeR = 2 * (r - nearRange) / LIGHTSPEED * Fs;
+		long timeID = floorf(timeR);
+		//if (timeID < 0 || timeID >= FreqIDmax) {
+		//	timeID = 0;
+		//	amp = 0;
+		//}
+
+		cuComplex echo = make_cuComplex(echophiexp.x, echophiexp.y);
+		echoArr[idx] = echo;
+		FreqID[idx] = timeID;
+	}
+}
+
+
+__global__ void CUDA_AntPatternInterpGain(float* anttheta, float* antphi, float* gain,
+	float* antpattern, float starttheta, float startphi, float dtheta, float dphi, int thetapoints, int phipoints, long len) {
+	int  idx = blockIdx.x * blockDim.x + threadIdx.x;
+
+	if (idx < len) {
+
+		float temptheta = anttheta[idx];
+		float tempphi = antphi[idx];
+		float antPatternGain = GPU_BillerInterpAntPattern(antpattern,
+			starttheta, startphi, dtheta, dphi, thetapoints, phipoints,
+			temptheta, tempphi);
+		gain[idx] = antPatternGain;
+	}
+}
+
+
+__global__ void CUDA_InterpSigma(
+	long* demcls, float* sigmaAmp, float* localanglearr, long len,
+	CUDASigmaParam* sigma0Paramslist, long sigmaparamslistlen) {
+	long  idx = blockIdx.x * blockDim.x + threadIdx.x;
+	if (idx < len) {
+		long clsid = demcls[idx];
+		float localangle = localanglearr[idx];
+		CUDASigmaParam tempsigma = sigma0Paramslist[clsid];
+		if (localangle < 0 || localangle >= LAMP_CUDA_PI / 2) {
+			sigmaAmp[idx] = 0;
+		}
+		else {}
+
+		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
+			) {
+			sigmaAmp[idx] = 0;
+		}
+		else {
+			float sigma = GPU_getSigma0dB(tempsigma, localangle);
+			sigma = powf(10.0, sigma / 10.0);// 后向散射系数
+			//printf("cls:%d;localangle=%f;sigma0=%f;\n", clsid, localangle, sigma);
+			sigmaAmp[idx] = sigma;
+		}
+	}
+}
+
+
+
+
+
+__global__ void CUDAKernel_RFPC_Caluation_R_Gain(
+	float antX, float antY, float antZ, // 天线的坐标
+	float* targetX, float* targetY, float* targetZ, long len, // 地面坐标
+	long* demCls,
+	float* demSlopeX, float* demSlopeY, float* demSlopeZ, // 地表坡度矢量
+	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,// 天线的指向
+	float Pt,// 发射能量
+	double 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 NearR, float FarR, // 距离范围
+	CUDASigmaParam* sigma0Paramslist, long sigmaparamslistlen,// 插值图
+	float* factorj, long freqnum,
+	double* outR,  // 输出距离
+	float* outAmp // 输出增益
+) {
+	long  idx = blockIdx.x * blockDim.x + threadIdx.x;
+	if (idx < len) {
+		double tx = targetX[idx];
+		double ty = targetY[idx];
+		double tz = targetZ[idx];
+		double RstX = antX - tx; // 计算坐标矢量
+		double RstY = antY - ty;
+		double RstZ = antZ - tz;
+
+		float slopeX = demSlopeX[idx];
+		float slopeY = demSlopeY[idx];
+		float slopeZ = demSlopeZ[idx];
+
+		double RstR2 = RstX * RstX + RstY * RstY + RstZ * RstZ;
+		double RstR = sqrt(RstR2); // 矢量距离
+
+		//printf("antX=%f;antY=%f;antZ=%f;targetX=%f;targetY=%f;targetZ=%f;RstR=%.6f;diffR=%.6f;\n",antX,antY,antZ,targetX,targetY,targetZ,RstR, RstR - 9.010858499003178e+05);
+
+		if (RstR<NearR || RstR>FarR) {
+			outAmp[idx] = 0;
+			outR[idx] = 0;
+		}
+		else {
+			// 求解坡度
+			float slopR = sqrtf(slopeX * slopeX + slopeY * slopeY + slopeZ * slopeZ); //  
+			float dotAB = RstX * slopeX + RstY * slopeY + RstZ * slopeZ;
+			float localangle = acosf(dotAB / (RstR * slopR)); // 局地入射角
+			float ampGain = 0;
+			// 求解天线方向图指向
+			CUDAVectorEllipsoidal antVector = GPU_SatelliteAntDirectNormal(
+				RstX, RstY, RstZ,
+				antXaxisX, antXaxisY, antXaxisZ,
+				antYaxisX, antYaxisY, antYaxisZ,
+				antZaxisX, antZaxisY, antZaxisZ,
+				antDirectX, antDirectY, antDirectZ
+			);
+			if (antVector.Rho > 0) {
+				// 发射方向图
+				float temptheta = antVector.theta * r2d;
+				float tempphi = antVector.phi * r2d;
+				float TansantPatternGain =
+					GPU_BillerInterpAntPattern(
+						TransAntpattern,
+						Transtarttheta, Transstartphi, Transdtheta, Transdphi, Transthetapoints, Transphipoints,
+						temptheta, tempphi);
+
+				// 接收方向图
+				float antPatternGain = GPU_BillerInterpAntPattern(
+					ReceiveAntpattern,
+					Receivestarttheta, Receivestartphi, Receivedtheta, Receivedphi, Receivethetapoints, Receivephipoints,
+					temptheta, tempphi);
+
+				// 计算
+				float sigma0 = 0;
+				{
+					long clsid = demCls[idx];
+					//printf("clsid=%d\n", clsid);
+					CUDASigmaParam tempsigma = sigma0Paramslist[clsid];
+					if (localangle < 0 || localangle >= LAMP_CUDA_PI / 2) {
+						sigma0 = 0;
+					}
+					else {}
+
+					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(tempsigma, localangle);
+						sigma0 = powf(10.0, sigma / 10.0);// 后向散射系数
+					}
+				}
+
+				ampGain = TansantPatternGain * antPatternGain;
+				ampGain = ampGain / (powf(4 * LAMP_CUDA_PI, 2) * powf(RstR, 4)); // 反射强度
+				outAmp[idx] = ampGain * Pt * sigma0;
+				outR[idx] = RstR - refPhaseRange;
+			}
+			else {
+				outAmp[idx] = 0;
+				outR[idx] = 0;
+			}
+		}
+	}
+}
+
+__global__ void CUDAKernel_PRF_CalFreqEcho(
+	double* Rarr, float* ampArr, long pixelcount,
+	float* factorj, long freqnum,
+	double dx, double nearR,
+	cuComplex* PRFEcho, long prfid) {
+	long  idx = blockIdx.x * blockDim.x + threadIdx.x;
+	if (idx < freqnum) {
+		float fatorj = factorj[idx];
+		float phi = 0;
+		float amptemp = 0;
+		cuComplex tempfreqEcho = PRFEcho[prfid * freqnum + idx];
+		for (long i = 0; i < pixelcount; i++) { // 区域积分
+			//phi = (R = R - (floor(R / lamda) - 1) * lamda)* fatorj; // 相位
+			float phi = Rarr[i] * factorj[idx]; // 相位
+			amptemp = ampArr[i];
+			//printf("amp=%f\n", amptemp);
+			// Eular; exp(ix)=cos(x)+isin(x)
+			tempfreqEcho.x = tempfreqEcho.x + amptemp * cos(phi); // 实部
+			tempfreqEcho.y = tempfreqEcho.y + amptemp * sin(phi); // 虚部
+			//printf("freqid=%d;fatorj=%.12f;d_R=%.10f;phi=%.10f;echo=complex(%.5f,%.5f)\n", idx,  fatorj, Rarr[i], phi, tempfreqEcho.x, tempfreqEcho.y);
+		}
+		PRFEcho[prfid * freqnum + idx] = tempfreqEcho;
+	}
+}
+
+
+
+__global__ void CUDAKernel_PRFSumEcho_Rows(
+	double* Rarr,float* ampArr,long Rows,long Cols,
+	long startRid,
+	float* factorj, long freqnum,
+	cuComplex* freqRowsbuffer, long tempRows
+){
+	long  idx = blockIdx.x * blockDim.x + threadIdx.x;
+	if (idx < Rows) { // 按行汇总
+		double R = 0;
+		double tempamp = 0;
+		float phi = 0;
+		long rid = idx + startRid;
+		float factor = 0;
+
+		for (long jj = 0; jj < freqnum; jj++) {
+			tempamp = ampArr[rid * Cols + jj];
+			R = Rarr[rid * Cols + jj];
+			for (long ii = 0; ii < freqnum; ii++) {
+				phi = R * factorj[ii];
+				freqRowsbuffer[idx * freqnum + ii].x = freqRowsbuffer[idx * freqnum + ii].x + tempamp * cos(phi); // 实部
+				freqRowsbuffer[idx * freqnum + ii].y = freqRowsbuffer[idx * freqnum + ii].y + tempamp * sin(phi); // 虚部
+			}
+			//freqRowsbuffer[idx * freqnum + ii] =  tempfreqEcho;
+		}
+	}
+}
+
+
+__global__ void  CUDAKernel_PRFSumEcho_Freq(
+	cuComplex* freqRowsbuffer, long tempRows,long freqnum,
+	cuComplex* PRFEcho, long prfid
+) {
+	long  idx = blockIdx.x * blockDim.x + threadIdx.x;
+	if (idx < freqnum) { // 按行汇总
+		cuComplex tempfreqEcho = freqRowsbuffer[prfid * freqnum + idx];
+		cuComplex temp = tempfreqEcho;
+		for (long ii = 0; ii < tempRows; ii++) { // 求和汇总
+			temp = freqRowsbuffer[ii * freqnum + idx];
+			tempfreqEcho.x = tempfreqEcho.x + temp.x;
+			tempfreqEcho.y = tempfreqEcho.y + temp.y;
+		}
+		freqRowsbuffer[prfid * freqnum + idx] = tempfreqEcho;
+	}
+}
+
+
+
+
+#endif

GPUTool/GPURFPC.cu

@@ -430,7 +430,7 @@ __global__ void CUDAKernel_RFPC_Caluation_R_Gain(
 	float* factorj, long freqnum,
 	double* outR,  // 输出距离
 	//float* outAmp // 输出增益
-	float* PRFEcho_real, float* PRFEcho_imag, long prfid
+	double* PRFEcho_real, double* PRFEcho_imag, long prfid
 ) {
 	long  idx = blockIdx.x * blockDim.x + threadIdx.x;
 	if (idx < len) {
@@ -514,10 +514,14 @@ __global__ void CUDAKernel_RFPC_Caluation_R_Gain(
 				double outR =  RstR- refPhaseRange;
 				for (long ii = 0; ii < freqnum; ii++) {
 					float phi= outR * factorj[ii]; // 相位
-					// Eular; exp(ix)=cos(x)+isin(x)
+					 // Eular; exp(ix)=cos(x)+isin(x)
-					float real =  outAmp * cos(phi); // 实部
+					double real =  outAmp * cos(phi); // ʵ²¿
-					float imag =  outAmp * sin(phi); // 虚部
+					double imag =  outAmp * sin(phi); // Ð鲿
+					atomicAdd(&PRFEcho_real[prfid * freqnum+ ii], real);
+					atomicAdd(&PRFEcho_imag[prfid * freqnum+ ii], imag);
 				}
+
+
 			}
 			else {
 			}
@@ -668,7 +672,7 @@ extern "C" void CUDARFPC_Caluation_R_Gain(
 	float* factorj, long freqnum,
 	double* outR,  // 输出距离
 	//float* outAmp // 输出增益
-	float* PRFEcho_real, float* PRFEcho_imag, long prfid
+	double* PRFEcho_real, double* PRFEcho_imag, long prfid
 )
 {
 

GPUTool/GPURFPC.cuh

@@ -18,6 +18,100 @@ extern "C" struct  CUDASigmaParam {
 	float p6;
 };
 
+
+extern __device__ float GPU_getSigma0dB(CUDASigmaParam param, float theta);
+
+extern __device__ CUDAVectorEllipsoidal GPU_SatelliteAntDirectNormal(
+	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(float* antpattern,
+	float starttheta, float startphi, float dtheta, float dphi,
+	long thetapoints, long phipoints,
+	float searththeta, float searchphi);
+
+
+extern __global__ void CUDA_AntPatternInterpGain(float* anttheta, float* antphi, float* gain,
+	float* antpattern, float starttheta, float startphi, float dtheta, float dphi, int thetapoints, int phipoints, long len);
+
+
+extern __global__ void CUDA_InterpSigma(
+	long* demcls, float* sigmaAmp, float* localanglearr, long len,
+	CUDASigmaParam* sigma0Paramslist, long sigmaparamslistlen);
+
+
+extern __global__ void CUDAKernel_RFPC_Caluation_R_Gain(
+	float antX, float antY, float antZ, // 天线的坐标
+	float* targetX, float* targetY, float* targetZ, long len, // 地面坐标
+	long* demCls,
+	float* demSlopeX, float* demSlopeY, float* demSlopeZ, // 地表坡度矢量
+	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,// 天线的指向
+	float Pt,// 发射能量
+	double 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 NearR, float FarR, // 距离范围
+	CUDASigmaParam* sigma0Paramslist, long sigmaparamslistlen,// 插值图
+	float* factorj, long freqnum,
+	double* outR,  // 输出距离
+	float* outAmp // 输出增益
+);
+
+
+extern __global__ void CUDAKernel_PRF_CalFreqEcho(
+	double* Rarr, float* ampArr, long pixelcount,
+	float* factorj, long freqnum,
+	double dx, double nearR,
+	cuComplex* PRFEcho, long prfid);
+
+extern __device__ cuComplex  GPU_calculationEcho(float sigma0, float TransAnt, float ReciveAnt,
+	float localangle, float R, float slopeangle, float Pt, float lamda);
+
+
+extern __global__ void CUDA_calculationEcho(float* sigma0, float* TransAnt, float* ReciveAnt,
+	float* localangle, float* R, float* slopeangle,
+	float nearRange, float Fs, float Pt, float lamda, long FreqIDmax,
+	cuComplex* echoArr, long* FreqID,
+	long len);
+
+
+extern __global__ void CUDA_BillerInterpAntPattern(float* antpattern,
+	float starttheta, float startphi, float dtheta, float dphi,
+	long thetapoints, long phipoints,
+	float* searththeta, float* searchphi, float* searchantpattern,
+	long len);
+
+
+extern __global__ void CUDA_SatelliteAntDirectNormal(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,
+	float* thetaAnt, float* phiAnt
+	, long len);
+
+
+extern __global__ void CUDAKernel_PRFSumEcho_Rows(
+	double* Rarr, float* ampArr, long Rows, long Cols,
+	long startRid,
+	float* factorj, long freqnum,
+	cuComplex* freqRowsbuffer, long tempRows
+);
+
+
+extern __global__ void  CUDAKernel_PRFSumEcho_Freq(
+	cuComplex* freqRowsbuffer, long tempRows, long freqnum,
+	cuComplex* PRFEcho, long prfid
+);
+
+
 extern "C" void SatelliteAntDirectNormal(float* RstX, float* RstY, float* RstZ, 
 	float  antXaxisX, float  antXaxisY, float  antXaxisZ, 
 	float  antYaxisX, float  antYaxisY, float  antYaxisZ, 
@@ -43,9 +137,9 @@ extern "C" void CUDAInterpSigma(
 
 // 计算坐标的 距离、增益
 extern "C" void CUDARFPC_Caluation_R_Gain(
-	float antX,float antY,float antZ, // 天线的坐标
+	float antX, float antY, float antZ, // 天线的坐标
-	float* targetX,float* targetY, float* targetZ, long TargetPixelNumber, // 地面坐标
+	float* targetX, float* targetY, float* targetZ, long TargetPixelNumber, // 地面坐标
-	long*  demCls,
+	long* demCls,
 	float* demSlopeX, float* demSlopeY, float* demSlopeZ, // 地表坡度矢量
 	float  antXaxisX, float  antXaxisY, float  antXaxisZ, // 天线坐标系的X轴
 	float  antYaxisX, float  antYaxisY, float  antYaxisZ,// 天线坐标系的Y轴
@@ -53,24 +147,22 @@ extern "C" void CUDARFPC_Caluation_R_Gain(
 	float  antDirectX, float  antDirectY, float  antDirectZ,// 天线的指向
 	float Pt,// 发射能量
 	double refPhaseRange,
-	float* TransAntpattern,  float Transtarttheta, float Transstartphi, float Transdtheta, float Transdphi, int Transthetapoints, int Transphipoints, // 发射天线方向图
+	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* ReceiveAntpattern, float Receivestarttheta, float Receivestartphi, float Receivedtheta, float Receivedphi, int Receivethetapoints, int Receivephipoints,//接收天线方向图
-	float NearR,float FarR, // 距离范围
+	float NearR, float FarR, // 距离范围
 	CUDASigmaParam* sigma0Paramslist, long sigmaparamslistlen,// 插值图
 	float* factorj, long freqnum,
 	double* outR,  // 输出距离
-	float* PRFEcho_real,float* PRFEcho_imag, long prfid
+	//float* outAmp // 输出增益
+	double* PRFEcho_real, double* PRFEcho_imag, long prfid
 );
 
 
-extern "C" void CUDA_PRF_CalFreqEcho(
+extern "C"  void CUDA_PRF_CalFreqEcho(
-	double* Rarr, float* amp, long pixelcount,// 
+	double* Rarr, float* ampArr, long pixelcount,
 	float* factorj, long freqnum,
 	double dx, double nearR,
-	cuComplex* PRFEcho, long prfid
+	cuComplex* PRFEcho, long prfid);
-);
-
-
 
 
 #endif

RasterProcessTool.vcxproj

@@ -84,6 +84,9 @@
       <LinkTimeCodeGeneration>Default</LinkTimeCodeGeneration>
       <LargeAddressAware>true</LargeAddressAware>
     </Link>
+    <CudaCompile>
+      <CodeGeneration>compute_86,sm_86</CodeGeneration>
+    </CudaCompile>
   </ItemDefinitionGroup>
   <ItemDefinitionGroup Condition="'$(Configuration)|$(Platform)' == 'Debug|x64'" Label="Configuration">
     <ClCompile>
@@ -222,12 +225,15 @@
   <ItemGroup>
     <CudaCompile Include="GPUTool\GPURFPC.cu">
       <GenerateRelocatableDeviceCode Condition="'$(Configuration)|$(Platform)'=='Release|x64'">true</GenerateRelocatableDeviceCode>
+      <CodeGeneration Condition="'$(Configuration)|$(Platform)'=='Release|x64'">compute_86,sm_86;%(CodeGeneration)</CodeGeneration>
     </CudaCompile>
     <CudaCompile Include="GPUTool\GPUTBPImage.cu">
       <GenerateRelocatableDeviceCode Condition="'$(Configuration)|$(Platform)'=='Release|x64'">true</GenerateRelocatableDeviceCode>
+      <CodeGeneration Condition="'$(Configuration)|$(Platform)'=='Release|x64'">compute_86,sm_86</CodeGeneration>
     </CudaCompile>
     <CudaCompile Include="GPUTool\GPUTool.cu">
       <GenerateRelocatableDeviceCode Condition="'$(Configuration)|$(Platform)'=='Release|x64'">true</GenerateRelocatableDeviceCode>
+      <CodeGeneration Condition="'$(Configuration)|$(Platform)'=='Release|x64'">compute_86,sm_86</CodeGeneration>
     </CudaCompile>
   </ItemGroup>
   <Import Project="$(VCTargetsPath)\Microsoft.Cpp.targets" />

SimulationSAR/RFPCProcessCls.cpp

@@ -141,7 +141,7 @@ ErrorCode RFPCProcessCls::Process(long num_thread)
 
 	//stateCode = this->RFPCMainProcess(num_thread);
 		// 初始化回波
-	this->EchoSimulationData->initEchoArr(std::complex<float>(0, 0));
+	this->EchoSimulationData->initEchoArr(std::complex<double>(0, 0));
 	stateCode = this->RFPCMainProcess_GPU();
 
 	if (stateCode != ErrorCode::SUCCESS) {
@@ -575,10 +575,10 @@ ErrorCode RFPCProcessCls::RFPCMainProcess_GPU()
 	long echoblockline = Memory1GB / 8 / 2 / PlusePoint * 2;
 	 
 	// 每一行的脉冲
-	float* h_PRFEcho_real = (float*)mallocCUDAHost(sizeof(float) * echoblockline * PlusePoint);
+	double* h_PRFEcho_real = (double*)mallocCUDAHost(sizeof(double) * echoblockline * PlusePoint);
-	float* h_PRFEcho_imag = (float*)mallocCUDAHost(sizeof(float) * echoblockline * PlusePoint);
+	double* h_PRFEcho_imag = (double*)mallocCUDAHost(sizeof(double) * echoblockline * PlusePoint);
-	float* d_PRFEcho_real = (float*)mallocCUDADevice(sizeof(float) * echoblockline * PlusePoint);
+	double* d_PRFEcho_real = (double*)mallocCUDADevice(sizeof(double) * echoblockline * PlusePoint);
-	float* d_PRFEcho_imag = (float*)mallocCUDADevice(sizeof(float) * echoblockline * PlusePoint);
+	double* d_PRFEcho_imag = (double*)mallocCUDADevice(sizeof(double) * echoblockline * PlusePoint);
 
 	float* h_factorj = (float*)mallocCUDAHost(sizeof(float) * freqlist.size());
 	float* d_factorj = (float*)mallocCUDADevice(sizeof(float) * freqlist.size());
@@ -607,6 +607,15 @@ ErrorCode RFPCProcessCls::RFPCMainProcess_GPU()
 
 		landcover = demlandcls.getData(startline, 0, newblokline, blockwidth, 1);
 
+		long calpluseFreqBufferLen = Memory1GB / 8 / 2 / PlusePoint * 2;
+		if (calpluseFreqBufferLen < 1000) {
+			qDebug() << "frequency point has morn than 50000";
+			QMessageBox::warning(nullptr,u8"frequency point has morn than 50000",u8"frequency point has morn than 50000");
+		}
+
+		cuComplex* pluseFreqPointBuffer = (cuComplex*)mallocCUDADevice(sizeof(cuComplex)*calpluseFreqBufferLen*PlusePoint);
+
+
 		if (bloklineflag) {
 			FreeCUDAHost(h_dem_x);						FreeCUDADevice(d_dem_x);
 			FreeCUDAHost(h_dem_y);						FreeCUDADevice(d_dem_y);
@@ -683,15 +692,15 @@ ErrorCode RFPCProcessCls::RFPCMainProcess_GPU()
 		long startprfid = 0;
 		for (startprfid = 0; startprfid < pluseCount; startprfid = startprfid + echoblockline) {
 			long templine = startprfid + echoblockline < PluseCount ? echoblockline : PluseCount - startprfid;
-			std::shared_ptr<std::complex<float>> echotemp = this->EchoSimulationData->getEchoArr(startprfid, templine);
+			std::shared_ptr<std::complex<double>> echotemp = this->EchoSimulationData->getEchoArr(startprfid, templine);
 			for (long tempprfid = 0; tempprfid < templine; tempprfid++) {
 				for (long freqid = 0; freqid < PlusePoint; freqid++) {
 					h_PRFEcho_real[tempprfid * PlusePoint + freqid] = echotemp.get()[tempprfid * PlusePoint + freqid].real();
 					h_PRFEcho_imag[tempprfid * PlusePoint + freqid] = echotemp.get()[tempprfid * PlusePoint + freqid].imag();
 				}
 			}
-			HostToDevice(h_PRFEcho_real, d_PRFEcho_real, sizeof(float) * echoblockline * PlusePoint);
+			HostToDevice(h_PRFEcho_real, d_PRFEcho_real, sizeof(double) * echoblockline * PlusePoint);
-			HostToDevice(h_PRFEcho_imag, d_PRFEcho_imag, sizeof(float) * echoblockline * PlusePoint);
+			HostToDevice(h_PRFEcho_imag, d_PRFEcho_imag, sizeof(double) * echoblockline * PlusePoint);
 			for (long tempprfid = 0; tempprfid < templine; tempprfid++) {
 				{// 计算
 					long prfid = tempprfid + startprfid;
@@ -717,8 +726,8 @@ ErrorCode RFPCProcessCls::RFPCMainProcess_GPU()
 					float antZaxisX = sateOirbtNodes[prfid].AntZaxisX;
 					float antZaxisY = sateOirbtNodes[prfid].AntZaxisY;
 					float antZaxisZ = sateOirbtNodes[prfid].AntZaxisZ;//18
-						// 计算距离、局地入射角、增益
+					// 计算距离、局地入射角、增益
-						CUDARFPC_Caluation_R_Gain(
+					CUDARFPC_Caluation_R_Gain(
 							antpx, antpy, antpz, // 天线的坐标
 							d_dem_x, d_dem_y, d_dem_z, pixelcount, // 地面坐标
 							d_demcls,
@@ -751,11 +760,10 @@ ErrorCode RFPCProcessCls::RFPCMainProcess_GPU()
 					h_temp_R.get()[prfid] = h_R[0];
 #endif
 				}
-
 			}
 
-			DeviceToHost(h_PRFEcho_real, d_PRFEcho_real, sizeof(float) * echoblockline * PlusePoint);
+			DeviceToHost(h_PRFEcho_real, d_PRFEcho_real, sizeof(double) * echoblockline * PlusePoint);
-			DeviceToHost(h_PRFEcho_imag, d_PRFEcho_imag, sizeof(float) * echoblockline * PlusePoint);
+			DeviceToHost(h_PRFEcho_imag, d_PRFEcho_imag, sizeof(double) * echoblockline * PlusePoint);
 			for (long tempprfid = 0; tempprfid < templine; tempprfid++) {
 				for (long freqid = 0; freqid < PlusePoint; freqid++) {
 					echotemp.get()[tempprfid * PlusePoint + freqid].real(h_PRFEcho_real[tempprfid * PlusePoint + freqid]);

SimulationSAR/TBPImageAlgCls.cpp

@@ -161,11 +161,11 @@ ErrorCode TBPImageAlgCls::Process(long num_thread)
 		}
 		std::cout << "\r[" << QDateTime::currentDateTime().toString("yyyy-MM-dd hh:mm:ss.zzz").toStdString() << "]  imgxyz :\t" << startline << "\t-\t" << startline + templine << " / " << imageheight << std::endl;
 
-		std::shared_ptr<std::complex<float>> imageRaster = this->L1ds->getImageRaster(startline, templine);
+		std::shared_ptr<std::complex<double>> imageRaster = this->L1ds->getImageRaster(startline, templine);
 
 		for (long i = 0; i < templine; i++) {
 			for (long j = 0; j < imagewidth; j++) {
-				imageRaster.get()[i * imagewidth + j] = std::complex<float>(0,0);
+				imageRaster.get()[i * imagewidth + j] = std::complex<double>(0,0);
 			}
 		}
 		this->L1ds->saveImageRaster(imageRaster, startline,templine);
@@ -253,7 +253,7 @@ ErrorCode TBPImageAlgCls::ProcessGPU()
 		std::shared_ptr<float> img_y = readDataArr<float>(imageXYZ,startimgrowid,0,tempimgBlockline,colCount,2,GDALREADARRCOPYMETHOD::VARIABLEMETHOD);
 		std::shared_ptr<float> img_z = readDataArr<float>(imageXYZ,startimgrowid,0,tempimgBlockline,colCount,3,GDALREADARRCOPYMETHOD::VARIABLEMETHOD);
 
-		std::shared_ptr<std::complex<float>> imgArr = this->L1ds->getImageRaster(startimgrowid, tempimgBlockline);
+		std::shared_ptr<std::complex<double>> imgArr = this->L1ds->getImageRaster(startimgrowid, tempimgBlockline);
 		// »ñÈ¡»Ø²¨
 		long startechoid = 0;
 		for (long startechoid = 0; startechoid < PRFCount; startechoid = startechoid + echoBlockline) {
@@ -261,7 +261,7 @@ ErrorCode TBPImageAlgCls::ProcessGPU()
 			if (startechoid + tempechoBlockline >= PRFCount) {
 				tempechoBlockline = PRFCount - startechoid;
 			}
-			std::shared_ptr<std::complex<float>> echoArr = this->L0ds->getEchoArr(startechoid, tempechoBlockline);
+			std::shared_ptr<std::complex<double>> echoArr = this->L0ds->getEchoArr(startechoid, tempechoBlockline);
 			std::shared_ptr<float> antpx(new float[tempechoBlockline*PlusePoints]);
 			std::shared_ptr<float> antpy(new float[tempechoBlockline* PlusePoints]);
 			std::shared_ptr<float> antpz(new float[tempechoBlockline* PlusePoints]);
@@ -293,8 +293,8 @@ ErrorCode TBPImageAlgCls::ProcessGPU()
 
 void TBPImageGPUAlg(std::shared_ptr<float> antPx, std::shared_ptr<float> antPy, std::shared_ptr<float> antPz,
 	std::shared_ptr<float> imgx, std::shared_ptr<float> imgy, std::shared_ptr<float> imgz,
-	std::shared_ptr<std::complex<float>> echoArr,
+	std::shared_ptr<std::complex<double>> echoArr,
-	std::shared_ptr<std::complex<float>> imgArr,
+	std::shared_ptr<std::complex<double>> imgArr,
 	float freq, float fs, float Rnear, float Rfar,
 	long rowcount, long colcount,
 	long prfcount, long freqcount,
@@ -456,7 +456,7 @@ void TBPImageGPUAlg(std::shared_ptr<float> antPx, std::shared_ptr<float> antPy,
 
 	for (long i = 0; i < rowcount; i++) {
 		for (long j = 0; j < colcount; j++) {
-			imgArr.get()[i * colcount + j] = std::complex<float>(h_imgArr[i * colcount + j].x,h_imgArr[i * colcount + j].y);
+			imgArr.get()[i * colcount + j] = std::complex<double>(h_imgArr[i * colcount + j].x,h_imgArr[i * colcount + j].y);
 		}
 	}
 

SimulationSAR/TBPImageAlgCls.h

@@ -67,8 +67,8 @@ void CreatePixelXYZ(std::shared_ptr<EchoL0Dataset> echoL0ds,QString outPixelXYZP
 void TBPImageProcess(QString echofile,QString outImageFolder,QString imagePlanePath,long num_thread);
 void TBPImageGPUAlg(std::shared_ptr<float> antPx, std::shared_ptr<float> antPy, std::shared_ptr<float> antPz,  
 	std::shared_ptr<float> img_x, std::shared_ptr<float> img_y, std::shared_ptr<float> img_z,  
-	std::shared_ptr<std::complex<float>> echoArr, 
+	std::shared_ptr<std::complex<double>> echoArr, 
-	std::shared_ptr<std::complex<float>> img_arr,
+	std::shared_ptr<std::complex<double>> img_arr,
 	float freq, float fs, float Rnear, float Rfar,
 	long rowcount, long colcount,
 	long prfcount,long freqcount,