#include <stdio.h>
#include <stdlib.h>
#include<time.h>
#include<math.h>

#include "TankPrediction.h" 

#pragma diag_suppress 1035,177,550,940

//cqr 20230224修改

//算法系数罐箱参数组分参数不在云端下发，直接设置参数 jzw 20230324
void SetTankPara(struct TankPara *tank)
{
	/////第2类：一次性下发的算法固定常数，按后面的注释值赋值，所有常数变量在苏醒后需要读入
	tank->DL0 = 5; //		0.05  =====液位差判断基准==4G输入给出，
	tank->Dis_0 = 100; //  1000 
	tank->DT_0 = 1; //采样时间周期ΔT  整型常数变量 1 小时 	
	tank->T_0 = 150; //		150  =====温度判断基准==4G输入给出，
	tank->Stable_Day = 1; //		1 ====静置天数==4G输入给出，
	tank->Lamda_0 = 0.1; //   0.1 
	tank->L_0 = 0.1; //充装状态判断条件Φ0 0.1 =====液位判断基准==4G输入给出， 
	tank->Eta_V = 0.95; //		0.95  =====容积安全系数==4G输入给出，
	tank->Pa_0 = 0.1013; //   0.1013 
	tank->Tc_0 = 303.15; //  =273.15+30 = 303.15 
	tank->R_0 = 8.314; //   8.314 
	tank->MTime_0 = 300; //   300
	tank->Fai_bz0 = 95; //   0.95 
	tank->Nday_0 = 20; // 整型常数 20  平均日漏热量计算周期CT
	tank->Tz_0 = 24; // 内能计算间隔周期tz  整型常数 24小时
	tank->Fai_0 = 75; //		0.75  =====充装率判断基准==4G输入给出，
	tank->Mavg_0 = 16; //   16
	tank->Pcr_0 = 4.59; // 常数  4.59
	tank->QHQR_0 = 511.12; //	511.12  =====标态LNG气化潜热==4G输入给出，
	tank->RhoLng_0 = 434.15; //   434.15 
	tank->K_1 = 0.9553; //		0.6583  =====标态LNG和液氮气化潜热转换系数=4G输入给出，

	tank->Default_CH4 = 96.299; //	0.95 ====默认甲烷成分==4G输入给出，
	tank->VN_0 = 111.91; //常数  0.01191
	tank->Nx_0 = 0.35; //   0.35
	tank->Ny_0 = 0.5; //   0.5
	tank->Nz_0 = 0.15; //   0.15
	tank->Kymax = 1.1; //常数 1.1
	tank->Kymin = 0.8; //常数 0.8

	tank->Da_0 = 20.2; //   0.00202 
	tank->Db_0 = 230; //   -0.023
	tank->Dc_0 = 125; //   0.0125 
	tank->Dd_0 = 309; //   -0.0309

	//预留
	tank->Default_C2 = 2.585; // 默认C2成分比
	tank->Default_C3 = 0.489; // 默认C3成分比
	tank->Default_C4 = 0.227; // 默认C4成分比
	tank->Default_N2 = 0.4; // 默认N2成分比
	tank->TimeYJ = 5; // 安全阀起跳预警值
	tank->RJJ = 0; //结果相关性修正系数—截距
	tank->RXL = 1; // 结果相关性修正系数—斜率
	tank->RQP = 0; // 结果相关性修正区间起始压力
	tank->QZP = 0.651; // 结果相关性修正区间终止压力
	tank->ZQQ = 0.03024; // 重罐轻罐权重系数-气相导热系数
	tank->ZQY = 0.22; // 重罐轻罐权重系数-液相导热系数
	tank->ZQFai0 = 5; // 重罐轻罐权重修正起始充装率
	tank->ZQFaiN = 95; // 重罐轻罐权重修正终止充装率
	tank->SFKZ1 = 2; // 算法控制参数1(是否启用区间结果相关性修正）
	tank->SFKZ2 = 1; // 算法控制参数2（重罐-轻罐权重系数）
	tank->SFKZ3 = 1; // 算法控制参数3（临界参数你和修正）
	tank->SFKZ4 = 1; // 算法控制参数4（温度系数K2修正）
	tank->SFKZ5 = 1; // 算法控制参数5（日漏热量修正）
	tank->SFKZ6 = 1; // 算法控制参数6（罐箱用途）

	tank->Tcr = 89; //临界温度
	tank->dTmax = 5; //罐内最大温差
	tank->dTP = 0.5; //安全阀起跳状态判断参数

	tank->csY3 = 0; // 常数预留变量
	tank->csY4 = 0;
	tank->csY5 = 0;
	tank->csY6 = 0;
	tank->csY7 = 0;
	tank->csY8 = 0;
	tank->csY9 = 0;
	tank->csY10 = 0;

	/////

	//////////////////////////////////////////////////////////////////////////

	////第3类 罐箱固定参数
	tank->P75 = 0.75; //满罐安全阀起跳压力，文件输入
	tank->a20 = 0.167; //液氮静态蒸发率，文件输入
	tank->TankVolYX = 39.24; // 有效容积，文件输入
	tank->TankVol = 43.6; // 总容积，文件输入
	tank->D = 2.21; // 罐子直径
	tank->fai = 0.82; // 额定充装率，文件输入	
	tank->GT = 1; //罐箱用途 运输罐箱为01，燃料罐箱为02
	tank->TankType = 1; // 罐箱绝热型式，文件输入
	tank->ZC = 2; //罐箱支撑型式  01代表两点支撑，02代表8点支撑，02为默认值
	tank->ZKD = 0.003; //真空度
	tank->Phz = 0.551; // 安全阀回座压力
	tank->minFai = 0.001001; //最小可识别充装率
	tank->RL = 11.744; // 罐箱圆柱段长度

	tank->gxY3 = 0; // 预留参数
	tank->gxY4 = 0;
	tank->gxY5 = 0;
	tank->gxY6 = 0;
	tank->gxY7 = 0;
	tank->gxY8 = 0;
	tank->gxY9 = 0;
	tank->gxY10 = 0;

	////////// 第5类 组分航次参数，营运输入，每航次变更
	tank->ch4 = 96.299; // 甲烷成分比
	tank->C2 = 2.585; // C2摩尔百分比（体积比）
	tank->C3 = 0.489; // C3摩尔百分比（体积比）
	tank->C4 = 0.227; // C4摩尔百分比（体积比）
	tank->N2 = 0.4; // N2摩尔百分比（体积比）
	tank->hcY1 = 0; //航次预留系数
	tank->hcY2 = 0;
	tank->hcY3 = 0;
	tank->hcY4 = 0;
	tank->hcY5 = 0;
	tank->hcY6 = 0;
	tank->hcY7 = 0;
	tank->hcY8 = 0;
	tank->hcY9 = 0;
	tank->hcY10 = 0;

	tank->sfbm=230401; //算法版本编码 00010001 0001算法版本 0001系数版本
}

//初始化算法基本常量参数
//在程序最初开始运行时调用此函数，进行算法基本常量的赋值，如遇程序异常重启，导致存储计算数据清零，也需要重新调用此函数进行初始化
void InitConstantPara(struct TankPara *tank)
{		
	// 应对以下参数赋初值
	tank->DT_0=0;

	tank->Nday_0=0;

	tank->Pcr_0=0;

	tank->VN_0=0;

	tank->Kymax=0;

	tank->Kymin=0;
	
	tank->Tz_0=0;
	
	tank->T_0=0;
	
	tank->L_0=0;
	
	tank->Fai_0=0;
	
	tank->DL0=0;
	
	tank->Eta_V=0;
	
	tank->QHQR_0=0;
	
	tank->K_1=0;
	
	tank->Stable_Day=0;
	
	tank->Default_CH4=0;
	
	tank->Dis_0=0;
	
	tank->Lamda_0=0;
	
	tank->Pa_0=0;
	
	tank->MTime_0=0;
	
	tank->R_0=0;
	
	tank->Mavg_0=0;
	
	tank->Tc_0=0;
	
	tank->RhoLng_0=0;
	
	tank->Fai_bz0=0;
	
	tank->Da_0=0;
	
	tank->Db_0=0;
	
	tank->Dc_0=0;
	
	tank->Dd_0=0;
	
	tank->Nx_0=0;
	
	tank->Ny_0=0;
	
	tank->Nz_0=0;

	tank->Default_C2=0;

	tank->Default_C3 = 0;

	tank->Default_C4 = 0;

	tank->Default_N2 = 0;

	tank->TimeYJ=0;

	tank->ZQQ=0;

	tank->ZQY=0;

	tank->Tcr=0;

	tank->dTmax=0;

	tank->dTP=0;
	
	tank->SFKZ1 = 0;
	
	tank->RQP = 0;
	
	tank->QZP = 0;

	tank->RJJ = 0;
	
	tank->RXL = 0;
	
	tank->ZQFai0 = 0;
	
	tank->ZQFaiN = 0;
	
//
//#define  T_0 150  // 温度判断基准
//#define  L_0 0.75  // 液位判断基准
//#define  Fai_0 0.75  //充装率判断基准
//#define  DL0 5  //充装率差基准，单位百分比，判断是否充装或使用
//#define  Eta_V 0.95  //容积安全系数，用于计算最小液体密度时的容积折减
//#define  QHQR_0 511.12  //标态LNG气化潜热
//#define  K_1 0.6583  //标态LNG和液氮气化潜热的转换系数
//#define  Stable_Day 1 // 静置天数
//#define  Default_CH4 95 // 默认甲烷成分比
//#define  Dis_0 1000 // 默认移动距离
//#define  Lamda_0 0.5 // 默认安全系数
//#define  Pa_0 0.1013 // 标准大气压
//
//#define  MTime_0 300 //安全预报时间
//#define  R_0 8.314 //理想气体常数
//#define  Mavg_0 16 //摩尔质量
//#define  Tc_0 293.15 //常温
//#define  RhoLng_0 434.15 //标准LNG密度
//#define  Fai_bz0 95 //标准充装率
//
//#define  Da_0 20.2 //运动模型参数Da
//#define  Db_0 230 //运动模型参数Da
//#define  Dc_0 125 //运动模型参数Da
//#define  Dd_0 309 //运动模型参数Da
//
//#define  Nx_0 0.35 //长度方向加权系数
//#define  Ny_0 0.5 //宽度方向加权系数
//#define  Nz_0 0.15 //高度方向加权系数

}

//初始化罐箱固有参数
void InitTankPara(struct TankPara *tank)
{
	
	  //应根据罐箱参数赋值
	tank->TankType=0;  //罐箱类型
	
	tank->ZC=0; //支撑型式
	
	tank->TankVol=0; // 总容积
	
	tank->TankVolYX=0; //有效容积
	
	tank->D=0; //直径
	
	tank->a20=0;//液氮静态蒸发率
	
	tank->P75=0;//起跳压力
	
	tank->GT=0; //罐箱用途

	tank->QDay_0=0; //默认日漏热量
	
	tank->Phz=0;//安全阀回座压力

	tank->minFai=0; // 最低可识别充装率

	tank->RL=0;// 罐箱圆柱段长度































}

//初始化罐箱每次充装参数
//每次罐箱重新充装后，充装率或甲烷成分比变化，需要调用此函数重新初始化
void InitLoadPara(struct TankPara *tank)
{
	///////// 应根据每次充装数据修改
	tank->fai=0; //充装率

	tank->ch4=0;//甲烷成分比

	//////////////////////////////////

	if (tank->ch4<0)
	{
	
	}

	//计算初始化，以下不用改
	tank->isStart=0;
	tank->MaintainTime=0;

	tank->ZeroDate = -1;
	tank->ZeroL =0; 
	tank->ForeCast_Sec = -1;
	tank->ForeCast_Date = -1; 
	tank->nData=0;
	tank->nday=0;  

	tank->iday=0;
	tank->ihour=0;
	tank->iQday=0;
	tank->QDay_0=0;

	tank->GT=1;
	tank->CY=0;
	
	tank->breakDate = -1;
	tank->isCal = 0;
	tank->isChange = 0;
	tank->isBreak = 0;
	tank->iMbreak = 0;
	tank->iMhour = 0;
	////////////////

}


//维持时间计算主函数，每次采集到罐箱状态（液位、压力、温度等）新数据后，调用此函数进行计算
//目前接口：Mdata存储每次采集的罐箱状态数据，tank存储算法和罐箱的基本参数，以及计算结果
//tank数据每次计算完之后需要保存，并且再给到下一次计算
void PredictTime( struct MonitoringData * MData, struct TankPara * tank ) 
{
	short int i,k,isStart;
	double ZeroDate;
	float MaintainTime,ZeroL, ForeCast_Sec, ForeCast_Date, breakSec;
	double Record_Date=MData->Record_Date;
	
	struct tm ZeroDate_Str;
	struct tm NowDate_Str;

	struct tm BreakDate_Str;

	struct tm PreDate_Str;

	float StablePeriod = 1 * 24* 3600; 
	struct PhysicochemicalPara  PCi[1];

	float chongzhuangbhl,dt,weizhibhl,lam;
	float totalM = 0,mi=0;  
	float m0Max=0,m0Min=0; 
	

	SetTankPara(tank);
	StablePeriod = tank->Stable_Day * 24* 3600; 
	tank->hcY1=1;
	if (abs(MData->WE)>0.0001&&abs(MData->SN)>0.0001)
	{
		//异常
		tank->hcY1=0;
	}
	

	//jzw 20230227 采集温度为℃，转换为k 测试数据单位为k
	MData->Record_OutT =tank->Tc_0;// cqr-20230420  安捷汇不采集环境温度    MData->Record_OutT + 273.15;  //cqr
	MData->Record_Tl = MData->Record_Tl + 273.15;
	if(MData->Record_Tg < -200 )
	{
		MData->Record_Tg = MData->Record_Tl;
		MData->Record_T = MData->Record_Tl;
	}
	else
	{
		MData->Record_Tg = MData->Record_Tg + 273.15;
		MData->Record_T = 0.5*(MData->Record_Tg + MData->Record_Tl);
	}

	MData->Record_L = MData->Record_L * 0.01;
	tank->SQdLng2 = tank->QDay_0;
	tank->fai=MData->Record_L; 
	//for(i=0;i<50;i++)
	{
	//	tank->RZ[i]=0;
	}


	tank->ztbm = 0;
	tank->ztbm += 1 * pow(10, 10);

	if (tank->ZeroDate<1) 
	{ 
		tank->ZeroDate =Record_Date;
		tank->ZeroL=MData->Record_L;

		ConvertDateFormat( Record_Date, &ZeroDate_Str);
		tank->ForeCast_Sec = 0;
		tank->ForeCast_Date = 0; 
		//tank->DisplayInfo = 1;
		//tank->CalModel=-1; 
		tank->Dis=0;

		tank->nday=0;  
		tank->iday=0;
		tank->ihour=0;
		tank->iQday=0;

		if (tank->hcY1 >0.5 )
		{
			tank->WE0=MData->WE; 
			tank->SN0=MData->SN;
		}
		

		tank->Pre_fai=MData->Record_L; 
		tank->preDate=Record_Date;

		tank->isChange=0; 
		tank->isCal=1;  
		tank->iMhour = 0; 
		tank->iMbreak = 0;		
		tank->isBreak = 0;
		tank->CZ = 0; //充装状态

	}
	else
	{
		
		ConvertDateFormat(tank->preDate , &PreDate_Str);
		ConvertDateFormat( Record_Date, &NowDate_Str);
		
		dt = (mktime(&NowDate_Str) - mktime(&PreDate_Str)) / 60; 

		if (tank->hcY1 >0.5 )
		{
			tank->Dis=getDistanceFromGPS(tank->WE0,tank->SN0,MData->WE,MData->SN);
			tank->WE0=MData->WE; 
			tank->SN0=MData->SN;
			
		}
		
		if (dt>5)  //5分钟
		{
			chongzhuangbhl=fabs(MData->Record_L - tank->Pre_fai )/(dt/60);
			if (tank->hcY1 >0.5 )
				weizhibhl=tank->Dis/(dt/60);

		}
		else
		{
			chongzhuangbhl=0;
			weizhibhl=0;
		}

		

		tank->Pre_fai=MData->Record_L; 
		tank->preDate=Record_Date;
		
		if (tank->GT==1)
		{
			if (chongzhuangbhl< tank->DL0*0.01)
			{
				if (weizhibhl < tank->Dis_0)
				{
					tank->CY=1;
				}
				else
					tank->CY=2;				
			}
			else
			{
				if (MData->Record_L > tank->Pre_fai )
				{
					tank->CY=3;	
				}
				else if (MData->Record_L < tank->Pre_fai && weizhibhl < tank->Dis_0 )
				{
					tank->CY=4;	
				}
				else
				{
					tank->CY=5;
					tank->ztbm += 3 * pow(10, 10); 
					tank->GJ = 1;
				}

				tank->isChange=0; 
				tank->isCal=1; 
				tank->iMhour = 0; 
			}
		}

		tank->ztbm += tank->CY * pow(10, 13); 
		
		if (tank->CY>2) 
		{
			ConvertDateFormat(tank->ZeroDate , &ZeroDate_Str);

			ConvertDateFormat( Record_Date, &NowDate_Str);
			
			
			tank->ForeCast_Sec = (mktime(&NowDate_Str) - mktime(&ZeroDate_Str));
			tank->ForeCast_Date = tank->ForeCast_Sec / 86400;

			
			tank->ZeroDate =-1;
			tank->ForeCast_Date = -1;
			tank->MaintainTime=-1;
			//tank->DisplayInfo = 1;
			tank->isStart=0;
			//tank->CalModel=-1; 


			
			tank->nday=0;  
			tank->iday=0;
			tank->ihour=0;
			tank->iQday=0;

			for (i=0;i<20;i++)
			{
				tank->Qday[i]=0;
				tank->day[i]=0;

				tank->Qhour[i]=0;
				tank->hour[i]=0;
			}
			for (i=20;i<24;i++)
			{
				tank->Qhour[i]=0;
				tank->hour[i]=0;
			}
			

			tank->Dis=0;

			
			
			
			
			
			

			tank->ZeroL = MData->Record_L;  
		}
		else   
		{
			
			ConvertDateFormat(tank->ZeroDate , &ZeroDate_Str);

			ConvertDateFormat( Record_Date, &NowDate_Str);
			
			
			tank->ForeCast_Sec = (mktime(&NowDate_Str) - mktime(&ZeroDate_Str));
			tank->ForeCast_Date = tank->ForeCast_Sec / 86400;

			

			
			

			
			
			
			
			
			
		}

		
	}
	
	if(tank->isChange == 0)
	{
		tank->m0 = Calculate_m(tank,MData);

		
		if(tank->ForeCast_Sec > StablePeriod) 
			tank->isChange =1;
	}


	if(tank->isCal == 1 && tank->CY == 1) 
	{
		
		tank->isChange =1;

		//if(tank->ForeCast_Sec < 3600 || tank->iMhour < 4) //1h内或数组数量不够方差计算时，进行质量实时计算
		if(tank->ForeCast_Sec <= StablePeriod)// cqr-20230420 基准静置时间内
		{
			mi=Calculate_m(tank,MData);
			if(tank->iMhour >= 12)
			{
				for (i=0;i<tank->iMhour-1;i++)
				{
					tank->Mhour[i]=tank->Mhour[i+1];
				}
				tank->Mhour[tank->iMhour-1] = mi;
			}
			else
			{
				tank->Mhour[tank->iMhour] =  mi;
				tank->iMhour = tank->iMhour + 1;
			}
			
			 
			if (tank->iMhour >= 4)
			{
			
				totalM = 0; 
				m0Max=0;
				m0Min= tank->Mhour[0]; 
				for (i=0;i<tank->iMhour;i++)
				{
					totalM = totalM + tank->Mhour[i];
					if(m0Max < tank->Mhour[i])
					{
						m0Max = tank->Mhour[i];
					}
					if(m0Min > tank->Mhour[i])
					{
						m0Min = tank->Mhour[i];
					}
				}
				tank->m0 = (totalM-m0Max-m0Min)/(tank->iMhour-2);
			}
			else
			{
				tank->m0 = tank->Mhour[tank->iMhour-1];
			}
		
			
		}
		else 
		{
			if(tank->iMhour >= 4) 
			{
				tank->isCal=0; 
				tank->iMhour = 0;
			}
			else
			{
				tank->Mhour[tank->iMhour] = Calculate_m(tank,MData);
				tank->iMhour = tank->iMhour + 1;
				
				if (tank->iMhour >= 4)
				{
					totalM = 0; mi=0; 
					m0Max=0;
					m0Min= tank->Mhour[0]; 
					for (i=0;i<tank->iMhour;i++)
					{
						totalM = totalM + tank->Mhour[i];
						if(m0Max < tank->Mhour[i])
						{
							m0Max = tank->Mhour[i];
						}
						if(m0Min > tank->Mhour[i])
						{
							m0Min = tank->Mhour[i];
						}
					}
					tank->m0 = (totalM-m0Max-m0Min)/(tank->iMhour-2);
				}
				else
				{
					tank->m0 = tank->Mhour[tank->iMhour-1];
				}
			}

		
		
		}
	
	}

	if(tank->isBreak == 1)
	{
		ConvertDateFormat(tank->breakDate , &BreakDate_Str);
		ConvertDateFormat( Record_Date, &NowDate_Str);
		breakSec = (mktime(&NowDate_Str) - mktime(&BreakDate_Str));
	}

	if(tank->isBreak == 1)
	{
		if(breakSec <= 3600 || tank->iMbreak < 4) 
		{
			tank->Mbreak[tank->iMbreak] = Calculate_m(tank,MData);
			tank->iMbreak = tank->iMbreak + 1;
		}
		else 
		{
			if(breakSec <= 3600*4) 
			{
				totalM = 0;mi = 0;
				m0Max=0;
				m0Min= tank->Mbreak[0]; 
				for (i=0;i<tank->iMbreak;i++)
				{
					totalM = totalM + tank->Mbreak[i];
					if(m0Max < tank->Mbreak[i])
					{
						m0Max = tank->Mbreak[i];
					}
					if(m0Min > tank->Mbreak[i])
					{
						m0Min = tank->Mbreak[i];
					}
				}
				mi = (totalM-m0Max-m0Min)/(tank->iMbreak-2);
				if(tank->iMbreak >= 12)
				{
					for (i=0;i<tank->iMbreak-1;i++) 
					{
						tank->Mbreak[i]=tank->Mbreak[i+1];
					}
					tank->Mhour[tank->iMbreak-1] = mi;
				}
				else
				{
					tank->Mbreak[tank->iMbreak] = mi;
					tank->iMbreak = tank->iMbreak + 1;
				}
				
			}
		}
		tank->m0 = tank->Mbreak[tank->iMbreak-1];
	}

	if(tank->isBreak == 1 && breakSec > 4*3600)
	{
		tank->iMbreak = 0;
		for (i=0;i<12;i++)
		{
			tank->Mbreak[i]=0;
		}
		tank->isBreak = 0;
	}
	if((tank->Pre_P - MData->Record_P) >= (tank->P75-tank->Phz)*tank->dTP) //Pi-1 - Pi >= (Pn-Pn1)*dtP 时判断起跳
	{
		tank->breakDate = Record_Date;
		tank->isBreak = 1;
		tank->isChange = 1;
	}
	tank->Pre_P = MData->Record_P; 

	
	



	if (MData->Record_T< tank->T_0) 
	{
		tank->ztbm += 1 * pow(10, 12); 
		
		if (tank->isStart <= 1 && tank->CY <= 2 && tank->CY > 0)
		{
			tank->isStart=1;
			TankInit(1, tank,MData);
		}

	}
	else 
	{
		tank->ztbm += 2 * pow(10, 12); 
		
		if(tank->CY <= 2 && (tank->CZ ==20 ||tank->CZ ==21))
		{
			if (tank->isStart <= 1 ) 
			{
				tank->isStart=1;
				TankInit(1, tank,MData);
			}
		}
		else
		{
			tank->MaintainTime=-1;
			tank->isStart=0;
			tank->GJ=2;
		}
		
	}
	StablePeriod =tank->Stable_Day * 24* 3600; 
	if (tank->ForeCast_Sec >= StablePeriod && (tank->isStart==1))
	{
		
		{
			
			
			
			

			tank->ztbm += 1 * pow(10, 9); 
			tank->isStart=2;
			tank->lamda=1;
		}
		
		//tank->DisplayInfo = 3;
	}
	else
	{
		tank->ztbm += 2 * pow(10, 9); 
		if (tank->isStart==1)
		{
			//tank->DisplayInfo = 2;
			lam=tank->ForeCast_Sec / StablePeriod + 0.2;
			tank->lamda=pow(lam, tank->Lamda_0);

			if(tank->ForeCast_Sec > StablePeriod *(1-0.2))
			{
				tank->lamda=1;
			}


			tank->YJ=1;
		}
	}

	tank->RhoLn = CalculateLiquidRho(tank->ch4, tank->Pn); 
	tank->RhoL = CalculateLiquidRho(tank->ch4, MData->Record_P); 
	tank->faiMax = tank->Eta_V * tank->RhoLn / tank->RhoL;
	tank->mMin = tank->Mavg_0 * (tank->Pn + tank->Pa_0) * tank->TankVolYX * 1000 / (tank->R_0 * tank->Tc_0);

	if (tank->fai > tank->L_0)
	{
		if (tank->fai > tank->faiMax)
		{
			tank->CZ = 12;
			tank->ztbm += 1 * pow(10, 11);
		}
		else
		{
			tank->CZ = 11;
			tank->ztbm += 2 * pow(10, 11);
		}
	}
	else
	{
		if (tank->m0 >= tank->mMin)
		{
			tank->CZ = 21;
			tank->ztbm += 3 * pow(10, 11);
		}
		else
		{
			tank->CZ = 22;
			tank->ztbm += 4 * pow(10, 11);
		}
	}

	if (tank->isStart>=1) 
	{
		

		tank->VLf=tank->TankVolYX*MData->Record_L;
		tank->VGf=tank->TankVolYX*(1-MData->Record_L);
		CalculatePhyChemPara(tank, MData->Record_P, PCi );	
		tank->TavgL=(MData->Record_T+tank->TLp)*0.5;  
		

		if (tank->CZ == 21||tank->CZ == 20) 
		{
	
			if (tank->fai < tank->minFai || tank->rhoG >= tank->m0 / tank->TankVolYX || MData->Record_Tl > tank->TLp +
				tank->dTmax)
			{
				//tank->XT = 1;
				tank->CZ = 20; 
				tank->MaintainTime = EmptyTankAlgorithm(tank, MData); // *PCi,  cqr-20230213  //cqr-1028修改
			}
			else
			{
				tank->CZ = 21; 
				SaveDailyHeatLoss(tank);
				tank->MaintainTime = FullTankAlgorithm(tank, MData);
			}
		}
		if (tank->CZ == 11 || tank->CZ == 12) 
		{
			SaveDailyHeatLoss(tank);
			tank->MaintainTime = FullTankAlgorithm(tank, MData);
		}
		if (tank->CZ == 22)
		{
			tank->MaintainTime = tank->MTime_0;
		}
	
		if (tank->MaintainTime > tank->MTime_0)
		{
			tank->MaintainTime = tank->MTime_0;
		}

	
		if (tank->SFKZ1 == 2 && MData->Record_P >= tank->RQP && MData->Record_P <= tank->QZP)
		{
			tank->MaintainTime = (tank->MaintainTime - tank->RJJ) / tank->RXL;
		}		
	}

	 
	if (tank->MaintainTime<0)
	{
		tank->MaintainTime=0;
	}

	if (tank->MaintainTime < 5)
	{
		tank->ztbm += 2 * pow(10, 10); 
	}

	
	tank->ztbm += tank->SFKZ1 * pow(10, 7);
	tank->ztbm += tank->SFKZ2 * pow(10, 6);
	tank->ztbm += tank->SFKZ3 * pow(10, 5);
	tank->ztbm += tank->SFKZ4 * pow(10, 4);
	tank->ztbm += tank->SFKZ5 * pow(10, 3);
	tank->ztbm += tank->SFKZ6 * pow(10, 2);
	//k=0; //cqr-202302323
	/*for(i=1;i<50;i++)
	{
		if (tank->RZ[i]>=1)
		{
			k=k+1;
			tank->RZ[k]=tank->RZ[i];
			if (i>k)
			{
				tank->RZ[i]=0;
			}
			
		}
		
	}
	tank->RZ[0]=k;*/
	tank->QDay_0 = tank->SQdLng2;

	
}

void TankInit(int InitTime, struct TankPara *tank, struct MonitoringData *data)
{
	float V0,ch4,fai0;//,p,k1;
	float RhoGn,Pn, P2,m0,RhoL,RhoG,VL,VG,HL,HG,fai_m;	

	float Factor_Rho=1.0; 
	float Factor_P=1.0;  



	// k1=tank->K_1;

	tank->P0=data->Record_P; 
	tank->T0=data->Record_T;

	//p=tank->P0;  
	V0=tank->TankVolYX;
	ch4=tank->ch4;
	fai0=tank->ZeroL;
	tank->aLNG20=tank->a20* tank->K_1;

	if (ch4<70) 
	{
		
	}
	if (ch4>98)
	{
		
	}

	tank->QHQR0=tank->QHQR_0; 
	tank->RhoLNG0=tank->RhoLng_0;
	tank->QdLNG20=V0* (tank->Fai_bz0 *0.01) * tank->RhoLNG0 *( tank->aLNG20*0.01) * tank->QHQR0;

	RhoL=CalculateLiquidRho(ch4,tank->P0); 
	RhoG=CalculateGasRho(ch4,tank->P0); 
	
	VL=V0*fai0;
	VG=V0-VL;//V0*(1-fai0);
	m0 = tank->m0;

	/*if (V0<0.01) 
	{
		
	}
	if (p<0) 
	{
		
	}
	if (p>1) 
	{
		
	}*/

	tank->RhoLimit=m0/(V0*tank->Eta_V);   
	//tank->m0=m0;

	HL=CalculateLiquid_H(ch4,tank->P0); 
	HG=CalculateGasH(ch4,tank->P0); 

	tank->Q0=RhoL*VL*HL+RhoG*VG*HG;


	Pn=tank->P75;
	tank->RhoLn=CalculateLiquidRho(ch4,Pn); 

	/*if (tank->RhoLimit<1.4) 
	{

	}
	if (tank->RhoLimit>500)  
	{

	}

	if (RhoL<350) 
	{

	}
	if (RhoL>500) 
	{

	}*/

	P2= Factor_P*CalculateLiquidPressure(ch4,tank->RhoLimit*Factor_Rho);

	fai_m=tank->Eta_V*tank->RhoLn/RhoL; 

	tank->faiMax=fai_m;

	
	////tank->CalModel=2; 
	if(data->Record_L>=fai_m)
	{
		Pn=P2;
		//tank->CalModel=3; 
		tank->CZ=12;
	}
	tank->Pn=Pn;
	tank->Tn=CalculateLiquidTemperature(ch4,Pn);
	tank->TAvg=(tank->T0 + tank->Tn)*0.5; 
	
	RhoGn=CalculateGasRho(ch4,Pn); 
	

	//if (tank->R_0*tank->Tc_0<0.1) 
	//{
	//	
	//}
	//
	
	tank->mMin=tank->Mavg_0*(Pn+tank->Pa_0)*V0*1000/(tank->R_0*tank->Tc_0);  
	
	HL=CalculateLiquid_H(ch4,Pn); 
	HG=CalculateGasH(ch4,Pn); 

	RhoL=CalculateLiquidRho(ch4,Pn); 
	RhoG=CalculateGasRho(ch4,Pn); 

	
	tank->rhoG_n=RhoG;
	tank->rhoL_n=RhoL;


	VL=(m0-V0*RhoG)/(RhoL-RhoG);
	VG=V0-VL;

	tank->Qn=RhoL*VL*HL+RhoG*VG*HG;


	if (InitTime==0)
	{	
		//tank->DisplayInfo=1;
		tank->isStart=0;
		tank->MaintainTime=0;

		tank->ZeroDate = -1;
		tank->ZeroL =0; 
		tank->ForeCast_Sec = 0;
		tank->ForeCast_Date = 0; 

		tank->nday=0;  
		tank->nData=0;
		tank->iday=0;

		tank->ihour=0;
		tank->iQday=0;
	}
	
}



void ConvertDateFormat( double Record_Date, struct tm *Date_Str ) 
{
	char Str_Temp[15];
	
	sprintf(Str_Temp, "%.0lf", Record_Date);
	
	Date_Str->tm_year = (Str_Temp[0] - '0') * 1000 + (Str_Temp[1] - '0') * 100 + (Str_Temp[2] - '0') * 10 + (Str_Temp[3] - '0') - 1900;
	Date_Str->tm_mon = (Str_Temp[4] - '0') * 10 + (Str_Temp[5] - '0') - 1;
	Date_Str->tm_mday = (Str_Temp[6] - '0') * 10 + (Str_Temp[7] - '0');
	Date_Str->tm_hour = (Str_Temp[8] - '0') * 10 + (Str_Temp[9] - '0');
	Date_Str->tm_min = (Str_Temp[10] - '0') * 10 + (Str_Temp[11] - '0');
	
	Date_Str->tm_sec = 0;
	Date_Str->tm_wday = 0;
	
}


void SaveDailyHeatLoss( struct TankPara * tank)//, struct PhysicochemicalPara * PCi )
{
	short int k,i=0;
	float tQ=0;
	float ForeCast_Date = tank->ForeCast_Date;
	float Thour=ForeCast_Date*24; 
	short int Nd=floor(ForeCast_Date);
	short int n=tank->Tz_0;         
	tank->DayNow=ForeCast_Date; 

	if (Thour<0)
	{
		Thour=0;
	}

	k=tank->iday;
	if (tank->ihour< n)  
	{
		tank->Qhour[tank->ihour]=tank->Qi; //cqr-20230230  PCi->Q; 
		tank->hour[tank->ihour]=Thour;
		tank->ihour=tank->ihour+1;
		tank->iQday=0;

	}
	else   
	{
		tank->dQi[tank->iQday]= 24* (tank->Qi - tank->Qhour[0]) / (Thour- tank->hour[0]);   //cqr-20230230  PCi->Q; 
		tank->iQday=tank->iQday+1;

		for (i=0;i<tank->ihour-1;i++) 
		{
			tank->Qhour[i]=tank->Qhour[i+1];
			tank->hour[i]=tank->hour[i+1];

		}
		tank->Qhour[tank->ihour-1]=tank->Qi; //cqr-20230230  PCi->Q; 
		tank->hour[tank->ihour-1]=Thour ;


		if (tank->iQday>0)  
		{
			tQ=0;
			for (i=0;i<tank->iQday;i++) 
				tQ+=tank->dQi[i];

			if (tank->iday >= tank->Nday_0 ) 
			{				
				for (i=0;i< tank->iday -1;i++) 
				{
					tank->Qday[i]=tank->Qday[i+1];
					tank->day[i]=tank->day[i+1];
				}
				tank->iday =tank->iday -1;
				k=tank->iday;
			}
			
			tank->Qday[k]=tQ/tank->iQday;  
			tank->day[k]=Nd;
		}


		if (tank->iQday==n) 
		{
			for (i=0;i<tank->iQday;i++) 
				tank->dQi[i]=0;
			tank->iQday=0;
			tank->iday=k+1;
			k=k+1;
		}
	}

	

//	tank->Qi=PCi->Q;

}


float Calculate_m(struct TankPara* tank, struct MonitoringData* data)
{
	float mi, RhoL, RhoG, VL, VG;
	RhoL = CalculateLiquidRho(tank->ch4, data->Record_P); 
	RhoG = CalculateGasRho(tank->ch4, data->Record_P); 
	VL = tank->TankVolYX * data->Record_L;
	VG = tank->TankVolYX - VL;
	mi = RhoL * VL + RhoG * VG;
	return mi;
}


void CalculatePhyChemPara(struct TankPara *tank, float p, struct PhysicochemicalPara *PCi )
{
	float m0=tank->m0;
	float V0=tank->TankVolYX;
	float ch4=tank->ch4;

	(*PCi).P=p;
	(*PCi).HL=CalculateLiquid_H(ch4,p); 
	(*PCi).HG=CalculateGasH(ch4,p); 

	(*PCi).RhoL=CalculateLiquidRho(ch4,p); 
	(*PCi).RhoG=CalculateGasRho(ch4,p); 

	//if ((*PCi).RhoL<350)
	{
//		tank->RZ[7]=61;
	}
//	if ((*PCi).RhoL>500)
	{
	//	tank->RZ[7]=62;
	}

	(*PCi).VL=(m0-V0*(*PCi).RhoG)/((*PCi).RhoL-(*PCi).RhoG);
	(*PCi).VG=V0-(*PCi).VL;

	tank->VLp=(*PCi).VL;
	tank->VGp=(*PCi).VG;
	tank->TLp=CalculateLiquidTemperature(ch4,p);

	tank->fai_p=tank->VLp/V0;
	tank->VLavg=(tank->VLp+tank->VLf)*0.5;
	tank->VGavg=(tank->VGp+tank->VGf)*0.5;


	
	tank->rhoG=(*PCi).RhoG;
	tank->rhoL=(*PCi).RhoL;
	tank->mL=tank->VLavg*(*PCi).RhoL;
	tank->mG=tank->VGavg*(*PCi).RhoG;
	

	PCi->QL=(*PCi).RhoL*(tank->VLp+tank->VLf)*0.5*(*PCi).HL; 
	PCi->QG=(*PCi).RhoG*(tank->VGp+tank->VGf)*0.5*(*PCi).HG; 
	PCi->Q=(PCi->QL)+(PCi->QG); 

	
	tank->Qi=PCi->Q; //cqr-20230230

//	if (PCi->QL <0.01)
	{
//		tank->RZ[8]=71;
	}

}


float RussiaCorrect(struct  TankPara *tank ,struct  MonitoringData *data, float F ) 
{	
	float LngFillRate=tank->fai;
	float Pn=tank->Pn;  
	float Pi=data->Record_P;

	



	float Pcr=tank->Pcr_0 - tank->Pa_0;
	float ZeroDimension_P = (Pn - Pi) / (Pcr - Pi);

	float s1=(2.2125 - 1.66*ZeroDimension_P + 0.83*ZeroDimension_P*ZeroDimension_P);
	float s2=(-4.38*LngFillRate*LngFillRate*LngFillRate + 8.861*LngFillRate*LngFillRate - 5.539*LngFillRate + 1.959);
	float s3=(-204 * F*F + 20.48*F + 0.8);

	


	float FR = s1*s2*s3;
	tank->Kr=FR;
	return FR;
}


float Calculate_k2(struct  TankPara *tank ,struct  MonitoringData *data) 
{
	float RegressionQHQR(float x);
	float t1=110.56; 
	float t2=293.15; 
	float k2=1;
	float QHQR0=tank->QHQR0; 
	float TAvg=tank->TAvg;

	float Pi=data->Record_P;
	float OutT=data->Record_OutT;

	float LT=(data->Record_T+tank->TLp)*0.5;
	float QHQRi=RegressionQHQR(Pi);
	//int tankType=tank->TankType;

	

	//
	//switch (tankType)
	//{
	//	case 1: 
	//		k2 = 0.7*(t2 - t1) / (OutT - LT) + 0.3* (t2*t2*t2*t2 - t1*t1*t1*t1) /(OutT*OutT*OutT*OutT - LT* LT* LT* LT);
	//		break;
	//	case 2:   
	//		k2 =(t2 - t1) / (OutT - LT);
	//		break;
	//	case 3:  
	//		k2 = (t2*t2*t2*t2 - t1*t1*t1*t1) /(OutT*OutT*OutT*OutT - LT* LT* LT* LT);
	//		break;
	//	default:
	//		k2 = 0.7*(t2 - t1) / (OutT - LT) + 0.3* (t2*t2*t2*t2 - t1*t1*t1*t1) /(OutT*OutT*OutT*OutT - LT* LT* LT* LT);
	//		break;
	//}

	switch (tank->TankType)
	{
	case 1: 
		k2 =127.813 / (OutT - LT) + 2170722109 /(OutT*OutT*OutT*OutT - LT* LT* LT* LT);
		break;
	case 2:   
		k2 =182.59/ (OutT - LT);
		break;
	case 3:  
		k2 = 7235740364 /(OutT*OutT*OutT*OutT - LT* LT* LT* LT);
		break;
	default:
		k2 = 127.813 / (OutT - LT) + 2170722109 /(OutT*OutT*OutT*OutT - LT* LT* LT* LT);
		break;
	}

	//if (QHQRi<0.001)
	{
//		tank->RZ[10]=91;
	}
	
	k2=k2*QHQRi/QHQR0;
	tank->K2=k2;
	return k2;
}


float Calculate_kzc(struct  TankPara *tank ) 
{
	float k;
	int zc=tank->ZC;
	switch (zc)
	{
	case 1: 
		k = 0.8;
		break;
	case 2:  
		k =1 ;
		break;	
	default:
		k = 1;
		break;
	}	
	tank->Kf=k;
	return k;
}


float Calculate_ky(struct  TankPara *tank ,struct  MonitoringData *data) 
{
	float k;
	float Re2,Re3;
	
	
	float Af=tank->Nx_0*data->Record_Amp_x*data->Record_fr_x+tank->Ny_0*data->Record_Amp_y*data->Record_fr_y+tank->Nz_0*data->Record_Amp_z*data->Record_fr_z;
	float a=0;
	float b=0;
	float c=1;
	float l=0;
	float m=0;
	float n=1;


	float Re=40000*Af*tank->D/tank->VN_0; 
	Re2=Re*Re;
	Re3=Re2*Re;
	k=tank->Da_0*Re3-tank->Db_0*Re2-tank->Dc_0*Re-tank->Dd_0;  

	k=k*0.0001; 

	if (k< tank->Kymin)
	{
		k=tank->Kymin;
	}
	if (k> tank->Kymax)
	{
		k=tank->Kymax;
	}

	
	{
		k=1;
	}
	
	tank->Ky=k;
	return k;
}



void CalSQdLng2( struct TankPara *tank,struct MonitoringData *data ) 
{
	short int i;
	float tQ=0;
	float k2 = Calculate_k2(tank,data);
	float Thour=tank->ForeCast_Date*24; 

	tank->LQdLng1=tank->QdLNG20 / k2;

	tQ=0;
	for (i=0;i<tank->iday;i++)  
	{
		tQ+=tank->Qday[i];
	}
	if (Thour < 12*tank->Nday_0) //24/2
	{		
		if (tank->SQdLng2<0.001) 
		{
			tank->LR=1;  
			tank->QdAvg=(tQ + tank->LQdLng1 *(tank->Nday_0-tank->iday))/tank->Nday_0; 
		}
		else  
		{
			tank->LR=2;  
			
			tank->SQdLng2=(tQ+tank->SQdLng2*(tank->Nday_0-tank->iday))/tank->Nday_0;
			tank->QdAvg=tank->SQdLng2;
		}
	}
	else if(Thour < 24*tank->Nday_0)
	{
		if (tank->SQdLng2<0.001)
		{
			tank->LR=3;
			tank->QdAvg=(tQ+tank->LQdLng1*(tank->Nday_0-tank->iday))/tank->Nday_0;
		}
		else
		{
			tank->LR=2;
			tank->SQdLng2=(tQ+tank->SQdLng2*(tank->Nday_0-tank->iday))/tank->Nday_0;
			tank->QdAvg=tank->SQdLng2;
		}
	}
	else
	{
		tank->LR=2;

		tank->SQdLng2=tQ/tank->iday;
		tank->QdAvg=tank->SQdLng2;
	}
	tank->ztbm += tank->LR * pow(10, 8);
}



float t1ByTheoretical(struct TankPara *tank,struct MonitoringData *data)
{
	float Calculate_k2(struct  TankPara *tank ,struct  MonitoringData *data);
	float RussiaCorrect(struct  TankPara *tank ,struct  MonitoringData *data, float F );
	
	

	float QHQR0=tank->QHQR0; 
	float QdLNG20=tank->QdLNG20;

	float Qn=tank->Qn; 
	float Qi=tank->Qi;
	
	float k2 = Calculate_k2(tank,data);

	

	float QdLNG1=tank->QdLNG20 / k2;
	float t0=(Qn-Qi)/QdLNG1;


	float ci=CalculateLiquidC(tank->ch4,data->Record_P); 

	float F=QdLNG1/(tank->m0*ci*tank->TavgL)/24; 

	
	//float Factor_Stable = 0.95;

	float Factor_Russia = RussiaCorrect(tank,data, F);
	float Factor_Kf = Calculate_kzc(tank) ;


	float Factor_y =Calculate_ky(tank,data) ;

	float MaintainTime = t0*Factor_Kf* Factor_y *tank->lamda/ Factor_Russia;



//	if (Factor_Russia<0.7)
	{
//		tank->RZ[9]=81;
	}
//	if (Factor_Russia>2)
	{
	//	tank->RZ[9]=82;
	}
	return MaintainTime;
}



float t2ByMeasuredData(struct TankPara *tank,struct MonitoringData *data) 
{
	float Calculate_k2(struct  TankPara *tank ,struct  MonitoringData *data);
	float RussiaCorrect(struct  TankPara *tank ,struct  MonitoringData *data, float F );
	float Qn=tank->Qn; 
	float Qi=tank->Qi;
	float tQ=Qn-Qi;

	float k2=Calculate_k2(tank,data);
	float QdLNG2=tank->QdAvg;

	float t0=tQ/QdLNG2; 

	
	//float Factor_Stable = 0.95;

	float ci=CalculateLiquidC(tank->ch4,data->Record_P); 
	float F=QdLNG2/(tank->m0*ci*tank->TavgL)/24;
	


	float Factor_Russia = RussiaCorrect(tank,data, F);

	float Factor_y =Calculate_ky(tank,data) ;

	

	float MaintainTime = t0*Factor_y / Factor_Russia*tank->lamda;

	
	
	return MaintainTime;
}


float FullTankAlgorithm(struct TankPara *tank,struct MonitoringData *data)
{
	float MT1,MT2;
	float MaintainTime=200;
	float Nday=tank->Nday_0;
	short int Nd2=Nday/2;
	float day = tank->ForeCast_Date;
	float Lamda1=1.0;

	CalSQdLng2(tank,data);

	if (tank->LR==1)
	{
		if (tank->CZ == 21)
		{
			MT1 = (tank->Qn - tank->Qi) / tank->QdLNG20;
		}
		else
		{
			Lamda1 = 1;
			MT1 = t1ByTheoretical(tank, data) * Lamda1; // pci,  cqr-20230213
		}

		MaintainTime = floor(MT1);
	}
	else if (tank->LR == 2)
	{
	
		if (tank->CZ == 21)
		{
			MT2 = (tank->Qn - tank->Qi) / tank->QdAvg;
		}
		else
		{
			MT2 = t2ByMeasuredData(tank, data); // pci,  cqr-20230213
		}
		MaintainTime = floor(MT2);
	}
	else 
	{
		
		if (tank->CZ == 21)
		{
			MT1 = (tank->Qn - tank->Qi) / tank->QdLNG20;
			MT2 = (tank->Qn - tank->Qi) / tank->QdAvg;
		}
		else
		{
			MT1 = t1ByTheoretical(tank, data); 
			MT2 = t2ByMeasuredData(tank, data); 
		}
		MaintainTime = (Nday - day) * MT1 / Nd2 + (day - Nd2) * MT2 / Nd2;
		MaintainTime = floor(MaintainTime);
	}
	return MaintainTime;
}


float UniteAlgorithm(struct TankPara *tank,struct MonitoringData *data)
{
	float t0,QdLNG2, F, Factor_Russia, Afai;
	short int MTk,MTm;
	float MT,MT1,MT2;
	short int MaintainTime=tank->MTime_0;
	short int  Nday=tank->Nday_0;
	short int Nd2=Nday/2;
	float day = tank->ForeCast_Date;


	float k2 = Calculate_k2(tank,data);

	float ci=CalculateLiquidC(tank->ch4,data->Record_P); 

	float Factor_Kf = Calculate_kzc(tank) ;

	float Factor_y =Calculate_ky(tank,data) ;

	CalSQdLng2(tank,data);

	if (tank->LR==1)
	{
		//////
		t0=(tank->Qn - tank->Qi) *k2 /tank->QdLNG20 ;

		F=tank->QdLNG20/ k2 /(tank->m0*ci*tank->TavgL)/24; 
		Factor_Russia = RussiaCorrect(tank,data, F);

		MT = t0*Factor_Kf* Factor_y *tank->lamda/ Factor_Russia;

		MTm= floor(MT);

		MT=(tank->Qn_kg - tank->Qi_kg) *k2 /tank->QdLNG20 ;
		MTk= floor(MT);


	}
	else if(tank->LR==2)
	{

		t0=(tank->Qn - tank->Qi)/tank->QdAvg; 

		F=tank->QdAvg/(tank->m0*ci*tank->TavgL)/24; // QdLNG2=tank->QdAvg

		Factor_Russia = RussiaCorrect(tank,data, F);

		MT = t0* Factor_y / Factor_Russia*tank->lamda;
		MTm= floor(MT);

		MT=(tank->Qn_kg - tank->Qi_kg)/tank->QdAvg; 
		MTk= floor(MT);
	}
	else  
	{
		t0=(tank->Qn - tank->Qi) *k2 /tank->QdLNG20 ;
		F=tank->QdLNG20/ k2 /(tank->m0*ci*tank->TavgL)/24; 
		Factor_Russia = RussiaCorrect(tank,data, F);
		MT1 = t0*Factor_Kf* Factor_y *tank->lamda/ Factor_Russia;

		t0=(tank->Qn - tank->Qi)/tank->QdAvg; 
		F=tank->QdAvg/(tank->m0*ci*tank->TavgL)/24; // QdLNG2=tank->QdAvg
		Factor_Russia = RussiaCorrect(tank,data, F);
		MT2 = t0* Factor_y / Factor_Russia*tank->lamda;

		MT=(Nday-day)*MT1/Nd2+(day-Nd2)*MT2/Nd2;
		MTm= floor(MT);
		
		MT1=(tank->Qn_kg - tank->Qi_kg) *k2 /tank->QdLNG20 ;
		MT2=(tank->Qn_kg - tank->Qi_kg)/tank->QdAvg; 
		MT=(Nday-day)*MT1/Nd2+(day-Nd2)*MT2/Nd2;
		MTk= floor(MT);

	}

	Afai=acos(1 - data->Yewei/(0.5*tank->D))/3.14159;
	MT=( tank->ZQY*Afai*MTm + tank->ZQQ*(1- Afai)*MTk ) / (tank->ZQY*Afai + tank->ZQQ*(1- Afai));  //ZQY =k2 ,ZQQ=k1

	MaintainTime= floor(MT);

	return MaintainTime;
}


float ExpandTankAlgorithm(struct PhysicochemicalPara *PCn)
{
	float Pn=0.75;
	float SafetyFactor=0.95;

	float RhoSafe=SafetyFactor*(*PCn).RhoL;
}




float EmptyTankAlgorithm(struct TankPara *tank,struct MonitoringData *data) 
{

	float MT1, MT2;
	float MaintainTime = tank->MTime_0;
	float Nday = tank->Nday_0;
	short int Nd2 = Nday / 2;
	float day = tank->ForeCast_Date;
	float v, Hi, Hn; 

	v = tank->R_0 * data->Record_Tg / ((data->Record_P + 0.1) * tank->Mavg_0 * 1000);
	Hi = CalculateGasH(tank->ch4, data->Record_P);
	Hn = CalculateGasH(tank->ch4, tank->Pn);
	tank->Qi = tank->TankVolYX * Hi / v;
	tank->Qn = tank->TankVolYX * Hn / v;

	SaveDailyHeatLoss(tank);

	CalSQdLng2(tank, data);

	if (tank->LR == 1)
	{
		MT1 = (tank->Qn - tank->Qi) / tank->QdLNG20;
		MaintainTime = floor(MT1);
	}
	else if (tank->LR == 2)
	{
		MT2 = (tank->Qn - tank->Qi) / tank->QdAvg;
		MaintainTime = floor(MT2);
	}
	else
	{
		MT1 = (tank->Qn - tank->Qi) / tank->QdLNG20;
		MT2 = (tank->Qn - tank->Qi) / tank->QdAvg;

		MaintainTime = (Nday - day) * MT1 / Nd2 + (day - Nd2) * MT2 / Nd2;
		MaintainTime = floor(MaintainTime);
	}

	return MaintainTime;

}



float RegressionQHQR(float x)
{
	float result = 0; 
	float x2, x3;
	x2 = x * x;
	x3 = x2 * x;

	if (x <= 0.5)
	{
		result = -178.95 * x3 + 226.49 * x2 - 192.60 * x + 510.65;
	}
	else
	{
		result = -12.14 * x3 + 42.16 * x2 - 123.63 * x + 501.63;
	}
	return result;
}


float CalculateLiquidTemperature(float ch4,float x)
{
	float RegressionFormulaOfLT(short int CH4,float x);
	float y1,y2, y =0;
	short int CH4,CH5;

	CH4 = 98 - floor((98-ch4)/4)*4;
	y1 = RegressionFormulaOfLT(CH4, x);

	CH5 = 98 - ceil((98-ch4)/4)*4;
	y2 = RegressionFormulaOfLT(CH5, x);

	if (CH4==CH5)
	{
		y=y1;
	}
	else
		y=y1+(y2-y1)*(ch4-CH4)/(CH5-CH4); 
	return y;

}

float RegressionFormulaOfLT(short int CH4,float x)
{
	float y = 0; 
	float x2, x3; 
	x2 = x * x;
	x3 = x2 * x;
	
	if (x <= 0.5)
	{
		switch (CH4)
		{
		case 98: y = 122.33 * x3 - 154.12 * x2 + 100.88 * x + 112.15;
			break;
		case 94: y = 123.45 * x3 - 155.49 * x2 + 102.03 * x + 112.64;
			break;
		case 90: y = 124.55 * x3 - 156.88 * x2 + 103.19 * x + 113.14;
			break;
		case 86: y = 125.46 * x3 - 158.15 * x2 + 104.35 * x + 113.63;
			break;
		case 82: y = 126.70 * x3 - 159.66 * x2 + 105.58 * x + 114.13;
			break;
		case 78: y = 127.80 * x3 - 161.10 * x2 + 106.85 * x + 114.65;
			break;
		case 74: y = 128.92 * x3 - 162.60 * x2 + 108.19 * x + 115.18;
			break;
		case 70: y = 130.42 * x3 - 164.42 * x2 + 109.65 * x + 115.73;
			break;

		}
	}
	else
	{
		switch (CH4)
		{
		case 98: y = 7.60 * x3 - 27.91 * x2 + 53.86 * x + 118.27;
			break;
		case 94: y = 7.67 * x3 - 28.14 * x2 + 54.57 * x + 118.83;
			break;
		case 90: y = 7.70 * x3 - 28.32 * x2 + 55.26 * x + 119.38;
			break;
		case 86: y = 7.83 * x3 - 28.71 * x2 + 56.12 * x + 119.91;
			break;
		case 82: y = 7.89 * x3 - 28.95 * x2 + 56.87 * x + 120.48;
			break;
		case 78: y = 7.92 * x3 - 29.13 * x2 + 57.62 * x + 121.07;
			break;
		case 74: y = 7.99 * x3 - 29.41 * x2 + 58.46 * x + 121.67;
			break;
		case 70: y = 8.13 * x3 - 29.85 * x2 + 59.47 * x + 122.27;
			break;

		}
	}


	return y;
}



float CalculateLiquidRhoAtP0(float ch4)
{
	float y =0;

	y= -0.0101343* ch4*ch4 - 0.748486*ch4 + 598.8173 ;
	return y;

}


float CalculateLiquidRho(float ch4,float x)
{
	float RegressionFormulaOfLRho(short int CH4, float x);
	float y1, y2, y = 0;
	short int CH4, CH5;

	CH4 = 98 - floor((98-ch4)/4)*4;
	y1 = RegressionFormulaOfLRho(CH4, x);

	CH5 = 98 - ceil((98-ch4)/4)*4;
	y2 = RegressionFormulaOfLRho(CH5, x);

	if (CH4==CH5)
	{
		y=y1;
	}
	else
		y=y1+(y2-y1)*(ch4-CH4)/(CH5-CH4); 
	return y;

}

float RegressionFormulaOfLRho(short int CH4,float x)
{
	float y = 0; 
	float x2, x3, x4, x5, x6;
	x2 = x * x;
	x3 = x2 * x;
	if (x <= 0.5)
	{
		switch (CH4)
		{
		case 98: y = -164.03 * x3 + 206.76 * x2 - 148.13 * x + 427.60;
			break;
		case 94: y = -163.15 * x3 + 205.76 * x2 - 146.79 * x + 438.41;
			break;
		case 90: y = -162.53 * x3 + 205.02 * x2 - 145.72 * x + 448.86;
			break;
		case 86: y = -162.12 * x3 + 204.50 * x2 - 144.89 * x + 458.99;
			break;
		case 82: y = -161.84 * x3 + 204.14 * x2 - 144.28 * x + 468.78;
			break;
		case 78: y = -161.74 * x3 + 203.97 * x2 - 143.89 * x + 478.26;
			break;
		case 74: y = -161.75 * x3 + 203.97 * x2 - 143.72 * x + 487.42;
			break;
		case 70: y = -161.96 * x3 + 204.20 * x2 - 143.79 * x + 496.27;
			break;

		}
	}
	else
	{
		switch (CH4)
		{
		case 98: y = -10.61 * x3 + 37.76 * x2 - 85.04 * x + 419.37;
			break;
		case 94: y = -10.37 * x3 + 37.30 * x2 - 83.78 * x + 430.16;
			break;
		case 90: y = -10.40 * x3 + 37.39 * x2 - 83.10 * x + 440.68;
			break;
		case 86: y = -10.28 * x3 + 37.16 * x2 - 82.36 * x + 450.81;
			break;
		case 82: y = -10.27 * x3 + 37.13 * x2 - 81.90 * x + 460.63;
			break;
		case 78: y = -10.22 * x3 + 37.03 * x2 - 81.54 * x + 470.11;
			break;
		case 74: y = -10.33 * x3 + 37.28 * x2 - 81.56 * x + 479.32;
			break;
		case 70: y = -10.24 * x3 + 37.10 * x2 - 81.41 * x + 488.12;
			break;

		}
	}

	return y;
}


float CalculateLiquidPressure(float ch4,float x)
{
	float RegressionFormulaOfLP(short int CH4, float x);
	float y1, y2, y = 0;
	short int CH4, CH5;

	CH4 = 98 - floor((98-ch4)/4)*4;
	y1 = RegressionFormulaOfLP(CH4, x);

	CH5 = 98 - ceil((98-ch4)/4)*4;
	y2 = RegressionFormulaOfLP(CH5, x);

	if (CH4==CH5)
	{
		y=y1;
	}
	else
		y=y1+(y2-y1)*(ch4-CH4)/(CH5-CH4); 

	return y;

}

float RegressionFormulaOfLP(short int CH4,float x)
{
	float y = 0; 
	float x2, x3,x4;


	x=x/1000.0;

	x2 = x * x;
	x3 = x2 * x;
	x4 = x3 * x;

	if (CH4 == 98) { y = -5705.7232 * x4 + 9012.0160 * x3 - 5185.3260 * x2 + 1269.7312 * x - 108.6789; }
	else if (CH4 == 94) { y = -5807.2057 * x4 + 9386.1696 * x3 - 5530.0392 * x2 + 1388.1830 * x - 122.0828; }
	else if (CH4 == 90) { y = -5957.1338 * x4 + 9848.0073 * x3 - 5940.2872 * x2 + 1529.3348 * x - 138.4160; }
	else if (CH4 == 86) { y = -5987.9029 * x4 + 10114.9401 * x3 - 6237.8634 * x2 + 1643.5345 * x - 152.5390; }
	else if (CH4 == 82) { y = -6065.9838 * x4 + 10465.5030 * x3 - 6597.6234 * x2 + 1779.7678 * x - 169.6392; }
	else if (CH4 == 78) { y = -6061.2901 * x4 + 10671.7238 * x3 - 6869.9006 * x2 + 1894.5804 * x - 185.0293; }
	else if (CH4 == 74) { y = -6111.6234 * x4 + 11082.4332 * x3 - 7358.2200 * x2 + 2098.0316 * x - 212.8264; }
	else { y = -6106.1082 * x4 + 10966.0253 * x3 - 7207.7210 * x2 + 2032.8882 * x - 203.6900; }


	return y;
}


float CalculateLiquid_H(float ch4,float x)
{
	float RegressionFormulaOfLH(short int CH4, float x);
	float y1, y2, y = 0;
	short int CH4, CH5;

	CH4 = 98 - floor((98-ch4)/4)*4;
	y1 = RegressionFormulaOfLH(CH4, x);

	CH5 = 98 - ceil((98-ch4)/4)*4;
	y2 = RegressionFormulaOfLH(CH5, x);

	if (CH4==CH5)
	{
		y=y1;
	}
	else
		y=y1+(y2-y1)*(ch4-CH4)/(CH5-CH4); 
	return y;

}

float RegressionFormulaOfLH(short int CH4,float x)
{
	float y = 0; 
	float x2, x3, x4, x5, x6;
	x2 = x * x;
	x3 = x2 * x;
	if (x <= 0.5)
	{
		switch (CH4)
		{
		case 98: y = 408.6621 * x3 - 514.7035 * x2 + 349.7897 * x + 0.7425;
			break;
		case 94: y = 401.8792 * x3 - 506.2911 * x2 + 344.0044 * x + 0.7301;
			break;
		case 90: y = 395.8130 * x3 - 498.7525 * x2 + 338.9289 * x + 0.7188;
			break;
		case 86: y = 390.4063 * x3 - 492.0258 * x2 + 334.5137 * x + 0.7087;
			break;
		case 82: y = 385.6236 * x3 - 486.0729 * x2 + 330.7247 * x + 0.6996;
			break;
		case 78: y = 381.4272 * x3 - 480.8513 * x2 + 327.5332 * x + 0.6913;
			break;
		case 74: y = 377.8063 * x3 - 476.3438 * x2 + 324.9221 * x + 0.6840;
			break;
		case 70: y = 374.7320 * x3 - 472.5257 * x2 + 322.8822 * x + 0.6773;
			break;

		}
	}
	else
	{
		switch (CH4)
		{
		case 98: y = 25.6124 * x3 - 92.9922 * x2 + 192.4485 * x + 21.2673;
			break;
		case 94: y = 25.1048 * x3 - 91.4094 * x2 + 189.1625 * x + 20.9381;
			break;
		case 90: y = 24.6776 * x3 - 90.0344 * x2 + 186.3589 * x + 20.6360;
			break;
		case 86: y = 24.4165 * x3 - 89.0386 * x2 + 184.1263 * x + 20.3316;
			break;
		case 82: y = 24.0721 * x3 - 87.9078 * x2 + 182.0935 * x + 20.1025;
			break;
		case 78: y = 23.9066 * x3 - 87.2079 * x2 + 180.6453 * x + 19.8552;
			break;
		case 74: y = 23.6234 * x3 - 86.3052 * x2 + 179.3353 * x + 19.6866;
			break;
		case 70: y = 23.4863 * x3 - 85.7562 * x2 + 178.5390 * x + 19.5130;
			break;

		}
	}

	return y;
}


float CalculateGasRho(float ch4,float x)
{
	float RegressionFormulaOfGasRho(short int CH4, float x);
	float y1, y2, y = 0;
	short int CH4, CH5;
	CH4 = 98 - floor((98-ch4)/4)*4;
	y1 = RegressionFormulaOfGasRho(CH4, x);

	CH5 = 98 - ceil((98-ch4)/4)*4;
	y2 = RegressionFormulaOfGasRho(CH5, x);

	if (CH4==CH5)
	{
		y=y1;
	}
	else
		y=y1+(y2-y1)*(ch4-CH4)/(CH5-CH4); 
	return y;

}

float RegressionFormulaOfGasRho(short int CH4,float x)
{
	float y = 0; 
	float x2, x3; //,x4,x5,x6;
	x2 = x * x;
	x3 = x2 * x;

	if (x <= 0.5)
	{
		switch (CH4)
		{
		case 98: y = 1.5037 * x3 - 1.2647 * x2 + 13.9240 * x + 1.4903;
			break;
		case 94: y = 1.5800 * x3 - 1.6362 * x2 + 13.2188 * x + 1.4277;
			break;
		case 90: y = 1.6362 * x3 - 1.8383 * x2 + 13.0824 * x + 1.4197;
			break;
		case 86: y = 1.6862 * x3 - 1.9767 * x2 + 13.1128 * x + 1.4276;
			break;
		case 82: y = 1.7397 * x3 - 2.0879 * x2 + 13.2236 * x + 1.4431;
			break;
		case 78: y = 1.7912 * x3 - 2.1799 * x2 + 13.3805 * x + 1.4629;
			break;
		case 74: y = 1.8404 * x3 - 2.2581 * x2 + 13.5671 * x + 1.4856;
			break;
		case 70: y = 1.8925 * x3 - 2.3297 * x2 + 13.7752 * x + 1.5102;
			break;


		}
	}
	else
	{
		switch (CH4)
		{
		case 98: y = 0.2936 * x3 + 0.1331 * x2 + 13.3682 * x + 1.5685;
			break;
		case 94: y = 0.2931 * x3 - 0.1545 * x2 + 12.6338 * x + 1.5089;
			break;
		case 90: y = 0.2746 * x3 - 0.2611 * x2 + 12.4535 * x + 1.5083;
			break;
		case 86: y = 0.2685 * x3 - 0.3282 * x2 + 12.4520 * x + 1.5214;
			break;
		case 82: y = 0.2749 * x3 - 0.3904 * x2 + 12.5471 * x + 1.5382;
			break;
		case 78: y = 0.2719 * x3 - 0.4115 * x2 + 12.6704 * x + 1.5639;
			break;
		case 74: y = 0.2788 * x3 - 0.4425 * x2 + 12.8398 * x + 1.5886;
			break;
		case 70: y = 0.2885 * x3 - 0.4690 * x2 + 13.0324 * x + 1.6149;
			break;

		}
	}

	return y;
}


float CalculateGasH(float ch4,float x)
{
	float RegressionFormulaOfGasH(short int CH4, float x);
	float y1, y2, y = 0;
	short int CH4, CH5;

	CH4 = 98 - floor((98-ch4)/4)*4;
	y1 = RegressionFormulaOfGasH(CH4, x);

	CH5 = 98 - ceil((98-ch4)/4)*4;
	y2 = RegressionFormulaOfGasH(CH5, x);

	if (CH4==CH5)
	{
		y=y1;
	}
	else
		y=y1+(y2-y1)*(ch4-CH4)/(CH5-CH4); 
	return y;

}

float RegressionFormulaOfGasH(short int CH4,float x)
{
	float y = 0; 
	float x2, x3; //,x4,x5,x6;
	x2 = x * x;
	x3 = x2 * x;

	if (x <= 0.5)
	{
		switch (CH4)
		{
		case 98: y = 198.79 * x3 - 245.75 * x2 + 123.85 * x + 561.67;
			break;
		case 94: y = 225.06 * x3 - 279.39 * x2 + 146.15 * x + 582.90;
			break;
		case 90: y = 236.53 * x3 - 294.02 * x2 + 157.09 * x + 593.51;
			break;
		case 86: y = 243.06 * x3 - 302.41 * x2 + 163.76 * x + 600.52;
			break;
		case 82: y = 246.96 * x3 - 307.52 * x2 + 168.13 * x + 605.66;
			break;
		case 78: y = 249.25 * x3 - 310.63 * x2 + 171.09 * x + 609.62;
			break;
		case 74: y = 250.75 * x3 - 312.64 * x2 + 173.14 * x + 612.75;
			break;
		case 70: y = 251.28 * x3 - 313.53 * x2 + 174.48 * x + 615.26;
			break;

		}
	}
	else
	{
		switch (CH4)
		{
		case 98: y = 10.64 * x3 - 39.70 * x2 + 47.50 * x + 571.55;
			break;
		case 94: y = 12.69 * x3 - 46.95 * x2 + 60.13 * x + 594.01;
			break;
		case 90: y = 13.47 * x3 - 49.74 * x2 + 66.62 * x + 605.20;
			break;
		case 86: y = 13.87 * x3 - 51.24 * x2 + 70.62 * x + 612.58;
			break;
		case 82: y = 14.20 * x3 - 52.38 * x2 + 73.51 * x + 617.91;
			break;
		case 78: y = 14.38 * x3 - 53.07 * x2 + 75.51 * x + 622.00;
			break;
		case 74: y = 14.49 * x3 - 53.52 * x2 + 76.96 * x + 625.22;
			break;
		case 70: y = 14.66 * x3 - 54.02 * x2 + 78.17 * x + 627.74;
			break;


		}
	}

	return y;
}


float CalculateLiquidC(float ch4,float x)
{
	float RegressionFormulaOfLC(short int CH4, float x);
	float y1, y2, y = 0;
	short int CH4, CH5;

	CH4 = 98 - floor((98-ch4)/4)*4;
	y1 = RegressionFormulaOfLC(CH4, x);

	CH5 = 98 - ceil((98-ch4)/4)*4;
	y2 = RegressionFormulaOfLC(CH5, x);

	if (CH4==CH5)
	{
		y=y1;
	}
	else
		y=y1+(y2-y1)*(ch4-CH4)/(CH5-CH4); 
	return y;

}

float RegressionFormulaOfLC(short int CH4,float x)
{
	float y = 0; 
	float x2, x3; //,x4,x5,x6;
	x2 = x * x;
	x3 = x2 * x;

	if (x <= 0.5)
	{
		switch (CH4)
		{
		case 98: y = 0.5366 * x3 - 0.5947 * x2 + 0.7390 * x + 3.4317;
			break;
		case 94: y = 0.4788 * x3 - 0.5426 * x2 + 0.6556 * x + 3.3389;
			break;
		case 90: y = 0.4280 * x3 - 0.4934 * x2 + 0.5850 * x + 3.2543;
			break;
		case 86: y = 0.3810 * x3 - 0.4460 * x2 + 0.5245 * x + 3.1768;
			break;
		case 82: y = 0.3383 * x3 - 0.4012 * x2 + 0.4724 * x + 3.1057;
			break;
		case 78: y = 0.3001 * x3 - 0.3599 * x2 + 0.4273 * x + 3.0400;
			break;
		case 74: y = 0.2651 * x3 - 0.3211 * x2 + 0.3881 * x + 2.9792;
			break;
		case 70: y = 0.2316 * x3 - 0.2839 * x2 + 0.3537 * x + 2.9228;
			break;

		}
	}
	else
	{
		switch (CH4)
		{
		case 98: y = 0.0575 * x3 - 0.0691 * x2 + 0.5441 * x + 3.4569;
			break;
		case 94: y = 0.0461 * x3 - 0.0662 * x2 + 0.4780 * x + 3.3620;
			break;
		case 90: y = 0.0371 * x3 - 0.0604 * x2 + 0.4221 * x + 3.2758;
			break;
		case 86: y = 0.0330 * x3 - 0.0605 * x2 + 0.3796 * x + 3.1959;
			break;
		case 82: y = 0.0285 * x3 - 0.0562 * x2 + 0.3418 * x + 3.1230;
			break;
		case 78: y = 0.0258 * x3 - 0.0538 * x2 + 0.3113 * x + 3.0554;
			break;
		case 74: y = 0.0237 * x3 - 0.0512 * x2 + 0.2856 * x + 2.9929;
			break;
		case 70: y = 0.0207 * x3 - 0.0458 * x2 + 0.2621 * x + 2.9352;
			break;

		}
	}

	return y;
}



double getDistanceFromGPS(double lat_a, double lng_a, double lat_b, double lng_b)
{
	double pk = 3.1415926 /180 ;

	double a1 = lat_a * pk;
	double a2 = lng_a * pk;
	double b1 = lat_b * pk;
	double b2 = lng_b * pk;


	double t1 = cos(a1) * cos(a2) * cos(b1) * cos(b2);
	double t2 = cos(a1) * sin(a2) * cos(b1) * sin(b2);
	double t3 = sin(a1) * sin(b1);
	double d = acos(t1 + t2 + t3);

	d = d * 6378138; 

	return  d;
}