haifeng.wang
/
MS-DTU
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity
MS-DTU/MS-DTU-V1/User/TankPrediction.c

2375 lines
47 KiB
C
Raw Permalink Normal View History

首次提交
2025-04-03 14:18:58 +08:00
#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;
}