/[ascend]/trunk/models/johnpye/fprops/helmholtz.c
ViewVC logotype

Annotation of /trunk/models/johnpye/fprops/helmholtz.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2275 - (hide annotations) (download) (as text)
Wed Aug 11 13:58:25 2010 UTC (13 years, 11 months ago) by jpye
File MIME type: text/x-csrc
File size: 38732 byte(s)
Fixing error return from asc_helmholtz function.
The bug with 'defaultall' functionality is still here... see boiler_simple_test mode in models/johnpye/fprops/rankine_fprops.a4c.
1 jpye 1810 /* ASCEND modelling environment
2 jpye 2113 Copyright (C) 2008-2009 Carnegie Mellon University
3 jpye 1810
4     This program is free software; you can redistribute it and/or modify
5     it under the terms of the GNU General Public License as published by
6     the Free Software Foundation; either version 2, or (at your option)
7     any later version.
8    
9     This program is distributed in the hope that it will be useful,
10     but WITHOUT ANY WARRANTY; without even the implied warranty of
11     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12     GNU General Public License for more details.
13    
14     You should have received a copy of the GNU General Public License
15     along with this program; if not, write to the Free Software
16     Foundation, Inc., 59 Temple Place - Suite 330,
17     Boston, MA 02111-1307, USA.
18     *//** @file
19     Implementation of the reduced molar Helmholtz free energy equation of state.
20    
21 jpye 1822 For nomenclature see Tillner-Roth, Harms-Watzenberg and Baehr, Eine neue
22     Fundamentalgleichung f端r Ammoniak.
23    
24 jpye 1810 John Pye, 29 Jul 2008.
25     */
26    
27 jpye 1822 #include <math.h>
28 jpye 1810
29 jpye 1822 #include "helmholtz.h"
30 jpye 1847 #include "ideal_impl.h"
31 jpye 2265 #include "sat.h"
32 jpye 1822
33 jpye 2271 #if 0
34     # include <assert.h>
35     #else
36     # define assert(ARGS...)
37     #endif
38    
39 jpye 1826 #include <stdlib.h>
40     #include <stdio.h>
41    
42 jpye 1996 /* macros and forward decls */
43    
44 jpye 1905 #define SQ(X) ((X)*(X))
45    
46 jpye 1985 #include "helmholtz_impl.h"
47 jpye 1822
48 jpye 1996 /* calculate tau and delta using a macro -- is used in most functions */
49     #define DEFINE_TD \
50     double tau = data->T_star / T; \
51     double delta = rho / data->rho_star
52    
53 jpye 2236 double helmholtz_p(double T, double rho, const HelmholtzData *d){
54     double p, rho_f, rho_g;
55 jpye 2270 #if 0
56 jpye 2260 if(T < d->T_t){
57     fprintf(stderr,"%s: Unable to calculate pressure, T = %e is below triple point.\n", __func__, T);
58     return d->p_t;
59     }
60 jpye 2270 /* but what if we're in the sublimation region?? */
61     #endif
62 jpye 2246 if(T < d->T_c){
63 jpye 2236 int res = fprops_sat_T(T, &p, &rho_f, &rho_g, d);
64     if(res){
65     //fprintf(stderr,"ERROR: got error % from saturation calc in %s",res,__func__);
66     return p;
67     }
68     if(rho < rho_f && rho > rho_g){
69     return p;
70     }
71     }
72     return helmholtz_p_raw(T,rho,d);
73     }
74    
75 jpye 2252 double helmholtz_h(double T, double rho, const HelmholtzData *d){
76     double p, rho_f, rho_g;
77     if(T < d->T_c){
78     int res = fprops_sat_T(T, &p, &rho_f, &rho_g, d);
79     if(res){
80     //fprintf(stderr,"ERROR: got error % from saturation calc in %s",res,__func__);
81     return d->rho_c;
82     }
83     if(rho < rho_f && rho > rho_g){
84     double x = rho_g*(rho_f/rho - 1)/(rho_f - rho_g);
85     return x*helmholtz_h_raw(T,rho_g,d) + (1.-x)*helmholtz_h_raw(T,rho_f,d);
86     }
87     }
88     return helmholtz_h_raw(T,rho,d);
89     }
90    
91 jpye 2255 double helmholtz_s(double T, double rho, const HelmholtzData *d){
92     double p, rho_f, rho_g;
93     if(T < d->T_c){
94     int res = fprops_sat_T(T, &p, &rho_f, &rho_g, d);
95     if(res){
96     //fprintf(stderr,"ERROR: got error % from saturation calc in %s",res,__func__);
97     return d->rho_c;
98     }
99     if(rho < rho_f && rho > rho_g){
100     double x = rho_g*(rho_f/rho - 1)/(rho_f - rho_g);
101     return x*helmholtz_s_raw(T,rho_g,d) + (1.-x)*helmholtz_s_raw(T,rho_f,d);
102     }
103     }
104     return helmholtz_s_raw(T,rho,d);
105     }
106 jpye 2252
107 jpye 1810 /**
108     Function to calculate pressure from Helmholtz free energy EOS, given temperature
109 jpye 1824 and mass density.
110 jpye 1810
111     @param T temperature in K
112 jpye 1824 @param rho mass density in kg/m続
113     @return pressure in Pa???
114 jpye 1810 */
115 jpye 2236 double helmholtz_p_raw(double T, double rho, const HelmholtzData *data){
116 jpye 1996 DEFINE_TD;
117 jpye 1810
118 jpye 1838 assert(data->rho_star!=0);
119     assert(T!=0);
120 jpye 2270 assert(!__isnan(T));
121     assert(!__isnan(rho));
122     assert(!__isnan(data->R));
123 jpye 1839
124 jpye 1840 //fprintf(stderr,"p calc: T = %f\n",T);
125     //fprintf(stderr,"p calc: tau = %f\n",tau);
126     //fprintf(stderr,"p calc: rho = %f\n",rho);
127     //fprintf(stderr,"p calc: delta = %f\n",delta);
128     //fprintf(stderr,"p calc: R*T*rho = %f\n",data->R * T * rho);
129 jpye 1841
130     //fprintf(stderr,"T = %f\n", T);
131     //fprintf(stderr,"rhob = %f, rhob* = %f, delta = %f\n", rho/data->M, data->rho_star/data->M, delta);
132 jpye 1839
133 jpye 2270 double p = data->R * T * rho * (1 + delta * helm_resid_del(tau,delta,data));
134     if(isnan(p)){
135     fprintf(stderr,"T = %.12e, rho = %.12e\n",T,rho);
136     }
137     assert(!__isnan(p));
138     return p;
139 jpye 1822 }
140 jpye 1810
141 jpye 1824 /**
142     Function to calculate internal energy from Helmholtz free energy EOS, given
143     temperature and mass density.
144    
145     @param T temperature in K
146     @param rho mass density in kg/m続
147     @return internal energy in ???
148     */
149 jpye 1826 double helmholtz_u(double T, double rho, const HelmholtzData *data){
150 jpye 1996 DEFINE_TD;
151 jpye 1824
152 jpye 1832 #ifdef TEST
153 jpye 1838 assert(data->rho_star!=0);
154     assert(T!=0);
155     assert(!isnan(tau));
156     assert(!isnan(delta));
157     assert(!isnan(data->R));
158     #endif
159    
160 jpye 1857 #if 0
161 jpye 1835 fprintf(stderr,"ideal_tau = %f\n",helm_ideal_tau(tau,delta,data->ideal));
162 jpye 1827 fprintf(stderr,"resid_tau = %f\n",helm_resid_tau(tau,delta,data));
163     fprintf(stderr,"R T = %f\n",data->R * data->T_star);
164 jpye 1832 #endif
165 jpye 1827
166 jpye 1835 return data->R * data->T_star * (helm_ideal_tau(tau,delta,data->ideal) + helm_resid_tau(tau,delta,data));
167 jpye 1824 }
168    
169     /**
170     Function to calculate enthalpy from Helmholtz free energy EOS, given
171     temperature and mass density.
172    
173     @param T temperature in K
174     @param rho mass density in kg/m続
175 jpye 1829 @return enthalpy in J/kg
176 jpye 1824 */
177 jpye 2252 double helmholtz_h_raw(double T, double rho, const HelmholtzData *data){
178 jpye 1996 DEFINE_TD;
179 jpye 1824
180 jpye 2268 //#ifdef TEST
181 jpye 1838 assert(data->rho_star!=0);
182     assert(T!=0);
183     assert(!isnan(tau));
184     assert(!isnan(delta));
185     assert(!isnan(data->R));
186 jpye 2268 //#endif
187     double h = 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));
188     assert(!__isnan(h));
189     return h;
190 jpye 1824 }
191    
192 jpye 1829 /**
193     Function to calculate entropy from Helmholtz free energy EOS, given
194     temperature and mass density.
195    
196     @param T temperature in K
197     @param rho mass density in kg/m続
198     @return entropy in J/kgK
199     */
200 jpye 2255 double helmholtz_s_raw(double T, double rho, const HelmholtzData *data){
201 jpye 1996 DEFINE_TD;
202 jpye 1829
203 jpye 1870 #ifdef ENTROPY_DEBUG
204 jpye 1838 assert(data->rho_star!=0);
205     assert(T!=0);
206     assert(!isnan(tau));
207     assert(!isnan(delta));
208     assert(!isnan(data->R));
209    
210 jpye 1845 fprintf(stderr,"helm_ideal_tau = %f\n",helm_ideal_tau(tau,delta,data->ideal));
211     fprintf(stderr,"helm_resid_tau = %f\n",helm_resid_tau(tau,delta,data));
212     fprintf(stderr,"helm_ideal = %f\n",helm_ideal(tau,delta,data->ideal));
213     fprintf(stderr,"helm_resid = %f\n",helm_resid(tau,delta,data));
214 jpye 1847 #endif
215 jpye 1829 return data->R * (
216 jpye 1835 tau * (helm_ideal_tau(tau,delta,data->ideal) + helm_resid_tau(tau,delta,data))
217 jpye 1871 - (helm_ideal(tau,delta,data->ideal) + helm_resid(tau,delta,data))
218 jpye 1829 );
219     }
220    
221 jpye 1850 /**
222 jpye 1863 Function to calculate Helmholtz energy from the Helmholtz free energy EOS,
223     given temperature and mass density.
224    
225     @param T temperature in K
226     @param rho mass density in kg/m続
227     @return Helmholtz energy 'a', in J/kg
228     */
229     double helmholtz_a(double T, double rho, const HelmholtzData *data){
230 jpye 1996 DEFINE_TD;
231 jpye 1863
232     #ifdef TEST
233     assert(data->rho_star!=0);
234     assert(T!=0);
235     assert(!isnan(tau));
236     assert(!isnan(delta));
237     assert(!isnan(data->R));
238     #endif
239    
240 jpye 1868 #ifdef HELMHOLTZ_DEBUG
241 jpye 1863 fprintf(stderr,"helmholtz_a: T = %f, rho = %f\n",T,rho);
242 jpye 1865 fprintf(stderr,"multiplying by RT = %f\n",data->R*T);
243 jpye 1863 #endif
244    
245 jpye 1870 return data->R * T * (helm_ideal(tau,delta,data->ideal) + helm_resid(tau,delta,data));
246 jpye 1863 }
247    
248 jpye 1988 /**
249     Function to calculate isochoric heat capacity from the Helmholtz free energy
250     EOS given temperature and mass density.
251 jpye 1863
252 jpye 1988 @param T temperature in K
253     @param rho mass density in kg/m続
254 jpye 1989 @return Isochoric specific heat capacity in J/kg/K.
255 jpye 1988 */
256     double helmholtz_cv(double T, double rho, const HelmholtzData *data){
257 jpye 1996 DEFINE_TD;
258 jpye 1988
259 jpye 1993 return - data->R * SQ(tau) * (helm_ideal_tautau(tau,data->ideal) + helm_resid_tautau(tau,delta,data));
260 jpye 1988 }
261    
262 jpye 1863 /**
263 jpye 1989 Function to calculate isobaric heat capacity from the Helmholtz free energy
264     EOS given temperature and mass density.
265    
266     @param T temperature in K
267     @param rho mass density in kg/m続
268     @return Isobaric specific heat capacity in J/kg/K.
269     */
270     double helmholtz_cp(double T, double rho, const HelmholtzData *data){
271 jpye 1996 DEFINE_TD;
272 jpye 1989
273     double phir_d = helm_resid_del(tau,delta,data);
274     double phir_dd = helm_resid_deldel(tau,delta,data);
275     double phir_dt = helm_resid_deltau(tau,delta,data);
276    
277     /* note similarities with helmholtz_w */
278     double temp1 = 1 + 2*delta*phir_d + SQ(delta)*phir_dd;
279     double temp2 = 1 + delta*phir_d - delta*tau*phir_dt;
280 jpye 1996 double temp3 = -SQ(tau)*(helm_ideal_tautau(tau,data->ideal) + helm_resid_tautau(tau,delta,data));
281 jpye 1989
282 jpye 1996 return data->R * (temp3 + SQ(temp2)/temp1);
283 jpye 1989 }
284    
285    
286     /**
287     Function to calculate the speed of sound in a fluid from the Helmholtz free
288     energy EOS, given temperature and mass density.
289    
290     @param T temperature in K
291     @param rho mass density in kg/m続
292     @return Speed of sound in m/s.
293     */
294     double helmholtz_w(double T, double rho, const HelmholtzData *data){
295 jpye 1996 DEFINE_TD;
296 jpye 1989
297     double phir_d = helm_resid_del(tau,delta,data);
298     double phir_dd = helm_resid_deldel(tau,delta,data);
299     double phir_dt = helm_resid_deltau(tau,delta,data);
300    
301     /* note similarities with helmholtz_cp */
302 jpye 1994 double temp1 = 1. + 2.*delta*phir_d + SQ(delta)*phir_dd;
303     double temp2 = 1. + delta*phir_d - delta*tau*phir_dt;
304 jpye 1996 double temp3 = -SQ(tau)*(helm_ideal_tautau(tau,data->ideal) + helm_resid_tautau(tau,delta,data));
305 jpye 1989
306 jpye 1996 return sqrt(data->R * T * (temp1 + SQ(temp2)/temp3));
307 jpye 1989
308     }
309    
310     /**
311 jpye 2119 Function to calculate the Gibbs energy fluid from the Helmholtz free
312     energy EOS, given temperature and mass density.
313    
314     @param T temperature in K
315     @param rho mass density in kg/m続
316     @return Gibbs energy, in J/kg.
317     */
318     double helmholtz_g(double T, double rho, const HelmholtzData *data){
319     DEFINE_TD;
320    
321     double phir_d = helm_resid_del(tau,delta,data);
322     double phir = helm_resid(tau,delta,data);
323     double phi0 = helm_ideal(tau,delta,data->ideal);
324    
325 jpye 2250 return data->R * T * (phi0 + phir + 1. + delta * phir_d);
326 jpye 2119 }
327    
328     /**
329 jpye 1850 Calculation zero-pressure specific heat capacity
330     */
331     double helmholtz_cp0(double T, const HelmholtzData *data){
332     double val = helm_cp0(T,data->ideal);
333 jpye 1855 #if 0
334 jpye 1850 fprintf(stderr,"val = %f\n",val);
335 jpye 1851 #endif
336 jpye 1850 return val;
337     }
338    
339 jpye 2247 /**
340 jpye 2262 alpha_p function from IAPWS Advisory Note 3, used in calculation of
341 jpye 2247 partial property derivatives.
342     */
343     double helmholtz_alphap(double T, double rho, const HelmholtzData *data){
344     DEFINE_TD;
345     double phir_d = helm_resid_del(tau,delta,data);
346     double phir_dt = helm_resid_deltau(tau,delta,data);
347     return 1./T * (1. - delta*tau*phir_dt/(1 + delta*phir_d));
348     }
349    
350 jpye 1905 /**
351 jpye 2262 beta_p function from IAPWS Advisory Note 3 , used in calculation of partial
352     property derivatives.
353 jpye 2247 */
354     double helmholtz_betap(double T, double rho, const HelmholtzData *data){
355     DEFINE_TD;
356     double phir_d = helm_resid_del(tau,delta,data);
357     double phir_dd = helm_resid_deldel(tau,delta,data);
358 jpye 2262 return rho*(1. + (delta*phir_d + SQ(delta)*phir_dd)/(1+delta*phir_d));
359 jpye 2247 }
360    
361    
362     /**
363 jpye 1996 Solve density given temperature and pressure and bounds for density,
364     used in solving the saturation curve. We use a single-variable Newton
365     method to solve for the desired pressure by iteratively adjusting density.
366    
367     @param rho_l lower bound on density
368     @param rho_u upper bound on density
369     @param rho (returned) solved value of density
370     @return 0 on success (1 = didn't converge within permitted iterations)
371     */
372     static int helm_find_rho_tp(double T, double p, double *rho
373     , const double rho_l, const double rho_u, const HelmholtzData *data
374     ){
375     double rho_1 = 0.5 *(rho_l + rho_u);
376     double p_1, dpdrho;
377     double niter = 0, maxiter = 100;
378     int res = 1; /* 1 = exceeded iterations, 0 = converged */
379    
380 jpye 1997 #ifdef TEST
381 jpye 1996 fprintf(stderr,"\nHELM_FIND_RHO_TP: rho_l = %f, rho_u = %f (v_u = %f, v_l = %f)\n", rho_l, rho_u, 1./rho_l, 1./rho_u);
382 jpye 1997 #endif
383 jpye 1996
384     #if 0
385     double r;
386     for(r = rho_l; r <= rho_u; r += (rho_u-rho_l)/20.){
387     p_1 = helmholtz_p(T, r, data);
388     fprintf(stderr,"T = %f, rho = %f --> p = %f MPa\n",T, r, p_1/1e6);
389     }
390     #endif
391    
392     while(res){
393     if(++niter>maxiter)break;
394     p_1 = helmholtz_p(T, rho_1, data);
395     if(p_1 < p){
396     p_1 = 0.5*(p_1 + p);
397     }
398     //fprintf(stderr,"T = %f, rho = %f --> p = %f MPa, err = %f %%\n",T, rho_1, p_1/1e6, (p - p_1)/p *100);
399     dpdrho = helmholtz_dpdrho_T(T, rho_1, data);
400     if(dpdrho < 0){
401     if(rho_l < data->rho_star){
402     /* looking for gas solution */
403     rho_1 = 0.5*(*rho + rho_l);
404     }
405     }
406     rho_1 += (p - p_1)/dpdrho;
407     if(rho_1 > rho_u){
408     rho_1 = 0.5*(*rho + rho_u);
409     //fprintf(stderr,"UPPERBOUND ");
410     }
411     if(rho_1 < rho_l){
412     rho_1 = 0.5*(*rho + rho_l);
413     //fprintf(stderr,"LOWERBOUND ");
414     }
415     *rho = rho_1;
416    
417     if(fabs((p - p_1)/p) < 1e-7){
418 jpye 1997 #ifdef TEST
419 jpye 1996 fprintf(stderr,"Converged to p = %f MPa with rho = %f\n", p/1e6, rho_1);
420 jpye 1997 #endif
421 jpye 1996 res = 0;
422     }
423     }
424    
425     return res;
426     }
427    
428     /**
429     Calculation of saturation conditions given temperature
430    
431     @param T temperature [K]
432     @param rho_f (returned) saturated liquid density [kg/m続]
433     @param rho_g (returned) saturated gas density [kg/m続]
434     @param p pressure p_sat(T) [Pa]
435     @return 0 on success
436     */
437     int helmholtz_sat_t(double T, double *p, double *rho_f, double *rho_g, const HelmholtzData *data){
438     double tau = data->T_star / T;
439    
440     if(T >= data->T_star){
441 jpye 1997 #ifdef TEST
442 jpye 1996 fprintf(stderr,"ERROR: temperature exceeds critical temperature in helmholtz_sat_t.\n");
443 jpye 1997 #endif
444 jpye 1996 /* return some reasonable values */
445     *rho_f = data->rho_star;
446     *rho_g = data->rho_star;
447     *p = helmholtz_p(data->T_star, data->rho_star, data);
448     return 1; /* error status */
449     }
450    
451     /* get a first estimate of saturation pressure using acentric factor */
452    
453     /* critical pressure */
454 jpye 1997 #ifdef TEST
455 jpye 1996 fprintf(stderr,"T_c = %f, rho_c = %f\n",data->T_star, data->rho_star);
456 jpye 1997 #endif
457 jpye 1996 double p_c = helmholtz_p(data->T_star, data->rho_star, data);
458    
459 jpye 1997 #ifdef TEST
460 jpye 1996 fprintf(stderr,"Critical pressure = %f MPa\n",p_c/1.e6);
461     fprintf(stderr,"Acentric factor = %f\n",data->omega);
462 jpye 1997 #endif
463 jpye 1996
464 jpye 2114
465     /* maybe use http://dx.doi.org/10.1016/S0009-2509(02)00017-9 for this part? */
466 jpye 1996 /* FIXME need to cite this formula */
467     *p = p_c * pow(10.,(data->omega + 1.)*-7./3.*(tau - 1.));
468    
469 jpye 1997 #ifdef TEST
470 jpye 1996 fprintf(stderr,"Estimated p_sat(T=%f) = %f MPa\n",T,(*p)/1e6);
471 jpye 1997 #endif
472 jpye 1996
473     if(tau < 1.01){
474     /* close to critical point: need a different approach */
475 jpye 1997 #ifdef TEST
476 jpye 1996 fprintf(stderr,"ERROR: not implemented\n");
477 jpye 1997 #endif
478 jpye 1996 return 1;
479     }else{
480     double niter = 0;
481 jpye 1997 (void)niter;
482    
483 jpye 1996 int res = helm_find_rho_tp(T, *p, rho_f, data->rho_star,data->rho_star*10, data);
484     if(res){
485 jpye 1997 #ifdef TEST
486 jpye 1996 fprintf(stderr,"ERROR: failed to solve rho_f\n");
487 jpye 1997 #endif
488 jpye 1996 }
489    
490     res = helm_find_rho_tp(T, *p, rho_g, 0.001*data->rho_star, data->rho_star, data);
491     if(res){
492 jpye 1997 #ifdef TEST
493 jpye 1996 fprintf(stderr,"ERROR: failed to solve rho_g\n");
494 jpye 1997 #endif
495 jpye 1996 }
496    
497 jpye 1997 #ifdef TEST
498 jpye 1996 fprintf(stderr,"p = %f MPa: rho_f = %f, rho_g = %f\n", *p, *rho_f, *rho_g);
499 jpye 1997 #endif
500 jpye 1996
501     double LHS = *p/data->R/T*(1./(*rho_g) - 1./(*rho_f)) - log((*rho_f)/(*rho_g));
502     double delta_f = (*rho_f) / data->rho_star;
503     double delta_g = (*rho_g) / data->rho_star;
504     double RHS = helm_resid(delta_f,tau,data) - helm_resid(delta_g,tau,data);
505    
506     double err = LHS - RHS;
507 jpye 1997 (void)err;
508 jpye 1996
509 jpye 1997 #ifdef TEST
510 jpye 1996 fprintf(stderr,"LHS = %f, RHS = %f, err = %f\n",LHS, RHS, err);
511 jpye 1997 #endif
512 jpye 1996 /* away from critical point... */
513     *rho_f = data->rho_star;
514     *rho_g = data->rho_star;
515 jpye 1997 #ifdef TEST
516 jpye 1996 fprintf(stderr,"ERROR: not implemented\n");
517     exit(1);
518 jpye 1997 #endif
519 jpye 1996 return 1;
520     }
521     }
522    
523     /*----------------------------------------------------------------------------
524     PARTIAL DERIVATIVES
525     */
526    
527     /**
528 jpye 1905 Calculate partial derivative of p with respect to T, with rho constant
529     */
530     double helmholtz_dpdT_rho(double T, double rho, const HelmholtzData *data){
531 jpye 1996 DEFINE_TD;
532 jpye 1906
533     double phir_del = helm_resid_del(tau,delta,data);
534     double phir_deltau = helm_resid_deltau(tau,delta,data);
535 jpye 1909 #ifdef TEST
536 jpye 1906 assert(!isinf(phir_del));
537     assert(!isinf(phir_deltau));
538     assert(!isnan(phir_del));
539     assert(!isnan(phir_deltau));
540     assert(!isnan(data->R));
541     assert(!isnan(rho));
542     assert(!isnan(tau));
543 jpye 1909 #endif
544 jpye 1906
545     double res = data->R * rho * (1 + delta*phir_del - delta*tau*phir_deltau);
546    
547 jpye 1909 #ifdef TEST
548 jpye 1906 assert(!isnan(res));
549     assert(!isinf(res));
550 jpye 1909 #endif
551 jpye 1906 return res;
552 jpye 1905 }
553    
554     /**
555     Calculate partial derivative of p with respect to rho, with T constant
556     */
557     double helmholtz_dpdrho_T(double T, double rho, const HelmholtzData *data){
558 jpye 1996 DEFINE_TD;
559 jpye 1906
560     double phir_del = helm_resid_del(tau,delta,data);
561     double phir_deldel = helm_resid_deldel(tau,delta,data);
562 jpye 1909 #ifdef TEST
563 jpye 1906 assert(!isinf(phir_del));
564     assert(!isinf(phir_deldel));
565 jpye 1909 #endif
566 jpye 1993 return data->R * T * (1 + 2*delta*phir_del + SQ(delta)*phir_deldel);
567 jpye 1905 }
568    
569 jpye 2234
570     double helmholtz_d2pdrho2_T(double T, double rho, const HelmholtzData *data){
571     DEFINE_TD;
572    
573     double phir_del = helm_resid_del(tau,delta,data);
574     double phir_deldel = helm_resid_deldel(tau,delta,data);
575     double phir_deldeldel = helm_resid_deldeldel(tau,delta,data);
576     #ifdef TEST
577     assert(!isinf(phir_del));
578     assert(!isinf(phir_deldel));
579     assert(!isinf(phir_deldeldel));
580     #endif
581    
582 jpye 2235 return data->R * T / rho * delta * (2*phir_del + delta*(4*phir_deldel + delta*phir_deldeldel));
583 jpye 2234 }
584    
585 jpye 1915 /**
586     Calculate partial derivative of h with respect to T, with rho constant
587     */
588     double helmholtz_dhdT_rho(double T, double rho, const HelmholtzData *data){
589 jpye 1996 DEFINE_TD;
590 jpye 1915
591     double phir_del = helm_resid_del(tau,delta,data);
592     double phir_deltau = helm_resid_deltau(tau,delta,data);
593     double phir_tautau = helm_resid_tautau(tau,delta,data);
594     double phi0_tautau = helm_ideal_tautau(tau,data->ideal);
595    
596 jpye 1917 //fprintf(stderr,"phir_del = %f, phir_deltau = %f, phir_tautau = %f, phi0_tautau = %f\n",phir_del,phir_deltau,phir_tautau,phi0_tautau);
597    
598     //return (helmholtz_h(T+0.01,rho,data) - helmholtz_h(T,rho,data)) / 0.01;
599 jpye 1993 return data->R * (1. + delta*phir_del - SQ(tau)*(phi0_tautau + phir_tautau) - delta*tau*phir_deltau);
600 jpye 1915 }
601    
602 jpye 1920 /**
603     Calculate partial derivative of h with respect to rho, with T constant
604     */
605 jpye 1915 double helmholtz_dhdrho_T(double T, double rho, const HelmholtzData *data){
606 jpye 1996 DEFINE_TD;
607 jpye 1915
608 jpye 1919 double phir_del = helm_resid_del(tau,delta,data);
609     double phir_deltau = helm_resid_deltau(tau,delta,data);
610     double phir_deldel = helm_resid_deldel(tau,delta,data);
611    
612     return data->R * T / rho * (tau*delta*(0 + phir_deltau) + delta * phir_del + SQ(delta)*phir_deldel);
613 jpye 1915 }
614 jpye 1919
615 jpye 1920
616     /**
617     Calculate partial derivative of u with respect to T, with rho constant
618     */
619     double helmholtz_dudT_rho(double T, double rho, const HelmholtzData *data){
620 jpye 1996 DEFINE_TD;
621 jpye 1920
622     double phir_tautau = helm_resid_tautau(tau,delta,data);
623     double phi0_tautau = helm_ideal_tautau(tau,data->ideal);
624    
625     return -data->R * SQ(tau) * (phi0_tautau + phir_tautau);
626     }
627    
628     /**
629     Calculate partial derivative of u with respect to rho, with T constant
630     */
631     double helmholtz_dudrho_T(double T, double rho, const HelmholtzData *data){
632 jpye 1996 DEFINE_TD;
633 jpye 1920
634     double phir_deltau = helm_resid_deltau(tau,delta,data);
635    
636 jpye 1921 return data->R * T / rho * (tau * delta * phir_deltau);
637 jpye 1920 }
638    
639 jpye 1824 /*---------------------------------------------
640 jpye 1826 UTILITY FUNCTION(S)
641     */
642    
643     /* ipow: public domain by Mark Stephen with suggestions by Keiichi Nakasato */
644     static double ipow(double x, int n){
645     double t = 1.0;
646    
647 jpye 1829 if(!n)return 1.0; /* At the top. x^0 = 1 */
648 jpye 1826
649     if (n < 0){
650     n = -n;
651     x = 1.0/x; /* error if x == 0. Good */
652     } /* ZTC/SC returns inf, which is even better */
653    
654     if (x == 0.0)return 0.0;
655    
656     do{
657     if(n & 1)t *= x;
658     n /= 2; /* KN prefers if (n/=2) x*=x; This avoids an */
659     x *= x; /* unnecessary but benign multiplication on */
660     }while(n); /* the last pass, but the comparison is always
661     true _except_ on the last pass. */
662    
663     return t;
664     }
665    
666 jpye 1889 //#define RESID_DEBUG
667 jpye 1872
668 jpye 2232 /* maxima expressions:
669     Psi(delta) := exp(-C*(delta-1)^2 -D*(tau-1)^2);
670     theta(delta) := (1-tau) + A*((delta-1)^2)^(1/(2*beta));
671     Delta(delta):= theta(delta)^2 + B*((delta-1)^2)^a;
672     n*Delta(delta)^b*delta*Psi(delta);
673     diff(%,delta,3);
674     yikes, that's scary! break down into steps.
675     */
676    
677 jpye 1993 /*
678     We avoid duplication by using the following #defines for common code in
679     calculation of critical terms.
680     */
681     #define DEFINE_DELTA \
682     double d1 = delta - 1.; \
683     double t1 = tau - 1.; \
684     double d12 = SQ(d1); \
685     double theta = (1. - tau) + ct->A * pow(d12, 0.5/ct->beta); \
686 jpye 2232 double PSI = exp(-ct->C*d12 - ct->D*SQ(t1)); \
687 jpye 2262 double DELTA = SQ(theta) + ct->B* pow(d12, ct->a)
688 jpye 1993
689 jpye 2232 #define DEFINE_DELB \
690     double DELB = pow(DELTA,ct->b)
691    
692 jpye 1993 #define DEFINE_DPSIDDELTA \
693 jpye 2232 double dPSIddelta = -2. * ct->C * d1 * PSI
694 jpye 1993
695     #define DEFINE_DDELDDELTA \
696 jpye 2262 double dDELddelta = d1 * (ct->A * theta * 2./ct->beta * pow(d12, 0.5/ct->beta - 1) + 2* ct->B * ct->a * pow(d12, ct->a - 1))
697 jpye 1993
698     #define DEFINE_DDELBDTAU \
699 jpye 2268 double dDELbdtau = (DELTA == 0) ? 0 : -2. * theta * ct->b * (DELB/DELTA);\
700     assert(!__isnan(dDELbdtau))
701 jpye 1993
702     #define DEFINE_DPSIDTAU \
703 jpye 2232 double dPSIdtau = -2. * ct->D * t1 * PSI
704 jpye 1993
705     #define DEFINE_DDELBDDELTA \
706 jpye 2262 double dDELbddelta = (DELTA==0?0:ct->b * (DELB/DELTA) * dDELddelta)
707 jpye 1993
708 jpye 2232 #define DEFINE_D2DELDDELTA2 \
709     double powd12bm1 = pow(d12,0.5/ct->beta-1.); \
710 jpye 2262 double d2DELddelta2 = 1./d1*dDELddelta + d12*( \
711 jpye 2232 4.*ct->B*ct->a*(ct->a-1.)*pow(d12,ct->a-2.) \
712     + 2.*SQ(ct->A)*SQ(1./ct->beta)*SQ(powd12bm1) \
713     + ct->A*theta*4./ct->beta*(0.5/ct->beta-1.)*powd12bm1/d12 \
714 jpye 2262 )
715 jpye 2232
716     #define DEFINE_D2DELBDDELTA2 \
717 jpye 2262 double d2DELbddelta2 = ct->b * ( (DELB/DELTA)*d2DELddelta2 + (ct->b-1.)*(DELB/SQ(DELTA)*SQ(dDELddelta)))
718 jpye 2232
719     #define DEFINE_D2PSIDDELTA2 \
720     double d2PSIddelta2 = (2.*ct->C*d12 - 1.)*2.*ct->C * PSI
721    
722     #define DEFINE_D3PSIDDELTA3 \
723     double d3PSIddelta3 = -4. * d1 * SQ(ct->C) * (2.*d12*ct->C - 3.) * PSI
724    
725     #define DEFINE_D3DELDDELTA3 \
726 jpye 2233 double d3DELddelta3 = 1./(d1*d12*ct->beta*SQ(ct->beta)) * (\
727     4*ct->B*ct->a*(1.+ct->a*(2*ct->a-3))*SQ(ct->beta)*pow(d12,ct->a)\
728     + ct->A * (1.+ct->beta*(2*ct->beta-3))*pow(d12,0.5/ct->beta)\
729     )
730 jpye 2232
731     #define DEFINE_D3DELBDDELTA3 \
732 jpye 2262 double d3DELbddelta3 = ct->b / (DELTA*SQ(DELTA)) * ( \
733 jpye 2233 (2+ct->b*(ct->b-3))*dDELddelta*SQ(dDELddelta)*DELB \
734 jpye 2232 + DELB*SQ(DELTA)*d3DELddelta3 \
735 jpye 2233 + 3*(ct->b-1) * DELB * DELTA * dDELddelta * d2DELddelta2 \
736 jpye 2232 )
737    
738 jpye 1824 /**
739 jpye 1845 Residual part of helmholtz function.
740 jpye 1822 */
741 jpye 1826 double helm_resid(double tau, double delta, const HelmholtzData *data){
742 jpye 1867 double dell,ldell, term, sum, res = 0;
743 jpye 1845 unsigned n, i;
744     const HelmholtzPowTerm *pt;
745 jpye 1887 const HelmholtzGausTerm *gt;
746 jpye 1985 const HelmholtzCritTerm *ct;
747 jpye 1845
748     n = data->np;
749     pt = &(data->pt[0]);
750    
751 jpye 1868 #ifdef RESID_DEBUG
752 jpye 1867 fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
753     #endif
754    
755 jpye 1845 /* power terms */
756     sum = 0;
757 jpye 1865 dell = ipow(delta,pt->l);
758     ldell = pt->l * dell;
759 jpye 1845 unsigned oldl;
760     for(i=0; i<n; ++i){
761 jpye 1867 term = pt->a * pow(tau, pt->t) * ipow(delta, pt->d);
762     sum += term;
763 jpye 1868 #ifdef RESID_DEBUG
764 jpye 1867 fprintf(stderr,"i = %d, a=%e, t=%f, d=%d, term = %f, sum = %f",i,pt->a,pt->t,pt->d,term,sum);
765     if(pt->l==0){
766     fprintf(stderr,",row=%e\n",term);
767     }else{
768     fprintf(stderr,",row=%e\n,",term*exp(-dell));
769     }
770     #endif
771 jpye 1845 oldl = pt->l;
772     ++pt;
773     if(i+1==n || oldl != pt->l){
774     if(oldl == 0){
775 jpye 1868 #ifdef RESID_DEBUG
776 jpye 1865 fprintf(stderr,"linear ");
777 jpye 1867 #endif
778 jpye 1845 res += sum;
779     }else{
780 jpye 1868 #ifdef RESID_DEBUG
781 jpye 1865 fprintf(stderr,"exp dell=%f, exp(-dell)=%f sum=%f: ",dell,exp(-dell),sum);
782 jpye 1867 #endif
783 jpye 1865 res += sum * exp(-dell);
784 jpye 1845 }
785 jpye 1868 #ifdef RESID_DEBUG
786 jpye 1865 fprintf(stderr,"i = %d, res = %f\n",i,res);
787 jpye 1867 #endif
788 jpye 1845 sum = 0;
789 jpye 1865 dell = ipow(delta,pt->l);
790     ldell = pt->l*dell;
791 jpye 1845 }
792     }
793 jpye 2270 assert(!__isnan(res));
794 jpye 1845
795 jpye 1887 /* gaussian terms */
796     n = data->ng;
797 jpye 1889 //fprintf(stderr,"THERE ARE %d GAUSSIAN TERMS\n",n);
798 jpye 1887 gt = &(data->gt[0]);
799     for(i=0; i<n; ++i){
800 jpye 1868 #ifdef RESID_DEBUG
801 jpye 1887 fprintf(stderr,"i = %d, GAUSSIAN, n = %e, t = %f, d = %f, alpha = %f, beta = %f, gamma = %f, epsilon = %f\n",i+1, gt->n, gt->t, gt->d, gt->alpha, gt->beta, gt->gamma, gt->epsilon);
802     #endif
803     double d1 = delta - gt->epsilon;
804     double t1 = tau - gt->gamma;
805 jpye 1993 double e1 = -gt->alpha*SQ(d1) - gt->beta*SQ(t1);
806 jpye 1887 sum = gt->n * pow(tau,gt->t) * pow(delta,gt->d) * exp(e1);
807     //fprintf(stderr,"sum = %f\n",sum);
808     res += sum;
809     ++gt;
810     }
811 jpye 2270 assert(!__isnan(res));
812 jpye 1887
813 jpye 1985 /* critical terms */
814     n = data->nc;
815     ct = &(data->ct[0]);
816     for(i=0; i<n; ++i){
817 jpye 1887 #ifdef RESID_DEBUG
818 jpye 1985 fprintf(stderr,"i = %d, CRITICAL, n = %e, a = %f, b = %f, beta = %f, A = %f, B = %f, C = %f, D = %f\n",i+1, ct->n, ct->a, ct->b, ct->beta, ct->A, ct->B, ct->C, ct->D);
819 jpye 1864 #endif
820 jpye 1993
821     DEFINE_DELTA;
822 jpye 2232 DEFINE_DELB;
823 jpye 1993
824 jpye 2232 sum = ct->n * DELB * delta * PSI;
825 jpye 1985 res += sum;
826     ++ct;
827     }
828 jpye 2270 assert(!__isnan(res));
829 jpye 1985
830     #ifdef RESID_DEBUG
831     fprintf(stderr,"CALCULATED RESULT FOR phir = %f\n",res);
832     #endif
833 jpye 1845 return res;
834     }
835    
836 jpye 1915 /*=================== FIRST DERIVATIVES =======================*/
837    
838 jpye 1822 /**
839     Derivative of the helmholtz residual function with respect to
840     delta.
841     */
842 jpye 1826 double helm_resid_del(double tau,double delta, const HelmholtzData *data){
843 jpye 1905 double sum = 0, res = 0;
844 jpye 1844 double dell, ldell;
845 jpye 1838 unsigned n, i;
846     const HelmholtzPowTerm *pt;
847 jpye 1887 const HelmholtzGausTerm *gt;
848 jpye 1986 const HelmholtzCritTerm *ct;
849 jpye 1810
850 jpye 1898 #ifdef RESID_DEBUG
851     fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
852     #endif
853    
854 jpye 1905 /* power terms */
855     n = data->np;
856     pt = &(data->pt[0]);
857 jpye 1844 dell = ipow(delta,pt->l);
858     ldell = pt->l * dell;
859     unsigned oldl;
860 jpye 1841 for(i=0; i<n; ++i){
861 jpye 2262 sum += pt->a * pow(tau, pt->t) * ipow(delta, pt->d - 1) * (pt->d - ldell);
862 jpye 1844 oldl = pt->l;
863 jpye 1841 ++pt;
864 jpye 1844 if(i+1==n || oldl != pt->l){
865     if(oldl == 0){
866 jpye 1836 res += sum;
867     }else{
868 jpye 1844 res += sum * exp(-dell);
869 jpye 1836 }
870 jpye 1844 sum = 0;
871     dell = ipow(delta,pt->l);
872     ldell = pt->l*dell;
873 jpye 1836 }
874 jpye 1822 }
875    
876 jpye 1887 /* gaussian terms */
877     n = data->ng;
878     gt = &(data->gt[0]);
879     for(i=0; i<n; ++i){
880     #ifdef RESID_DEBUG
881 jpye 1891 fprintf(stderr,"i = %d, GAUSSIAN, n = %e, t = %f, d = %f, alpha = %f, beta = %f, gamma = %f, epsilon = %f\n",i+1, gt->n, gt->t, gt->d, gt->alpha, gt->beta, gt->gamma, gt->epsilon);
882 jpye 1840 #endif
883 jpye 1986 sum = - gt->n * pow(tau,gt->t) * pow(delta, -1. + gt->d)
884 jpye 1898 * (2. * gt->alpha * delta * (delta - gt->epsilon) - gt->d)
885     * exp(-(gt->alpha * SQ(delta-gt->epsilon) + gt->beta*SQ(tau-gt->gamma)));
886 jpye 1986 res += sum;
887 jpye 1891 #ifdef RESID_DEBUG
888 jpye 1986 fprintf(stderr,"sum = %f --> res = %f\n",sum,res);
889 jpye 1891 #endif
890 jpye 1887 ++gt;
891     }
892 jpye 1838
893 jpye 1986 /* critical terms */
894     n = data->nc;
895     ct = &(data->ct[0]);
896     for(i=0; i<n; ++i){
897     #ifdef RESID_DEBUG
898     fprintf(stderr,"i = %d, CRITICAL, n = %e, a = %f, b = %f, beta = %f, A = %f, B = %f, C = %f, D = %f\n",i+1, ct->n, ct->a, ct->b, ct->beta, ct->A, ct->B, ct->C, ct->D);
899     #endif
900 jpye 1993 DEFINE_DELTA;
901 jpye 2232 DEFINE_DELB;
902 jpye 1993 DEFINE_DPSIDDELTA;
903     DEFINE_DDELDDELTA;
904     DEFINE_DDELBDDELTA;
905 jpye 1985
906 jpye 1996 #if 0
907     if(fabs(dpsiddelta) ==0)fprintf(stderr,"WARNING: dpsiddelta == 0\n");
908     if(fabs(dpsiddelta) ==0)fprintf(stderr,"WARNING: dpsiddelta == 0\n");
909     fprintf(stderr,"psi = %f\n",psi);
910     fprintf(stderr,"DELTA = %f\n",DELTA);
911    
912     fprintf(stderr,"dDELddelta = %f\n",dDELddelta);
913     fprintf(stderr,"ct->b - 1. = %f\n",ct->b - 1.);
914     fprintf(stderr,"pow(DELTA,ct->b - 1.) = %f\n",pow(DELTA,ct->b - 1.));
915     assert(!isnan(pow(DELTA,ct->b - 1.)));
916     assert(!isnan(dDELddelta));
917     assert(!isnan(dDELbddelta));
918     //double dDELbddelta = ct->b * pow(DELTA,ct->b - 1.) * dDELddelta
919     fprintf(stderr,"sum = %f\n",sum);
920     if(isnan(sum))fprintf(stderr,"ERROR: sum isnan with i=%d at %d\n",i,__LINE__);
921     #endif
922 jpye 2232 sum = ct->n * (DELB * (PSI + delta * dPSIddelta) + dDELbddelta * delta * PSI);
923 jpye 1986 res += sum;
924 jpye 1996
925 jpye 1986 ++ct;
926     }
927    
928 jpye 1836 return res;
929 jpye 1822 }
930    
931 jpye 1824 /**
932     Derivative of the helmholtz residual function with respect to
933     tau.
934     */
935 jpye 1826 double helm_resid_tau(double tau,double delta,const HelmholtzData *data){
936 jpye 1824
937     double sum;
938 jpye 1827 double res = 0;
939 jpye 1832 double delX;
940     unsigned l;
941 jpye 1845 unsigned n, i;
942 jpye 1838 const HelmholtzPowTerm *pt;
943 jpye 1891 const HelmholtzGausTerm *gt;
944 jpye 1987 const HelmholtzCritTerm *ct;
945 jpye 1822
946 jpye 1845 n = data->np;
947 jpye 1838 pt = &(data->pt[0]);
948 jpye 1832
949     delX = 1;
950    
951     l = 0;
952     sum = 0;
953 jpye 1845 for(i=0; i<n; ++i){
954 jpye 1838 if(pt->t){
955     //fprintf(stderr,"i = %d, a = %e, t = %f, d = %d, l = %d\n",i+1, pt->a, pt->t, pt->d, pt->l);
956     sum += pt->a * pow(tau, pt->t - 1) * ipow(delta, pt->d) * pt->t;
957 jpye 1832 }
958 jpye 1838 ++pt;
959 jpye 1832 //fprintf(stderr,"l = %d\n",l);
960 jpye 1845 if(i+1==n || l != pt->l){
961 jpye 1832 if(l==0){
962     //fprintf(stderr,"Adding non-exp term\n");
963     res += sum;
964     }else{
965     //fprintf(stderr,"Adding exp term with l = %d, delX = %e\n",l,delX);
966     res += sum * exp(-delX);
967     }
968     /* set l to new value */
969 jpye 1845 if(i+1!=n){
970 jpye 1838 l = pt->l;
971 jpye 1832 //fprintf(stderr,"New l = %d\n",l);
972     delX = ipow(delta,l);
973     sum = 0;
974     }
975     }
976     }
977 jpye 2268 assert(!__isnan(res));
978 jpye 1832
979 jpye 1897 //#define RESID_DEBUG
980 jpye 1891 /* gaussian terms */
981     n = data->ng;
982     gt = &(data->gt[0]);
983     for(i=0; i<n; ++i){
984     #ifdef RESID_DEBUG
985     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);
986     #endif
987 jpye 1897
988     double val2;
989 jpye 1900 val2 = -gt->n * pow(tau,gt->t - 1.) * pow(delta, gt->d)
990     * (2. * gt->beta * tau * (tau - gt->gamma) - gt->t)
991 jpye 1898 * exp(-(gt->alpha * SQ(delta-gt->epsilon) + gt->beta*SQ(tau-gt->gamma)));
992     res += val2;
993 jpye 1891 #ifdef RESID_DEBUG
994 jpye 1897 fprintf(stderr,"res = %f\n",res);
995 jpye 1891 #endif
996 jpye 1897
997 jpye 1891 ++gt;
998     }
999 jpye 2268 assert(!__isnan(res));
1000 jpye 1891
1001 jpye 1987 /* critical terms */
1002     n = data->nc;
1003     ct = &(data->ct[0]);
1004     for(i=0; i<n; ++i){
1005     #ifdef RESID_DEBUG
1006     fprintf(stderr,"i = %d, CRITICAL, n = %e, a = %f, b = %f, beta = %f, A = %f, B = %f, C = %f, D = %f\n",i+1, ct->n, ct->a, ct->b, ct->beta, ct->A, ct->B, ct->C, ct->D);
1007     #endif
1008 jpye 1993 DEFINE_DELTA;
1009 jpye 2232 DEFINE_DELB;
1010 jpye 1993 DEFINE_DDELBDTAU;
1011     DEFINE_DPSIDTAU;
1012 jpye 1987
1013 jpye 2232 sum = ct->n * delta * (dDELbdtau * PSI + DELB * dPSIdtau);
1014 jpye 1987 res += sum;
1015     ++ct;
1016     }
1017 jpye 2268 assert(!__isnan(res));
1018 jpye 1987
1019 jpye 1827 return res;
1020 jpye 1824 }
1021    
1022    
1023 jpye 1915 /*=================== SECOND DERIVATIVES =======================*/
1024 jpye 1824
1025 jpye 1826 /**
1026     Mixed derivative of the helmholtz residual function with respect to
1027 jpye 1907 delta and tau.
1028 jpye 1868 */
1029 jpye 1826 double helm_resid_deltau(double tau,double delta,const HelmholtzData *data){
1030 jpye 2265 double dell,ldell, sum = 0, res = 0;
1031 jpye 1905 unsigned n, i;
1032     const HelmholtzPowTerm *pt;
1033     const HelmholtzGausTerm *gt;
1034 jpye 1990 const HelmholtzCritTerm *ct;
1035 jpye 1824
1036 jpye 1905 /* power terms */
1037     n = data->np;
1038     pt = &(data->pt[0]);
1039     dell = ipow(delta,pt->l);
1040     ldell = pt->l * dell;
1041     unsigned oldl;
1042     for(i=0; i<n; ++i){
1043     sum += pt->a * pt->t * pow(tau, pt->t - 1) * ipow(delta, pt->d - 1) * (pt->d - ldell);
1044     oldl = pt->l;
1045 jpye 1838 ++pt;
1046 jpye 1905 if(i+1==n || oldl != pt->l){
1047     if(oldl == 0){
1048     res += sum;
1049     }else{
1050     res += sum * exp(-dell);
1051     }
1052     sum = 0;
1053     dell = ipow(delta,pt->l);
1054     ldell = pt->l*dell;
1055     }
1056 jpye 1826 }
1057    
1058 jpye 1909 #ifdef TEST
1059 jpye 1906 assert(!isinf(res));
1060 jpye 1909 #endif
1061 jpye 1906
1062 jpye 1905 /* gaussian terms */
1063     n = data->ng;
1064     gt = &(data->gt[0]);
1065     for(i=0; i<n; ++i){
1066     #ifdef RESID_DEBUG
1067     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);
1068     #endif
1069     double d1 = delta - gt->epsilon;
1070     double t1 = tau - gt->gamma;
1071     double e1 = -gt->alpha*SQ(d1) - gt->beta*SQ(t1);
1072 jpye 1826
1073 jpye 1905 double f1 = gt->t - 2*gt->beta*tau*(tau - gt->gamma);
1074     double g1 = gt->d - 2*gt->alpha*delta*(delta - gt->epsilon);
1075 jpye 1826
1076 jpye 1905 sum = gt->n * f1 * pow(tau,gt->t-1) * g1 * pow(delta,gt->d-1) * exp(e1);
1077 jpye 1906
1078 jpye 1905 //fprintf(stderr,"sum = %f\n",sum);
1079     res += sum;
1080 jpye 1909 #ifdef TEST
1081 jpye 1906 assert(!isinf(res));
1082 jpye 1909 #endif
1083 jpye 1905 ++gt;
1084 jpye 1826 }
1085    
1086 jpye 1990 /* critical terms */
1087     n = data->nc;
1088     ct = &(data->ct[0]);
1089     for(i=0; i<n; ++i){
1090     #ifdef RESID_DEBUG
1091     fprintf(stderr,"i = %d, CRITICAL, n = %e, a = %f, b = %f, beta = %f, A = %f, B = %f, C = %f, D = %f\n",i+1, ct->n, ct->a, ct->b, ct->beta, ct->A, ct->B, ct->C, ct->D);
1092     #endif
1093 jpye 1993 DEFINE_DELTA;
1094 jpye 2232 DEFINE_DELB;
1095 jpye 1993 DEFINE_DPSIDDELTA;
1096     DEFINE_DDELBDTAU;
1097     DEFINE_DDELDDELTA;
1098 jpye 1985
1099 jpye 2262 double d2DELbddeldtau = -ct->A * ct->b * 2./ct->beta * (DELB/DELTA)*d1*pow(d12,0.5/ct->beta-1) \
1100     - 2. * theta * ct->b * (ct->b - 1) * (DELB/SQ(DELTA)) * dDELddelta;
1101 jpye 1990
1102 jpye 2232 double d2PSIddeldtau = 4. * ct->C*ct->D*d1*t1*PSI;
1103 jpye 1990
1104 jpye 1993 DEFINE_DPSIDTAU;
1105 jpye 1990
1106 jpye 2232 sum = ct->n * (DELB * (dPSIdtau + delta * d2PSIddeldtau) \
1107     + delta *dDELbdtau*dPSIdtau \
1108     + dDELbdtau*(PSI+delta*dPSIddelta) \
1109     + d2DELbddeldtau*delta*PSI
1110 jpye 1990 );
1111     res += sum;
1112     ++ct;
1113     }
1114    
1115 jpye 1905 #ifdef RESID_DEBUG
1116     fprintf(stderr,"phir = %f\n",res);
1117     #endif
1118 jpye 1906
1119 jpye 1909 #ifdef TEST
1120 jpye 1906 assert(!isnan(res));
1121     assert(!isinf(res));
1122 jpye 1909 #endif
1123 jpye 1905 return res;
1124 jpye 1826 }
1125    
1126     /**
1127     Second derivative of helmholtz residual function with respect to
1128     delta (twice).
1129     */
1130     double helm_resid_deldel(double tau,double delta,const HelmholtzData *data){
1131 jpye 1907 double sum = 0, res = 0;
1132     double dell, ldell;
1133     unsigned n, i;
1134     const HelmholtzPowTerm *pt;
1135     const HelmholtzGausTerm *gt;
1136 jpye 1990 const HelmholtzCritTerm *ct;
1137 jpye 2271
1138 jpye 1907 #ifdef RESID_DEBUG
1139     fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
1140     #endif
1141 jpye 1826
1142 jpye 1907 /* power terms */
1143     n = data->np;
1144     pt = &(data->pt[0]);
1145     dell = ipow(delta,pt->l);
1146     ldell = pt->l * dell;
1147     unsigned oldl;
1148     for(i=0; i<n; ++i){
1149     double lpart = pt->l ? SQ(ldell) + ldell*(1. - 2*pt->d - pt->l) : 0;
1150     sum += pt->a * pow(tau, pt->t) * ipow(delta, pt->d - 2) * (pt->d*(pt->d - 1) + lpart);
1151     oldl = pt->l;
1152 jpye 1838 ++pt;
1153 jpye 1907 if(i+1==n || oldl != pt->l){
1154     if(oldl == 0){
1155     res += sum;
1156     }else{
1157     res += sum * exp(-dell);
1158     }
1159     sum = 0;
1160     dell = ipow(delta,pt->l);
1161     ldell = pt->l*dell;
1162     }
1163 jpye 1826 }
1164 jpye 2270 if(__isnan(res)){
1165     fprintf(stderr,"tau = %.12e, del = %.12e\n",tau,delta);
1166     }
1167     assert(!__isnan(res));
1168 jpye 1826
1169 jpye 1907 /* gaussian terms */
1170     n = data->ng;
1171     //fprintf(stderr,"THERE ARE %d GAUSSIAN TERMS\n",n);
1172     gt = &(data->gt[0]);
1173     for(i=0; i<n; ++i){
1174     double s1 = SQ(delta - gt->epsilon);
1175 jpye 2262 double f1 = gt->d*(gt->d - 1)
1176     + 2.*gt->alpha*delta * (
1177     delta * (2. * gt->alpha * s1 - 1)
1178     - 2. * gt->d * (delta - gt->epsilon)
1179     );
1180 jpye 1908 res += gt->n * pow(tau,gt->t) * pow(delta, gt->d - 2.)
1181 jpye 1907 * f1
1182     * exp(-(gt->alpha * s1 + gt->beta*SQ(tau-gt->gamma)));
1183     ++gt;
1184 jpye 1826 }
1185 jpye 2270 assert(!__isnan(res));
1186 jpye 1826
1187 jpye 1990 /* critical terms */
1188     n = data->nc;
1189     ct = &(data->ct[0]);
1190     for(i=0; i<n; ++i){
1191     #ifdef RESID_DEBUG
1192     fprintf(stderr,"i = %d, CRITICAL, n = %e, a = %f, b = %f, beta = %f, A = %f, B = %f, C = %f, D = %f\n",i+1, ct->n, ct->a, ct->b, ct->beta, ct->A, ct->B, ct->C, ct->D);
1193     #endif
1194 jpye 1985
1195 jpye 1993 DEFINE_DELTA;
1196 jpye 2232 DEFINE_DELB;
1197 jpye 1993 DEFINE_DPSIDDELTA;
1198     DEFINE_DDELDDELTA;
1199     DEFINE_DDELBDDELTA;
1200 jpye 1990
1201 jpye 2232 DEFINE_D2DELDDELTA2;
1202     DEFINE_D2DELBDDELTA2;
1203 jpye 1990
1204 jpye 2232 DEFINE_D2PSIDDELTA2;
1205 jpye 1990
1206 jpye 2232 sum = ct->n * (DELB*(2.*dPSIddelta + delta*d2PSIddelta2) + 2.*dDELbddelta*(PSI+delta*dPSIddelta) + d2DELbddelta2*delta*PSI);
1207 jpye 1990
1208     res += sum;
1209     ++ct;
1210     }
1211 jpye 2270 assert(!__isnan(res));
1212 jpye 1990
1213 jpye 1907 return res;
1214 jpye 1826 }
1215    
1216 jpye 1915
1217    
1218     /**
1219     Residual part of helmholtz function.
1220     */
1221     double helm_resid_tautau(double tau, double delta, const HelmholtzData *data){
1222     double dell,ldell, term, sum, res = 0;
1223     unsigned n, i;
1224     const HelmholtzPowTerm *pt;
1225     const HelmholtzGausTerm *gt;
1226 jpye 1988 const HelmholtzCritTerm *ct;
1227 jpye 1915
1228     n = data->np;
1229     pt = &(data->pt[0]);
1230    
1231     #ifdef RESID_DEBUG
1232     fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
1233     #endif
1234    
1235     /* power terms */
1236     sum = 0;
1237     dell = ipow(delta,pt->l);
1238     ldell = pt->l * dell;
1239     unsigned oldl;
1240     for(i=0; i<n; ++i){
1241 jpye 1918 term = pt->a * pt->t * (pt->t - 1) * pow(tau, pt->t - 2) * ipow(delta, pt->d);
1242 jpye 1915 sum += term;
1243     #ifdef RESID_DEBUG
1244     fprintf(stderr,"i = %d, a=%e, t=%f, d=%d, term = %f, sum = %f",i,pt->a,pt->t,pt->d,term,sum);
1245     if(pt->l==0){
1246     fprintf(stderr,",row=%e\n",term);
1247     }else{
1248     fprintf(stderr,",row=%e\n,",term*exp(-dell));
1249     }
1250     #endif
1251     oldl = pt->l;
1252     ++pt;
1253     if(i+1==n || oldl != pt->l){
1254     if(oldl == 0){
1255     #ifdef RESID_DEBUG
1256     fprintf(stderr,"linear ");
1257     #endif
1258     res += sum;
1259     }else{
1260     #ifdef RESID_DEBUG
1261     fprintf(stderr,"exp dell=%f, exp(-dell)=%f sum=%f: ",dell,exp(-dell),sum);
1262     #endif
1263     res += sum * exp(-dell);
1264     }
1265     #ifdef RESID_DEBUG
1266     fprintf(stderr,"i = %d, res = %f\n",i,res);
1267     #endif
1268     sum = 0;
1269     dell = ipow(delta,pt->l);
1270     ldell = pt->l*dell;
1271     }
1272     }
1273    
1274     /* gaussian terms */
1275     n = data->ng;
1276     //fprintf(stderr,"THERE ARE %d GAUSSIAN TERMS\n",n);
1277     gt = &(data->gt[0]);
1278     for(i=0; i<n; ++i){
1279     #ifdef RESID_DEBUG
1280     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);
1281     #endif
1282     double d1 = delta - gt->epsilon;
1283     double t1 = tau - gt->gamma;
1284     double f1 = gt->t*(gt->t - 1) + 4. * gt->beta * tau * (tau * (gt->beta*SQ(t1) - 0.5) - t1*gt->t);
1285     double e1 = -gt->alpha*SQ(d1) - gt->beta*SQ(t1);
1286 jpye 1920 sum = gt->n * f1 * pow(tau,gt->t - 2) * pow(delta,gt->d) * exp(e1);
1287 jpye 1915 //fprintf(stderr,"sum = %f\n",sum);
1288     res += sum;
1289     ++gt;
1290     }
1291    
1292 jpye 1988 /* critical terms */
1293     n = data->nc;
1294     ct = &(data->ct[0]);
1295     for(i=0; i<n; ++i){
1296     #ifdef RESID_DEBUG
1297     fprintf(stderr,"i = %d, CRITICAL, n = %e, a = %f, b = %f, beta = %f, A = %f, B = %f, C = %f, D = %f\n",i+1, ct->n, ct->a, ct->b, ct->beta, ct->A, ct->B, ct->C, ct->D);
1298     #endif
1299 jpye 1993 DEFINE_DELTA;
1300 jpye 2232 DEFINE_DELB;
1301 jpye 1993 DEFINE_DDELBDTAU;
1302     DEFINE_DPSIDTAU;
1303 jpye 1988
1304 jpye 2262 double d2DELbdtau2 = 2. * ct->b * (DELB/DELTA) + 4. * SQ(theta) * ct->b * (ct->b - 1) * (DELB/SQ(DELTA));
1305 jpye 1988
1306 jpye 2232 double d2PSIdtau2 = 2. * ct->D * PSI * (2. * ct->D * SQ(t1) -1.);
1307 jpye 1993
1308 jpye 2232 sum = ct->n * delta * (d2DELbdtau2 * PSI + 2 * dDELbdtau*dPSIdtau + DELB * d2PSIdtau2);
1309 jpye 1988 res += sum;
1310     ++ct;
1311     }
1312    
1313 jpye 1915 #ifdef RESID_DEBUG
1314     fprintf(stderr,"phir_tautau = %f\n",res);
1315     #endif
1316     return res;
1317     }
1318    
1319 jpye 2231 /* === THIRD DERIVATIVES (this is getting boring now) === */
1320    
1321 jpye 2275 #ifdef INCLUDE_THIRD_DERIV_CODE
1322 jpye 2231 /**
1323     Third derivative of helmholtz residual function, with respect to
1324     delta (thrice).
1325     */
1326 jpye 2232 double helm_resid_deldeldel(double tau,double delta,const HelmholtzData *data){
1327 jpye 2231 double sum = 0, res = 0;
1328 jpye 2235 double D,L;
1329 jpye 2231 unsigned n, i;
1330     const HelmholtzPowTerm *pt;
1331     const HelmholtzGausTerm *gt;
1332     const HelmholtzCritTerm *ct;
1333    
1334     #ifdef RESID_DEBUG
1335     fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
1336     #endif
1337    
1338 jpye 2235 #if 1
1339     /* major shortcut, but not efficient */
1340     double ddel = 0.0000000001;
1341     return (helm_resid_deldel(tau,delta+ddel,data) - helm_resid_deldel(tau,delta,data))/ddel;
1342     #endif
1343    
1344     /* seem to be errors in the following, still haven't tracked them all down. */
1345    
1346     #if 0
1347     /* wxmaxima code:
1348     a*delta^d*%e^(-delta^l)*tau^t
1349     diff(%,delta,3);
1350     */
1351 jpye 2231 /* power terms */
1352     n = data->np;
1353     pt = &(data->pt[0]);
1354 jpye 2235 D = ipow(delta,pt->l);
1355 jpye 2231 unsigned oldl;
1356     for(i=0; i<n; ++i){
1357 jpye 2235 double d = pt->d;
1358     double l = pt->l;
1359 jpye 2232 double lpart = pt->l
1360 jpye 2235 ? D*((D-1)*(D-2)-1) * l*SQ(l)
1361     + 3*D*(1-d)*(D-1) * SQ(l)
1362     + D*(3*SQ(d-1)-1) * l
1363 jpye 2232 : 0;
1364 jpye 2235 sum += pt->a * pow(tau,pt->t) * ipow(delta, d-3) * (d*(d-1)*(d-2) + lpart);
1365 jpye 2231 oldl = pt->l;
1366     ++pt;
1367     if(i+1==n || oldl != pt->l){
1368     if(oldl == 0){
1369 jpye 2232 res += sum; // note special meaning of l==0 case: no exponential
1370 jpye 2231 }else{
1371 jpye 2235 res += sum * exp(-D);
1372 jpye 2231 }
1373     sum = 0;
1374 jpye 2235 D = ipow(delta,pt->l);
1375 jpye 2231 }
1376     }
1377    
1378 jpye 2235 //fprintf(stderr,"DELDELDEL fiff = %f, sum = %f ",fdiff, res);
1379     #endif
1380    
1381     #if 0
1382 jpye 2232 /* gaussian terms */
1383     n = data->ng;
1384     //fprintf(stderr,"THERE ARE %d GAUSSIAN TERMS\n",n);
1385     gt = &(data->gt[0]);
1386     for(i=0; i<n; ++i){
1387     double D = delta - gt->epsilon;
1388     double D2 = SQ(D);
1389     double T2 = SQ(tau - gt->gamma);
1390     double A = gt->alpha * delta;
1391     double A2 = SQ(A);
1392     double d = gt->d;
1393     double d2 = SQ(d);
1394 jpye 2231
1395 jpye 2232 // this expression calculated from wxMaxima using subsitutions for
1396     // D=delta-epsilon and A=alpha*delta.
1397     double f1 =
1398     - (8*A*A2) * D*D2
1399     + (12*d * A2) * D2
1400     + (12 * delta * A2 + (6*d - 6*d2)*A) * D
1401     - (6 * d * delta * A + d*d2 - 3*d2 + 2*d);
1402 jpye 2231
1403 jpye 2232 res += gt->n * pow(tau,gt->t) * pow(delta, d - 3.)
1404     * f1
1405     * exp(-(gt->alpha * D2 + gt->beta * T2));
1406     ++gt;
1407     }
1408 jpye 2235 #endif
1409 jpye 2231
1410 jpye 2235 #if 0
1411 jpye 2232 /* critical terms */
1412     n = data->nc;
1413     ct = &(data->ct[0]);
1414     for(i=0; i<n; ++i){
1415     #ifdef RESID_DEBUG
1416     fprintf(stderr,"i = %d, CRITICAL, n = %e, a = %f, b = %f, beta = %f, A = %f, B = %f, C = %f, D = %f\n",i+1, ct->n, ct->a, ct->b, ct->beta, ct->A, ct->B, ct->C, ct->D);
1417     #endif
1418 jpye 2231
1419 jpye 2232 DEFINE_DELTA;
1420     DEFINE_DELB;
1421     DEFINE_DPSIDDELTA;
1422     DEFINE_DDELDDELTA;
1423     DEFINE_DDELBDDELTA;
1424 jpye 2231
1425 jpye 2232 DEFINE_D2PSIDDELTA2;
1426     DEFINE_D2DELDDELTA2;
1427     DEFINE_D2DELBDDELTA2;
1428    
1429     DEFINE_D3PSIDDELTA3;
1430     DEFINE_D3DELDDELTA3;
1431     DEFINE_D3DELBDDELTA3;
1432    
1433     sum = ct->n * (
1434     delta * (DELB*d3PSIddelta3 + 3 * dDELbddelta * d2PSIddelta2 + 3 * d2DELbddelta2 * dPSIddelta + PSI * d3DELbddelta3)
1435     + 3 * (DELB*d2PSIddelta2 + 2 * dDELbddelta * dPSIddelta + d2DELbddelta2 * PSI)
1436     );
1437    
1438     res += sum;
1439     ++ct;
1440     }
1441 jpye 2231 #endif
1442    
1443 jpye 2232 return res;
1444     }
1445 jpye 2231
1446 jpye 2275 #endif
1447 jpye 2231
1448    

john.pye@anu.edu.au
ViewVC Help
Powered by ViewVC 1.1.22