ZNY_Pakistan/Anjiehui7_ZNY/User/TankPrediction.c

2375 lines
47 KiB
C
Raw Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

#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;
}