/*************************************************************************
* File:	aga10.cpp	
* Description:	Manages overall process of calculating speed of sound
*	or C*; creates and uses objects based on Detail and Therm classes
*	Contains the following functions:
*	AGA10_Init(), AGA10_UnInit(), SOS(), Crit(), Cperf(), CRi()
* Version:	ver 1.7	2002.11.17
* Author:	W.B. Peterson
*Revisions:
*Copyright (c) 2002  American Gas Association

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

#include "aga10.h" #include "therm.h" #include "detail.h"  #include "math.h"


//Create file-scope pointers to objects we will need; one of Therm class

//and one of Detail class.

static Therm *ptTherm ; static Detail *ptDetail ;

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

*	Function	:	AGA10_Init()
*	Arguments	:	void
*	Returns	:	int
*	Purpose	:	Initializes library; creates required objects
*	Revisions	:	
**************************************************************************/

 int AGA10_Init(void)

{
//create object for calculating density if (NULL == (ptDetail = new Detail))
{
return MEMORY_ALLOCATION_ERROR ;
}

//create object for calculating thermodynamic properties
if (NULL == (ptTherm = new Therm))

{
return MEMORY_ALLOCATION_ERROR ;
}

return AGA10_INITIALIZED ;

}// AGA10_Init


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

*	Function	:	AGA10_UnInit()
*	Arguments	:	void
*	Returns	:	int
*	Purpose	:	Un-initializes library; deletes objects
*	Revisions	:	
**************************************************************************/

 int AGA10_UnInit(void)

{
// delete the objects (if they exist) if (ptDetail) delete ptDetail ;
if (ptTherm) delete ptTherm ; return 0 ;

}// AGA10_UnInit


/**************************************************************************
*	Function	:	SOS()
*	Arguments	:	Pointers to external AGA10 data struct
*	Returns	:	double
*	Purpose	:	calculates speed of sound and other parameters
*	Revisions	:	
**************************************************************************/

 double SOS(AGA10STRUCT *ptAGA10)

{
// check if library is ready; initialize if necessary if (NULL == ptDetail || NULL == ptTherm)

{
AGA10_UnInit() ;
AGA10_Init() ;

}

//Call function to calculate densities and thermodynamic properties ptTherm->Run(ptAGA10, ptDetail) ;

//the basic sound speed calculation doesn't calculate C*; initialize to zero ptAGA10->dCstar = 0.0 ;

//return the speed of sound to caller

return ptAGA10->dSOS ;

}// VOS()


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

*	Function	:	Crit()
*	Arguments	:	Pointers to external AGA10 data struct, Detail and Therm
*			objects and a double precision float (gas velocity in plenum)
*	Returns	:	double
*	Purpose	:	calculates C*
*	Revisions	:	
**************************************************************************/

 double Crit(AGA10STRUCT *ptAGA10, double dPlenumVelocity)

{
//variables local to function double DH, DDH, S, H;
double tolerance = 1.0 ; double R, P, T, Z ;

int i ;

//check objects for readiness; try to initialize if not if (NULL == ptDetail || NULL == ptTherm)

{
AGA10_UnInit() ;

if (AGA10_INITIALIZED != AGA10_Init())
{

ptAGA10->lStatus = MEMORY_ALLOCATION_ERROR ; return 0.0 ;

}
}

//begin by calculating densities and thermodynamic properties ptTherm->Run(ptAGA10, ptDetail) ;

//DH is enthalpy change from plenum to throat; this is our initial guess

DH = (ptAGA10->dSOS * ptAGA10->dSOS - dPlenumVelocity * dPlenumVelocity)	/ 2.0 ;

//trap plenum conditions before we alter the data stucture's contents S = ptAGA10->dS ;

H = ptAGA10->dH ;

R = ptAGA10->dRhof ; P = ptAGA10->dPf ; Z = ptAGA10->dZf ; T = ptAGA10->dTf ;

//initialize delta of DH to an arbitrary value outside of

//convergence tolerance
DDH = 10.0 ;

//Via simple repetition, search for a pressure, temperature and sound speed

//at a nozzle throat which provide constant enthalpy, given the entropy known
//at the plenum. Abort if loop executes more than 100 times without convergence. for (i = 1; i < MAX_NUM_OF_ITERATIONS; i++)

{

// calculate P and T to satisfy H and S ptTherm->HS_Mode(ptAGA10, ptDetail, H - DH, S, true) ;

//calculate new thermo, including SOS ptTherm->Run(ptAGA10, ptDetail) ;

//hold DH for tolerance check

DDH = DH ;

// recalculate DH

DH = (ptAGA10->dSOS * ptAGA10->dSOS - dPlenumVelocity * dPlenumVelocity) / 2.0 ;

// end loop if tolerance reached

if (fabs(DDH - DH) < tolerance) break ;
}

//C* is the real gas critical flow constant (not to be confused with Cperf or CRi) ptAGA10->dCstar = (ptAGA10->dRhof * ptAGA10->dSOS) / sqrt(R * P * Z) ;

//put the original plenum pressure and temperature back

ptAGA10->dPf = P ; ptAGA10->dTf = T ;

//restore fluid props to plenum conditions ptTherm->Run(ptAGA10, ptDetail) ;

//return the critical flow function to caller return ptAGA10->dCstar ;

}// Crit()



/**************************************************************************
*	Function	:	Cperf()
*	Arguments	:	pointer to external AGA10 data struct
*	Returns	:	double
*	Purpose	:	calculates isentropic perfect gas critical flow function
*	Revisions	:	
**************************************************************************/

double Cperf(AGA10STRUCT *ptAGA10)

{

double k, root, exponent ;

k = ptAGA10->dKappa ; root = 2.0 / (k + 1.0) ;
exponent = (k + 1.0) / (k - 1.0) ;

// isentropic perfect gas critical flow function C*i return(sqrt(k * pow(root, exponent))) ;

}// Cperf
/**************************************************************************
*	Function	:	CRi()
*	Arguments	:	pointer to external AGA10 data struct
*	Returns	:	double
*	Purpose	:	calculates isentropic real gas critical flow function CRi
*	Revisions	:	
**************************************************************************/


double CRi(AGA10STRUCT *ptAGA10)

{

return (Cperf(ptAGA10) / sqrt(ptAGA10->dZf)) ;
}// CRi()