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

/* 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;
}

/*---------------------------------------------
  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;
        }
    }

#if 1
    /* 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;
    }
#endif

#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, res = 0;
    double dell, ldell;
    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

    sum = 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;
        }
    }

#if 1
    /* 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
#define SQ(X) ((X)*(X))
        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
#undef SQ
        ++gt;
    }
#endif

    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

#define SQ(X) ((X)*(X))
        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
#undef SQ
            
        ++gt;
    }

    return res;
}   


/**
    Mixed derivative of the helmholtz residual function with respect to
    delta and tau

    FIXME this function is WRONG.
*/
double helm_resid_deltau(double tau,double delta,const HelmholtzData *data){
    
    double sum;
    double phir = 0;
    unsigned i;
    double XdelX;

    const HelmholtzPowTerm *pt = &(data->pt[0]);
    
    for(i=0; i<5; ++i){
        phir += pt->a * pow(tau, pt->t - 1) * ipow(delta, pt->d - 1) * pt->d * pt->t;
        ++pt;
    }

    sum = 0;
    XdelX = delta;
    for(i=5; i<10; ++i){
        sum += pt->a * pow(tau, pt->t - 1) * ipow(delta, pt->d - 1) * pt->t *(pt->d - XdelX);
        ++pt;
    }
    phir += exp(-delta) * sum;

    sum = 0; 
    XdelX = 2*delta*delta;
    for(i=10; i<17; ++i){
        sum += pt->a * pow(tau, pt->t - 1) * ipow(delta, pt->d - 1) * pt->t *(pt->d - XdelX);
        ++pt;
    }
    phir += exp(-delta*delta) * sum;

    sum = 0;
    XdelX = 3*delta*delta*delta;
    for(i=17; i<21; ++i){
        sum += pt->a * pow(tau, pt->t - 1) * ipow(delta, pt->d - 1) * pt->t *(pt->d - XdelX);
        ++pt;
    }
    phir += exp(-delta*delta*delta) * sum;

    return phir;
}

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

/**
    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;
    double phir = 0;
    unsigned i;
    unsigned X;
    double XdelX;

    const HelmholtzPowTerm *pt = &(data->pt[0]);
    
    for(i=0; i<5; ++i){
        phir += pt->a * pow(tau, pt->t) * ipow(delta, pt->d - 2) * (SQ(pt->d) - X);
        ++pt;
    }

    sum = 0;
    X = 1;
    XdelX = delta;
    for(i=5; i<10; ++i){
        sum += pt->a * pow(tau, pt->t) * ipow(delta, pt->d - 2) * (SQ(XdelX) - X*XdelX - 2*pt->d*XdelX + XdelX + SQ(pt->d) - pt->d);
        ++pt;
    }
    phir += exp(-delta) * sum;

    sum = 0; 
    X = 2;
    XdelX = 2*delta*delta;
    for(i=10; i<17; ++i){
        sum += pt->a * pow(tau, pt->t) * ipow(delta, pt->d - 2) * (SQ(XdelX) - X*XdelX - 2*pt->d*XdelX + XdelX + SQ(pt->d) - pt->d);
        ++pt;
    }
    phir += exp(-delta*delta) * sum;

    sum = 0;
    X = 3;
    XdelX = 3*delta*delta*delta;
    for(i=17; i<21; ++i){
        sum += pt->a * pow(tau, pt->t) * ipow(delta, pt->d - 2) * (SQ(XdelX) - X*XdelX - 2*pt->d*XdelX + XdelX + SQ(pt->d) - pt->d);
        ++pt;
    }
    phir += exp(-delta*delta*delta) * sum;

    return phir;
}

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	/* forward decls */
39
40	static double helm_resid(double tau, double delta, const HelmholtzData *data);
41	static double helm_resid_del(double tau, double delta, const HelmholtzData *data);
42	static double helm_resid_tau(double tau, double delta, const HelmholtzData *data);
43	static double helm_resid_deltau(double tau, double delta, const HelmholtzData *data);
44	static double helm_resid_deldel(double tau, double delta, const HelmholtzData *data);
45
46	/**
47	Function to calculate pressure from Helmholtz free energy EOS, given temperature
48	and mass density.
49
50	@param T temperature in K
51	@param rho mass density in kg/m続
52	@return pressure in Pa???
53	*/
54	double helmholtz_p(double T, double rho, const HelmholtzData *data){
55
56	double tau = data->T_star / T;
57	double delta = rho / data->rho_star;
58
59	#ifdef TEST
60	assert(data->rho_star!=0);
61	assert(T!=0);
62	assert(!isnan(tau));
63	assert(!isnan(delta));
64	assert(!isnan(data->R));
65
66	//fprintf(stderr,"p calc: T = %f\n",T);
67	//fprintf(stderr,"p calc: tau = %f\n",tau);
68	//fprintf(stderr,"p calc: rho = %f\n",rho);
69	//fprintf(stderr,"p calc: delta = %f\n",delta);
70	//fprintf(stderr,"p calc: RTrho = %f\n",data->R * T * rho);
71
72	//fprintf(stderr,"T = %f\n", T);
73	//fprintf(stderr,"rhob = %f, rhob* = %f, delta = %f\n", rho/data->M, data->rho_star/data->M, delta);
74	#endif
75
76	return data->R * T * rho * (1 + delta * helm_resid_del(tau,delta,data));
77	}
78
79	/**
80	Function to calculate internal energy from Helmholtz free energy EOS, given
81	temperature and mass density.
82
83	@param T temperature in K
84	@param rho mass density in kg/m続
85	@return internal energy in ???
86	*/
87	double helmholtz_u(double T, double rho, const HelmholtzData *data){
88
89	double tau = data->T_star / T;
90	double delta = rho / data->rho_star;
91
92	#ifdef TEST
93	assert(data->rho_star!=0);
94	assert(T!=0);
95	assert(!isnan(tau));
96	assert(!isnan(delta));
97	assert(!isnan(data->R));
98	#endif
99
100	#if 0
101	fprintf(stderr,"ideal_tau = %f\n",helm_ideal_tau(tau,delta,data->ideal));
102	fprintf(stderr,"resid_tau = %f\n",helm_resid_tau(tau,delta,data));
103	fprintf(stderr,"R T = %f\n",data->R * data->T_star);
104	#endif
105
106	return data->R * data->T_star * (helm_ideal_tau(tau,delta,data->ideal) + helm_resid_tau(tau,delta,data));
107	}
108
109	/**
110	Function to calculate enthalpy from Helmholtz free energy EOS, given
111	temperature and mass density.
112
113	@param T temperature in K
114	@param rho mass density in kg/m続
115	@return enthalpy in J/kg
116	*/
117	double helmholtz_h(double T, double rho, const HelmholtzData *data){
118
119	double tau = data->T_star / T;
120	double delta = rho / data->rho_star;
121
122	#ifdef TEST
123	assert(data->rho_star!=0);
124	assert(T!=0);
125	assert(!isnan(tau));
126	assert(!isnan(delta));
127	assert(!isnan(data->R));
128	#endif
129
130	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));
131	}
132
133	/**
134	Function to calculate entropy from Helmholtz free energy EOS, given
135	temperature and mass density.
136
137	@param T temperature in K
138	@param rho mass density in kg/m続
139	@return entropy in J/kgK
140	*/
141	double helmholtz_s(double T, double rho, const HelmholtzData *data){
142
143	double tau = data->T_star / T;
144	double delta = rho / data->rho_star;
145
146	#ifdef ENTROPY_DEBUG
147	assert(data->rho_star!=0);
148	assert(T!=0);
149	assert(!isnan(tau));
150	assert(!isnan(delta));
151	assert(!isnan(data->R));
152
153	fprintf(stderr,"helm_ideal_tau = %f\n",helm_ideal_tau(tau,delta,data->ideal));
154	fprintf(stderr,"helm_resid_tau = %f\n",helm_resid_tau(tau,delta,data));
155	fprintf(stderr,"helm_ideal = %f\n",helm_ideal(tau,delta,data->ideal));
156	fprintf(stderr,"helm_resid = %f\n",helm_resid(tau,delta,data));
157	#endif
158	return data->R * (
159	tau * (helm_ideal_tau(tau,delta,data->ideal) + helm_resid_tau(tau,delta,data))
160	- (helm_ideal(tau,delta,data->ideal) + helm_resid(tau,delta,data))
161	);
162	}
163
164	/**
165	Function to calculate Helmholtz energy from the Helmholtz free energy EOS,
166	given temperature and mass density.
167
168	@param T temperature in K
169	@param rho mass density in kg/m続
170	@return Helmholtz energy 'a', in J/kg
171	*/
172	double helmholtz_a(double T, double rho, const HelmholtzData *data){
173
174	double tau = data->T_star / T;
175	double delta = rho / data->rho_star;
176
177	#ifdef TEST
178	assert(data->rho_star!=0);
179	assert(T!=0);
180	assert(!isnan(tau));
181	assert(!isnan(delta));
182	assert(!isnan(data->R));
183	#endif
184
185	#ifdef HELMHOLTZ_DEBUG
186	fprintf(stderr,"helmholtz_a: T = %f, rho = %f\n",T,rho);
187	fprintf(stderr,"multiplying by RT = %f\n",data->R*T);
188	#endif
189
190	return data->R * T * (helm_ideal(tau,delta,data->ideal) + helm_resid(tau,delta,data));
191	}
192
193
194	/**
195	Calculation zero-pressure specific heat capacity
196	*/
197	double helmholtz_cp0(double T, const HelmholtzData *data){
198	double val = helm_cp0(T,data->ideal);
199	#if 0
200	fprintf(stderr,"val = %f\n",val);
201	#endif
202	return val;
203	}
204
205	/*---------------------------------------------
206	UTILITY FUNCTION(S)
207	*/
208
209	/* ipow: public domain by Mark Stephen with suggestions by Keiichi Nakasato */
210	static double ipow(double x, int n){
211	double t = 1.0;
212
213	if(!n)return 1.0; /* At the top. x^0 = 1 */
214
215	if (n < 0){
216	n = -n;
217	x = 1.0/x; /* error if x == 0. Good */
218	} /* ZTC/SC returns inf, which is even better */
219
220	if (x == 0.0)return 0.0;
221
222	do{
223	if(n & 1)t *= x;
224	n /= 2; /* KN prefers if (n/=2) x=x; This avoids an /
225	x = x; / unnecessary but benign multiplication on */
226	}while(n); /* the last pass, but the comparison is always
227	true _except_ on the last pass. */
228
229	return t;
230	}
231
232	//#define RESID_DEBUG
233
234	/**
235	Residual part of helmholtz function.
236	*/
237	double helm_resid(double tau, double delta, const HelmholtzData *data){
238	double dell,ldell, term, sum, res = 0;
239	unsigned n, i;
240	const HelmholtzPowTerm *pt;
241	const HelmholtzGausTerm *gt;
242
243	n = data->np;
244	pt = &(data->pt[0]);
245
246	#ifdef RESID_DEBUG
247	fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
248	#endif
249
250	/* power terms */
251	sum = 0;
252	dell = ipow(delta,pt->l);
253	ldell = pt->l * dell;
254	unsigned oldl;
255	for(i=0; i<n; ++i){
256	term = pt->a * pow(tau, pt->t) * ipow(delta, pt->d);
257	sum += term;
258	#ifdef RESID_DEBUG
259	fprintf(stderr,"i = %d, a=%e, t=%f, d=%d, term = %f, sum = %f",i,pt->a,pt->t,pt->d,term,sum);
260	if(pt->l==0){
261	fprintf(stderr,",row=%e\n",term);
262	}else{
263	fprintf(stderr,",row=%e\n,",term*exp(-dell));
264	}
265	#endif
266	oldl = pt->l;
267	++pt;
268	if(i+1==n \|\| oldl != pt->l){
269	if(oldl == 0){
270	#ifdef RESID_DEBUG
271	fprintf(stderr,"linear ");
272	#endif
273	res += sum;
274	}else{
275	#ifdef RESID_DEBUG
276	fprintf(stderr,"exp dell=%f, exp(-dell)=%f sum=%f: ",dell,exp(-dell),sum);
277	#endif
278	res += sum * exp(-dell);
279	}
280	#ifdef RESID_DEBUG
281	fprintf(stderr,"i = %d, res = %f\n",i,res);
282	#endif
283	sum = 0;
284	dell = ipow(delta,pt->l);
285	ldell = pt->l*dell;
286	}
287	}
288
289	#if 1
290	/* gaussian terms */
291	n = data->ng;
292	//fprintf(stderr,"THERE ARE %d GAUSSIAN TERMS\n",n);
293	gt = &(data->gt[0]);
294	for(i=0; i<n; ++i){
295	#ifdef RESID_DEBUG
296	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);
297	#endif
298	double d1 = delta - gt->epsilon;
299	double t1 = tau - gt->gamma;
300	double e1 = -gt->alphad1d1 - gt->betat1t1;
301	sum = gt->n * pow(tau,gt->t) * pow(delta,gt->d) * exp(e1);
302	//fprintf(stderr,"sum = %f\n",sum);
303	res += sum;
304	++gt;
305	}
306	#endif
307
308	#ifdef RESID_DEBUG
309	fprintf(stderr,"phir = %f\n",res);
310	#endif
311	return res;
312	}
313
314	/**
315	Derivative of the helmholtz residual function with respect to
316	delta.
317	*/
318	double helm_resid_del(double tau,double delta, const HelmholtzData *data){
319	double sum, res = 0;
320	double dell, ldell;
321	unsigned n, i;
322	const HelmholtzPowTerm *pt;
323	const HelmholtzGausTerm *gt;
324
325	n = data->np;
326	pt = &(data->pt[0]);
327
328	#ifdef RESID_DEBUG
329	fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
330	#endif
331
332	sum = 0;
333	dell = ipow(delta,pt->l);
334	ldell = pt->l * dell;
335	unsigned oldl;
336	for(i=0; i<n; ++i){
337	sum += pt->a * pow(tau, pt->t) * ipow(delta, pt->d - 1) * (pt->d - ldell);
338	oldl = pt->l;
339	++pt;
340	if(i+1==n \|\| oldl != pt->l){
341	if(oldl == 0){
342	res += sum;
343	}else{
344	res += sum * exp(-dell);
345	}
346	sum = 0;
347	dell = ipow(delta,pt->l);
348	ldell = pt->l*dell;
349	}
350	}
351
352	#if 1
353	/* gaussian terms */
354	n = data->ng;
355	//fprintf(stderr,"THERE ARE %d GAUSSIAN TERMS\n",n);
356	gt = &(data->gt[0]);
357	for(i=0; i<n; ++i){
358	#ifdef RESID_DEBUG
359	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);
360	#endif
361	#define SQ(X) ((X)*(X))
362	double val2;
363	val2 = - gt->n * pow(tau,gt->t) * pow(delta, -1. + gt->d)
364	* (2. * gt->alpha * delta * (delta - gt->epsilon) - gt->d)
365	* exp(-(gt->alpha * SQ(delta-gt->epsilon) + gt->beta*SQ(tau-gt->gamma)));
366	res += val2;
367	#ifdef RESID_DEBUG
368	fprintf(stderr,"val2 = %f --> res = %f\n",val2,res);
369	#endif
370	#undef SQ
371	++gt;
372	}
373	#endif
374
375	return res;
376	}
377
378	/**
379	Derivative of the helmholtz residual function with respect to
380	tau.
381	*/
382	double helm_resid_tau(double tau,double delta,const HelmholtzData *data){
383
384	double sum;
385	double res = 0;
386	double delX;
387	unsigned l;
388	unsigned n, i;
389	const HelmholtzPowTerm *pt;
390	const HelmholtzGausTerm *gt;
391
392	n = data->np;
393	pt = &(data->pt[0]);
394
395	delX = 1;
396
397	l = 0;
398	sum = 0;
399	for(i=0; i<n; ++i){
400	if(pt->t){
401	//fprintf(stderr,"i = %d, a = %e, t = %f, d = %d, l = %d\n",i+1, pt->a, pt->t, pt->d, pt->l);
402	sum += pt->a * pow(tau, pt->t - 1) * ipow(delta, pt->d) * pt->t;
403	}
404	++pt;
405	//fprintf(stderr,"l = %d\n",l);
406	if(i+1==n \|\| l != pt->l){
407	if(l==0){
408	//fprintf(stderr,"Adding non-exp term\n");
409	res += sum;
410	}else{
411	//fprintf(stderr,"Adding exp term with l = %d, delX = %e\n",l,delX);
412	res += sum * exp(-delX);
413	}
414	/* set l to new value */
415	if(i+1!=n){
416	l = pt->l;
417	//fprintf(stderr,"New l = %d\n",l);
418	delX = ipow(delta,l);
419	sum = 0;
420	}
421	}
422	}
423
424	//#define RESID_DEBUG
425	/* gaussian terms */
426	n = data->ng;
427	gt = &(data->gt[0]);
428	for(i=0; i<n; ++i){
429	#ifdef RESID_DEBUG
430	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);
431	#endif
432
433	#define SQ(X) ((X)*(X))
434	double val2;
435	val2 = -gt->n * pow(tau,gt->t - 1.) * pow(delta, gt->d)
436	* (2. * gt->beta * tau * (tau - gt->gamma) - gt->t)
437	* exp(-(gt->alpha * SQ(delta-gt->epsilon) + gt->beta*SQ(tau-gt->gamma)));
438	res += val2;
439	#ifdef RESID_DEBUG
440	fprintf(stderr,"res = %f\n",res);
441	#endif
442	#undef SQ
443
444	++gt;
445	}
446
447	return res;
448	}
449
450
451
452	/**
453	Mixed derivative of the helmholtz residual function with respect to
454	delta and tau
455
456	FIXME this function is WRONG.
457	*/
458	double helm_resid_deltau(double tau,double delta,const HelmholtzData *data){
459
460	double sum;
461	double phir = 0;
462	unsigned i;
463	double XdelX;
464
465	const HelmholtzPowTerm *pt = &(data->pt[0]);
466
467	for(i=0; i<5; ++i){
468	phir += pt->a * pow(tau, pt->t - 1) * ipow(delta, pt->d - 1) * pt->d * pt->t;
469	++pt;
470	}
471
472	sum = 0;
473	XdelX = delta;
474	for(i=5; i<10; ++i){
475	sum += pt->a * pow(tau, pt->t - 1) * ipow(delta, pt->d - 1) * pt->t *(pt->d - XdelX);
476	++pt;
477	}
478	phir += exp(-delta) * sum;
479
480	sum = 0;
481	XdelX = 2deltadelta;
482	for(i=10; i<17; ++i){
483	sum += pt->a * pow(tau, pt->t - 1) * ipow(delta, pt->d - 1) * pt->t *(pt->d - XdelX);
484	++pt;
485	}
486	phir += exp(-deltadelta) sum;
487
488	sum = 0;
489	XdelX = 3deltadelta*delta;
490	for(i=17; i<21; ++i){
491	sum += pt->a * pow(tau, pt->t - 1) * ipow(delta, pt->d - 1) * pt->t *(pt->d - XdelX);
492	++pt;
493	}
494	phir += exp(-deltadeltadelta) * sum;
495
496	return phir;
497	}
498
499	#define SQ(X) ((X)*(X))
500
501	/**
502	Second derivative of helmholtz residual function with respect to
503	delta (twice).
504
505	FIXME this function is WRONG.
506	*/
507	double helm_resid_deldel(double tau,double delta,const HelmholtzData *data){
508
509	double sum;
510	double phir = 0;
511	unsigned i;
512	unsigned X;
513	double XdelX;
514
515	const HelmholtzPowTerm *pt = &(data->pt[0]);
516
517	for(i=0; i<5; ++i){
518	phir += pt->a * pow(tau, pt->t) * ipow(delta, pt->d - 2) * (SQ(pt->d) - X);
519	++pt;
520	}
521
522	sum = 0;
523	X = 1;
524	XdelX = delta;
525	for(i=5; i<10; ++i){
526	sum += pt->a * pow(tau, pt->t) * ipow(delta, pt->d - 2) * (SQ(XdelX) - XXdelX - 2pt->d*XdelX + XdelX + SQ(pt->d) - pt->d);
527	++pt;
528	}
529	phir += exp(-delta) * sum;
530
531	sum = 0;
532	X = 2;
533	XdelX = 2deltadelta;
534	for(i=10; i<17; ++i){
535	sum += pt->a * pow(tau, pt->t) * ipow(delta, pt->d - 2) * (SQ(XdelX) - XXdelX - 2pt->d*XdelX + XdelX + SQ(pt->d) - pt->d);
536	++pt;
537	}
538	phir += exp(-deltadelta) sum;
539
540	sum = 0;
541	X = 3;
542	XdelX = 3deltadelta*delta;
543	for(i=17; i<21; ++i){
544	sum += pt->a * pow(tau, pt->t) * ipow(delta, pt->d - 2) * (SQ(XdelX) - XXdelX - 2pt->d*XdelX + XdelX + SQ(pt->d) - pt->d);
545	++pt;
546	}
547	phir += exp(-deltadeltadelta) * sum;
548
549	return phir;
550	}
551