/*************************************************************************
* 文件:	aga10.c	
* 描述:	管理计算声速或C*的整体流程
*       创建并使用基于Detail和Therm类的对象
*       包含以下函数:
*       AGA10_Init(), AGA10_UnInit(), SOS(), Crit(), Cperf(), CRi()
* 版本:	ver 1.7	2002.11.17
* 作者:	W.B. Peterson
* 修订:
* 版权 (c) 2002  美国天然气协会

**************************************************************************/

#include "NGCal.h" 
#include "Therm.h" 
#include "Detail.h"

// 创建文件作用域的对象指针：一个Therm类对象和一个Detail类对象
static Therm *ptTherm;
static Detail *ptDetail;

/**************************************************************************
*	函数	:	AGA10_Init()
*	参数	:	void
*	返回值	:	int
*	功能	:	初始化库；创建所需对象
*	修订	:	
**************************************************************************/

 int NGCal_Init(void)
{
	ptDetail = (Detail *)malloc(sizeof(Detail));
    // 创建用于计算密度的对象
    if (NULL == ptDetail )
    {
        return MEMORY_ALLOCATION_ERROR;
    }

    // 创建用于计算热力学性质的对象
		 ptTherm = (Therm *)malloc(sizeof(Therm));
    if (NULL == ptTherm )
    {
        return MEMORY_ALLOCATION_ERROR;
    }

    return AGA10_INITIALIZED;
}

/**************************************************************************
*	函数	:	AGA10_UnInit()
*	参数	:	void
*	返回值	:	int
*	功能	:	反初始化库；删除对象
*	修订	:	
**************************************************************************/

 int NGCal_UnInit(void)
{
    // 删除对象（如果存在）
    if (ptDetail) free( ptDetail);
    if (ptTherm) free (ptTherm);
    return 0;
}

/**************************************************************************
*	函数	:	SOS()
*	参数	:	指向外部AGA10数据结构的指针
*	返回值	:	double
*	功能	:	计算声速和其他参数
*	修订	:	
**************************************************************************/

 double SOS(AGA10STRUCT *ptAGA10)
{
    // 检查库是否就绪；必要时初始化
    if (NULL == ptDetail || NULL == ptTherm)
    {
        NGCal_UnInit();
        NGCal_Init();
    }

    // 调用函数计算密度和热力学性质
    ptTherm->Run(ptTherm,ptAGA10, ptDetail);

    // 基本声速计算不计算C*；初始化为零
    ptAGA10->dCstar = 0.0;

    // 返回声速给调用者
    return ptAGA10->dSOS;
}

/**************************************************************************
*	函数	:	Crit()
*	参数	:	指向外部AGA10数据结构的指针、Detail和Therm对象以及一个双精度浮点数（管道中的气体速度）
*	返回值	:	double
*	功能	:	计算C*
*	修订	:	
**************************************************************************/

 double Crit(AGA10STRUCT *ptAGA10, double dPlenumVelocity)
{
    // 函数局部变量
    double DH, DDH, S, H;
    double tolerance = 1.0;
    double R, P, T, Z;
    int i;

    // 检查对象是否就绪；如果未就绪则尝试初始化
    if (NULL == ptDetail || NULL == ptTherm)
    {
        NGCal_UnInit();
        if (AGA10_INITIALIZED != NGCal_Init())
        {
            ptAGA10->lStatus = MEMORY_ALLOCATION_ERROR;
            return 0.0;
        }
    }

    // 首先计算密度和热力学性质
    ptTherm->Run(ptTherm,ptAGA10, ptDetail);

    // DH是从管道到喉部的焓变；这是我们的初始猜测值
    DH = (ptAGA10->dSOS * ptAGA10->dSOS - dPlenumVelocity * dPlenumVelocity) / 2.0;

    // 在我们改变数据结构内容之前捕获管道条件
    S = ptAGA10->dS;
    H = ptAGA10->dH;
    R = ptAGA10->dRhof;
    P = ptAGA10->dPf;
    Z = ptAGA10->dZf;
    T = ptAGA10->dTf;

    // 将DH的增量初始化为超出收敛容差的任意值
    DDH = 10.0;

    // 通过简单重复，搜索喷嘴喉部的压力、温度和声速，
    // 在给定管道熵的情况下提供恒定焓。
    // 如果循环执行超过100次仍未收敛，则中止。
    for (i = 1; i < MAX_NUM_OF_ITERATIONS; i++)
    {
        // 计算满足H和S的P和T
        ptTherm->HS_Mode(ptTherm,ptAGA10, ptDetail, H - DH, S, true);

        // 计算新的热力学性质，包括SOS
        ptTherm->Run(ptTherm,ptAGA10, ptDetail);

        // 保存DH用于容差检查
        DDH = DH;

        // 重新计算DH
        DH = (ptAGA10->dSOS * ptAGA10->dSOS - dPlenumVelocity * dPlenumVelocity) / 2.0;

        // 如果达到容差则结束循环
        if (fabs(DDH - DH) < tolerance) break;
    }

    // C*是真实气体临界流常数（不要与Cperf或CRi混淆）
    ptAGA10->dCstar = (ptAGA10->dRhof * ptAGA10->dSOS) / sqrt(R * P * Z);

    // 恢复原始管道压力和温度
    ptAGA10->dPf = P;
    ptAGA10->dTf = T;

    // 将流体属性恢复为管道条件
    ptTherm->Run(ptTherm,ptAGA10, ptDetail);

    // 返回临界流函数给调用者
    return ptAGA10->dCstar;
}

/**************************************************************************
*	函数	:	Cperf()
*	参数	:	指向外部AGA10数据结构的指针
*	返回值	:	double
*	功能	:	计算等熵完美气体临界流函数
*	修订	:	
**************************************************************************/

double Cperf(AGA10STRUCT *ptAGA10)
{
    double k, root, exponent;
    k = ptAGA10->dKappa;
    root = 2.0 / (k + 1.0);
    exponent = (k + 1.0) / (k - 1.0);

    // 等熵完美气体临界流函数C*i
    return(sqrt(k * pow(root, exponent)));
}

/**************************************************************************
*	函数	:	CRi()
*	参数	:	指向外部AGA10数据结构的指针
*	返回值	:	double
*	功能	:	计算等熵真实气体临界流函数CRi
*	修订	:	
**************************************************************************/

double CRi(AGA10STRUCT *ptAGA10)
{
    return (Cperf(ptAGA10) / sqrt(ptAGA10->dZf));
}