#include "stdafx.h"
#include "RTPCProcessCls.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>

 
#ifdef DEBUGSHOWDIALOG
#include "ImageShowDialogClass.h"


#endif

RTPCProcessCls::RTPCProcessCls()
{
	this->PluseCount = 0;
	this->PlusePoint = 0;
	this->TaskSetting = nullptr;
	this->EchoSimulationData = nullptr;
	this->DEMTiffPath="";
	this->LandCoverPath="";
	this->HHSigmaPath="";
	this->HVSigmaPath="";
	this->VHSigmaPath="";
	this->VVSigmaPath = "";
	this->OutEchoPath = "";

	this->DEMTiffPath.clear();
	this->LandCoverPath.clear();
	this->HHSigmaPath.clear();
	this->HVSigmaPath.clear();
	this->VHSigmaPath.clear();
	this->VVSigmaPath.clear();
	this->OutEchoPath.clear();

	this->SigmaDatabasePtr=std::shared_ptr<SigmaDatabase>(new SigmaDatabase);
 
}

RTPCProcessCls::~RTPCProcessCls()
{
}

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

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

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

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

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

void RTPCProcessCls::setHHSigmaPath(QString HHSigmaPath)
{
	this->HHSigmaPath = HHSigmaPath;
}

void RTPCProcessCls::setHVSigmaPath(QString HVSigmaPath)
{
	this->HVSigmaPath = HVSigmaPath;
}

void RTPCProcessCls::setVHSigmaPath(QString VHSigmaPath)
{
	this->VHSigmaPath = VHSigmaPath;
}

void RTPCProcessCls::setVVSigmaPath(QString VVSigmaPath)
{
	this->VVSigmaPath = VVSigmaPath;
}

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



ErrorCode RTPCProcessCls::Process(long num_thread)
{
	// RTPC 算法
	qDebug() << u8"params init ....";
	ErrorCode stateCode = this->InitParams();
	if (stateCode != ErrorCode::SUCCESS){
		return stateCode;
	}
	else {}
	qDebug() << "DEMMainProcess";
	stateCode = this->DEMPreprocess();
	if (stateCode != ErrorCode::SUCCESS) {
		return stateCode;
	}
	else {}
	qDebug() << "RTPCMainProcess";


	stateCode = this->RTPCMainProcess(num_thread);
	
	
	
	if (stateCode != ErrorCode::SUCCESS) {
		return stateCode;
	}else{}
 

	return ErrorCode::SUCCESS;
}

ErrorCode RTPCProcessCls::InitParams()
{
	if (nullptr==this->TaskSetting||this->DEMTiffPath.isEmpty() ||
		this->LandCoverPath.isEmpty() || this->HHSigmaPath.isEmpty() ||
		this->HVSigmaPath.isEmpty() || this->VHSigmaPath.isEmpty() || 
		this->VVSigmaPath.isEmpty()) {
		return ErrorCode::RTPC_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->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 RTPCProcessCls::DEMPreprocess()
{

	this->demxyzPath = QDir(tmpfolderPath).filePath("demxyz.tif");
	gdalImage demds(this->DEMTiffPath);
	gdalImage demxyz = CreategdalImage(demxyzPath, demds.height, demds.width, 3, demds.gt, demds.projection, true, true);// X,Y,Z

	// 分块计算并转换为XYZ

	Eigen::MatrixXd demArr = demds.getData(0, 0, demds.height, demds.width, 1);
	Eigen::MatrixXd demR = demArr;
	Landpoint LandP{ 0,0,0 };
	Point3 GERpoint{ 0,0,0 };
	double R = 0;
	double dem_row = 0, dem_col = 0, dem_alt = 0;

	long line_invert = 1000;

	double rowidx = 0;
	double colidx = 0;
	for (int max_rows_ids = 0; max_rows_ids < demds.height; max_rows_ids = max_rows_ids + line_invert) {
			Eigen::MatrixXd demdata = demds.getData(max_rows_ids, 0, line_invert, demds.width, 1);
			Eigen::MatrixXd xyzdata_x = demdata.array() * 0;
			Eigen::MatrixXd xyzdata_y = demdata.array() * 0;
			Eigen::MatrixXd xyzdata_z = demdata.array() * 0;

			int datarows = demdata.rows();
			int datacols = demdata.cols();

			for (int i = 0; i < datarows; i++) {
				for (int j = 0; j < datacols; j++) {
					rowidx = i + max_rows_ids;
					colidx = j;
					demds.getLandPoint(rowidx, colidx, demdata(i, j), LandP); // 获取地理坐标
					LLA2XYZ(LandP, GERpoint); // 经纬度转换为地心坐标系
					xyzdata_x(i, j) = GERpoint.x;
					xyzdata_y(i, j) = GERpoint.y;
					xyzdata_z(i, j) = GERpoint.z;
				}
			}
			demxyz.saveImage(xyzdata_x, max_rows_ids, 0, 1);
			demxyz.saveImage(xyzdata_y, max_rows_ids, 0, 2);
			demxyz.saveImage(xyzdata_z, max_rows_ids, 0, 3);
	}


	// 计算坡向角

	this->demsloperPath=QDir(tmpfolderPath).filePath("demsloper.tif");
	this->demmaskPath = QDir(tmpfolderPath).filePath("demmask.tif");

	gdalImage demsloperxyz = CreategdalImage(this->demsloperPath, demds.height, demds.width, 4, demds.gt, demds.projection, true, true);// X,Y,Z,cosangle
	gdalImage demmask = CreategdalImage(this->demmaskPath, demxyz.height, demxyz.width, 1, demxyz.gt, demxyz.projection, true, true);// X,Y,Z

	line_invert = 1000;
	long start_ids = 0;
	long dem_rows = 0, dem_cols = 0;

	for (start_ids = 1; start_ids < demds.height; start_ids = start_ids + line_invert) {
		Eigen::MatrixXd demdata = demds.getData(start_ids - 1, 0, line_invert + 2, demxyz.width, 1);
		long startlineid = start_ids;
		Eigen::MatrixXd maskdata = demmask.getData(start_ids - 1, 0, line_invert + 2, demxyz.width, 1);
		Eigen::MatrixXd demsloper_x = demsloperxyz.getData(start_ids - 1, 0, line_invert + 2, demxyz.width, 1);
		Eigen::MatrixXd demsloper_y = demsloperxyz.getData(start_ids - 1, 0, line_invert + 2, demxyz.width, 2);
		Eigen::MatrixXd demsloper_z = demsloperxyz.getData(start_ids - 1, 0, line_invert + 2, demxyz.width, 3);
		Eigen::MatrixXd demsloper_angle = demsloperxyz.getData(start_ids - 1, 0, line_invert + 2, demxyz.width, 4);

		maskdata = maskdata.array() * 0;
		Landpoint p0, p1, p2, p3, p4,pslopeVector,pp;
		Vector3D  slopeVector;

		dem_rows = maskdata.rows();
		dem_cols = maskdata.cols();
		double sloperAngle = 0;
		Vector3D Zaxis = { 0,0,1 };

		double rowidx = 0, colidx = 0;

		for (long i = 1; i < dem_rows - 1; i++) {
			for (long j = 1; j < dem_cols - 1; j++) {
				rowidx = i + startlineid;
				colidx = j;
				demds.getLandPoint(rowidx, colidx, demdata(i, j), p0);
				demds.getLandPoint(rowidx-1, colidx, demdata(i-1, j), p1);
				demds.getLandPoint(rowidx, colidx-1, demdata(i, j-1), p2);
				demds.getLandPoint(rowidx+1, colidx, demdata(i+1, j), p3);
				demds.getLandPoint(rowidx, colidx+1, demdata(i, j+1), p4);
 
				pslopeVector = getSlopeVector(p0, p1, p2, p3, p4); // 地面坡向矢量
				slopeVector = { pslopeVector.lon,pslopeVector.lat,pslopeVector.ati };
				pp = LLA2XYZ(p0);
				Zaxis.x = pp.lon;
				Zaxis.y = pp.lat;
				Zaxis.z = pp.ati;
				sloperAngle = getCosAngle(slopeVector, Zaxis) ; // 地面坡向角

				demsloper_x(i, j) = slopeVector.x;
				demsloper_y(i, j) = slopeVector.y;
				demsloper_z(i, j) = slopeVector.z;
				demsloper_angle(i, j) = sloperAngle;
				maskdata(i, j)++;
			}
		}
		demmask.saveImage(maskdata, start_ids - 1, 0, 1);
		demsloperxyz.saveImage(demsloper_x, start_ids - 1, 0, 1);
		demsloperxyz.saveImage(demsloper_y, start_ids - 1, 0, 2);
		demsloperxyz.saveImage(demsloper_z, start_ids - 1, 0, 3);
		demsloperxyz.saveImage(demsloper_angle, start_ids - 1, 0, 4);
	}

 
	return ErrorCode::SUCCESS;
}

 
ErrorCode RTPCProcessCls::RTPCMainProcess(long num_thread)
{
	omp_set_num_threads(num_thread);// 设置openmp 线程数量
	double widthSpace = LIGHTSPEED/2/this->TaskSetting->getFs();
 
	double prf_time = 0;
	double dt = 1 / this->TaskSetting->getPRF();// 获取每次脉冲的时间间隔
	bool antflag = true; // 计算天线方向图
	Landpoint LandP{ 0,0,0 };
	Point3 GERpoint{ 0,0,0 };
	double R = 0;
	double dem_row = 0 , dem_col=0, dem_alt=0;

	long double imageStarttime = 0;
	imageStarttime=this->TaskSetting->getSARImageStartTime();
	//std::vector<SatelliteOribtNode> sateOirbtNodes(this->PluseCount);
	std::shared_ptr<SatelliteOribtNode[]> sateOirbtNodes(new SatelliteOribtNode[this->PluseCount], delArrPtr);
	{ // 姿态计算不同
		// 计算姿态
		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 * 16 + 0] = prf_time + imageStarttime;
			antpos.get()[prf_id * 16 + 1] = sateOirbtNode.Px;
			antpos.get()[prf_id * 16 + 2] = sateOirbtNode.Py;
			antpos.get()[prf_id * 16 + 3] = sateOirbtNode.Pz;
			antpos.get()[prf_id * 16 + 4] = sateOirbtNode.Vx;
			antpos.get()[prf_id * 16 + 5] = sateOirbtNode.Vy;
			antpos.get()[prf_id * 16 + 6] = sateOirbtNode.Vz;
			antpos.get()[prf_id * 16 + 7] = sateOirbtNode.AntDirecX;
			antpos.get()[prf_id * 16 + 8] = sateOirbtNode.AntDirecY;
			antpos.get()[prf_id * 16 + 9] = sateOirbtNode.AntDirecZ;
			antpos.get()[prf_id * 16 + 10] = sateOirbtNode.AVx;
			antpos.get()[prf_id * 16 + 11] = sateOirbtNode.AVy;
			antpos.get()[prf_id * 16 + 12] = sateOirbtNode.AVz;
			antpos.get()[prf_id * 16 + 13] = outP.lon;
			antpos.get()[prf_id * 16 + 14] = outP.lat;
			antpos.get()[prf_id * 16 + 15] = outP.ati;

			sateOirbtNodes[prf_id] = sateOirbtNode;
		}
		this->EchoSimulationData->saveAntPos(antpos);
		antpos.reset();
		qDebug() << "Ant position finished sucessfully !!!";
	}
	
	
	
	// 回波
	long echoIdx = 0;

	gdalImage demxyz (demxyzPath );// 地面点坐标
	gdalImage demlandcls(this->LandCoverPath);// 地表覆盖类型
	gdalImage demsloperxyz(this->demsloperPath);// 地面坡向
 
	double NearRange = this->EchoSimulationData->getNearRange(); // 近斜据
	double FarRange = this->EchoSimulationData->getFarRange();

	double TimgNearRange = 2 * NearRange / LIGHTSPEED;
	double TimgFarRange = 2 * FarRange / LIGHTSPEED;

	double Fs = this->TaskSetting->getFs(); // 距离向采样率
	double Pt =  this->TaskSetting->getPt()*this->TaskSetting->getGri();// 发射电压 1v
	//double GainAntLen = -3;// -3dB 为天线半径
	long pluseCount = this->PluseCount;

	// 天线方向图
	std::shared_ptr<AbstractRadiationPattern> TransformPattern = this->TaskSetting->getTransformRadiationPattern(); // 发射天线方向图
	std::shared_ptr< AbstractRadiationPattern> ReceivePattern = this->TaskSetting->getReceiveRadiationPattern();   // 接收天线方向图
	std::shared_ptr<std::complex<double>> echo = this->EchoSimulationData->getEchoArr();
	long PlusePoint = this->EchoSimulationData->getPlusePoints();
	 
	// 初始化 为 0
	for (long i = 0; i < pluseCount * PlusePoint; i++) {
		echo.get()[i] =  std::complex<double>(0,0);
	}
	this->EchoSimulationData->saveEchoArr(echo, 0, PluseCount);

	POLARTYPEENUM polartype = this->TaskSetting->getPolarType();


	
	long demrowStep = std::ceil(Memory1MB*10 / 8 / demxyz.width) ;
	long line_invert = demrowStep > 3 ? demrowStep : 3;

	double lamda = this->TaskSetting->getCenterLamda(); // 波长

	omp_lock_t lock; // 定义锁
	omp_init_lock(&lock); // 初始化锁

	// DEM计算
	long start_ids = 1250;
	for (start_ids = 1; start_ids < demxyz.height; start_ids = start_ids + line_invert) {  // 8+ 17 + 0.3 MB
		QDateTime current = QDateTime::currentDateTime();
		long pluseStep = Memory1MB * 100 / 3 / PlusePoint;
		if (pluseStep * num_thread *3 > this->PluseCount) {
			pluseStep = this->PluseCount / num_thread /3;
		}
		pluseStep = pluseStep > 50 ? pluseStep : 50;


		qDebug()<< current.toString("yyyy-MM-dd HH:mm:ss.zzz") <<" \tstart \t "<< start_ids << " - "<< start_ids + line_invert <<"\t" << demxyz.height<<"\t pluseCount:\t"<< pluseStep;
		// 文件读取
		Eigen::MatrixXd dem_x = demxyz.getData(start_ids - 1, 0, line_invert + 1, demxyz.width, 1); //       
		Eigen::MatrixXd dem_y = demxyz.getData(start_ids - 1, 0, line_invert + 1, demxyz.width, 2); //     
		Eigen::MatrixXd dem_z = demxyz.getData(start_ids - 1, 0, line_invert + 1, demxyz.width, 3); //   

		// 地表覆盖
		std::shared_ptr<long> dem_landcls = readDataArr<long>(demlandcls, start_ids - 1, 0, line_invert + 1, demxyz.width, 1, GDALREADARRCOPYMETHOD::VARIABLEMETHOD); // 地表覆盖类型
		long* dem_landcls_ptr = dem_landcls.get();
		double localAngle = 30.0;

		bool sigmaNoZeroFlag = true;
		for (long ii = 0; ii < dem_x.rows(); ii++) {
			for (long jj = 0; jj < dem_y.cols(); jj++) {
				if (0 != this->SigmaDatabasePtr->getAmp(dem_landcls_ptr[dem_x.cols() * ii + jj], localAngle, polartype)) {
					sigmaNoZeroFlag = false;
					break;
				}
			}
			if (!sigmaNoZeroFlag) {
				break;
			}
		}
		if (sigmaNoZeroFlag) {
			continue;
		}

//#ifdef DEBUGSHOWDIALOG
//		dialog->load_double_MatrixX_data(dem_z, "dem_z");
//#endif


		Eigen::MatrixXd demsloper_x = demsloperxyz.getData(start_ids - 1, 0, line_invert + 1, demxyz.width, 1); //  
		Eigen::MatrixXd demsloper_y = demsloperxyz.getData(start_ids - 1, 0, line_invert + 1, demxyz.width, 2); //  
		Eigen::MatrixXd demsloper_z = demsloperxyz.getData(start_ids - 1, 0, line_invert + 1, demxyz.width, 3); //  
		Eigen::MatrixXd sloperAngle = demsloperxyz.getData(start_ids - 1, 0, line_invert + 1, demxyz.width, 4); //  
		sloperAngle = sloperAngle.array() * T180_PI;

		
		long dem_rows = dem_x.rows();
		long dem_cols = dem_x.cols();
	
		long freqidx = 0;// 



#ifdef DEBUGSHOWDIALOG
		ImageShowDialogClass* dialog = new ImageShowDialogClass(nullptr);
		dialog->show();
		
		Eigen::MatrixXd landaArr = Eigen::MatrixXd::Zero(dem_rows, dem_cols);
		for (long i = 0; i < dem_rows; i++) {
			for (long j = 0; j < dem_cols; j++) {
				landaArr(i, j) = dem_landcls.get()[i * dem_cols + j];
			}
		}
		dialog->load_double_MatrixX_data(landaArr, "landCover");
		
#endif
		//qDebug() << " pluse bolck size :\t " << pluseStep << " all size:\t" << this->PluseCount;
		

		#pragma omp parallel for
		for (long startprfidx = 0; startprfidx < this->PluseCount; startprfidx=startprfidx+ pluseStep) { // 17 + 0.3 MB
			long prfcount_step = startprfidx + pluseStep < this->PluseCount ? pluseStep : this->PluseCount- startprfidx;
			Eigen::MatrixXcd echoPluse = Eigen::MatrixXcd::Zero(prfcount_step, PlusePoint); // 当前脉冲的回波积分情况
			// 内存预分配
			Eigen::MatrixXd Rst_x = Eigen::MatrixXd::Zero(dem_x.rows(), dem_x.cols());   
			Eigen::MatrixXd Rst_y = Eigen::MatrixXd::Zero(dem_x.rows(), dem_x.cols());
			Eigen::MatrixXd Rst_z = Eigen::MatrixXd::Zero(dem_x.rows(), dem_x.cols());
			Eigen::MatrixXd R = Eigen::MatrixXd::Zero(dem_x.rows(), dem_x.cols());
			Eigen::MatrixXd localangle = Eigen::MatrixXd::Zero(dem_x.rows(), dem_x.cols());
			Eigen::MatrixXd Vst_x = Eigen::MatrixXd::Zero(dem_x.rows(), dem_x.cols());
			Eigen::MatrixXd Vst_y = Eigen::MatrixXd::Zero(dem_x.rows(), dem_x.cols());
			Eigen::MatrixXd Vst_z = Eigen::MatrixXd::Zero(dem_x.rows(), dem_x.cols());
			Eigen::MatrixXd fde = Eigen::MatrixXd::Zero(dem_x.rows(), dem_x.cols());
			Eigen::MatrixXd fr = Eigen::MatrixXd::Zero(dem_x.rows(), dem_x.cols());
			Eigen::MatrixXd Rx = Eigen::MatrixXd::Zero(dem_x.rows(), dem_x.cols());
			Eigen::MatrixXd sigam = Eigen::MatrixXd::Zero(dem_x.rows(), dem_x.cols());
			Eigen::MatrixXd echoAmp = Eigen::MatrixXd::Zero(dem_x.rows(), dem_x.cols()).array() + Pt;
			Eigen::MatrixXd Rphi = Eigen::MatrixXd::Zero(dem_x.rows(), dem_x.cols());
			Eigen::MatrixXd TimeRange = Eigen::MatrixXd::Zero(dem_x.rows(), dem_x.cols());
			Eigen::MatrixXd TransAnt = Eigen::MatrixXd::Zero(dem_x.rows(), dem_x.cols());
			Eigen::MatrixXd ReciveAnt = Eigen::MatrixXd::Zero(dem_x.rows(), dem_x.cols());

			Eigen::MatrixXd AntTheta = Eigen::MatrixXd::Zero(dem_x.rows(), dem_x.cols());
			Eigen::MatrixXd AntPhi = Eigen::MatrixXd::Zero(dem_x.rows(), dem_x.cols());

			double minR = 0, maxR = 0;
			double minLocalAngle = 0, maxLocalAngle = 0;

			Vector3D Rt = { 0,0,0 };
			SatelliteOribtNode oRs = SatelliteOribtNode{0};;

			Vector3D	p0 = {}, slopeVector = {}, sateAntDirect = {};
			Vector3D	Rs = {}, Vs = {}, Ast = {};
			SatelliteAntDirect antdirectNode = {};
			std::complex<double> echofreq;
			std::complex<double> Imag1(0, 1);
			double TAntPattern = 1;	// 发射天线方向图
			double RAntPanttern = 1;// 接收天线方向图
			double maxechoAmp = 1;
			double tempAmp = 1;
			for (long prfidx = 0; prfidx < prfcount_step; prfidx++)
			{
				oRs = sateOirbtNodes[prfidx + startprfidx];
				// 计算天线方向图
				for (long jj = 1; jj < dem_cols - 1; jj++) {
					for (long ii = 1; ii < dem_rows - 1; ii++) {
						p0.x = dem_x(ii, jj);
						p0.y = dem_y(ii, jj);
						p0.z = dem_z(ii, jj);
						this->TaskSetting->getSatelliteAntDirectNormal(oRs, p0, antdirectNode);
						//antdirectNode.ThetaAnt = antdirectNode.ThetaAnt * r2d;
						//antdirectNode.PhiAnt = antdirectNode.PhiAnt * r2d;
						AntTheta(ii, jj) = antdirectNode.ThetaAnt * r2d;
						AntPhi(ii, jj) = antdirectNode.PhiAnt * r2d;
					}
				}

				// 计算发射天线方向图
				for (long jj = 1; jj < dem_cols - 1; jj++) {
					for (long ii = 1; ii < dem_rows - 1; ii++) {
						TransformPattern->getGainLinear(AntTheta(ii, jj), AntPhi(ii, jj), TransAnt(ii, jj));
						//TransAnt(ii, jj) = TAntPattern;
					}
				}
 
				// 计算接收天线方向图
				for (long jj = 1; jj < dem_cols - 1; jj++) {
					for (long ii = 1; ii < dem_rows - 1; ii++) {
						TransformPattern->getGainLinear(AntTheta(ii, jj), AntPhi(ii, jj), ReciveAnt(ii, jj));
						//ReciveAnt(ii, jj) = RAntPanttern;
					}
				}


				// 计算经过增益的能量
				echoAmp = Pt * TransAnt.array() * ReciveAnt.array();

				maxechoAmp = echoAmp.maxCoeff();
				if (std::abs(maxechoAmp) < PRECISIONTOLERANCE) { // 这种情况下，不在合成孔径范围中
					continue;
				}

				Rs.x = sateOirbtNodes[prfidx+startprfidx].Px; // 卫星位置
				Rs.y = sateOirbtNodes[prfidx+startprfidx].Py;
				Rs.z = sateOirbtNodes[prfidx+startprfidx].Pz;

				Vs.x = sateOirbtNodes[prfidx+startprfidx].Vx; // 卫星速度
				Vs.y = sateOirbtNodes[prfidx+startprfidx].Vy;
				Vs.z = sateOirbtNodes[prfidx+startprfidx].Vz;

				Ast.x = sateOirbtNodes[prfidx+startprfidx].AVx;// 卫星加速度
				Ast.y = sateOirbtNodes[prfidx+startprfidx].AVy;
				Ast.z = sateOirbtNodes[prfidx+startprfidx].AVz;

				Rst_x = Rs.x - dem_x.array(); //  Rst = Rs - Rt;
				Rst_y = Rs.y - dem_y.array();
				Rst_z = Rs.z - dem_z.array();
				R	  = (Rst_x.array().pow(2) + Rst_y.array().pow(2) + Rst_z.array().pow(2)).array().sqrt(); // R

				minR = R.minCoeff();
				maxR = R.maxCoeff();
				//qDebug() << "minR:\t" << minR << "  maxR:\t" << maxR;
				if (maxR<NearRange || minR>FarRange) {
					continue;
				}
				else {}

				// getCosAngle
				// double c = dot(a, b) / (getlength(a) * getlength(b));
				// return acos(c > 1 ? 1 : c < -1 ? -1 : c) * r2d;
				// localangle = getCosAngle(Rst, slopeVector) * T180_PI; // 注意这个只能实时计算，因为非实时计算代价太大
				localangle = (Rst_x.array() * demsloper_x.array() + Rst_y.array() * demsloper_y.array() + Rst_z.array() * demsloper_z.array()).array(); // dot(a, b)
				localangle = localangle.array() / R.array();
				localangle = localangle.array() / (demsloper_x.array().pow(2) + demsloper_y.array().pow(2) + demsloper_z.array().pow(2)).array().sqrt().array();
				localangle = localangle.array().acos(); // 弧度值 
	 
				minLocalAngle = localangle.minCoeff();
				maxLocalAngle = localangle.maxCoeff();
				if (maxLocalAngle<0 || minLocalAngle>PI/2) {
					continue;
				}
				else {}

				//Vst_x = Vs.x + 1 * earthRoute * dem_y.array(); // Vst = Vs - Vt;
				//Vst_y = Vs.y - 1 * earthRoute * dem_x.array();
				//Vst_z = Vs.z - Eigen::MatrixXd::Zero(dem_x.rows(), dem_y.cols()).array();

				//// 计算多普勒中心频率 Rst, Vst : ( - 2 / lamda) * dot(Rs - Rt, Vs - Vt) / R; // 星载合成孔径雷达原始回波数据模拟研究 3.18
				//fde = (-2 / lamda) * (Rst_x.array() * Vst_x.array() + Rst_y.array() * Vst_y.array() + Rst_z.array() * Vst_z.array()).array() / (R.array());
				//// 计算多普勒频率斜率   // 星载合成孔径雷达原始回波数据模拟研究 3.19
				//// -(2/lamda)*(  dot(Vs - Vt, Vs - Vt)/R  +   dot(Ast, Rs - Rt)/R  - std::pow(dot(Vs - Vt, Rs - Rt),2 )/std::pow(R,3));
				//fr = (-2 / lamda) *
				//	(Vst_x.array() * Vst_x.array() + Vst_y.array() * Vst_y.array() + Vst_z.array() * Vst_z.array()).array() / (R.array()) +
				//	(-2 / lamda) *
				//	(Ast.x * Rst_x.array() + Ast.y * Rst_y.array() + Ast.z * Rst_z.array()).array() / (R.array()) -
				//	(-2 / lamda) *
				//	(Vst_x.array() * Rst_x.array() + Vst_y.array() * Rst_y.array() + Vst_z.array() * Rst_z.array()).array().pow(2) / (R.array().pow(3));
				// 计算回波
				Rx = R;//  -(lamda / 2) * (fde * TRx + 0.5 * fr * TRx * TRx); // 斜距历程值 

				
				// 逐点计算 this->SigmaDatabasePtr->getAmp(covercls, localangle, polartype); // 后向散射系数 HH
			 
				for (long ii = 0; ii < dem_x.rows(); ii++) {
					for (long jj = 0; jj < dem_y.cols(); jj++) {
						sigam(ii, jj) = this->SigmaDatabasePtr->getAmp(dem_landcls_ptr[dem_x.cols() * ii + jj], localangle(ii, jj)*r2d, polartype);
					}
				}

				if (sigam.maxCoeff() > 0) {}
				else {
					continue;
				}
				// projArea = 1 / std::cos(sloperAngle) * std::sin(localangle); // 投影面积系数，单位投影面积 1m x 1m --注意这里是假设，后期再补充
				// echoAmp =  projArea*TAntPattern * RAntPanttern * sigam  / (4 * PI * R * R);
				
				echoAmp = echoAmp.array()*sigam.array() * (1 / sloperAngle.array().cos() * localangle.array().sin());  // 反射强度
				echoAmp = echoAmp.array() / (4 * PI * R.array().pow(2));				// 距离衰减
							
				Rphi = -4 * PI / lamda * Rx.array();// 距离徙动相位
				// 积分
				TimeRange = ((2 * R.array() / LIGHTSPEED - TimgNearRange).array() * Fs).array();
				double localAnglepoint = -1;
				long prf_freq_id = 0;
 
				
				for (long jj = 1; jj < dem_cols - 1; jj++) {
					for (long ii = 1; ii < dem_rows - 1; ii++) {
						prf_freq_id =std::floor(TimeRange(ii, jj));
						if (prf_freq_id < 0 || prf_freq_id >= PlusePoint|| localangle(ii, jj) <0 || localangle(ii, jj) >PI / 2|| echoAmp(ii, jj)==0) {
							continue;
						}
						echofreq = echoAmp(ii, jj) * std::exp( Rphi(ii, jj)* Imag1);
						echoPluse(prfidx, prf_freq_id) = echoPluse(prfidx, prf_freq_id) + echofreq;
					}
				}

#ifdef DEBUGSHOWDIALOG
				ImageShowDialogClass* localangledialog = new ImageShowDialogClass(dialog);
				localangledialog->show();
				localangledialog->load_double_MatrixX_data(localangle.array()*r2d, "localangle");


				ImageShowDialogClass* sigamdialog = new ImageShowDialogClass(dialog);
				sigamdialog->show();
				sigamdialog->load_double_MatrixX_data(TimeRange, "TimeRange");


				ImageShowDialogClass* ampdialog = new ImageShowDialogClass(dialog);
				ampdialog->show();
				ampdialog->load_double_MatrixX_data(echoAmp, "echoAmp");

				Eigen::MatrixXd echoPluseamp = echoPluse.array().abs().cast<double>().array();
				ImageShowDialogClass* echoampdialog = new ImageShowDialogClass(dialog);
				echoampdialog->show();
				echoampdialog->load_double_MatrixX_data(echoPluseamp, "echoPluseamp");


				dialog->exec();
#endif
				//qDebug() << QDateTime::currentDateTime().toString("yyyy-MM-dd HH:mm:ss.zzz") << " end  " << prfidx;
			}
			//qDebug() << QDateTime::currentDateTime().toString("yyyy-MM-dd HH:mm:ss.zzz")<<" step "<< prfcount_step;
			
			omp_set_lock(&lock); // 回波整体赋值处理
			for (long prfidx = 0; prfidx < prfcount_step; prfidx++) {
				for (long freqidx = 0; freqidx < PlusePoint; freqidx++)
				{
					//qDebug() << prfidx << "  " << freqidx << "  " << echoPluse(prfidx, freqidx).real() << " + " << echoPluse(prfidx, freqidx).imag() << "  j";
					echo.get()[(prfidx + startprfidx) * PlusePoint + freqidx] = echo.get()[(prfidx + startprfidx) * PlusePoint + freqidx] + echoPluse(prfidx, freqidx);
				}
			}
			//this->EchoSimulationData->saveEchoArr(echo, 0, PluseCount);
			omp_unset_lock(&lock); // 解锁
			//qDebug() << QDateTime::currentDateTime().toString("yyyy-MM-dd HH:mm:ss.zzz") << " step 2" << prfcount_step;
		}
		
		omp_set_lock(&lock); // 保存文件
		this->EchoSimulationData->saveEchoArr(echo ,0, PluseCount);
		omp_unset_lock(&lock); // 解锁

		qDebug() << QDateTime::currentDateTime().toString("yyyy-MM-dd HH:mm:ss.zzz") << "  \t " << start_ids << "\t--\t " << start_ids + line_invert<<"\t/\t" << demxyz.height;
	}
	omp_destroy_lock(&lock); // 销毁锁
	this->EchoSimulationData->saveEchoArr(echo, 0, PluseCount);
	return ErrorCode::SUCCESS;
}

void RTPCProcessMain(long num_thread,QString TansformPatternFilePath, QString ReceivePatternFilePath, QString simulationtaskName,QString OutEchoPath,QString GPSXmlPath, QString TaskXmlPath, QString demTiffPath   ,	QString LandCoverPath ,	QString HHSigmaPath	  ,	QString HVSigmaPath	  ,	QString VHSigmaPath	  ,	QString VVSigmaPath)
{
	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);

	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() << "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 设置天线方向图
	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() << "--------------  RTPC init ---------------------------------------";
	RTPCProcessCls rtpc;
	rtpc.setTaskSetting(task);						//qDebug() << "setTaskSetting";
	rtpc.setTaskFileName(simulationtaskName);		//qDebug() << "setTaskFileName";
	rtpc.setDEMTiffPath(demTiffPath);				//qDebug() << "setDEMTiffPath";
	rtpc.setLandCoverPath( LandCoverPath);			//qDebug() << "setLandCoverPath";
	rtpc.setHHSigmaPath(   HHSigmaPath  );			//qDebug() << "setHHSigmaPath";
	rtpc.setHVSigmaPath(   HVSigmaPath  );			//qDebug() << "setHVSigmaPath";
	rtpc.setVHSigmaPath(   VHSigmaPath  );			//qDebug() << "setVHSigmaPath";
	rtpc.setVVSigmaPath(   VVSigmaPath   );			//qDebug() << "setVVSigmaPath";
	rtpc.setOutEchoPath(OutEchoPath);				//qDebug() << "setOutEchoPath";
	qDebug() << "--------------  RTPC start---------------------------------------";
	rtpc.Process(num_thread); // 处理程序
	qDebug() << "--------------  RTPC end---------------------------------------";
}