johnpye/fprops/helmholtz.c

/*  ASCEND modelling environment
    Copyright (C) 2008 Carnegie Mellon University

    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, 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.
*//** @file
    Implementation of the reduced molar Helmholtz free energy equation of state.

    For nomenclature see Tillner-Roth, Harms-Watzenberg and Baehr, Eine neue
    Fundamentalgleichung f端r Ammoniak.

    John Pye, 29 Jul 2008.
*/

#include <math.h>

#include "helmholtz.h"
#include "ideal_impl.h"

#ifdef TEST
#include <assert.h>
#include <stdlib.h>
#include <stdio.h>
#endif

#define SQ(X) ((X)*(X))

/* forward decls */

static double helm_resid(double tau, double delta, const HelmholtzData *data);
static double helm_resid_del(double tau, double delta, const HelmholtzData *data);
static double helm_resid_tau(double tau, double delta, const HelmholtzData *data);
static double helm_resid_deltau(double tau, double delta, const HelmholtzData *data);
static double helm_resid_deldel(double tau, double delta, const HelmholtzData *data);

/**
    Function to calculate pressure from Helmholtz free energy EOS, given temperature
    and mass density.

    @param T temperature in K
    @param rho mass density in kg/m続
    @return pressure in Pa???
*/
double helmholtz_p(double T, double rho, const HelmholtzData *data){
    
    double tau = data->T_star / T;
    double delta = rho / data->rho_star;

#ifdef TEST
    assert(data->rho_star!=0);
    assert(T!=0);
    assert(!isnan(tau));
    assert(!isnan(delta));
    assert(!isnan(data->R));

    //fprintf(stderr,"p calc: T = %f\n",T);
    //fprintf(stderr,"p calc: tau = %f\n",tau);
    //fprintf(stderr,"p calc: rho = %f\n",rho);
    //fprintf(stderr,"p calc: delta = %f\n",delta);
    //fprintf(stderr,"p calc: R*T*rho = %f\n",data->R * T * rho);

    //fprintf(stderr,"T = %f\n", T);
    //fprintf(stderr,"rhob = %f, rhob* = %f, delta = %f\n", rho/data->M, data->rho_star/data->M, delta);
#endif
    
    return data->R * T * rho * (1 + delta * helm_resid_del(tau,delta,data));
}

/**
    Function to calculate internal energy from Helmholtz free energy EOS, given
    temperature and mass density.

    @param T temperature in K
    @param rho mass density in kg/m続
    @return internal energy in ???
*/
double helmholtz_u(double T, double rho, const HelmholtzData *data){
    
    double tau = data->T_star / T;
    double delta = rho / data->rho_star;

#ifdef TEST
    assert(data->rho_star!=0);
    assert(T!=0);
    assert(!isnan(tau));
    assert(!isnan(delta));
    assert(!isnan(data->R));
#endif

#if 0
    fprintf(stderr,"ideal_tau = %f\n",helm_ideal_tau(tau,delta,data->ideal));
    fprintf(stderr,"resid_tau = %f\n",helm_resid_tau(tau,delta,data));
    fprintf(stderr,"R T = %f\n",data->R * data->T_star);
#endif

    return data->R * data->T_star * (helm_ideal_tau(tau,delta,data->ideal) + helm_resid_tau(tau,delta,data));
}

/**
    Function to calculate enthalpy from Helmholtz free energy EOS, given
    temperature and mass density.

    @param T temperature in K
    @param rho mass density in kg/m続
    @return enthalpy in J/kg
*/
double helmholtz_h(double T, double rho, const HelmholtzData *data){
    
    double tau = data->T_star / T;
    double delta = rho / data->rho_star;

#ifdef TEST
    assert(data->rho_star!=0);
    assert(T!=0);
    assert(!isnan(tau));
    assert(!isnan(delta));
    assert(!isnan(data->R));
#endif

    return data->R * T * (1 + tau * (helm_ideal_tau(tau,delta,data->ideal) + helm_resid_tau(tau,delta,data)) + delta*helm_resid_del(tau,delta,data));
}

/**
    Function to calculate entropy from Helmholtz free energy EOS, given
    temperature and mass density.

    @param T temperature in K
    @param rho mass density in kg/m続
    @return entropy in J/kgK
*/
double helmholtz_s(double T, double rho, const HelmholtzData *data){
    
    double tau = data->T_star / T;
    double delta = rho / data->rho_star;

#ifdef ENTROPY_DEBUG
    assert(data->rho_star!=0);
    assert(T!=0);
    assert(!isnan(tau));
    assert(!isnan(delta));
    assert(!isnan(data->R));

    fprintf(stderr,"helm_ideal_tau = %f\n",helm_ideal_tau(tau,delta,data->ideal));
    fprintf(stderr,"helm_resid_tau = %f\n",helm_resid_tau(tau,delta,data));
    fprintf(stderr,"helm_ideal = %f\n",helm_ideal(tau,delta,data->ideal));
    fprintf(stderr,"helm_resid = %f\n",helm_resid(tau,delta,data));
#endif
    return data->R * (
        tau * (helm_ideal_tau(tau,delta,data->ideal) + helm_resid_tau(tau,delta,data))
        - (helm_ideal(tau,delta,data->ideal) + helm_resid(tau,delta,data))
    );
}

/**
    Function to calculate Helmholtz energy from the Helmholtz free energy EOS,
    given temperature and mass density.

    @param T temperature in K
    @param rho mass density in kg/m続
    @return Helmholtz energy 'a', in J/kg
*/
double helmholtz_a(double T, double rho, const HelmholtzData *data){

    double tau = data->T_star / T;
    double delta = rho / data->rho_star;

#ifdef TEST
    assert(data->rho_star!=0);
    assert(T!=0);
    assert(!isnan(tau));
    assert(!isnan(delta));
    assert(!isnan(data->R));
#endif

#ifdef HELMHOLTZ_DEBUG
    fprintf(stderr,"helmholtz_a: T = %f, rho = %f\n",T,rho);
    fprintf(stderr,"multiplying by RT = %f\n",data->R*T);
#endif

    return data->R * T * (helm_ideal(tau,delta,data->ideal) + helm_resid(tau,delta,data));
}


/**
    Calculation zero-pressure specific heat capacity
*/
double helmholtz_cp0(double T, const HelmholtzData *data){
    double val = helm_cp0(T,data->ideal);
#if 0
    fprintf(stderr,"val = %f\n",val);
#endif
    return val;
}

/**
    Calculate partial derivative of p with respect to T, with rho constant
*/
double helmholtz_dpdT_rho(double T, double rho, const HelmholtzData *data){
    double tau = data->T_star / T;
    double delta = rho / data->rho_star;

    double phir_del = helm_resid_del(tau,delta,data);
    double phir_deltau = helm_resid_deltau(tau,delta,data);
    assert(!isinf(phir_del));
    assert(!isinf(phir_deltau));
    assert(!isnan(phir_del));
    assert(!isnan(phir_deltau));
    assert(!isnan(data->R));
    assert(!isnan(rho));
    assert(!isnan(tau));

    double res =  data->R * rho * (1 + delta*phir_del - delta*tau*phir_deltau);

    assert(!isnan(res));
    assert(!isinf(res));
    return res;
}

/**
    Calculate partial derivative of p with respect to rho, with T constant
*/
double helmholtz_dpdrho_T(double T, double rho, const HelmholtzData *data){
    double tau = data->T_star / T;
    double delta = rho / data->rho_star;

    double phir_del = helm_resid_del(tau,delta,data);
    double phir_deldel = helm_resid_deldel(tau,delta,data);
    assert(!isinf(phir_del));
    assert(!isinf(phir_deldel));
    
    return data->R * T * (1 + 2*delta*phir_del + delta*delta* phir_deldel);
}

/*---------------------------------------------
  UTILITY FUNCTION(S)
*/

/* ipow:  public domain by Mark Stephen with suggestions by Keiichi Nakasato */
static double ipow(double x, int n){
    double t = 1.0;

    if(!n)return 1.0;    /* At the top. x^0 = 1 */

    if (n < 0){
        n = -n;
        x = 1.0/x;  /* error if x == 0. Good                        */
    }                 /* ZTC/SC returns inf, which is even better     */

    if (x == 0.0)return 0.0;

    do{
        if(n & 1)t *= x;
        n /= 2;     /* KN prefers if (n/=2) x*=x; This avoids an    */
        x *= x;     /* unnecessary but benign multiplication on     */
    }while(n);      /* the last pass, but the comparison is always
                       true _except_ on the last pass. */

    return t; 
}

//#define RESID_DEBUG

/**
    Residual part of helmholtz function.
*/
double helm_resid(double tau, double delta, const HelmholtzData *data){
    double dell,ldell, term, sum, res = 0;
    unsigned n, i;
    const HelmholtzPowTerm *pt;
    const HelmholtzGausTerm *gt;

    n = data->np;
    pt = &(data->pt[0]);

#ifdef RESID_DEBUG
        fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
#endif

    /* power terms */
    sum = 0;
    dell = ipow(delta,pt->l);
    ldell = pt->l * dell;
    unsigned oldl;
    for(i=0; i<n; ++i){
        term = pt->a * pow(tau, pt->t) * ipow(delta, pt->d);
        sum += term;
#ifdef RESID_DEBUG
        fprintf(stderr,"i = %d,               a=%e, t=%f, d=%d, term = %f, sum = %f",i,pt->a,pt->t,pt->d,term,sum);
        if(pt->l==0){
            fprintf(stderr,",row=%e\n",term);
        }else{
            fprintf(stderr,",row=%e\n,",term*exp(-dell));
        }
#endif
        oldl = pt->l;
        ++pt;
        if(i+1==n || oldl != pt->l){
            if(oldl == 0){
#ifdef RESID_DEBUG
                fprintf(stderr,"linear ");
#endif
                res += sum;
            }else{
#ifdef RESID_DEBUG
                fprintf(stderr,"exp dell=%f, exp(-dell)=%f sum=%f: ",dell,exp(-dell),sum);
#endif
                res += sum * exp(-dell);
            }
#ifdef RESID_DEBUG
            fprintf(stderr,"i = %d, res = %f\n",i,res);
#endif
            sum = 0;
            dell = ipow(delta,pt->l);
            ldell = pt->l*dell;
        }
    }

    /* gaussian terms */
    n = data->ng;
    //fprintf(stderr,"THERE ARE %d GAUSSIAN TERMS\n",n);
    gt = &(data->gt[0]);
    for(i=0; i<n; ++i){
#ifdef RESID_DEBUG
        fprintf(stderr,"i = %d, GAUSSIAN, n = %e, t = %f, d = %f, alpha = %f, beta = %f, gamma = %f, epsilon = %f\n",i+1, gt->n, gt->t, gt->d, gt->alpha, gt->beta, gt->gamma, gt->epsilon);
#endif
        double d1 = delta - gt->epsilon;
        double t1 = tau - gt->gamma;
        double e1 = -gt->alpha*d1*d1 - gt->beta*t1*t1;
        sum = gt->n * pow(tau,gt->t) * pow(delta,gt->d) * exp(e1);
        //fprintf(stderr,"sum = %f\n",sum);
        res += sum;
        ++gt;
    }

#ifdef RESID_DEBUG
    fprintf(stderr,"phir = %f\n",res);
#endif
    return res;
}

/**
    Derivative of the helmholtz residual function with respect to
    delta.
*/  
double helm_resid_del(double tau,double delta, const HelmholtzData *data){
    double sum = 0, res = 0;
    double dell, ldell;
    unsigned n, i;
    const HelmholtzPowTerm *pt;
    const HelmholtzGausTerm *gt;


#ifdef RESID_DEBUG
        fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
#endif

    /* power terms */
    n = data->np;
    pt = &(data->pt[0]);
    dell = ipow(delta,pt->l);
    ldell = pt->l * dell;
    unsigned oldl;
    for(i=0; i<n; ++i){
        sum += pt->a * pow(tau, pt->t) * ipow(delta, pt->d - 1) * (pt->d - ldell);
        oldl = pt->l;
        ++pt;
        if(i+1==n || oldl != pt->l){
            if(oldl == 0){
                res += sum;
            }else{
                res += sum * exp(-dell);
            }
            sum = 0;
            dell = ipow(delta,pt->l);
            ldell = pt->l*dell;
        }
    }

    /* gaussian terms */
    n = data->ng;
    //fprintf(stderr,"THERE ARE %d GAUSSIAN TERMS\n",n);
    gt = &(data->gt[0]);
    for(i=0; i<n; ++i){
#ifdef RESID_DEBUG
        fprintf(stderr,"i = %d, GAUSSIAN, n = %e, t = %f, d = %f, alpha = %f, beta = %f, gamma = %f, epsilon = %f\n",i+1, gt->n, gt->t, gt->d, gt->alpha, gt->beta, gt->gamma, gt->epsilon);
#endif
        double val2;
        val2 = - gt->n * pow(tau,gt->t) * pow(delta, -1. + gt->d)
            * (2. * gt->alpha * delta * (delta - gt->epsilon) - gt->d)
            * exp(-(gt->alpha * SQ(delta-gt->epsilon) + gt->beta*SQ(tau-gt->gamma)));
        res += val2;
#ifdef RESID_DEBUG
        fprintf(stderr,"val2 = %f --> res = %f\n",val2,res);
#endif
        ++gt;
    }

    return res;
}

/**
    Derivative of the helmholtz residual function with respect to
    tau.
*/          
double helm_resid_tau(double tau,double delta,const HelmholtzData *data){
    
    double sum;
    double res = 0;
    double delX;
    unsigned l;
    unsigned n, i;
    const HelmholtzPowTerm *pt;
    const HelmholtzGausTerm *gt;

    n = data->np;
    pt = &(data->pt[0]);

    delX = 1;

    l = 0;
    sum = 0;
    for(i=0; i<n; ++i){
        if(pt->t){
            //fprintf(stderr,"i = %d, a = %e, t = %f, d = %d, l = %d\n",i+1, pt->a, pt->t, pt->d, pt->l);
            sum += pt->a * pow(tau, pt->t - 1) * ipow(delta, pt->d) * pt->t;
        }
        ++pt;
        //fprintf(stderr,"l = %d\n",l);
        if(i+1==n || l != pt->l){
            if(l==0){
                //fprintf(stderr,"Adding non-exp term\n");
                res += sum;
            }else{
                //fprintf(stderr,"Adding exp term with l = %d, delX = %e\n",l,delX);
                res += sum * exp(-delX);
            }
            /* set l to new value */
            if(i+1!=n){
                l = pt->l;
                //fprintf(stderr,"New l = %d\n",l);
                delX = ipow(delta,l);
                sum = 0;
            }
        }
    }

//#define RESID_DEBUG
    /* gaussian terms */
    n = data->ng;
    gt = &(data->gt[0]);
    for(i=0; i<n; ++i){
#ifdef RESID_DEBUG
        fprintf(stderr,"i = %d, GAUSSIAN, n = %e, t = %f, d = %f, alpha = %f, beta = %f, gamma = %f, epsilon = %f\n",i+1, gt->n, gt->t, gt->d, gt->alpha, gt->beta, gt->gamma, gt->epsilon);
#endif

        double val2;
        val2 = -gt->n * pow(tau,gt->t - 1.) * pow(delta, gt->d)
            * (2. * gt->beta * tau * (tau - gt->gamma) - gt->t)
            * exp(-(gt->alpha * SQ(delta-gt->epsilon) + gt->beta*SQ(tau-gt->gamma)));
        res += val2;
#ifdef RESID_DEBUG
        fprintf(stderr,"res = %f\n",res);
#endif
            
        ++gt;
    }

    return res;
}   


/**
    Mixed derivative of the helmholtz residual function with respect to
    delta and tau.
*/
double helm_resid_deltau(double tau,double delta,const HelmholtzData *data){
    double dell,ldell, term, sum = 0, res = 0;
    unsigned n, i;
    const HelmholtzPowTerm *pt;
    const HelmholtzGausTerm *gt;

    /* power terms */
    n = data->np;
    pt = &(data->pt[0]);
    dell = ipow(delta,pt->l);
    ldell = pt->l * dell;
    unsigned oldl;
    for(i=0; i<n; ++i){
        sum += pt->a * pt->t * pow(tau, pt->t - 1) * ipow(delta, pt->d - 1) * (pt->d - ldell);
        oldl = pt->l;
        ++pt;
        if(i+1==n || oldl != pt->l){
            if(oldl == 0){
                res += sum;
            }else{
                res += sum * exp(-dell);
            }
            sum = 0;
            dell = ipow(delta,pt->l);
            ldell = pt->l*dell;
        }
    }

    assert(!isinf(res));

    /* gaussian terms */
    n = data->ng;
    gt = &(data->gt[0]);
    for(i=0; i<n; ++i){
#ifdef RESID_DEBUG
        fprintf(stderr,"i = %d, GAUSSIAN, n = %e, t = %f, d = %f, alpha = %f, beta = %f, gamma = %f, epsilon = %f\n",i+1, gt->n, gt->t, gt->d, gt->alpha, gt->beta, gt->gamma, gt->epsilon);
#endif
        double d1 = delta - gt->epsilon;
        double t1 = tau - gt->gamma;
        double e1 = -gt->alpha*SQ(d1) - gt->beta*SQ(t1);

        double f1 = gt->t - 2*gt->beta*tau*(tau - gt->gamma);
        double g1 = gt->d - 2*gt->alpha*delta*(delta - gt->epsilon);

        sum = gt->n * f1 * pow(tau,gt->t-1) * g1 * pow(delta,gt->d-1) * exp(e1);

        //fprintf(stderr,"sum = %f\n",sum);
        res += sum;
        assert(!isinf(res));

        ++gt;
    }

#ifdef RESID_DEBUG
    fprintf(stderr,"phir = %f\n",res);
#endif

    assert(!isnan(res));
    assert(!isinf(res));
    return res;
}

/**
    Second derivative of helmholtz residual function with respect to
    delta (twice).

    FIXME this function is WRONG.
*/
double helm_resid_deldel(double tau,double delta,const HelmholtzData *data){
    double sum = 0, res = 0;
    double dell, ldell;
    unsigned n, i;
    const HelmholtzPowTerm *pt;
    const HelmholtzGausTerm *gt;


#ifdef RESID_DEBUG
        fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
#endif

    /* power terms */
    n = data->np;
    pt = &(data->pt[0]);
    dell = ipow(delta,pt->l);
    ldell = pt->l * dell;
    unsigned oldl;
    for(i=0; i<n; ++i){
        double lpart = pt->l ? SQ(ldell) + ldell*(1. - 2*pt->d - pt->l) : 0;
        sum += pt->a * pow(tau, pt->t) * ipow(delta, pt->d - 2) * (pt->d*(pt->d - 1) + lpart);
        oldl = pt->l;
        ++pt;
        if(i+1==n || oldl != pt->l){
            if(oldl == 0){
                res += sum;
            }else{
                res += sum * exp(-dell);
            }
            sum = 0;
            dell = ipow(delta,pt->l);
            ldell = pt->l*dell;
        }
    }

    /* gaussian terms */
    n = data->ng;
    //fprintf(stderr,"THERE ARE %d GAUSSIAN TERMS\n",n);
    gt = &(data->gt[0]);
    for(i=0; i<n; ++i){
        double s1 = SQ(delta - gt->epsilon);
        double f1 = 2*delta*gt->alpha *(2*gt->d*gt->epsilon
            - delta * (2*gt->d + 1 - 2 * gt->alpha * s1));
        res += - gt->n * pow(tau,gt->t) * pow(delta, -1. + gt->d)
            * f1
            * exp(-(gt->alpha * s1 + gt->beta*SQ(tau-gt->gamma)));
        ++gt;
    }

    return res;
}

1	jpye	1810	/* ASCEND modelling environment
2	jpye	1831	Copyright (C) 2008 Carnegie Mellon University
3	jpye	1810
4			This program is free software; you can redistribute it and/or modify
5			it under the terms of the GNU General Public License as published by
6			the Free Software Foundation; either version 2, or (at your option)
7			any later version.
8
9			This program is distributed in the hope that it will be useful,
10			but WITHOUT ANY WARRANTY; without even the implied warranty of
11			MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12			GNU General Public License for more details.
13
14			You should have received a copy of the GNU General Public License
15			along with this program; if not, write to the Free Software
16			Foundation, Inc., 59 Temple Place - Suite 330,
17			Boston, MA 02111-1307, USA.
18			//* @file
19			Implementation of the reduced molar Helmholtz free energy equation of state.
20
21	jpye	1822	For nomenclature see Tillner-Roth, Harms-Watzenberg and Baehr, Eine neue
22			Fundamentalgleichung f端r Ammoniak.
23
24	jpye	1810	John Pye, 29 Jul 2008.
25			*/
26
27	jpye	1822	#include <math.h>
28	jpye	1810
29	jpye	1822	#include "helmholtz.h"
30	jpye	1847	#include "ideal_impl.h"
31	jpye	1822
32	jpye	1826	#ifdef TEST
33			#include <assert.h>
34			#include <stdlib.h>
35			#include <stdio.h>
36			#endif
37
38	jpye	1905	#define SQ(X) ((X)*(X))
39
40	jpye	1822	/* forward decls */
41
42	jpye	1826	static double helm_resid(double tau, double delta, const HelmholtzData *data);
43			static double helm_resid_del(double tau, double delta, const HelmholtzData *data);
44			static double helm_resid_tau(double tau, double delta, const HelmholtzData *data);
45			static double helm_resid_deltau(double tau, double delta, const HelmholtzData *data);
46			static double helm_resid_deldel(double tau, double delta, const HelmholtzData *data);
47	jpye	1822
48	jpye	1810	/**
49			Function to calculate pressure from Helmholtz free energy EOS, given temperature
50	jpye	1824	and mass density.
51	jpye	1810
52			@param T temperature in K
53	jpye	1824	@param rho mass density in kg/m続
54			@return pressure in Pa???
55	jpye	1810	*/
56	jpye	1826	double helmholtz_p(double T, double rho, const HelmholtzData *data){
57	jpye	1822
58			double tau = data->T_star / T;
59			double delta = rho / data->rho_star;
60	jpye	1810
61	jpye	1838	#ifdef TEST
62			assert(data->rho_star!=0);
63			assert(T!=0);
64			assert(!isnan(tau));
65			assert(!isnan(delta));
66			assert(!isnan(data->R));
67	jpye	1839
68	jpye	1840	//fprintf(stderr,"p calc: T = %f\n",T);
69			//fprintf(stderr,"p calc: tau = %f\n",tau);
70			//fprintf(stderr,"p calc: rho = %f\n",rho);
71			//fprintf(stderr,"p calc: delta = %f\n",delta);
72			//fprintf(stderr,"p calc: RTrho = %f\n",data->R * T * rho);
73	jpye	1841
74			//fprintf(stderr,"T = %f\n", T);
75			//fprintf(stderr,"rhob = %f, rhob* = %f, delta = %f\n", rho/data->M, data->rho_star/data->M, delta);
76	jpye	1838	#endif
77	jpye	1839
78	jpye	1841	return data->R * T * rho * (1 + delta * helm_resid_del(tau,delta,data));
79	jpye	1822	}
80	jpye	1810
81	jpye	1824	/**
82			Function to calculate internal energy from Helmholtz free energy EOS, given
83			temperature and mass density.
84
85			@param T temperature in K
86			@param rho mass density in kg/m続
87			@return internal energy in ???
88			*/
89	jpye	1826	double helmholtz_u(double T, double rho, const HelmholtzData *data){
90	jpye	1824
91			double tau = data->T_star / T;
92			double delta = rho / data->rho_star;
93
94	jpye	1832	#ifdef TEST
95	jpye	1838	assert(data->rho_star!=0);
96			assert(T!=0);
97			assert(!isnan(tau));
98			assert(!isnan(delta));
99			assert(!isnan(data->R));
100			#endif
101
102	jpye	1857	#if 0
103	jpye	1835	fprintf(stderr,"ideal_tau = %f\n",helm_ideal_tau(tau,delta,data->ideal));
104	jpye	1827	fprintf(stderr,"resid_tau = %f\n",helm_resid_tau(tau,delta,data));
105			fprintf(stderr,"R T = %f\n",data->R * data->T_star);
106	jpye	1832	#endif
107	jpye	1827
108	jpye	1835	return data->R * data->T_star * (helm_ideal_tau(tau,delta,data->ideal) + helm_resid_tau(tau,delta,data));
109	jpye	1824	}
110
111			/**
112			Function to calculate enthalpy from Helmholtz free energy EOS, given
113			temperature and mass density.
114
115			@param T temperature in K
116			@param rho mass density in kg/m続
117	jpye	1829	@return enthalpy in J/kg
118	jpye	1824	*/
119	jpye	1826	double helmholtz_h(double T, double rho, const HelmholtzData *data){
120	jpye	1824
121			double tau = data->T_star / T;
122			double delta = rho / data->rho_star;
123
124	jpye	1838	#ifdef TEST
125			assert(data->rho_star!=0);
126			assert(T!=0);
127			assert(!isnan(tau));
128			assert(!isnan(delta));
129			assert(!isnan(data->R));
130			#endif
131
132	jpye	1835	return data->R * T * (1 + tau * (helm_ideal_tau(tau,delta,data->ideal) + helm_resid_tau(tau,delta,data)) + delta*helm_resid_del(tau,delta,data));
133	jpye	1824	}
134
135	jpye	1829	/**
136			Function to calculate entropy from Helmholtz free energy EOS, given
137			temperature and mass density.
138
139			@param T temperature in K
140			@param rho mass density in kg/m続
141			@return entropy in J/kgK
142			*/
143			double helmholtz_s(double T, double rho, const HelmholtzData *data){
144
145			double tau = data->T_star / T;
146			double delta = rho / data->rho_star;
147
148	jpye	1870	#ifdef ENTROPY_DEBUG
149	jpye	1838	assert(data->rho_star!=0);
150			assert(T!=0);
151			assert(!isnan(tau));
152			assert(!isnan(delta));
153			assert(!isnan(data->R));
154
155	jpye	1845	fprintf(stderr,"helm_ideal_tau = %f\n",helm_ideal_tau(tau,delta,data->ideal));
156			fprintf(stderr,"helm_resid_tau = %f\n",helm_resid_tau(tau,delta,data));
157			fprintf(stderr,"helm_ideal = %f\n",helm_ideal(tau,delta,data->ideal));
158			fprintf(stderr,"helm_resid = %f\n",helm_resid(tau,delta,data));
159	jpye	1847	#endif
160	jpye	1829	return data->R * (
161	jpye	1835	tau * (helm_ideal_tau(tau,delta,data->ideal) + helm_resid_tau(tau,delta,data))
162	jpye	1871	- (helm_ideal(tau,delta,data->ideal) + helm_resid(tau,delta,data))
163	jpye	1829	);
164			}
165
166	jpye	1850	/**
167	jpye	1863	Function to calculate Helmholtz energy from the Helmholtz free energy EOS,
168			given temperature and mass density.
169
170			@param T temperature in K
171			@param rho mass density in kg/m続
172			@return Helmholtz energy 'a', in J/kg
173			*/
174			double helmholtz_a(double T, double rho, const HelmholtzData *data){
175
176			double tau = data->T_star / T;
177			double delta = rho / data->rho_star;
178
179			#ifdef TEST
180			assert(data->rho_star!=0);
181			assert(T!=0);
182			assert(!isnan(tau));
183			assert(!isnan(delta));
184			assert(!isnan(data->R));
185			#endif
186
187	jpye	1868	#ifdef HELMHOLTZ_DEBUG
188	jpye	1863	fprintf(stderr,"helmholtz_a: T = %f, rho = %f\n",T,rho);
189	jpye	1865	fprintf(stderr,"multiplying by RT = %f\n",data->R*T);
190	jpye	1863	#endif
191
192	jpye	1870	return data->R * T * (helm_ideal(tau,delta,data->ideal) + helm_resid(tau,delta,data));
193	jpye	1863	}
194
195
196			/**
197	jpye	1850	Calculation zero-pressure specific heat capacity
198			*/
199			double helmholtz_cp0(double T, const HelmholtzData *data){
200			double val = helm_cp0(T,data->ideal);
201	jpye	1855	#if 0
202	jpye	1850	fprintf(stderr,"val = %f\n",val);
203	jpye	1851	#endif
204	jpye	1850	return val;
205			}
206
207	jpye	1905	/**
208			Calculate partial derivative of p with respect to T, with rho constant
209			*/
210			double helmholtz_dpdT_rho(double T, double rho, const HelmholtzData *data){
211			double tau = data->T_star / T;
212			double delta = rho / data->rho_star;
213	jpye	1906
214			double phir_del = helm_resid_del(tau,delta,data);
215			double phir_deltau = helm_resid_deltau(tau,delta,data);
216			assert(!isinf(phir_del));
217			assert(!isinf(phir_deltau));
218			assert(!isnan(phir_del));
219			assert(!isnan(phir_deltau));
220			assert(!isnan(data->R));
221			assert(!isnan(rho));
222			assert(!isnan(tau));
223
224			double res = data->R * rho * (1 + deltaphir_del - deltatau*phir_deltau);
225
226			assert(!isnan(res));
227			assert(!isinf(res));
228			return res;
229	jpye	1905	}
230
231			/**
232			Calculate partial derivative of p with respect to rho, with T constant
233			*/
234			double helmholtz_dpdrho_T(double T, double rho, const HelmholtzData *data){
235			double tau = data->T_star / T;
236			double delta = rho / data->rho_star;
237	jpye	1906
238			double phir_del = helm_resid_del(tau,delta,data);
239			double phir_deldel = helm_resid_deldel(tau,delta,data);
240			assert(!isinf(phir_del));
241			assert(!isinf(phir_deldel));
242	jpye	1905
243	jpye	1906	return data->R * T * (1 + 2deltaphir_del + deltadelta phir_deldel);
244	jpye	1905	}
245
246	jpye	1824	/*---------------------------------------------
247	jpye	1826	UTILITY FUNCTION(S)
248			*/
249
250			/* ipow: public domain by Mark Stephen with suggestions by Keiichi Nakasato */
251			static double ipow(double x, int n){
252			double t = 1.0;
253
254	jpye	1829	if(!n)return 1.0; /* At the top. x^0 = 1 */
255	jpye	1826
256			if (n < 0){
257			n = -n;
258			x = 1.0/x; /* error if x == 0. Good */
259			} /* ZTC/SC returns inf, which is even better */
260
261			if (x == 0.0)return 0.0;
262
263			do{
264			if(n & 1)t *= x;
265			n /= 2; /* KN prefers if (n/=2) x=x; This avoids an /
266			x = x; / unnecessary but benign multiplication on */
267			}while(n); /* the last pass, but the comparison is always
268			true _except_ on the last pass. */
269
270			return t;
271			}
272
273	jpye	1889	//#define RESID_DEBUG
274	jpye	1872
275	jpye	1824	/**
276	jpye	1845	Residual part of helmholtz function.
277	jpye	1822	*/
278	jpye	1826	double helm_resid(double tau, double delta, const HelmholtzData *data){
279	jpye	1867	double dell,ldell, term, sum, res = 0;
280	jpye	1845	unsigned n, i;
281			const HelmholtzPowTerm *pt;
282	jpye	1887	const HelmholtzGausTerm *gt;
283	jpye	1845
284			n = data->np;
285			pt = &(data->pt[0]);
286
287	jpye	1868	#ifdef RESID_DEBUG
288	jpye	1867	fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
289			#endif
290
291	jpye	1845	/* power terms */
292			sum = 0;
293	jpye	1865	dell = ipow(delta,pt->l);
294			ldell = pt->l * dell;
295	jpye	1845	unsigned oldl;
296			for(i=0; i<n; ++i){
297	jpye	1867	term = pt->a * pow(tau, pt->t) * ipow(delta, pt->d);
298			sum += term;
299	jpye	1868	#ifdef RESID_DEBUG
300	jpye	1867	fprintf(stderr,"i = %d, a=%e, t=%f, d=%d, term = %f, sum = %f",i,pt->a,pt->t,pt->d,term,sum);
301			if(pt->l==0){
302			fprintf(stderr,",row=%e\n",term);
303			}else{
304			fprintf(stderr,",row=%e\n,",term*exp(-dell));
305			}
306			#endif
307	jpye	1845	oldl = pt->l;
308			++pt;
309			if(i+1==n \|\| oldl != pt->l){
310			if(oldl == 0){
311	jpye	1868	#ifdef RESID_DEBUG
312	jpye	1865	fprintf(stderr,"linear ");
313	jpye	1867	#endif
314	jpye	1845	res += sum;
315			}else{
316	jpye	1868	#ifdef RESID_DEBUG
317	jpye	1865	fprintf(stderr,"exp dell=%f, exp(-dell)=%f sum=%f: ",dell,exp(-dell),sum);
318	jpye	1867	#endif
319	jpye	1865	res += sum * exp(-dell);
320	jpye	1845	}
321	jpye	1868	#ifdef RESID_DEBUG
322	jpye	1865	fprintf(stderr,"i = %d, res = %f\n",i,res);
323	jpye	1867	#endif
324	jpye	1845	sum = 0;
325	jpye	1865	dell = ipow(delta,pt->l);
326			ldell = pt->l*dell;
327	jpye	1845	}
328			}
329
330	jpye	1887	/* gaussian terms */
331			n = data->ng;
332	jpye	1889	//fprintf(stderr,"THERE ARE %d GAUSSIAN TERMS\n",n);
333	jpye	1887	gt = &(data->gt[0]);
334			for(i=0; i<n; ++i){
335	jpye	1868	#ifdef RESID_DEBUG
336	jpye	1887	fprintf(stderr,"i = %d, GAUSSIAN, n = %e, t = %f, d = %f, alpha = %f, beta = %f, gamma = %f, epsilon = %f\n",i+1, gt->n, gt->t, gt->d, gt->alpha, gt->beta, gt->gamma, gt->epsilon);
337			#endif
338			double d1 = delta - gt->epsilon;
339			double t1 = tau - gt->gamma;
340			double e1 = -gt->alphad1d1 - gt->betat1t1;
341			sum = gt->n * pow(tau,gt->t) * pow(delta,gt->d) * exp(e1);
342			//fprintf(stderr,"sum = %f\n",sum);
343			res += sum;
344			++gt;
345			}
346
347			#ifdef RESID_DEBUG
348	jpye	1864	fprintf(stderr,"phir = %f\n",res);
349			#endif
350	jpye	1845	return res;
351			}
352
353	jpye	1822	/**
354			Derivative of the helmholtz residual function with respect to
355			delta.
356			*/
357	jpye	1826	double helm_resid_del(double tau,double delta, const HelmholtzData *data){
358	jpye	1905	double sum = 0, res = 0;
359	jpye	1844	double dell, ldell;
360	jpye	1838	unsigned n, i;
361			const HelmholtzPowTerm *pt;
362	jpye	1887	const HelmholtzGausTerm *gt;
363	jpye	1810
364	jpye	1822
365	jpye	1898	#ifdef RESID_DEBUG
366			fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
367			#endif
368
369	jpye	1905	/* power terms */
370			n = data->np;
371			pt = &(data->pt[0]);
372	jpye	1844	dell = ipow(delta,pt->l);
373			ldell = pt->l * dell;
374			unsigned oldl;
375	jpye	1841	for(i=0; i<n; ++i){
376	jpye	1844	sum += pt->a * pow(tau, pt->t) * ipow(delta, pt->d - 1) * (pt->d - ldell);
377			oldl = pt->l;
378	jpye	1841	++pt;
379	jpye	1844	if(i+1==n \|\| oldl != pt->l){
380			if(oldl == 0){
381	jpye	1836	res += sum;
382			}else{
383	jpye	1844	res += sum * exp(-dell);
384	jpye	1836	}
385	jpye	1844	sum = 0;
386			dell = ipow(delta,pt->l);
387			ldell = pt->l*dell;
388	jpye	1836	}
389	jpye	1822	}
390
391	jpye	1887	/* gaussian terms */
392			n = data->ng;
393	jpye	1889	//fprintf(stderr,"THERE ARE %d GAUSSIAN TERMS\n",n);
394	jpye	1887	gt = &(data->gt[0]);
395			for(i=0; i<n; ++i){
396			#ifdef RESID_DEBUG
397	jpye	1891	fprintf(stderr,"i = %d, GAUSSIAN, n = %e, t = %f, d = %f, alpha = %f, beta = %f, gamma = %f, epsilon = %f\n",i+1, gt->n, gt->t, gt->d, gt->alpha, gt->beta, gt->gamma, gt->epsilon);
398	jpye	1840	#endif
399	jpye	1898	double val2;
400			val2 = - gt->n * pow(tau,gt->t) * pow(delta, -1. + gt->d)
401			* (2. * gt->alpha * delta * (delta - gt->epsilon) - gt->d)
402			* exp(-(gt->alpha * SQ(delta-gt->epsilon) + gt->beta*SQ(tau-gt->gamma)));
403			res += val2;
404	jpye	1891	#ifdef RESID_DEBUG
405	jpye	1898	fprintf(stderr,"val2 = %f --> res = %f\n",val2,res);
406	jpye	1891	#endif
407	jpye	1887	++gt;
408			}
409	jpye	1838
410	jpye	1836	return res;
411	jpye	1822	}
412
413	jpye	1824	/**
414			Derivative of the helmholtz residual function with respect to
415			tau.
416			*/
417	jpye	1826	double helm_resid_tau(double tau,double delta,const HelmholtzData *data){
418	jpye	1824
419			double sum;
420	jpye	1827	double res = 0;
421	jpye	1832	double delX;
422			unsigned l;
423	jpye	1845	unsigned n, i;
424	jpye	1838	const HelmholtzPowTerm *pt;
425	jpye	1891	const HelmholtzGausTerm *gt;
426	jpye	1822
427	jpye	1845	n = data->np;
428	jpye	1838	pt = &(data->pt[0]);
429	jpye	1832
430			delX = 1;
431
432			l = 0;
433			sum = 0;
434	jpye	1845	for(i=0; i<n; ++i){
435	jpye	1838	if(pt->t){
436			//fprintf(stderr,"i = %d, a = %e, t = %f, d = %d, l = %d\n",i+1, pt->a, pt->t, pt->d, pt->l);
437			sum += pt->a * pow(tau, pt->t - 1) * ipow(delta, pt->d) * pt->t;
438	jpye	1832	}
439	jpye	1838	++pt;
440	jpye	1832	//fprintf(stderr,"l = %d\n",l);
441	jpye	1845	if(i+1==n \|\| l != pt->l){
442	jpye	1832	if(l==0){
443			//fprintf(stderr,"Adding non-exp term\n");
444			res += sum;
445			}else{
446			//fprintf(stderr,"Adding exp term with l = %d, delX = %e\n",l,delX);
447			res += sum * exp(-delX);
448			}
449			/* set l to new value */
450	jpye	1845	if(i+1!=n){
451	jpye	1838	l = pt->l;
452	jpye	1832	//fprintf(stderr,"New l = %d\n",l);
453			delX = ipow(delta,l);
454			sum = 0;
455			}
456			}
457			}
458
459	jpye	1897	//#define RESID_DEBUG
460	jpye	1891	/* gaussian terms */
461			n = data->ng;
462			gt = &(data->gt[0]);
463			for(i=0; i<n; ++i){
464			#ifdef RESID_DEBUG
465			fprintf(stderr,"i = %d, GAUSSIAN, n = %e, t = %f, d = %f, alpha = %f, beta = %f, gamma = %f, epsilon = %f\n",i+1, gt->n, gt->t, gt->d, gt->alpha, gt->beta, gt->gamma, gt->epsilon);
466			#endif
467	jpye	1897
468			double val2;
469	jpye	1900	val2 = -gt->n * pow(tau,gt->t - 1.) * pow(delta, gt->d)
470			* (2. * gt->beta * tau * (tau - gt->gamma) - gt->t)
471	jpye	1898	* exp(-(gt->alpha * SQ(delta-gt->epsilon) + gt->beta*SQ(tau-gt->gamma)));
472			res += val2;
473	jpye	1891	#ifdef RESID_DEBUG
474	jpye	1897	fprintf(stderr,"res = %f\n",res);
475	jpye	1891	#endif
476	jpye	1897
477	jpye	1891	++gt;
478			}
479
480	jpye	1827	return res;
481	jpye	1824	}
482
483
484
485	jpye	1826	/**
486			Mixed derivative of the helmholtz residual function with respect to
487	jpye	1907	delta and tau.
488	jpye	1868	*/
489	jpye	1826	double helm_resid_deltau(double tau,double delta,const HelmholtzData *data){
490	jpye	1906	double dell,ldell, term, sum = 0, res = 0;
491	jpye	1905	unsigned n, i;
492			const HelmholtzPowTerm *pt;
493			const HelmholtzGausTerm *gt;
494	jpye	1824
495	jpye	1905	/* power terms */
496			n = data->np;
497			pt = &(data->pt[0]);
498			dell = ipow(delta,pt->l);
499			ldell = pt->l * dell;
500			unsigned oldl;
501			for(i=0; i<n; ++i){
502			sum += pt->a * pt->t * pow(tau, pt->t - 1) * ipow(delta, pt->d - 1) * (pt->d - ldell);
503			oldl = pt->l;
504	jpye	1838	++pt;
505	jpye	1905	if(i+1==n \|\| oldl != pt->l){
506			if(oldl == 0){
507			res += sum;
508			}else{
509			res += sum * exp(-dell);
510			}
511			sum = 0;
512			dell = ipow(delta,pt->l);
513			ldell = pt->l*dell;
514			}
515	jpye	1826	}
516
517	jpye	1906	assert(!isinf(res));
518
519	jpye	1905	/* gaussian terms */
520			n = data->ng;
521			gt = &(data->gt[0]);
522			for(i=0; i<n; ++i){
523			#ifdef RESID_DEBUG
524			fprintf(stderr,"i = %d, GAUSSIAN, n = %e, t = %f, d = %f, alpha = %f, beta = %f, gamma = %f, epsilon = %f\n",i+1, gt->n, gt->t, gt->d, gt->alpha, gt->beta, gt->gamma, gt->epsilon);
525			#endif
526			double d1 = delta - gt->epsilon;
527			double t1 = tau - gt->gamma;
528			double e1 = -gt->alphaSQ(d1) - gt->betaSQ(t1);
529	jpye	1826
530	jpye	1905	double f1 = gt->t - 2gt->betatau*(tau - gt->gamma);
531			double g1 = gt->d - 2gt->alphadelta*(delta - gt->epsilon);
532	jpye	1826
533	jpye	1905	sum = gt->n * f1 * pow(tau,gt->t-1) * g1 * pow(delta,gt->d-1) * exp(e1);
534	jpye	1906
535	jpye	1905	//fprintf(stderr,"sum = %f\n",sum);
536			res += sum;
537	jpye	1906	assert(!isinf(res));
538
539	jpye	1905	++gt;
540	jpye	1826	}
541
542	jpye	1905	#ifdef RESID_DEBUG
543			fprintf(stderr,"phir = %f\n",res);
544			#endif
545	jpye	1906
546			assert(!isnan(res));
547			assert(!isinf(res));
548	jpye	1905	return res;
549	jpye	1826	}
550
551			/**
552			Second derivative of helmholtz residual function with respect to
553			delta (twice).
554	jpye	1890
555			FIXME this function is WRONG.
556	jpye	1826	*/
557			double helm_resid_deldel(double tau,double delta,const HelmholtzData *data){
558	jpye	1907	double sum = 0, res = 0;
559			double dell, ldell;
560			unsigned n, i;
561			const HelmholtzPowTerm *pt;
562			const HelmholtzGausTerm *gt;
563	jpye	1826
564
565	jpye	1907	#ifdef RESID_DEBUG
566			fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
567			#endif
568	jpye	1826
569	jpye	1907	/* power terms */
570			n = data->np;
571			pt = &(data->pt[0]);
572			dell = ipow(delta,pt->l);
573			ldell = pt->l * dell;
574			unsigned oldl;
575			for(i=0; i<n; ++i){
576			double lpart = pt->l ? SQ(ldell) + ldell(1. - 2pt->d - pt->l) : 0;
577			sum += pt->a * pow(tau, pt->t) * ipow(delta, pt->d - 2) * (pt->d*(pt->d - 1) + lpart);
578			oldl = pt->l;
579	jpye	1838	++pt;
580	jpye	1907	if(i+1==n \|\| oldl != pt->l){
581			if(oldl == 0){
582			res += sum;
583			}else{
584			res += sum * exp(-dell);
585			}
586			sum = 0;
587			dell = ipow(delta,pt->l);
588			ldell = pt->l*dell;
589			}
590	jpye	1826	}
591
592	jpye	1907	/* gaussian terms */
593			n = data->ng;
594			//fprintf(stderr,"THERE ARE %d GAUSSIAN TERMS\n",n);
595			gt = &(data->gt[0]);
596			for(i=0; i<n; ++i){
597			double s1 = SQ(delta - gt->epsilon);
598			double f1 = 2deltagt->alpha (2gt->d*gt->epsilon
599			- delta * (2gt->d + 1 - 2 gt->alpha * s1));
600			res += - gt->n * pow(tau,gt->t) * pow(delta, -1. + gt->d)
601			* f1
602			* exp(-(gt->alpha * s1 + gt->beta*SQ(tau-gt->gamma)));
603			++gt;
604	jpye	1826	}
605
606	jpye	1907	return res;
607	jpye	1826	}
608