Initial commit from Polytechnique Montreal

1 files changed, 278 insertions, 0 deletions

diff --git a/Utilib/src/region2.c b/Utilib/src/region2.c
new file mode 100644
index 0000000..8ae1c33
--- /dev/null
+++ b/Utilib/src/region2.c
@@ -0,0 +1,278 @@
+/*
+freesteam - IAPWS-IF97 steam tables library
+Copyright (C) 2004-2009  John Pye
+
+This program is free software; you can redistribute it and/or
+modify it under the terms of the GNU General Public License
+as published by the Free Software Foundation; either version 2
+of the License, or (at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with this program; if not, write to the Free Software
+Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
+*/
+
+#define FREESTEAM_BUILDING_LIB
+#include "region.h"
+
+#define GAM0(PI,TAU) gam0(PI,TAU)
+#define GAM0PI(PI,TAU) (1./PI)
+#define GAM0PIPI(PI,TAU) (-1./SQ(PI))
+#define GAM0TAU(PI,TAU) gam0tau(TAU)
+#define GAM0TAUTAU(PI,TAU) gam0tautau(TAU)
+#define GAM0PITAU(PI,TAU) (0)
+
+#define gam(PI,TAU) (GAM0(PI,TAU) + gamr(PI,TAU))
+#define gampi(PI,TAU) (GAM0PI(PI,TAU) + gamrpi(PI,TAU))
+#define gampipi(PI,TAU) (GAM0PIPI(PI,TAU) + gamrpipi(PI,TAU))
+#define gamtau(PI,TAU) (GAM0TAU(PI,TAU) + gamrtau(PI,TAU))
+#define gamtautau(PI,TAU) (GAM0TAUTAU(PI,TAU) + gamrtautau(PI,TAU))
+#define gampitau(PI,TAU) (GAM0PITAU(PI,TAU) + gamrpitau(PI,TAU))
+
+static double gamr(double pi, double tau);
+static double gamrpi(double pi, double tau);
+static double gamrpipi(double pi, double tau);
+static double gamrtau(double pi, double tau);
+static double gamrtautau(double pi, double tau);
+static double gamrpitau(double pi, double tau);
+
+static double gam0(double pi, double tau);
+static double gam0tau(double tau);
+static double gam0tautau(double tau);
+
+#include <math.h>
+#include "common.h"
+
+#define R IAPWS97_R
+
+#define REGION2_GPT_PSTAR 1.e6 /* Pa */
+#define REGION2_GPT_TSTAR 540. /* K */
+
+#define DEFINE_PITAU(P,T) \
+	double pi = p / REGION2_GPT_PSTAR; \
+	double tau = REGION2_GPT_TSTAR / T
+
+
+double freesteam_region2_v_pT(double p, double T){
+	DEFINE_PITAU(p,T);
+	return (R * T / p) * pi * gampi(pi,tau);
+}
+
+double freesteam_region2_u_pT(double p, double T){
+	DEFINE_PITAU(p,T);
+	return (R * T) * (tau * gamtau(pi,tau) - pi * gampi(pi,tau));
+}
+
+double freesteam_region2_s_pT(double p, double T){
+	DEFINE_PITAU(p,T);
+	return R * (tau * gamtau(pi,tau) - gam(pi,tau));
+}
+
+double freesteam_region2_h_pT(double p, double T){
+	DEFINE_PITAU(p,T);
+	return R * T * (tau * gamtau(pi,tau));
+}
+
+double freesteam_region2_cp_pT(double p, double T){
+	DEFINE_PITAU(p,T);
+	return R * (-SQ(tau) * gamtautau(pi,tau));
+}
+
+double freesteam_region2_cv_pT(double p, double T){
+	DEFINE_PITAU(p,T);
+	return R * (-SQ(tau) * gamtautau(pi,tau) + SQ(gampi(pi,tau) -
+		tau * gampitau(pi,tau)) / gampipi(pi,tau)
+	);
+}
+
+double freesteam_region2_w_pT(double p, double T){
+        double gp;
+	DEFINE_PITAU(p,T);
+	gp = gamrpi(pi,tau);
+	return sqrt(R * T * (1. + 2.*pi*gp+SQ(pi*gp))/
+		((1. - SQ(pi)*gamrpipi(pi,tau))	+ SQ(1. + pi*gp - tau*pi*gamrpitau(pi,tau))/SQ(tau)/gamtautau(pi,tau))
+	);
+}
+
+double freesteam_region2_g_pT(double p, double T){
+	DEFINE_PITAU(p,T);
+	return R * T * gam(pi,tau);
+}	
+
+double freesteam_region2_a_pT(double p, double T){
+	DEFINE_PITAU(p,T);
+	return R * T * (gam(pi,tau) - gampi(pi,tau) * pi);
+}
+
+double freesteam_region2_alphav_pT(double p, double T){
+        double pigamrpi, alphav;
+	DEFINE_PITAU(p,T);
+	pigamrpi = pi*gamrpi(pi,tau);
+	alphav = 1./T * (1. + pigamrpi - tau*pi*gamrpitau(pi,tau))/(1. + pigamrpi);
+	return alphav;
+}
+
+double freesteam_region2_kappaT_pT(double p, double T){
+        double kappaT;
+	DEFINE_PITAU(p,T);
+	kappaT = 1./p * (1.-SQ(pi)*gamrpipi(pi,tau)) / (1.+pi*gamrpi(pi,tau));
+	return kappaT;
+}
+
+/*------------------------------------------------------------------------------
+  REGION 2 IDEAL PART - GAM0(PI,TAU)
+*/
+
+typedef struct{
+	int J;
+	double n;
+} JNData;
+
+const JNData REGION2_GPT_IDEAL_DATA[] = {
+	{0,	-0.96927686500217E+01}
+	,{1,	0.10086655968018E+02}
+	,{-5,	-0.56087911283020E-02}
+	,{-4,	0.71452738081455E-01}
+	,{-3,	-0.40710498223928E+00}
+	,{-2,	0.14240819171444E+01}
+	,{-1,	-0.43839511319450E+01}
+	,{2,	-0.28408632460772E+00}
+	,{3,	0.21268463753307E-01}
+};
+
+const unsigned REGION2_GPT_IDEAL_MAX = sizeof(REGION2_GPT_IDEAL_DATA)/sizeof(JNData);
+
+#define REGION2_GPT_IDEAL_LOOP \
+	double sum = 0; \
+	const JNData *d, *e = REGION2_GPT_IDEAL_DATA + REGION2_GPT_IDEAL_MAX; \
+	for(d = REGION2_GPT_IDEAL_DATA; d < e; ++d)
+
+double gam0(double pi, double tau){
+	REGION2_GPT_IDEAL_LOOP{
+		sum += d->n * ipow(tau, d->J);
+	}
+	return log(pi) + sum;
+}
+
+double gam0tau(double tau){
+	REGION2_GPT_IDEAL_LOOP{
+		sum += d->n * d->J * ipow(tau, d->J - 1);
+	}
+	return sum;
+}
+
+double gam0tautau(double tau){
+	REGION2_GPT_IDEAL_LOOP{
+		sum += d->n * d->J * (d->J - 1) * ipow(tau, d->J - 2);
+	}
+	return sum;
+}
+
+/*------------------------------------------------------------------------------
+  REGION 2 RESIDUAL PART - GAMR(PI,TAU)
+*/
+
+typedef struct{
+	int I, J;
+	double n;
+} IJNData;
+
+const IJNData REGION2_GPT_RESID_DATA[] = {
+	{1,	0,	-0.17731742473213E-02}
+	,{1,	1,	-0.17834862292358E-01}
+	,{1,	2,	-0.45996013696365E-01}
+	,{1,	3,	-0.57581259083432E-01}
+	,{1,	6,	-0.50325278727930E-01}
+	,{2,	1,	-0.33032641670203E-04}
+	,{2,	2,	-0.18948987516315E-03}
+	,{2,	4,	-0.39392777243355E-02}
+	,{2,	7,	-0.43797295650573E-01}
+	,{2,	36,	-0.26674547914087E-04}
+	,{3,	0,	0.20481737692309E-07}
+	,{3,	1,	0.43870667284435E-06}
+	,{3,	3,	-0.32277677238570E-04}
+	,{3,	6,	-0.15033924542148E-02}
+	,{3,	35,	-0.40668253562649E-01}
+	,{4,	1,	-0.78847309559367E-09}
+	,{4,	2,	0.12790717852285E-07}
+	,{4,	3,	0.48225372718507E-06}
+	,{5,	7,	0.22922076337661E-05}
+	,{6,	3,	-0.16714766451061E-10}
+	,{6,	16,	-0.21171472321355E-02}
+	,{6,	35,	-0.23895741934104E+02}
+	,{7,	0,	-0.59059564324270E-17}
+	,{7,	11,	-0.12621808899101E-05}
+	,{7,	25,	-0.38946842435739E-01}
+	,{8,	8,	0.11256211360459E-10}
+	,{8,	36,	-0.82311340897998E+01}
+	,{9,	13,	0.19809712802088E-07}
+	,{10,	4,	0.10406965210174E-18}
+	,{10,	10,	-0.10234747095929E-12}
+	,{10,	14,	-0.10018179379511E-08}
+	,{16,	29,	-0.80882908646985E-10}
+	,{16,	50,	0.10693031879409E+00}
+	,{18,	57,	-0.33662250574171E+00}
+	,{20,	20,	0.89185845355421E-24}
+	,{20,	35,	0.30629316876232E-12}
+	,{20,	48,	-0.42002467698208E-05}
+	,{21,	21,	-0.59056029685639E-25}
+	,{22,	53,	0.37826947613457E-05}
+	,{23,	39,	-0.12768608934681E-14}
+	,{24,	26,	0.73087610595061E-28}
+	,{24,	40,	0.55414715350778E-16}
+	,{24,	58,	-0.94369707241210E-06}
+};
+
+const unsigned REGION2_GPT_RESID_MAX = sizeof(REGION2_GPT_RESID_DATA)/sizeof(IJNData);
+
+#define REGION2_GPT_RESID_LOOP \
+	double sum = 0; \
+	const IJNData *d, *e = REGION2_GPT_RESID_DATA + REGION2_GPT_RESID_MAX; \
+	for(d = REGION2_GPT_RESID_DATA; d < e; ++d)
+
+double gamr(double pi, double tau){
+	REGION2_GPT_RESID_LOOP{
+		sum += d->n * ipow(pi,d->I) * ipow(tau - 0.5,d->J);
+	}
+	return sum;
+}
+
+double gamrpi(double pi, double tau){
+	REGION2_GPT_RESID_LOOP{
+		sum +=  d->n * d->I * ipow(pi,d->I -1) * ipow(tau - 0.5,d->J);
+	}
+	return sum;
+}
+
+double gamrpipi(double pi, double tau){
+	REGION2_GPT_RESID_LOOP{
+		sum += d->n * d->I * (d->I - 1) * ipow(pi, d->I - 2) * ipow(tau - 0.5, d->J);
+    }
+	return sum;
+}
+
+double gamrtau(double pi, double tau){
+	REGION2_GPT_RESID_LOOP{
+		sum += d->n * ipow(pi, d->I) * d->J * ipow(tau - 0.5, d->J - 1);
+	}
+	return sum;
+}
+
+double gamrtautau(double pi, double tau){
+	REGION2_GPT_RESID_LOOP{
+		sum += d->n * ipow(pi, d->I) * d->J * (d->J - 1) * ipow(tau - 0.5, d->J - 2);
+	}
+	return sum;
+}
+
+double gamrpitau(double pi, double tau){
+	REGION2_GPT_RESID_LOOP{
+		sum += d->n * d->I * ipow(pi, d->I - 1) * d->J * ipow(tau - 0.5, d->J - 1);
+	}
+	return sum;
+}