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);
#ifdef TEST
    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));
#endif

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

#ifdef TEST
    assert(!isnan(res));
    assert(!isinf(res));
#endif
    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);
#ifdef TEST
    assert(!isinf(phir_del));
    assert(!isinf(phir_deldel));
#endif  
    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;
        }
    }

#ifdef TEST
    assert(!isinf(res));
#endif

    /* 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;
#ifdef TEST
        assert(!isinf(res));
#endif
        ++gt;
    }

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

#ifdef TEST
    assert(!isnan(res));
    assert(!isinf(res));
#endif
    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 = gt->d*(gt->d - 1) 
            + 2.*gt->alpha*delta * (
                delta * (2. * gt->alpha * s1 - 1) 
                - 2. * gt->d * (delta - gt->epsilon)
            );
        res += gt->n * pow(tau,gt->t) * pow(delta, gt->d - 2.)
            * f1
            * exp(-(gt->alpha * s1 + gt->beta*SQ(tau-gt->gamma)));
        ++gt;
    }

    return res;
}

1	/* ASCEND modelling environment
2	Copyright (C) 2008 Carnegie Mellon University
3
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	For nomenclature see Tillner-Roth, Harms-Watzenberg and Baehr, Eine neue
22	Fundamentalgleichung f端r Ammoniak.
23
24	John Pye, 29 Jul 2008.
25	*/
26
27	#include <math.h>
28
29	#include "helmholtz.h"
30	#include "ideal_impl.h"
31
32	#ifdef TEST
33	#include <assert.h>
34	#include <stdlib.h>
35	#include <stdio.h>
36	#endif
37
38	#define SQ(X) ((X)*(X))
39
40	/* forward decls */
41
42	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
48	/**
49	Function to calculate pressure from Helmholtz free energy EOS, given temperature
50	and mass density.
51
52	@param T temperature in K
53	@param rho mass density in kg/m続
54	@return pressure in Pa???
55	*/
56	double helmholtz_p(double T, double rho, const HelmholtzData *data){
57
58	double tau = data->T_star / T;
59	double delta = rho / data->rho_star;
60
61	#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
68	//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
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	#endif
77
78	return data->R * T * rho * (1 + delta * helm_resid_del(tau,delta,data));
79	}
80
81	/**
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	double helmholtz_u(double T, double rho, const HelmholtzData *data){
90
91	double tau = data->T_star / T;
92	double delta = rho / data->rho_star;
93
94	#ifdef TEST
95	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	#if 0
103	fprintf(stderr,"ideal_tau = %f\n",helm_ideal_tau(tau,delta,data->ideal));
104	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	#endif
107
108	return data->R * data->T_star * (helm_ideal_tau(tau,delta,data->ideal) + helm_resid_tau(tau,delta,data));
109	}
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	@return enthalpy in J/kg
118	*/
119	double helmholtz_h(double T, double rho, const HelmholtzData *data){
120
121	double tau = data->T_star / T;
122	double delta = rho / data->rho_star;
123
124	#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	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	}
134
135	/**
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	#ifdef ENTROPY_DEBUG
149	assert(data->rho_star!=0);
150	assert(T!=0);
151	assert(!isnan(tau));
152	assert(!isnan(delta));
153	assert(!isnan(data->R));
154
155	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	#endif
160	return data->R * (
161	tau * (helm_ideal_tau(tau,delta,data->ideal) + helm_resid_tau(tau,delta,data))
162	- (helm_ideal(tau,delta,data->ideal) + helm_resid(tau,delta,data))
163	);
164	}
165
166	/**
167	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	#ifdef HELMHOLTZ_DEBUG
188	fprintf(stderr,"helmholtz_a: T = %f, rho = %f\n",T,rho);
189	fprintf(stderr,"multiplying by RT = %f\n",data->R*T);
190	#endif
191
192	return data->R * T * (helm_ideal(tau,delta,data->ideal) + helm_resid(tau,delta,data));
193	}
194
195
196	/**
197	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	#if 0
202	fprintf(stderr,"val = %f\n",val);
203	#endif
204	return val;
205	}
206
207	/**
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
214	double phir_del = helm_resid_del(tau,delta,data);
215	double phir_deltau = helm_resid_deltau(tau,delta,data);
216	#ifdef TEST
217	assert(!isinf(phir_del));
218	assert(!isinf(phir_deltau));
219	assert(!isnan(phir_del));
220	assert(!isnan(phir_deltau));
221	assert(!isnan(data->R));
222	assert(!isnan(rho));
223	assert(!isnan(tau));
224	#endif
225
226	double res = data->R * rho * (1 + deltaphir_del - deltatau*phir_deltau);
227
228	#ifdef TEST
229	assert(!isnan(res));
230	assert(!isinf(res));
231	#endif
232	return res;
233	}
234
235	/**
236	Calculate partial derivative of p with respect to rho, with T constant
237	*/
238	double helmholtz_dpdrho_T(double T, double rho, const HelmholtzData *data){
239	double tau = data->T_star / T;
240	double delta = rho / data->rho_star;
241
242	double phir_del = helm_resid_del(tau,delta,data);
243	double phir_deldel = helm_resid_deldel(tau,delta,data);
244	#ifdef TEST
245	assert(!isinf(phir_del));
246	assert(!isinf(phir_deldel));
247	#endif
248	return data->R * T * (1 + 2deltaphir_del + deltadelta phir_deldel);
249	}
250
251	/*---------------------------------------------
252	UTILITY FUNCTION(S)
253	*/
254
255	/* ipow: public domain by Mark Stephen with suggestions by Keiichi Nakasato */
256	static double ipow(double x, int n){
257	double t = 1.0;
258
259	if(!n)return 1.0; /* At the top. x^0 = 1 */
260
261	if (n < 0){
262	n = -n;
263	x = 1.0/x; /* error if x == 0. Good */
264	} /* ZTC/SC returns inf, which is even better */
265
266	if (x == 0.0)return 0.0;
267
268	do{
269	if(n & 1)t *= x;
270	n /= 2; /* KN prefers if (n/=2) x=x; This avoids an /
271	x = x; / unnecessary but benign multiplication on */
272	}while(n); /* the last pass, but the comparison is always
273	true _except_ on the last pass. */
274
275	return t;
276	}
277
278	//#define RESID_DEBUG
279
280	/**
281	Residual part of helmholtz function.
282	*/
283	double helm_resid(double tau, double delta, const HelmholtzData *data){
284	double dell,ldell, term, sum, res = 0;
285	unsigned n, i;
286	const HelmholtzPowTerm *pt;
287	const HelmholtzGausTerm *gt;
288
289	n = data->np;
290	pt = &(data->pt[0]);
291
292	#ifdef RESID_DEBUG
293	fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
294	#endif
295
296	/* power terms */
297	sum = 0;
298	dell = ipow(delta,pt->l);
299	ldell = pt->l * dell;
300	unsigned oldl;
301	for(i=0; i<n; ++i){
302	term = pt->a * pow(tau, pt->t) * ipow(delta, pt->d);
303	sum += term;
304	#ifdef RESID_DEBUG
305	fprintf(stderr,"i = %d, a=%e, t=%f, d=%d, term = %f, sum = %f",i,pt->a,pt->t,pt->d,term,sum);
306	if(pt->l==0){
307	fprintf(stderr,",row=%e\n",term);
308	}else{
309	fprintf(stderr,",row=%e\n,",term*exp(-dell));
310	}
311	#endif
312	oldl = pt->l;
313	++pt;
314	if(i+1==n \|\| oldl != pt->l){
315	if(oldl == 0){
316	#ifdef RESID_DEBUG
317	fprintf(stderr,"linear ");
318	#endif
319	res += sum;
320	}else{
321	#ifdef RESID_DEBUG
322	fprintf(stderr,"exp dell=%f, exp(-dell)=%f sum=%f: ",dell,exp(-dell),sum);
323	#endif
324	res += sum * exp(-dell);
325	}
326	#ifdef RESID_DEBUG
327	fprintf(stderr,"i = %d, res = %f\n",i,res);
328	#endif
329	sum = 0;
330	dell = ipow(delta,pt->l);
331	ldell = pt->l*dell;
332	}
333	}
334
335	/* gaussian terms */
336	n = data->ng;
337	//fprintf(stderr,"THERE ARE %d GAUSSIAN TERMS\n",n);
338	gt = &(data->gt[0]);
339	for(i=0; i<n; ++i){
340	#ifdef RESID_DEBUG
341	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);
342	#endif
343	double d1 = delta - gt->epsilon;
344	double t1 = tau - gt->gamma;
345	double e1 = -gt->alphad1d1 - gt->betat1t1;
346	sum = gt->n * pow(tau,gt->t) * pow(delta,gt->d) * exp(e1);
347	//fprintf(stderr,"sum = %f\n",sum);
348	res += sum;
349	++gt;
350	}
351
352	#ifdef RESID_DEBUG
353	fprintf(stderr,"phir = %f\n",res);
354	#endif
355	return res;
356	}
357
358	/**
359	Derivative of the helmholtz residual function with respect to
360	delta.
361	*/
362	double helm_resid_del(double tau,double delta, const HelmholtzData *data){
363	double sum = 0, res = 0;
364	double dell, ldell;
365	unsigned n, i;
366	const HelmholtzPowTerm *pt;
367	const HelmholtzGausTerm *gt;
368
369
370	#ifdef RESID_DEBUG
371	fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
372	#endif
373
374	/* power terms */
375	n = data->np;
376	pt = &(data->pt[0]);
377	dell = ipow(delta,pt->l);
378	ldell = pt->l * dell;
379	unsigned oldl;
380	for(i=0; i<n; ++i){
381	sum += pt->a * pow(tau, pt->t) * ipow(delta, pt->d - 1) * (pt->d - ldell);
382	oldl = pt->l;
383	++pt;
384	if(i+1==n \|\| oldl != pt->l){
385	if(oldl == 0){
386	res += sum;
387	}else{
388	res += sum * exp(-dell);
389	}
390	sum = 0;
391	dell = ipow(delta,pt->l);
392	ldell = pt->l*dell;
393	}
394	}
395
396	/* gaussian terms */
397	n = data->ng;
398	//fprintf(stderr,"THERE ARE %d GAUSSIAN TERMS\n",n);
399	gt = &(data->gt[0]);
400	for(i=0; i<n; ++i){
401	#ifdef RESID_DEBUG
402	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);
403	#endif
404	double val2;
405	val2 = - gt->n * pow(tau,gt->t) * pow(delta, -1. + gt->d)
406	* (2. * gt->alpha * delta * (delta - gt->epsilon) - gt->d)
407	* exp(-(gt->alpha * SQ(delta-gt->epsilon) + gt->beta*SQ(tau-gt->gamma)));
408	res += val2;
409	#ifdef RESID_DEBUG
410	fprintf(stderr,"val2 = %f --> res = %f\n",val2,res);
411	#endif
412	++gt;
413	}
414
415	return res;
416	}
417
418	/**
419	Derivative of the helmholtz residual function with respect to
420	tau.
421	*/
422	double helm_resid_tau(double tau,double delta,const HelmholtzData *data){
423
424	double sum;
425	double res = 0;
426	double delX;
427	unsigned l;
428	unsigned n, i;
429	const HelmholtzPowTerm *pt;
430	const HelmholtzGausTerm *gt;
431
432	n = data->np;
433	pt = &(data->pt[0]);
434
435	delX = 1;
436
437	l = 0;
438	sum = 0;
439	for(i=0; i<n; ++i){
440	if(pt->t){
441	//fprintf(stderr,"i = %d, a = %e, t = %f, d = %d, l = %d\n",i+1, pt->a, pt->t, pt->d, pt->l);
442	sum += pt->a * pow(tau, pt->t - 1) * ipow(delta, pt->d) * pt->t;
443	}
444	++pt;
445	//fprintf(stderr,"l = %d\n",l);
446	if(i+1==n \|\| l != pt->l){
447	if(l==0){
448	//fprintf(stderr,"Adding non-exp term\n");
449	res += sum;
450	}else{
451	//fprintf(stderr,"Adding exp term with l = %d, delX = %e\n",l,delX);
452	res += sum * exp(-delX);
453	}
454	/* set l to new value */
455	if(i+1!=n){
456	l = pt->l;
457	//fprintf(stderr,"New l = %d\n",l);
458	delX = ipow(delta,l);
459	sum = 0;
460	}
461	}
462	}
463
464	//#define RESID_DEBUG
465	/* gaussian terms */
466	n = data->ng;
467	gt = &(data->gt[0]);
468	for(i=0; i<n; ++i){
469	#ifdef RESID_DEBUG
470	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);
471	#endif
472
473	double val2;
474	val2 = -gt->n * pow(tau,gt->t - 1.) * pow(delta, gt->d)
475	* (2. * gt->beta * tau * (tau - gt->gamma) - gt->t)
476	* exp(-(gt->alpha * SQ(delta-gt->epsilon) + gt->beta*SQ(tau-gt->gamma)));
477	res += val2;
478	#ifdef RESID_DEBUG
479	fprintf(stderr,"res = %f\n",res);
480	#endif
481
482	++gt;
483	}
484
485	return res;
486	}
487
488
489
490	/**
491	Mixed derivative of the helmholtz residual function with respect to
492	delta and tau.
493	*/
494	double helm_resid_deltau(double tau,double delta,const HelmholtzData *data){
495	double dell,ldell, term, sum = 0, res = 0;
496	unsigned n, i;
497	const HelmholtzPowTerm *pt;
498	const HelmholtzGausTerm *gt;
499
500	/* power terms */
501	n = data->np;
502	pt = &(data->pt[0]);
503	dell = ipow(delta,pt->l);
504	ldell = pt->l * dell;
505	unsigned oldl;
506	for(i=0; i<n; ++i){
507	sum += pt->a * pt->t * pow(tau, pt->t - 1) * ipow(delta, pt->d - 1) * (pt->d - ldell);
508	oldl = pt->l;
509	++pt;
510	if(i+1==n \|\| oldl != pt->l){
511	if(oldl == 0){
512	res += sum;
513	}else{
514	res += sum * exp(-dell);
515	}
516	sum = 0;
517	dell = ipow(delta,pt->l);
518	ldell = pt->l*dell;
519	}
520	}
521
522	#ifdef TEST
523	assert(!isinf(res));
524	#endif
525
526	/* gaussian terms */
527	n = data->ng;
528	gt = &(data->gt[0]);
529	for(i=0; i<n; ++i){
530	#ifdef RESID_DEBUG
531	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);
532	#endif
533	double d1 = delta - gt->epsilon;
534	double t1 = tau - gt->gamma;
535	double e1 = -gt->alphaSQ(d1) - gt->betaSQ(t1);
536
537	double f1 = gt->t - 2gt->betatau*(tau - gt->gamma);
538	double g1 = gt->d - 2gt->alphadelta*(delta - gt->epsilon);
539
540	sum = gt->n * f1 * pow(tau,gt->t-1) * g1 * pow(delta,gt->d-1) * exp(e1);
541
542	//fprintf(stderr,"sum = %f\n",sum);
543	res += sum;
544	#ifdef TEST
545	assert(!isinf(res));
546	#endif
547	++gt;
548	}
549
550	#ifdef RESID_DEBUG
551	fprintf(stderr,"phir = %f\n",res);
552	#endif
553
554	#ifdef TEST
555	assert(!isnan(res));
556	assert(!isinf(res));
557	#endif
558	return res;
559	}
560
561	/**
562	Second derivative of helmholtz residual function with respect to
563	delta (twice).
564
565	FIXME this function is WRONG.
566	*/
567	double helm_resid_deldel(double tau,double delta,const HelmholtzData *data){
568	double sum = 0, res = 0;
569	double dell, ldell;
570	unsigned n, i;
571	const HelmholtzPowTerm *pt;
572	const HelmholtzGausTerm *gt;
573
574
575	#ifdef RESID_DEBUG
576	fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
577	#endif
578
579	/* power terms */
580	n = data->np;
581	pt = &(data->pt[0]);
582	dell = ipow(delta,pt->l);
583	ldell = pt->l * dell;
584	unsigned oldl;
585	for(i=0; i<n; ++i){
586	double lpart = pt->l ? SQ(ldell) + ldell(1. - 2pt->d - pt->l) : 0;
587	sum += pt->a * pow(tau, pt->t) * ipow(delta, pt->d - 2) * (pt->d*(pt->d - 1) + lpart);
588	oldl = pt->l;
589	++pt;
590	if(i+1==n \|\| oldl != pt->l){
591	if(oldl == 0){
592	res += sum;
593	}else{
594	res += sum * exp(-dell);
595	}
596	sum = 0;
597	dell = ipow(delta,pt->l);
598	ldell = pt->l*dell;
599	}
600	}
601
602	/* gaussian terms */
603	n = data->ng;
604	//fprintf(stderr,"THERE ARE %d GAUSSIAN TERMS\n",n);
605	gt = &(data->gt[0]);
606	for(i=0; i<n; ++i){
607	double s1 = SQ(delta - gt->epsilon);
608	double f1 = gt->d*(gt->d - 1)
609	+ 2.gt->alphadelta * (
610	delta * (2. * gt->alpha * s1 - 1)
611	- 2. * gt->d * (delta - gt->epsilon)
612	);
613	res += gt->n * pow(tau,gt->t) * pow(delta, gt->d - 2.)
614	* f1
615	* exp(-(gt->alpha * s1 + gt->beta*SQ(tau-gt->gamma)));
616	++gt;
617	}
618
619	return res;
620	}
621