#include "stdafx.h"
#include "RFPCProcessCls.h"
#include "BaseConstVariable.h"
#include "SARSatelliteSimulationAbstractCls.h"
#include "SARSimulationTaskSetting.h"
#include "SatelliteOribtModel.h"
#include <QDebug>
#include "ImageOperatorBase.h"
#include "GeoOperator.h"
#include "EchoDataFormat.h"
#include <QDir>
#include <QDatetime>
#include <omp.h>
#include <QProgressDialog>
#include <QMessageBox>

#ifdef DEBUGSHOWDIALOG
#include "ImageShowDialogClass.h"
#endif

#ifdef __CUDANVCC___
#include "GPUTool.cuh"
#include "GPURFPC.cuh"
#include <cuda_runtime.h>
#include <cublas_v2.h>
#endif // __CUDANVCC___
//#include <Imageshow/ImageShowDialogClass.h>







RFPCProcessCls::RFPCProcessCls()
{
	this->PluseCount = 0;
	this->PlusePoint = 0;
	this->TaskSetting = nullptr;
	this->EchoSimulationData = nullptr;
 	this->LandCoverPath = "";
 	this->OutEchoPath = "";

 
	this->LandCoverPath.clear();
 	this->OutEchoPath.clear();
	this->SigmaDatabasePtr = std::shared_ptr<SigmaDatabase>(new SigmaDatabase);

}

RFPCProcessCls::~RFPCProcessCls()
{
}

void RFPCProcessCls::setTaskSetting(std::shared_ptr < AbstractSARSatelliteModel> TaskSetting)
{
	this->TaskSetting = std::shared_ptr < AbstractSARSatelliteModel>(TaskSetting);
	qDebug() << "RFPCProcessCls::setTaskSetting";
}

void RFPCProcessCls::setEchoSimulationDataSetting(std::shared_ptr<EchoL0Dataset> EchoSimulationData)
{
	this->EchoSimulationData = std::shared_ptr<EchoL0Dataset>(EchoSimulationData);
	qDebug() << "RFPCProcessCls::setEchoSimulationDataSetting";
}

void RFPCProcessCls::setTaskFileName(QString EchoFileName)
{
	this->TaskFileName = EchoFileName;
}

void RFPCProcessCls::setDEMTiffPath(QString DEMTiffPath)
{
	this->demxyzPath = DEMTiffPath;
}

void RFPCProcessCls::setSloperPath(QString InSloperPath)
{
	this->demsloperPath = InSloperPath;
}

void RFPCProcessCls::setLandCoverPath(QString LandCoverPath)
{
	this->LandCoverPath = LandCoverPath;
}


void RFPCProcessCls::setOutEchoPath(QString OutEchoPath)
{
	this->OutEchoPath = OutEchoPath;
}



ErrorCode RFPCProcessCls::Process(long num_thread)
{
	// RFPC 算法
	qDebug() << u8"params init ....";
	ErrorCode stateCode = this->InitParams();
	if (stateCode != ErrorCode::SUCCESS) {
		return stateCode;
	}
	else {}
 
	qDebug() << "RFPCMainProcess";

	stateCode = this->InitEchoMaskArray();
	if (stateCode != ErrorCode::SUCCESS) {
		return stateCode;
	}
	else {}
	qDebug() << "InitEchoMaskArray";

	//stateCode = this->RFPCMainProcess(num_thread);
		// 初始化回波
	this->EchoSimulationData->initEchoArr(std::complex<double>(0, 0));
	stateCode = this->RFPCMainProcess_GPU();

	if (stateCode != ErrorCode::SUCCESS) {
		return stateCode;
	}
	else {}


	return ErrorCode::SUCCESS;
}

ErrorCode RFPCProcessCls::InitParams()
{
	if (nullptr == this->TaskSetting || this->demxyzPath.isEmpty() ||
		this->LandCoverPath.isEmpty() || this->demsloperPath.isEmpty()) {
		return ErrorCode::RFPC_PARAMSISEMPTY;
	}
	else {

	}

	// 归一化绝对路径

	this->OutEchoPath = QDir(this->OutEchoPath).absolutePath();
	// 回波大小
	double imgStart_end = this->TaskSetting->getSARImageEndTime() - this->TaskSetting->getSARImageStartTime();
	this->PluseCount = ceil(imgStart_end * this->TaskSetting->getPRF());

	double rangeTimeSample = (this->TaskSetting->getFarRange() - this->TaskSetting->getNearRange()) * 2.0 / LIGHTSPEED;
	this->PlusePoint = ceil(rangeTimeSample * this->TaskSetting->getFs());

	// 初始化回波存放位置
	qDebug() << "--------------Echo Data Setting ---------------------------------------";
	this->EchoSimulationData = std::shared_ptr<EchoL0Dataset>(new EchoL0Dataset);
	this->EchoSimulationData->setCenterFreq(this->TaskSetting->getCenterFreq());
	this->EchoSimulationData->setNearRange(this->TaskSetting->getNearRange());
	this->EchoSimulationData->setFarRange(this->TaskSetting->getFarRange());
	this->EchoSimulationData->setFs(this->TaskSetting->getFs());
	this->EchoSimulationData->setBandwidth(this->TaskSetting->getBandWidth());
	this->EchoSimulationData->setCenterAngle(this->TaskSetting->getCenterLookAngle());
	this->EchoSimulationData->setLookSide(this->TaskSetting->getIsRightLook() ? "R" : "L");
	this->EchoSimulationData->OpenOrNew(OutEchoPath, TaskFileName, PluseCount, PlusePoint);

	QString tmpfolderPath = QDir(OutEchoPath).filePath("tmp");
	if (QDir(tmpfolderPath).exists() == false) {
		QDir(OutEchoPath).mkpath(tmpfolderPath);
	}
	this->tmpfolderPath = tmpfolderPath;
	return ErrorCode::SUCCESS;
}
 
ErrorCode RFPCProcessCls::InitEchoMaskArray()
{
	QString name = this->EchoSimulationData->getSimulationTaskName();
	this->OutEchoMaskPath = JoinPath(this->OutEchoPath, name + "_echomask.bin");
	Eigen::MatrixXd gt(2, 3);
	gt(0, 0) = 0;
	gt(0, 1) = 1;
	gt(0, 2) = 0;
	gt(1, 0) = 0;
	gt(1, 1) = 0;
	gt(1, 2) = 1;
	gdalImage echomaskImg = CreategdalImage(this->OutEchoMaskPath,
		this->EchoSimulationData->getPluseCount(),
		this->EchoSimulationData->getPlusePoints(),
		1,
		gt, "",
		false, true, true);
	long cols = this->EchoSimulationData->getPlusePoints();
	long rows = this->EchoSimulationData->getPluseCount();
	long blocksize = Memory1GB / 8 / this->EchoSimulationData->getPlusePoints() * 4;
	for (long startid = 0; startid < rows; startid = startid + blocksize) {
		Eigen::MatrixXd data = echomaskImg.getData(startid, 0, blocksize, cols, 1);
		data = data.array() * 0;
		echomaskImg.saveImage(data, startid, 0, 1);
	}
	return ErrorCode::SUCCESS;
}


std::shared_ptr<SatelliteOribtNode[]> RFPCProcessCls::getSatelliteOribtNodes(double prf_time, double dt, bool antflag, long double imageStarttime)
{


	std::shared_ptr<SatelliteOribtNode[]> sateOirbtNodes(new SatelliteOribtNode[this->PluseCount], delArrPtr);
	{ // 姿态计算不同
		qDebug() << "Ant position finished started !!!";
		// 计算姿态
		std::shared_ptr<double> antpos = this->EchoSimulationData->getAntPos();
		double dAt = 1e-6;
		double prf_time_dt = 0;
		Landpoint InP{ 0,0,0 }, outP{ 0,0,0 };
		for (long prf_id = 0; prf_id < this->PluseCount; prf_id++) {
			prf_time = dt * prf_id;
			prf_time_dt = prf_time + dAt;
			SatelliteOribtNode sateOirbtNode;
			SatelliteOribtNode sateOirbtNode_dAt;
			this->TaskSetting->getSatelliteOribtNode(prf_time, sateOirbtNode, antflag);
			this->TaskSetting->getSatelliteOribtNode(prf_time_dt, sateOirbtNode_dAt, antflag);

			sateOirbtNode.AVx = (sateOirbtNode_dAt.Vx - sateOirbtNode.Vx) / dAt; // 加速度
			sateOirbtNode.AVy = (sateOirbtNode_dAt.Vy - sateOirbtNode.Vy) / dAt;
			sateOirbtNode.AVz = (sateOirbtNode_dAt.Vz - sateOirbtNode.Vz) / dAt;

			InP.lon = sateOirbtNode.Px;
			InP.lat = sateOirbtNode.Py;
			InP.ati = sateOirbtNode.Pz;

			outP = XYZ2LLA(InP);

			antpos.get()[prf_id * 19 + 0] = prf_time + imageStarttime;
			antpos.get()[prf_id * 19 + 1] = sateOirbtNode.Px;
			antpos.get()[prf_id * 19 + 2] = sateOirbtNode.Py;
			antpos.get()[prf_id * 19 + 3] = sateOirbtNode.Pz;
			antpos.get()[prf_id * 19 + 4] = sateOirbtNode.Vx;
			antpos.get()[prf_id * 19 + 5] = sateOirbtNode.Vy;
			antpos.get()[prf_id * 19 + 6] = sateOirbtNode.Vz;
			antpos.get()[prf_id * 19 + 7] = sateOirbtNode.AntDirecX;
			antpos.get()[prf_id * 19 + 8] = sateOirbtNode.AntDirecY;
			antpos.get()[prf_id * 19 + 9] = sateOirbtNode.AntDirecZ;
			antpos.get()[prf_id * 19 + 10] = sateOirbtNode.AVx;
			antpos.get()[prf_id * 19 + 11] = sateOirbtNode.AVy;
			antpos.get()[prf_id * 19 + 12] = sateOirbtNode.AVz;
			antpos.get()[prf_id * 19 + 13] = sateOirbtNode.zeroDopplerDirectX;
			antpos.get()[prf_id * 19 + 14] = sateOirbtNode.zeroDopplerDirectY;
			antpos.get()[prf_id * 19 + 15] = sateOirbtNode.zeroDopplerDirectZ;
			antpos.get()[prf_id * 19 + 16] = outP.lon;
			antpos.get()[prf_id * 19 + 17] = outP.lat;
			antpos.get()[prf_id * 19 + 18] = outP.ati;
			sateOirbtNodes[prf_id] = sateOirbtNode;
		}
		this->EchoSimulationData->saveAntPos(antpos);
		antpos.reset();
		qDebug() << "Ant position finished sucessfully !!!";
	}

	return sateOirbtNodes;
}


void RFPCProcessMain(long num_thread,
	QString TansformPatternFilePath, QString ReceivePatternFilePath, 
	QString simulationtaskName, QString OutEchoPath, 
	QString GPSXmlPath, QString TaskXmlPath,QString demTiffPath, QString sloperPath, QString LandCoverPath) 
{

	std::shared_ptr < AbstractSARSatelliteModel>  task = ReadSimulationSettingsXML(TaskXmlPath);

	if (nullptr == task)
	{
		return;
	}
	else {
		// 打印参数
		qDebug() << "--------------Task Seting ---------------------------------------";
		qDebug() << "SARImageStartTime: " << task->getSARImageStartTime();
		qDebug() << "SARImageEndTime: " << task->getSARImageEndTime();
		qDebug() << "BandWidth: " << task->getBandWidth();
		qDebug() << "CenterFreq: " << task->getCenterFreq();
		qDebug() << "PRF: " << task->getPRF();
		qDebug() << "Lamda: " << task->getCenterLamda();
		qDebug() << "Fs: " << task->getFs();
		qDebug() << "POLAR: " << task->getPolarType();
		qDebug() << "NearRange: " << task->getNearRange();
		qDebug() << "FarRange: " << task->getFarRange();
		qDebug() << (task->getFarRange() - task->getNearRange()) * 2 / LIGHTSPEED * task->getFs();
		qDebug() << "\n\n";
	}
	// 1.2 设置天线方向图
	std::vector<RadiationPatternGainPoint> TansformPatternGainpoints = ReadGainFile(TansformPatternFilePath);
	std::shared_ptr<AbstractRadiationPattern> TansformPatternGainPtr = CreateAbstractRadiationPattern(TansformPatternGainpoints);

	std::vector<RadiationPatternGainPoint> ReceivePatternGainpoints = ReadGainFile(ReceivePatternFilePath);
	std::shared_ptr<AbstractRadiationPattern> ReceivePatternGainPtr = CreateAbstractRadiationPattern(ReceivePatternGainpoints);

	task->setTransformRadiationPattern(TansformPatternGainPtr);
	task->setReceiveRadiationPattern(ReceivePatternGainPtr);

	//2. 读取GPS节点
	std::vector<SatelliteOribtNode> nodes;
	ErrorCode stateCode = ReadSateGPSPointsXML(GPSXmlPath, nodes);

	if (stateCode != ErrorCode::SUCCESS)
	{
		qWarning() << QString::fromStdString(errorCode2errInfo(stateCode));
		return;
	}
	else {}

	std::shared_ptr<AbstractSatelliteOribtModel> SatelliteOribtModel = CreataPolyfitSatelliteOribtModel(nodes, task->getSARImageStartTime(), 3); // 以成像开始时间作为 时间参考起点
	SatelliteOribtModel->setbeamAngle(task->getCenterLookAngle(), task->getIsRightLook()); // 设置天线方向图

	if (nullptr == SatelliteOribtModel)
	{
		return;
	}
	else {
		task->setSatelliteOribtModel(SatelliteOribtModel);
	}

	qDebug() << "--------------  RFPC init ---------------------------------------";
	RFPCProcessCls RFPC;
	RFPC.setTaskSetting(task);						//qDebug() << "setTaskSetting";
	RFPC.setTaskFileName(simulationtaskName);		//qDebug() << "setTaskFileName";
	RFPC.setDEMTiffPath(demTiffPath);				//qDebug() << "setDEMTiffPath";
	RFPC.setSloperPath(sloperPath);				//qDebug() << "setDEMTiffPath";
	RFPC.setLandCoverPath(LandCoverPath);			//qDebug() << "setLandCoverPath";
	RFPC.setOutEchoPath(OutEchoPath);				//qDebug() << "setOutEchoPath";
	qDebug() << "--------------  RFPC start---------------------------------------";
	RFPC.Process(num_thread); // 处理程序
	qDebug() << "--------------  RFPC end---------------------------------------";
}


ErrorCode RFPCProcessCls::RFPCMainProcess_GPU() {
	/** 内存分配***************************************************/
	long TargetMemoryMB = 500;


	/** 参数区域***************************************************/
	QVector<double> freqlist = this->TaskSetting->getFreqList();
	long freqnum = freqlist.count();
	float f0 = float(freqlist[0] / 1e9);
	float dfreq = float((freqlist[1] - freqlist[0]) / 1e9);
	long PRFCount = this->EchoSimulationData->getPluseCount();

	double NearRange = this->EchoSimulationData->getNearRange(); // 近斜距
	double FarRange = this->EchoSimulationData->getFarRange();

	double Pt = this->TaskSetting->getPt() * this->TaskSetting->getGri();// 发射电压 1v

	double lamda = this->TaskSetting->getCenterLamda(); // 波长
	double refphaseRange = this->TaskSetting->getRefphaseRange(); // 参考相位斜距

	double prf_time = 0;
	double dt = 1 / this->TaskSetting->getPRF();// 获取每次脉冲的时间间隔
	bool antflag = true; // 计算天线方向图
	long double  imageStarttime = this->TaskSetting->getSARImageStartTime();
 
	// 卫星
	double* h_antpx, * d_antpx;
	double* h_antpy, * d_antpy;
	double* h_antpz, * d_antpz;
	double* h_antvx, * d_antvx;
	double* h_antvy, * d_antvy;
	double* h_antvz, * d_antvz;
	double* h_antdirectx, * d_antdirectx;
	double* h_antdirecty, * d_antdirecty;
	double* h_antdirectz, * d_antdirectz;
	double* h_antXaxisX, * d_antXaxisX;
	double* h_antXaxisY, * d_antXaxisY;
	double* h_antXaxisZ, * d_antXaxisZ;
	double* h_antYaxisX, * d_antYaxisX;
	double* h_antYaxisY, * d_antYaxisY;
	double* h_antYaxisZ, * d_antYaxisZ;
	double* h_antZaxisX, * d_antZaxisX;
	double* h_antZaxisY, * d_antZaxisY;
	double* h_antZaxisZ, * d_antZaxisZ;
	
	
	{
		h_antpx			= (double*)mallocCUDAHost(sizeof(double) * PRFCount);
		h_antpy			= (double*)mallocCUDAHost(sizeof(double) * PRFCount);
		h_antpz			= (double*)mallocCUDAHost(sizeof(double) * PRFCount);
		h_antvx			= (double*)mallocCUDAHost(sizeof(double) * PRFCount);
		h_antvy			= (double*)mallocCUDAHost(sizeof(double) * PRFCount);
		h_antvz			= (double*)mallocCUDAHost(sizeof(double) * PRFCount);
		h_antdirectx	= (double*)mallocCUDAHost(sizeof(double) * PRFCount);
		h_antdirecty	= (double*)mallocCUDAHost(sizeof(double) * PRFCount);
		h_antdirectz	= (double*)mallocCUDAHost(sizeof(double) * PRFCount);
		h_antXaxisX		= (double*)mallocCUDAHost(sizeof(double) * PRFCount);
		h_antXaxisY		= (double*)mallocCUDAHost(sizeof(double) * PRFCount);
		h_antXaxisZ		= (double*)mallocCUDAHost(sizeof(double) * PRFCount);
		h_antYaxisX		= (double*)mallocCUDAHost(sizeof(double) * PRFCount);
		h_antYaxisY		= (double*)mallocCUDAHost(sizeof(double) * PRFCount);
		h_antYaxisZ		= (double*)mallocCUDAHost(sizeof(double) * PRFCount);
		h_antZaxisX		= (double*)mallocCUDAHost(sizeof(double) * PRFCount);
		h_antZaxisY		= (double*)mallocCUDAHost(sizeof(double) * PRFCount);
		h_antZaxisZ		= (double*)mallocCUDAHost(sizeof(double) * PRFCount);


		d_antpx			= (double*)mallocCUDADevice(sizeof(double) * PRFCount);
		d_antpy			= (double*)mallocCUDADevice(sizeof(double) * PRFCount);
		d_antpz			= (double*)mallocCUDADevice(sizeof(double) * PRFCount);
		d_antvx			= (double*)mallocCUDADevice(sizeof(double) * PRFCount);
		d_antvy			= (double*)mallocCUDADevice(sizeof(double) * PRFCount);
		d_antvz			= (double*)mallocCUDADevice(sizeof(double) * PRFCount);
		d_antdirectx	= (double*)mallocCUDADevice(sizeof(double) * PRFCount);
		d_antdirecty	= (double*)mallocCUDADevice(sizeof(double) * PRFCount);
		d_antdirectz	= (double*)mallocCUDADevice(sizeof(double) * PRFCount);
		d_antXaxisX		= (double*)mallocCUDADevice(sizeof(double) * PRFCount);
		d_antXaxisY		= (double*)mallocCUDADevice(sizeof(double) * PRFCount);
		d_antXaxisZ		= (double*)mallocCUDADevice(sizeof(double) * PRFCount);
		d_antYaxisX		= (double*)mallocCUDADevice(sizeof(double) * PRFCount);
		d_antYaxisY		= (double*)mallocCUDADevice(sizeof(double) * PRFCount);
		d_antYaxisZ		= (double*)mallocCUDADevice(sizeof(double) * PRFCount);
		d_antZaxisX		= (double*)mallocCUDADevice(sizeof(double) * PRFCount);
		d_antZaxisY		= (double*)mallocCUDADevice(sizeof(double) * PRFCount);
		d_antZaxisZ		= (double*)mallocCUDADevice(sizeof(double) * PRFCount);


		this->EchoSimulationData->getAntPos();
		std::shared_ptr<SatelliteOribtNode[]> sateOirbtNodes = this->getSatelliteOribtNodes(prf_time, dt, antflag, imageStarttime);
		for (long tempprfid = 0; tempprfid < PRFCount; tempprfid++) {
			long prfid = tempprfid;
			h_antpx[tempprfid] = sateOirbtNodes[prfid].Px;
			h_antpy[tempprfid] = sateOirbtNodes[prfid].Py;
			h_antpz[tempprfid] = sateOirbtNodes[prfid].Pz;
			h_antvx[tempprfid] = sateOirbtNodes[prfid].Vx;
			h_antvy[tempprfid] = sateOirbtNodes[prfid].Vy;
			h_antvz[tempprfid] = sateOirbtNodes[prfid].Vz; //6
			h_antdirectx[tempprfid] = sateOirbtNodes[prfid].AntDirecX;
			h_antdirecty[tempprfid] = sateOirbtNodes[prfid].AntDirecY;
			h_antdirectz[tempprfid] = sateOirbtNodes[prfid].AntDirecZ; // 9 天线指向
			h_antXaxisX[tempprfid] = sateOirbtNodes[prfid].AntXaxisX;
			h_antXaxisY[tempprfid] = sateOirbtNodes[prfid].AntXaxisY;
			h_antXaxisZ[tempprfid] = sateOirbtNodes[prfid].AntXaxisZ;//12 天线坐标系
			h_antYaxisX[tempprfid] = sateOirbtNodes[prfid].AntYaxisX;
			h_antYaxisY[tempprfid] = sateOirbtNodes[prfid].AntYaxisY;
			h_antYaxisZ[tempprfid] = sateOirbtNodes[prfid].AntYaxisZ;//15
			h_antZaxisX[tempprfid] = sateOirbtNodes[prfid].AntZaxisX;
			h_antZaxisY[tempprfid] = sateOirbtNodes[prfid].AntZaxisY;
			h_antZaxisZ[tempprfid] = sateOirbtNodes[prfid].AntZaxisZ;//18
		}

		DeviceToDevice(h_antpx, d_antpx, sizeof(double) * PRFCount);
		DeviceToDevice(h_antpy, d_antpy, sizeof(double) * PRFCount);
		DeviceToDevice(h_antpz, d_antpz, sizeof(double) * PRFCount);
		DeviceToDevice(h_antvx, d_antvx, sizeof(double) * PRFCount);
		DeviceToDevice(h_antvy, d_antvy, sizeof(double) * PRFCount);
		DeviceToDevice(h_antvz, d_antvz, sizeof(double) * PRFCount);
		DeviceToDevice(h_antdirectx, d_antdirectx, sizeof(double) * PRFCount);
		DeviceToDevice(h_antdirecty, d_antdirecty, sizeof(double) * PRFCount);
		DeviceToDevice(h_antdirectz, d_antdirectz, sizeof(double) * PRFCount);
		DeviceToDevice(h_antXaxisX, d_antXaxisX, sizeof(double) * PRFCount);
		DeviceToDevice(h_antXaxisY, d_antXaxisY, sizeof(double) * PRFCount);
		DeviceToDevice(h_antXaxisZ, d_antXaxisZ, sizeof(double) * PRFCount);
		DeviceToDevice(h_antYaxisX, d_antYaxisX, sizeof(double) * PRFCount);
		DeviceToDevice(h_antYaxisY, d_antYaxisY, sizeof(double) * PRFCount);
		DeviceToDevice(h_antYaxisZ, d_antYaxisZ, sizeof(double) * PRFCount);
		DeviceToDevice(h_antZaxisX, d_antZaxisX, sizeof(double) * PRFCount);
		DeviceToDevice(h_antZaxisY, d_antZaxisY, sizeof(double) * PRFCount);
		DeviceToDevice(h_antZaxisZ, d_antZaxisZ, sizeof(double) * PRFCount);

	}


	/** 天线方向图***************************************************/
	std::shared_ptr<AbstractRadiationPattern> TransformPattern = this->TaskSetting->getTransformRadiationPattern(); // 发射天线方向图

	std::shared_ptr<AbstractRadiationPattern> ReceivePattern = this->TaskSetting->getReceiveRadiationPattern();   // 接收天线方向图

	POLARTYPEENUM polartype = this->TaskSetting->getPolarType();
	PatternImageDesc TantPatternDesc = {};
	double* h_TantPattern = nullptr;
	double* d_TantPattern = nullptr;

	{
		// 处理发射天线方向图
		double Tminphi = TransformPattern->getMinPhi();
		double Tmaxphi = TransformPattern->getMaxPhi();
		double Tmintheta = TransformPattern->getMinTheta();
		double Tmaxtheta = TransformPattern->getMaxTheta();

		long Tphinum = TransformPattern->getPhis().size();
		long Tthetanum = TransformPattern->getThetas().size();

		double TstartTheta = Tmintheta;
		double TstartPhi = Tminphi;

		double Tdtheta = (Tmaxtheta - Tmintheta) / (Tthetanum - 1);
		double Tdphi = (Tmaxphi - Tminphi) / (Tphinum - 1);

		h_TantPattern = (double*)mallocCUDAHost(sizeof(double) * Tthetanum * Tphinum);
		d_TantPattern = (double*)mallocCUDADevice(sizeof(double) * Tthetanum * Tphinum);

		for (long i = 0; i < Tthetanum; i++) {
			for (long j = Tphinum - 1; j >= 0; j--) {
				//h_TantPattern[i * Tphinum + j] = TransformPattern->getGainLearThetaPhi(TstartTheta + i * Tdtheta, TstartPhi + j * Tdphi);
				h_TantPattern[i * Tphinum + j] = TransformPattern->getGain(TstartTheta + i * Tdtheta, TstartPhi + j * Tdphi);
			}
		}

		testOutAntPatternTrans("TransPattern.bin", h_TantPattern, TstartTheta, Tdtheta, TstartPhi, Tdphi, Tthetanum, Tphinum);
		for (long i = 0; i < Tthetanum; i++) {
			for (long j = 0; j < Tphinum; j++) {
				h_TantPattern[i * Tphinum + j] = powf(10.0, h_TantPattern[i * Tphinum + j] / 10);
			}
		}
		HostToDevice(h_TantPattern, d_TantPattern, sizeof(double) * Tthetanum * Tphinum);
		TantPatternDesc.startTheta = TstartTheta;
		TantPatternDesc.startPhi = TstartPhi;
		TantPatternDesc.dtheta = Tdtheta;
		TantPatternDesc.dphi = Tdphi;
		TantPatternDesc.phinum = Tphinum;
		TantPatternDesc.thetanum = Tthetanum;
	}

	PatternImageDesc RantPatternDesc = {};
	double* h_RantPattern = nullptr;
	double* d_RantPattern = nullptr;
	{
		// 处理接收天线方向图
		double Rminphi = ReceivePattern->getMinPhi();
		double Rmaxphi = ReceivePattern->getMaxPhi();
		double Rmintheta = ReceivePattern->getMinTheta();
		double Rmaxtheta = ReceivePattern->getMaxTheta();

		long Rphinum = ReceivePattern->getPhis().size();
		long Rthetanum = ReceivePattern->getThetas().size();

		double RstartTheta = Rmintheta;
		double RstartPhi = Rminphi;

		double Rdtheta = (Rmaxtheta - Rmintheta) / (Rthetanum - 1);
		double Rdphi = (Rmaxphi - Rminphi) / (Rphinum - 1);

		h_RantPattern = (double*)mallocCUDAHost(sizeof(double) * Rthetanum * Rphinum);
		d_RantPattern = (double*)mallocCUDADevice(sizeof(double) * Rthetanum * Rphinum);

		for (long i = 0; i < Rthetanum; i++) {
			for (long j = 0; j < Rphinum; j++) {
				//h_RantPattern[i * Rphinum + j] = ReceivePattern->getGainLearThetaPhi(RstartTheta + i * Rdtheta, RstartPhi + j * Rdphi);
				h_RantPattern[i * Rphinum + j] = ReceivePattern->getGain(RstartTheta + i * Rdtheta, RstartPhi + j * Rdphi);
			}
		}

		testOutAntPatternTrans("ReceivePattern.bin", h_RantPattern, Rmintheta, Rdtheta, RstartPhi, Rdphi, Rthetanum, Rphinum);
		for (long i = 0; i < Rthetanum; i++) {
			for (long j = 0; j < Rphinum; j++) {
				h_RantPattern[i * Rphinum + j] = powf(10.0, h_RantPattern[i * Rphinum + j] / 10);
			}
		}
		HostToDevice(h_RantPattern, d_RantPattern, sizeof(double) * Rthetanum * Rphinum);
		RantPatternDesc.startTheta = RstartTheta;
		RantPatternDesc.startPhi = RstartPhi;
		RantPatternDesc.dtheta = Rdtheta;
		RantPatternDesc.dphi = Rdphi;
		RantPatternDesc.phinum = Rphinum;
		RantPatternDesc.thetanum = Rthetanum;
	}



	/** 坐标区域点***************************************************/
	gdalImage demxyz(this->demxyzPath);// 地面点坐标
	gdalImage demlandcls(this->LandCoverPath);// 地表覆盖类型 
	gdalImage demsloperxyz(this->demsloperPath);// 地面坡向

	long demRow = demxyz.height;
	long demCol = demxyz.width;


	//处理地表覆盖
	QMap<long, long> clamap;
	long clamapid = 0;
	long startline = 0;

	{
		long blokline = getBlockRows(2e4, demCol, sizeof(double));
		for (startline = 0; startline < demRow; startline = startline + blokline) {
			Eigen::MatrixXd clsland = demlandcls.getData(startline, 0, blokline, demlandcls.width, 1);
			long clsrows = clsland.rows();
			long clscols = clsland.cols();
			long clsid = 0;
			for (long ii = 0; ii < clsrows; ii++) {
				for (long jj = 0; jj < clscols; jj++) {
					clsid = clsland(ii, jj);
					if (clamap.contains(clsid)) {}
					else {
						clamap.insert(clsid, clamapid);
						clamapid = clamapid + 1;
					}
				}
			}
		}

		std::cout << "class id recoding" << std::endl;
		for (long id : clamap.keys()) {
			std::cout << id << " -> " << clamap[id] << std::endl;
		}
	}

	CUDASigmaParam* h_clsSigmaParam = (CUDASigmaParam*)mallocCUDAHost(sizeof(CUDASigmaParam) * clamapid);
	CUDASigmaParam* d_clsSigmaParam = (CUDASigmaParam*)mallocCUDADevice(sizeof(CUDASigmaParam) * clamapid);

	{
		std::map<long, SigmaParam> tempSigmaParam = this->SigmaDatabasePtr->getsigmaParams(polartype);
		for (long id : clamap.keys()) {
			SigmaParam tempp = tempSigmaParam[id];
			h_clsSigmaParam[clamap[id]].p1 = tempp.p1;
			h_clsSigmaParam[clamap[id]].p2 = tempp.p2;
			h_clsSigmaParam[clamap[id]].p3 = tempp.p3;
			h_clsSigmaParam[clamap[id]].p4 = tempp.p4;
			h_clsSigmaParam[clamap[id]].p5 = tempp.p5;
			h_clsSigmaParam[clamap[id]].p6 = tempp.p6;
		}

		// 打印日志
		std::cout << "sigma params:" << std::endl;
		std::cout << "classid:\tp1\tp2\tp3\tp4\tp5\tp6" << std::endl;
		for (long ii = 0; ii < clamapid; ii++) {
			std::cout << ii << ":\t" << h_clsSigmaParam[ii].p1;
			std::cout << "\t" << h_clsSigmaParam[ii].p2;
			std::cout << "\t" << h_clsSigmaParam[ii].p3;
			std::cout << "\t" << h_clsSigmaParam[ii].p4;
			std::cout << "\t" << h_clsSigmaParam[ii].p5;
			std::cout << "\t" << h_clsSigmaParam[ii].p6 << std::endl;
		}
		std::cout << "";
	}

	HostToDevice(h_clsSigmaParam, d_clsSigmaParam, sizeof(CUDASigmaParam) * clamapid);

	// 处理地面坐标
	long blockline = getBlockRows(TargetMemoryMB, demCol, sizeof(double));
	double* h_dem_x = (double*)mallocCUDAHost(sizeof(double) * blockline * demCol);
	double* h_dem_y = (double*)mallocCUDAHost(sizeof(double) * blockline * demCol);
	double* h_dem_z = (double*)mallocCUDAHost(sizeof(double) * blockline * demCol);
	double* h_demsloper_x = (double*)mallocCUDAHost(sizeof(double) * blockline * demCol);
	double* h_demsloper_y = (double*)mallocCUDAHost(sizeof(double) * blockline * demCol);
	double* h_demsloper_z = (double*)mallocCUDAHost(sizeof(double) * blockline * demCol);
	long* h_demcls = (long*)mallocCUDAHost(sizeof(long) * blockline * demCol);


	/** 处理回波***************************************************/
	long echo_block_rows = getBlockRows(5000, freqnum, sizeof(float));
	float* h_echo_block_real = (float*)mallocCUDAHost(sizeof(float) * echo_block_rows * freqnum);
	float* h_echo_block_imag = (float*)mallocCUDAHost(sizeof(float) * echo_block_rows * freqnum);



	/** 主流程处理 ***************************************************/

	for (long sprfid = 0; sprfid < PRFCount; sprfid = sprfid + echo_block_rows) {
		long PRF_len = (sprfid + echo_block_rows) < PRFCount ? echo_block_rows : (PRFCount - sprfid);
		std::shared_ptr<std::complex<double>> echo_temp = this->EchoSimulationData->getEchoArr(sprfid, PRF_len);
		for (long ii = 0; ii < PRF_len; ii++) {
			for (long jj = 0; jj < freqnum; jj++) {
				h_echo_block_real[ii * freqnum + jj]=echo_temp.get()[ii * freqnum + jj].real();
				h_echo_block_imag[ii * freqnum + jj]=echo_temp.get()[ii * freqnum + jj].imag();
			}
		}

		for (long startline = 0; startline < demRow; startline = startline + blockline) {
			Eigen::MatrixXd dem_x = demxyz.getData(startline, 0, blockline, demCol, 1);  // 地面坐标
			Eigen::MatrixXd dem_y = demxyz.getData(startline, 0, blockline, demCol, 2);
			Eigen::MatrixXd dem_z = demxyz.getData(startline, 0, blockline, demCol, 3);
			Eigen::MatrixXd demsloper_x = demsloperxyz.getData(startline, 0, blockline, demCol, 1);
			Eigen::MatrixXd demsloper_y = demsloperxyz.getData(startline, 0, blockline, demCol, 2);
			Eigen::MatrixXd demsloper_z = demsloperxyz.getData(startline, 0, blockline, demCol, 3);
			Eigen::MatrixXd landcover = demlandcls.getData(startline, 0, blockline, demCol, 1);

			long temp_dem_row = dem_x.rows();
			long temp_dem_col = dem_x.cols();
			long temp_dem_count = dem_x.count();
			// 更新数据格式
			for (long i = 0; i < temp_dem_row; i++) {
				for (long j = 0; j < temp_dem_col; j++) {

#ifdef __PRFDEBUG__
					h_dem_x[i * temp_dem_col + j] = -2028380.6250000; double(dem_x(i, j));
					h_dem_y[i * temp_dem_col + j] = 4139373.250000; double(dem_y(i, j));
					h_dem_z[i * temp_dem_col + j] = 4393382.500000; double(dem_z(i, j));
					h_demsloper_x[i * temp_dem_col + j] = 4393382.500000; double(demsloper_x(i, j));
					h_demsloper_y[i * temp_dem_col + j] = 446.923950; double(demsloper_y(i, j));
					h_demsloper_z[i * temp_dem_col + j] = -219.002213; double(demsloper_z(i, j));
					h_demcls[i * temp_dem_col + j] = clamap[80];// clamap[long(landcover(i, j))];
#else
					h_dem_x[i * temp_dem_col + j] = double(dem_x(i, j));
					h_dem_y[i * temp_dem_col + j] = double(dem_y(i, j));
					h_dem_z[i * temp_dem_col + j] = double(dem_z(i, j));
					h_demsloper_x[i * temp_dem_col + j] = double(demsloper_x(i, j));
					h_demsloper_y[i * temp_dem_col + j] = double(demsloper_y(i, j));
					h_demsloper_z[i * temp_dem_col + j] = double(demsloper_z(i, j));
					h_demcls[i * temp_dem_col + j] = clamap[long(landcover(i, j))];
#endif			


				}
			}

			// 分块处理
			CUDA_RFPC_MainProcess(
				d_antpx,		d_antpy,		d_antpz,
				d_antXaxisX,	d_antXaxisY,	d_antXaxisZ, // 天线坐标系的X轴
				d_antYaxisX,	d_antYaxisY,	d_antYaxisZ,// 天线坐标系的Y轴
				d_antZaxisX,	d_antZaxisY,	d_antZaxisZ,// 天线坐标系的Z轴
				d_antdirectx,	d_antdirecty,	d_antdirectz,// 天线的指向
				PRFCount, freqnum,
				f0,dfreq,
				Pt, 
				refphaseRange,
				// 天线方向图
				d_TantPattern,
				TantPatternDesc.startTheta, TantPatternDesc.startPhi, TantPatternDesc.dtheta, TantPatternDesc.dphi, TantPatternDesc.thetanum, TantPatternDesc.phinum,
				d_RantPattern,
				RantPatternDesc.startTheta, RantPatternDesc.startPhi, RantPatternDesc.dtheta, RantPatternDesc.dphi, RantPatternDesc.thetanum, RantPatternDesc.phinum,
				
				NearRange, FarRange, // 近斜据
				
				h_dem_x, h_dem_y, h_dem_z, h_demcls, temp_dem_count, // 地面坐标				
				h_demsloper_x, h_demsloper_y, h_demsloper_z, // 地表坡度矢量
				d_clsSigmaParam, clamapid,

				h_echo_block_real, h_echo_block_imag// 输出回波
			);

			PRINT("dem : %d - %d / %d , echo: %d -%d / %d", startline, startline+ temp_dem_row, demRow, sprfid, sprfid+ PRF_len, PRFCount);
		}


		for (long ii = 0; ii < PRF_len; ii++) {
			for (long jj = 0; jj < freqnum; jj++) {
				echo_temp.get()[ii * freqnum + jj].real(h_echo_block_real[ii * freqnum + jj]);
				echo_temp.get()[ii * freqnum + jj].imag(h_echo_block_imag[ii * freqnum + jj]);
			}
		}
		this->EchoSimulationData->saveEchoArr(echo_temp, sprfid, PRF_len);
	}





	/** 内存释放***************************************************/
	FreeCUDAHost(h_TantPattern);
	FreeCUDAHost(h_RantPattern);
	FreeCUDADevice(d_TantPattern);
	FreeCUDADevice(d_RantPattern);

	FreeCUDAHost(h_dem_x);
	FreeCUDAHost(h_dem_y);
	FreeCUDAHost(h_dem_z);
	FreeCUDAHost(h_demsloper_x);
	FreeCUDAHost(h_demsloper_y);
	FreeCUDAHost(h_demsloper_z);
	FreeCUDAHost(h_demsloper_z);
	FreeCUDAHost(h_demcls);
	FreeCUDAHost(h_echo_block_real);
	FreeCUDAHost(h_echo_block_imag);

	FreeCUDAHost(h_antpx);
	FreeCUDAHost(h_antpy);
	FreeCUDAHost(h_antpz);
	FreeCUDAHost(h_antvx);
	FreeCUDAHost(h_antvy);
	FreeCUDAHost(h_antvz);
	FreeCUDAHost(h_antdirectx);
	FreeCUDAHost(h_antdirecty);
	FreeCUDAHost(h_antdirectz);
	FreeCUDAHost(h_antXaxisX);
	FreeCUDAHost(h_antXaxisY);
	FreeCUDAHost(h_antXaxisZ);
	FreeCUDAHost(h_antYaxisX);
	FreeCUDAHost(h_antYaxisY);
	FreeCUDAHost(h_antYaxisZ);
	FreeCUDAHost(h_antZaxisX);
	FreeCUDAHost(h_antZaxisY);
	FreeCUDAHost(h_antZaxisZ);


	FreeCUDADevice(d_antpx);
	FreeCUDADevice(d_antpy);
	FreeCUDADevice(d_antpz);
	FreeCUDADevice(d_antvx);
	FreeCUDADevice(d_antvy);
	FreeCUDADevice(d_antvz);
	FreeCUDADevice(d_antdirectx);
	FreeCUDADevice(d_antdirecty);
	FreeCUDADevice(d_antdirectz);
	FreeCUDADevice(d_antXaxisX);
	FreeCUDADevice(d_antXaxisY);
	FreeCUDADevice(d_antXaxisZ);
	FreeCUDADevice(d_antYaxisX);
	FreeCUDADevice(d_antYaxisY);
	FreeCUDADevice(d_antYaxisZ);
	FreeCUDADevice(d_antZaxisX);
	FreeCUDADevice(d_antZaxisY);
	FreeCUDADevice(d_antZaxisZ);


	return ErrorCode::SUCCESS;
}