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"

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

/* forward decls */

static double helm_ideal(double tau, double delta, const IdealData *data);
static double helm_ideal_tau(double tau, double delta, const IdealData *data);
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

#ifdef TEST
    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 TEST
    assert(data->rho_star!=0);
    assert(T!=0);
    assert(!isnan(tau));
    assert(!isnan(delta));
    assert(!isnan(data->R));
#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)
    );
}

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

/*---------------------------------------------
  IDEAL COMPONENT RELATIONS
*/

/**
    Ideal component of helmholtz function
*/  
double helm_ideal(double tau, double delta, const IdealData *data){

    const IdealPowTerm *pt;
    const IdealExpTerm *et;

    unsigned i;
    double sum = log(delta) + data->c + data->m * tau - log(tau);
    double term;

    //fprintf(stderr,"constant = %f, linear = %f", data->c, data->m);
    //fprintf(stderr,"initial terms = %f\n",sum);
    pt = &(data->pt[0]);
    for(i = 0; i<data->np; ++i, ++pt){
        term = pt->a0 * pow(tau, pt->t0);
        //fprintf(stderr,"i = %d: a0 = %f, t0 = %f, term = %f\n",i,pt->a0, pt->t0, term);
        sum += pt->a0 * pow(tau, pt->t0);
    }

#if 0
    et = &(data->et[0]);
    for(i=0; i<data->ne; ++i, ++et){
        sum += et->b * log( 1 - exp(- et->B * tau));
    }
#endif

    return sum;
}

/**
    Partial dervivative of ideal component of helmholtz residual function with 
    respect to tau.
*/  
double helm_ideal_tau(double tau, double delta, const IdealData *data){
    const IdealPowTerm *pt;
    const IdealExpTerm *et;

    unsigned i;
    double sum = -1./tau + data->m;

    pt = &(data->pt[0]);
    for(i = 0; i<data->np; ++i, ++pt){
        sum += pt->a0 * pt->t0 * pow(tau, pt->t0 - 1);
    }

#if 0
    et = &(data->et[0]);
    for(i=0; i<data->ne; ++i, ++et){
        sum += et->b * log( 1 - exp(- et->B * tau));
    }
#endif
    return sum;
}   

/**
    Residual part of helmholtz function. Note: we have NOT prematurely
    optimised here ;-)
*/
double helm_resid(double tau, double delta, const HelmholtzData *data){
    double sum;
    double res = 0;
    double delX;
    unsigned l;
    unsigned np, i;
    const HelmholtzPowTerm *pt;

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

    delX = 1;

    l = 0;
    sum = 0;
    for(i=0; i<np; ++i){
        //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) * ipow(delta, pt->d);
        ++pt;
        //fprintf(stderr,"l = %d\n",l);
        if(i+1==np || 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!=np){
                l = pt->l;
                //fprintf(stderr,"New l = %d\n",l);
                delX = ipow(delta,l);
                sum = 0;
            }
        }
    }

    return res;
}

/**
    Derivative of the helmholtz residual function with respect to
    delta.

    NOTE: POSSIBLY STILL AN ERROR IN THIS FUNCTION.
*/  
double helm_resid_del(double tau,double delta, const HelmholtzData *data){
    double term, x, sum;
    double res = 0;
    double delX, XdelX;
    unsigned l;
    unsigned n, i;
    const HelmholtzPowTerm *pt;
    const HelmholtzExpTerm *et;

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


    sum = 0;


    for(i=0; i<n; ++i){
        term = pt->a * pow(tau, pt->t) * pow(delta, pt->d - 1) * (pt->d - pt->l*pow(delta,pt->l));
        if(pt->l==0){   
            x = 1;
        }else{
            x = exp(-pow(delta,pt->l));
        }
        //fprintf(stderr,"i = %d, a = %e, t = %f, d = %d, l = %d --> pow = %f, exp = %f, term = %f\n",i+1, pt->a, pt->t, pt->d, pt->l, term, x, term*x);
        res += term * x;
        ++pt;
    }
#if 0
    delX = 1;
    l = 0;
    XdelX = 0;
    for(i=0; i<n; ++i){
        term = pt->a * pow(tau, pt->t) * ipow(delta, pt->d - 1) * (pt->d - XdelX);
        fprintf(stderr,"i = %d, a = %e, t = %f, d = %d, l = %d --> term = %f\n",i+1, pt->a, pt->t, pt->d, pt->l, term);
        sum += term;
        ++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);
            }
            fprintf(stderr,"sum = %f, mult = %f, res = %f\n",sum,exp(-delX),res);
            /* set l to new value */
            if(i+1!=n){
                l = pt->l;
                delX = ipow(delta,l);
                XdelX = l * delX;
                //fprintf(stderr,"New l = %d, XdelX = %f\n",l,XdelX);
                sum = 0;
                fprintf(stderr,"sum zero\n");
            }
        }
    }
#endif

#if 1
    /* now the exponential terms */
    n = data->ne;
    et = &(data->et[0]);
    for(i=0; i< n; ++i){
        //fprintf(stderr,"i = %d, a = %e, t = %f, d = %d, phi = %d, beta = %d, gamma = %f\n",i+1, et->a, et->t, et->d, et->phi, et->beta, et->gamma);
        
        double del2 = delta*delta;
        double tau2 = tau*tau;
        double gam2 = et->gamma * et->gamma;
        sum = -et->a * pow(tau,et->t) * ipow(delta,et->d-1)
            * (2 * et->phi * del2 - 2 * et->phi * delta - et->d)
            * exp(-et->phi * del2
                     + 2 * et->phi * delta
                     - et->beta * tau2
                     + 2 * et->beta * et->gamma * tau
                     - et->phi 
                     - et->beta * gam2
               );
        //fprintf(stderr,"sum = %f\n",sum);
        res += sum;
        ++et;
    }
#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 np, i;
    const HelmholtzPowTerm *pt;

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

    delX = 1;

    l = 0;
    sum = 0;
    for(i=0; i<np; ++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==np || 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!=np){
                l = pt->l;
                //fprintf(stderr,"New l = %d\n",l);
                delX = ipow(delta,l);
                sum = 0;
            }
        }
    }

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