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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2271 - (show annotations) (download) (as text)
Sun Aug 8 12:17:00 2010 UTC (13 years, 11 months ago) by jpye
File MIME type: text/x-csrc
File size: 38730 byte(s)
Removed p_t, p_c from HelmholtzData, added routines that calculate those values directly.
'satcvgc.py' now runs without errors for 100% of data points.
1 /* ASCEND modelling environment
2 Copyright (C) 2008-2009 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 #include "sat.h"
32
33 #if 0
34 # include <assert.h>
35 #else
36 # define assert(ARGS...)
37 #endif
38
39 #include <stdlib.h>
40 #include <stdio.h>
41
42 /* macros and forward decls */
43
44 #define SQ(X) ((X)*(X))
45
46 #include "helmholtz_impl.h"
47
48 /* 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 double helmholtz_p(double T, double rho, const HelmholtzData *d){
54 double p, rho_f, rho_g;
55 #if 0
56 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 /* but what if we're in the sublimation region?? */
61 #endif
62 if(T < d->T_c){
63 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 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 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
107 /**
108 Function to calculate pressure from Helmholtz free energy EOS, given temperature
109 and mass density.
110
111 @param T temperature in K
112 @param rho mass density in kg/m続
113 @return pressure in Pa???
114 */
115 double helmholtz_p_raw(double T, double rho, const HelmholtzData *data){
116 DEFINE_TD;
117
118 assert(data->rho_star!=0);
119 assert(T!=0);
120 assert(!__isnan(T));
121 assert(!__isnan(rho));
122 assert(!__isnan(data->R));
123
124 //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
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
133 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 }
140
141 /**
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 double helmholtz_u(double T, double rho, const HelmholtzData *data){
150 DEFINE_TD;
151
152 #ifdef TEST
153 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 #if 0
161 fprintf(stderr,"ideal_tau = %f\n",helm_ideal_tau(tau,delta,data->ideal));
162 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 #endif
165
166 return data->R * data->T_star * (helm_ideal_tau(tau,delta,data->ideal) + helm_resid_tau(tau,delta,data));
167 }
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 @return enthalpy in J/kg
176 */
177 double helmholtz_h_raw(double T, double rho, const HelmholtzData *data){
178 DEFINE_TD;
179
180 //#ifdef TEST
181 assert(data->rho_star!=0);
182 assert(T!=0);
183 assert(!isnan(tau));
184 assert(!isnan(delta));
185 assert(!isnan(data->R));
186 //#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 }
191
192 /**
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 double helmholtz_s_raw(double T, double rho, const HelmholtzData *data){
201 DEFINE_TD;
202
203 #ifdef ENTROPY_DEBUG
204 assert(data->rho_star!=0);
205 assert(T!=0);
206 assert(!isnan(tau));
207 assert(!isnan(delta));
208 assert(!isnan(data->R));
209
210 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 #endif
215 return data->R * (
216 tau * (helm_ideal_tau(tau,delta,data->ideal) + helm_resid_tau(tau,delta,data))
217 - (helm_ideal(tau,delta,data->ideal) + helm_resid(tau,delta,data))
218 );
219 }
220
221 /**
222 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 DEFINE_TD;
231
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 #ifdef HELMHOLTZ_DEBUG
241 fprintf(stderr,"helmholtz_a: T = %f, rho = %f\n",T,rho);
242 fprintf(stderr,"multiplying by RT = %f\n",data->R*T);
243 #endif
244
245 return data->R * T * (helm_ideal(tau,delta,data->ideal) + helm_resid(tau,delta,data));
246 }
247
248 /**
249 Function to calculate isochoric heat capacity from the Helmholtz free energy
250 EOS given temperature and mass density.
251
252 @param T temperature in K
253 @param rho mass density in kg/m続
254 @return Isochoric specific heat capacity in J/kg/K.
255 */
256 double helmholtz_cv(double T, double rho, const HelmholtzData *data){
257 DEFINE_TD;
258
259 return - data->R * SQ(tau) * (helm_ideal_tautau(tau,data->ideal) + helm_resid_tautau(tau,delta,data));
260 }
261
262 /**
263 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 DEFINE_TD;
272
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 double temp3 = -SQ(tau)*(helm_ideal_tautau(tau,data->ideal) + helm_resid_tautau(tau,delta,data));
281
282 return data->R * (temp3 + SQ(temp2)/temp1);
283 }
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 DEFINE_TD;
296
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 double temp1 = 1. + 2.*delta*phir_d + SQ(delta)*phir_dd;
303 double temp2 = 1. + delta*phir_d - delta*tau*phir_dt;
304 double temp3 = -SQ(tau)*(helm_ideal_tautau(tau,data->ideal) + helm_resid_tautau(tau,delta,data));
305
306 return sqrt(data->R * T * (temp1 + SQ(temp2)/temp3));
307
308 }
309
310 /**
311 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 return data->R * T * (phi0 + phir + 1. + delta * phir_d);
326 }
327
328 /**
329 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 #if 0
334 fprintf(stderr,"val = %f\n",val);
335 #endif
336 return val;
337 }
338
339 /**
340 alpha_p function from IAPWS Advisory Note 3, used in calculation of
341 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 /**
351 beta_p function from IAPWS Advisory Note 3 , used in calculation of partial
352 property derivatives.
353 */
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 return rho*(1. + (delta*phir_d + SQ(delta)*phir_dd)/(1+delta*phir_d));
359 }
360
361
362 /**
363 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 #ifdef TEST
381 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 #endif
383
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 #ifdef TEST
419 fprintf(stderr,"Converged to p = %f MPa with rho = %f\n", p/1e6, rho_1);
420 #endif
421 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 #ifdef TEST
442 fprintf(stderr,"ERROR: temperature exceeds critical temperature in helmholtz_sat_t.\n");
443 #endif
444 /* 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 #ifdef TEST
455 fprintf(stderr,"T_c = %f, rho_c = %f\n",data->T_star, data->rho_star);
456 #endif
457 double p_c = helmholtz_p(data->T_star, data->rho_star, data);
458
459 #ifdef TEST
460 fprintf(stderr,"Critical pressure = %f MPa\n",p_c/1.e6);
461 fprintf(stderr,"Acentric factor = %f\n",data->omega);
462 #endif
463
464
465 /* maybe use http://dx.doi.org/10.1016/S0009-2509(02)00017-9 for this part? */
466 /* FIXME need to cite this formula */
467 *p = p_c * pow(10.,(data->omega + 1.)*-7./3.*(tau - 1.));
468
469 #ifdef TEST
470 fprintf(stderr,"Estimated p_sat(T=%f) = %f MPa\n",T,(*p)/1e6);
471 #endif
472
473 if(tau < 1.01){
474 /* close to critical point: need a different approach */
475 #ifdef TEST
476 fprintf(stderr,"ERROR: not implemented\n");
477 #endif
478 return 1;
479 }else{
480 double niter = 0;
481 (void)niter;
482
483 int res = helm_find_rho_tp(T, *p, rho_f, data->rho_star,data->rho_star*10, data);
484 if(res){
485 #ifdef TEST
486 fprintf(stderr,"ERROR: failed to solve rho_f\n");
487 #endif
488 }
489
490 res = helm_find_rho_tp(T, *p, rho_g, 0.001*data->rho_star, data->rho_star, data);
491 if(res){
492 #ifdef TEST
493 fprintf(stderr,"ERROR: failed to solve rho_g\n");
494 #endif
495 }
496
497 #ifdef TEST
498 fprintf(stderr,"p = %f MPa: rho_f = %f, rho_g = %f\n", *p, *rho_f, *rho_g);
499 #endif
500
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 (void)err;
508
509 #ifdef TEST
510 fprintf(stderr,"LHS = %f, RHS = %f, err = %f\n",LHS, RHS, err);
511 #endif
512 /* away from critical point... */
513 *rho_f = data->rho_star;
514 *rho_g = data->rho_star;
515 #ifdef TEST
516 fprintf(stderr,"ERROR: not implemented\n");
517 exit(1);
518 #endif
519 return 1;
520 }
521 }
522
523 /*----------------------------------------------------------------------------
524 PARTIAL DERIVATIVES
525 */
526
527 /**
528 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 DEFINE_TD;
532
533 double phir_del = helm_resid_del(tau,delta,data);
534 double phir_deltau = helm_resid_deltau(tau,delta,data);
535 #ifdef TEST
536 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 #endif
544
545 double res = data->R * rho * (1 + delta*phir_del - delta*tau*phir_deltau);
546
547 #ifdef TEST
548 assert(!isnan(res));
549 assert(!isinf(res));
550 #endif
551 return res;
552 }
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 DEFINE_TD;
559
560 double phir_del = helm_resid_del(tau,delta,data);
561 double phir_deldel = helm_resid_deldel(tau,delta,data);
562 #ifdef TEST
563 assert(!isinf(phir_del));
564 assert(!isinf(phir_deldel));
565 #endif
566 return data->R * T * (1 + 2*delta*phir_del + SQ(delta)*phir_deldel);
567 }
568
569
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 return data->R * T / rho * delta * (2*phir_del + delta*(4*phir_deldel + delta*phir_deldeldel));
583 }
584
585 /**
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 DEFINE_TD;
590
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 //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 return data->R * (1. + delta*phir_del - SQ(tau)*(phi0_tautau + phir_tautau) - delta*tau*phir_deltau);
600 }
601
602 /**
603 Calculate partial derivative of h with respect to rho, with T constant
604 */
605 double helmholtz_dhdrho_T(double T, double rho, const HelmholtzData *data){
606 DEFINE_TD;
607
608 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 }
614
615
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 DEFINE_TD;
621
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 DEFINE_TD;
633
634 double phir_deltau = helm_resid_deltau(tau,delta,data);
635
636 return data->R * T / rho * (tau * delta * phir_deltau);
637 }
638
639 /*---------------------------------------------
640 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 if(!n)return 1.0; /* At the top. x^0 = 1 */
648
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 //#define RESID_DEBUG
667
668 /* 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 /*
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 double PSI = exp(-ct->C*d12 - ct->D*SQ(t1)); \
687 double DELTA = SQ(theta) + ct->B* pow(d12, ct->a)
688
689 #define DEFINE_DELB \
690 double DELB = pow(DELTA,ct->b)
691
692 #define DEFINE_DPSIDDELTA \
693 double dPSIddelta = -2. * ct->C * d1 * PSI
694
695 #define DEFINE_DDELDDELTA \
696 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
698 #define DEFINE_DDELBDTAU \
699 double dDELbdtau = (DELTA == 0) ? 0 : -2. * theta * ct->b * (DELB/DELTA);\
700 assert(!__isnan(dDELbdtau))
701
702 #define DEFINE_DPSIDTAU \
703 double dPSIdtau = -2. * ct->D * t1 * PSI
704
705 #define DEFINE_DDELBDDELTA \
706 double dDELbddelta = (DELTA==0?0:ct->b * (DELB/DELTA) * dDELddelta)
707
708 #define DEFINE_D2DELDDELTA2 \
709 double powd12bm1 = pow(d12,0.5/ct->beta-1.); \
710 double d2DELddelta2 = 1./d1*dDELddelta + d12*( \
711 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 )
715
716 #define DEFINE_D2DELBDDELTA2 \
717 double d2DELbddelta2 = ct->b * ( (DELB/DELTA)*d2DELddelta2 + (ct->b-1.)*(DELB/SQ(DELTA)*SQ(dDELddelta)))
718
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 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
731 #define DEFINE_D3DELBDDELTA3 \
732 double d3DELbddelta3 = ct->b / (DELTA*SQ(DELTA)) * ( \
733 (2+ct->b*(ct->b-3))*dDELddelta*SQ(dDELddelta)*DELB \
734 + DELB*SQ(DELTA)*d3DELddelta3 \
735 + 3*(ct->b-1) * DELB * DELTA * dDELddelta * d2DELddelta2 \
736 )
737
738 /**
739 Residual part of helmholtz function.
740 */
741 double helm_resid(double tau, double delta, const HelmholtzData *data){
742 double dell,ldell, term, sum, res = 0;
743 unsigned n, i;
744 const HelmholtzPowTerm *pt;
745 const HelmholtzGausTerm *gt;
746 const HelmholtzCritTerm *ct;
747
748 n = data->np;
749 pt = &(data->pt[0]);
750
751 #ifdef RESID_DEBUG
752 fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
753 #endif
754
755 /* power terms */
756 sum = 0;
757 dell = ipow(delta,pt->l);
758 ldell = pt->l * dell;
759 unsigned oldl;
760 for(i=0; i<n; ++i){
761 term = pt->a * pow(tau, pt->t) * ipow(delta, pt->d);
762 sum += term;
763 #ifdef RESID_DEBUG
764 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 oldl = pt->l;
772 ++pt;
773 if(i+1==n || oldl != pt->l){
774 if(oldl == 0){
775 #ifdef RESID_DEBUG
776 fprintf(stderr,"linear ");
777 #endif
778 res += sum;
779 }else{
780 #ifdef RESID_DEBUG
781 fprintf(stderr,"exp dell=%f, exp(-dell)=%f sum=%f: ",dell,exp(-dell),sum);
782 #endif
783 res += sum * exp(-dell);
784 }
785 #ifdef RESID_DEBUG
786 fprintf(stderr,"i = %d, res = %f\n",i,res);
787 #endif
788 sum = 0;
789 dell = ipow(delta,pt->l);
790 ldell = pt->l*dell;
791 }
792 }
793 assert(!__isnan(res));
794
795 /* gaussian terms */
796 n = data->ng;
797 //fprintf(stderr,"THERE ARE %d GAUSSIAN TERMS\n",n);
798 gt = &(data->gt[0]);
799 for(i=0; i<n; ++i){
800 #ifdef RESID_DEBUG
801 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 double e1 = -gt->alpha*SQ(d1) - gt->beta*SQ(t1);
806 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 assert(!__isnan(res));
812
813 /* critical terms */
814 n = data->nc;
815 ct = &(data->ct[0]);
816 for(i=0; i<n; ++i){
817 #ifdef RESID_DEBUG
818 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 #endif
820
821 DEFINE_DELTA;
822 DEFINE_DELB;
823
824 sum = ct->n * DELB * delta * PSI;
825 res += sum;
826 ++ct;
827 }
828 assert(!__isnan(res));
829
830 #ifdef RESID_DEBUG
831 fprintf(stderr,"CALCULATED RESULT FOR phir = %f\n",res);
832 #endif
833 return res;
834 }
835
836 /*=================== FIRST DERIVATIVES =======================*/
837
838 /**
839 Derivative of the helmholtz residual function with respect to
840 delta.
841 */
842 double helm_resid_del(double tau,double delta, const HelmholtzData *data){
843 double sum = 0, res = 0;
844 double dell, ldell;
845 unsigned n, i;
846 const HelmholtzPowTerm *pt;
847 const HelmholtzGausTerm *gt;
848 const HelmholtzCritTerm *ct;
849
850 #ifdef RESID_DEBUG
851 fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
852 #endif
853
854 /* power terms */
855 n = data->np;
856 pt = &(data->pt[0]);
857 dell = ipow(delta,pt->l);
858 ldell = pt->l * dell;
859 unsigned oldl;
860 for(i=0; i<n; ++i){
861 sum += pt->a * pow(tau, pt->t) * ipow(delta, pt->d - 1) * (pt->d - ldell);
862 oldl = pt->l;
863 ++pt;
864 if(i+1==n || oldl != pt->l){
865 if(oldl == 0){
866 res += sum;
867 }else{
868 res += sum * exp(-dell);
869 }
870 sum = 0;
871 dell = ipow(delta,pt->l);
872 ldell = pt->l*dell;
873 }
874 }
875
876 /* gaussian terms */
877 n = data->ng;
878 gt = &(data->gt[0]);
879 for(i=0; i<n; ++i){
880 #ifdef RESID_DEBUG
881 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 #endif
883 sum = - gt->n * pow(tau,gt->t) * pow(delta, -1. + gt->d)
884 * (2. * gt->alpha * delta * (delta - gt->epsilon) - gt->d)
885 * exp(-(gt->alpha * SQ(delta-gt->epsilon) + gt->beta*SQ(tau-gt->gamma)));
886 res += sum;
887 #ifdef RESID_DEBUG
888 fprintf(stderr,"sum = %f --> res = %f\n",sum,res);
889 #endif
890 ++gt;
891 }
892
893 /* 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 DEFINE_DELTA;
901 DEFINE_DELB;
902 DEFINE_DPSIDDELTA;
903 DEFINE_DDELDDELTA;
904 DEFINE_DDELBDDELTA;
905
906 #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 sum = ct->n * (DELB * (PSI + delta * dPSIddelta) + dDELbddelta * delta * PSI);
923 res += sum;
924
925 ++ct;
926 }
927
928 return res;
929 }
930
931 /**
932 Derivative of the helmholtz residual function with respect to
933 tau.
934 */
935 double helm_resid_tau(double tau,double delta,const HelmholtzData *data){
936
937 double sum;
938 double res = 0;
939 double delX;
940 unsigned l;
941 unsigned n, i;
942 const HelmholtzPowTerm *pt;
943 const HelmholtzGausTerm *gt;
944 const HelmholtzCritTerm *ct;
945
946 n = data->np;
947 pt = &(data->pt[0]);
948
949 delX = 1;
950
951 l = 0;
952 sum = 0;
953 for(i=0; i<n; ++i){
954 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 }
958 ++pt;
959 //fprintf(stderr,"l = %d\n",l);
960 if(i+1==n || l != pt->l){
961 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 if(i+1!=n){
970 l = pt->l;
971 //fprintf(stderr,"New l = %d\n",l);
972 delX = ipow(delta,l);
973 sum = 0;
974 }
975 }
976 }
977 assert(!__isnan(res));
978
979 //#define RESID_DEBUG
980 /* 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
988 double val2;
989 val2 = -gt->n * pow(tau,gt->t - 1.) * pow(delta, gt->d)
990 * (2. * gt->beta * tau * (tau - gt->gamma) - gt->t)
991 * exp(-(gt->alpha * SQ(delta-gt->epsilon) + gt->beta*SQ(tau-gt->gamma)));
992 res += val2;
993 #ifdef RESID_DEBUG
994 fprintf(stderr,"res = %f\n",res);
995 #endif
996
997 ++gt;
998 }
999 assert(!__isnan(res));
1000
1001 /* 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 DEFINE_DELTA;
1009 DEFINE_DELB;
1010 DEFINE_DDELBDTAU;
1011 DEFINE_DPSIDTAU;
1012
1013 sum = ct->n * delta * (dDELbdtau * PSI + DELB * dPSIdtau);
1014 res += sum;
1015 ++ct;
1016 }
1017 assert(!__isnan(res));
1018
1019 return res;
1020 }
1021
1022
1023 /*=================== SECOND DERIVATIVES =======================*/
1024
1025 /**
1026 Mixed derivative of the helmholtz residual function with respect to
1027 delta and tau.
1028 */
1029 double helm_resid_deltau(double tau,double delta,const HelmholtzData *data){
1030 double dell,ldell, sum = 0, res = 0;
1031 unsigned n, i;
1032 const HelmholtzPowTerm *pt;
1033 const HelmholtzGausTerm *gt;
1034 const HelmholtzCritTerm *ct;
1035
1036 /* 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 ++pt;
1046 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 }
1057
1058 #ifdef TEST
1059 assert(!isinf(res));
1060 #endif
1061
1062 /* 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
1073 double f1 = gt->t - 2*gt->beta*tau*(tau - gt->gamma);
1074 double g1 = gt->d - 2*gt->alpha*delta*(delta - gt->epsilon);
1075
1076 sum = gt->n * f1 * pow(tau,gt->t-1) * g1 * pow(delta,gt->d-1) * exp(e1);
1077
1078 //fprintf(stderr,"sum = %f\n",sum);
1079 res += sum;
1080 #ifdef TEST
1081 assert(!isinf(res));
1082 #endif
1083 ++gt;
1084 }
1085
1086 /* 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 DEFINE_DELTA;
1094 DEFINE_DELB;
1095 DEFINE_DPSIDDELTA;
1096 DEFINE_DDELBDTAU;
1097 DEFINE_DDELDDELTA;
1098
1099 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
1102 double d2PSIddeldtau = 4. * ct->C*ct->D*d1*t1*PSI;
1103
1104 DEFINE_DPSIDTAU;
1105
1106 sum = ct->n * (DELB * (dPSIdtau + delta * d2PSIddeldtau) \
1107 + delta *dDELbdtau*dPSIdtau \
1108 + dDELbdtau*(PSI+delta*dPSIddelta) \
1109 + d2DELbddeldtau*delta*PSI
1110 );
1111 res += sum;
1112 ++ct;
1113 }
1114
1115 #ifdef RESID_DEBUG
1116 fprintf(stderr,"phir = %f\n",res);
1117 #endif
1118
1119 #ifdef TEST
1120 assert(!isnan(res));
1121 assert(!isinf(res));
1122 #endif
1123 return res;
1124 }
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 double sum = 0, res = 0;
1132 double dell, ldell;
1133 unsigned n, i;
1134 const HelmholtzPowTerm *pt;
1135 const HelmholtzGausTerm *gt;
1136 const HelmholtzCritTerm *ct;
1137
1138 if(__isnan(delta)){
1139 abort();
1140 }
1141
1142 #ifdef RESID_DEBUG
1143 fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
1144 #endif
1145
1146 /* power terms */
1147 n = data->np;
1148 pt = &(data->pt[0]);
1149 dell = ipow(delta,pt->l);
1150 ldell = pt->l * dell;
1151 unsigned oldl;
1152 for(i=0; i<n; ++i){
1153 double lpart = pt->l ? SQ(ldell) + ldell*(1. - 2*pt->d - pt->l) : 0;
1154 sum += pt->a * pow(tau, pt->t) * ipow(delta, pt->d - 2) * (pt->d*(pt->d - 1) + lpart);
1155 oldl = pt->l;
1156 ++pt;
1157 if(i+1==n || oldl != pt->l){
1158 if(oldl == 0){
1159 res += sum;
1160 }else{
1161 res += sum * exp(-dell);
1162 }
1163 sum = 0;
1164 dell = ipow(delta,pt->l);
1165 ldell = pt->l*dell;
1166 }
1167 }
1168 if(__isnan(res)){
1169 fprintf(stderr,"tau = %.12e, del = %.12e\n",tau,delta);
1170 }
1171 assert(!__isnan(res));
1172
1173 /* gaussian terms */
1174 n = data->ng;
1175 //fprintf(stderr,"THERE ARE %d GAUSSIAN TERMS\n",n);
1176 gt = &(data->gt[0]);
1177 for(i=0; i<n; ++i){
1178 double s1 = SQ(delta - gt->epsilon);
1179 double f1 = gt->d*(gt->d - 1)
1180 + 2.*gt->alpha*delta * (
1181 delta * (2. * gt->alpha * s1 - 1)
1182 - 2. * gt->d * (delta - gt->epsilon)
1183 );
1184 res += gt->n * pow(tau,gt->t) * pow(delta, gt->d - 2.)
1185 * f1
1186 * exp(-(gt->alpha * s1 + gt->beta*SQ(tau-gt->gamma)));
1187 ++gt;
1188 }
1189 assert(!__isnan(res));
1190
1191 /* critical terms */
1192 n = data->nc;
1193 ct = &(data->ct[0]);
1194 for(i=0; i<n; ++i){
1195 #ifdef RESID_DEBUG
1196 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);
1197 #endif
1198
1199 DEFINE_DELTA;
1200 DEFINE_DELB;
1201 DEFINE_DPSIDDELTA;
1202 DEFINE_DDELDDELTA;
1203 DEFINE_DDELBDDELTA;
1204
1205 DEFINE_D2DELDDELTA2;
1206 DEFINE_D2DELBDDELTA2;
1207
1208 DEFINE_D2PSIDDELTA2;
1209
1210 sum = ct->n * (DELB*(2.*dPSIddelta + delta*d2PSIddelta2) + 2.*dDELbddelta*(PSI+delta*dPSIddelta) + d2DELbddelta2*delta*PSI);
1211
1212 res += sum;
1213 ++ct;
1214 }
1215 assert(!__isnan(res));
1216
1217 return res;
1218 }
1219
1220
1221
1222 /**
1223 Residual part of helmholtz function.
1224 */
1225 double helm_resid_tautau(double tau, double delta, const HelmholtzData *data){
1226 double dell,ldell, term, sum, res = 0;
1227 unsigned n, i;
1228 const HelmholtzPowTerm *pt;
1229 const HelmholtzGausTerm *gt;
1230 const HelmholtzCritTerm *ct;
1231
1232 n = data->np;
1233 pt = &(data->pt[0]);
1234
1235 #ifdef RESID_DEBUG
1236 fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
1237 #endif
1238
1239 /* power terms */
1240 sum = 0;
1241 dell = ipow(delta,pt->l);
1242 ldell = pt->l * dell;
1243 unsigned oldl;
1244 for(i=0; i<n; ++i){
1245 term = pt->a * pt->t * (pt->t - 1) * pow(tau, pt->t - 2) * ipow(delta, pt->d);
1246 sum += term;
1247 #ifdef RESID_DEBUG
1248 fprintf(stderr,"i = %d, a=%e, t=%f, d=%d, term = %f, sum = %f",i,pt->a,pt->t,pt->d,term,sum);
1249 if(pt->l==0){
1250 fprintf(stderr,",row=%e\n",term);
1251 }else{
1252 fprintf(stderr,",row=%e\n,",term*exp(-dell));
1253 }
1254 #endif
1255 oldl = pt->l;
1256 ++pt;
1257 if(i+1==n || oldl != pt->l){
1258 if(oldl == 0){
1259 #ifdef RESID_DEBUG
1260 fprintf(stderr,"linear ");
1261 #endif
1262 res += sum;
1263 }else{
1264 #ifdef RESID_DEBUG
1265 fprintf(stderr,"exp dell=%f, exp(-dell)=%f sum=%f: ",dell,exp(-dell),sum);
1266 #endif
1267 res += sum * exp(-dell);
1268 }
1269 #ifdef RESID_DEBUG
1270 fprintf(stderr,"i = %d, res = %f\n",i,res);
1271 #endif
1272 sum = 0;
1273 dell = ipow(delta,pt->l);
1274 ldell = pt->l*dell;
1275 }
1276 }
1277
1278 /* gaussian terms */
1279 n = data->ng;
1280 //fprintf(stderr,"THERE ARE %d GAUSSIAN TERMS\n",n);
1281 gt = &(data->gt[0]);
1282 for(i=0; i<n; ++i){
1283 #ifdef RESID_DEBUG
1284 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);
1285 #endif
1286 double d1 = delta - gt->epsilon;
1287 double t1 = tau - gt->gamma;
1288 double f1 = gt->t*(gt->t - 1) + 4. * gt->beta * tau * (tau * (gt->beta*SQ(t1) - 0.5) - t1*gt->t);
1289 double e1 = -gt->alpha*SQ(d1) - gt->beta*SQ(t1);
1290 sum = gt->n * f1 * pow(tau,gt->t - 2) * pow(delta,gt->d) * exp(e1);
1291 //fprintf(stderr,"sum = %f\n",sum);
1292 res += sum;
1293 ++gt;
1294 }
1295
1296 /* critical terms */
1297 n = data->nc;
1298 ct = &(data->ct[0]);
1299 for(i=0; i<n; ++i){
1300 #ifdef RESID_DEBUG
1301 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);
1302 #endif
1303 DEFINE_DELTA;
1304 DEFINE_DELB;
1305 DEFINE_DDELBDTAU;
1306 DEFINE_DPSIDTAU;
1307
1308 double d2DELbdtau2 = 2. * ct->b * (DELB/DELTA) + 4. * SQ(theta) * ct->b * (ct->b - 1) * (DELB/SQ(DELTA));
1309
1310 double d2PSIdtau2 = 2. * ct->D * PSI * (2. * ct->D * SQ(t1) -1.);
1311
1312 sum = ct->n * delta * (d2DELbdtau2 * PSI + 2 * dDELbdtau*dPSIdtau + DELB * d2PSIdtau2);
1313 res += sum;
1314 ++ct;
1315 }
1316
1317 #ifdef RESID_DEBUG
1318 fprintf(stderr,"phir_tautau = %f\n",res);
1319 #endif
1320 return res;
1321 }
1322
1323 /* === THIRD DERIVATIVES (this is getting boring now) === */
1324
1325 /**
1326 Third derivative of helmholtz residual function, with respect to
1327 delta (thrice).
1328 */
1329 double helm_resid_deldeldel(double tau,double delta,const HelmholtzData *data){
1330 double sum = 0, res = 0;
1331 double D,L;
1332 unsigned n, i;
1333 const HelmholtzPowTerm *pt;
1334 const HelmholtzGausTerm *gt;
1335 const HelmholtzCritTerm *ct;
1336
1337 #ifdef RESID_DEBUG
1338 fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
1339 #endif
1340
1341 #if 1
1342 /* major shortcut, but not efficient */
1343 double ddel = 0.0000000001;
1344 return (helm_resid_deldel(tau,delta+ddel,data) - helm_resid_deldel(tau,delta,data))/ddel;
1345 #endif
1346
1347 /* seem to be errors in the following, still haven't tracked them all down. */
1348
1349 #if 0
1350 /* wxmaxima code:
1351 a*delta^d*%e^(-delta^l)*tau^t
1352 diff(%,delta,3);
1353 */
1354 /* power terms */
1355 n = data->np;
1356 pt = &(data->pt[0]);
1357 D = ipow(delta,pt->l);
1358 unsigned oldl;
1359 for(i=0; i<n; ++i){
1360 double d = pt->d;
1361 double l = pt->l;
1362 double lpart = pt->l
1363 ? D*((D-1)*(D-2)-1) * l*SQ(l)
1364 + 3*D*(1-d)*(D-1) * SQ(l)
1365 + D*(3*SQ(d-1)-1) * l
1366 : 0;
1367 sum += pt->a * pow(tau,pt->t) * ipow(delta, d-3) * (d*(d-1)*(d-2) + lpart);
1368 oldl = pt->l;
1369 ++pt;
1370 if(i+1==n || oldl != pt->l){
1371 if(oldl == 0){
1372 res += sum; // note special meaning of l==0 case: no exponential
1373 }else{
1374 res += sum * exp(-D);
1375 }
1376 sum = 0;
1377 D = ipow(delta,pt->l);
1378 }
1379 }
1380
1381 //fprintf(stderr,"DELDELDEL fiff = %f, sum = %f ",fdiff, res);
1382 #endif
1383
1384 #if 0
1385 /* gaussian terms */
1386 n = data->ng;
1387 //fprintf(stderr,"THERE ARE %d GAUSSIAN TERMS\n",n);
1388 gt = &(data->gt[0]);
1389 for(i=0; i<n; ++i){
1390 double D = delta - gt->epsilon;
1391 double D2 = SQ(D);
1392 double T2 = SQ(tau - gt->gamma);
1393 double A = gt->alpha * delta;
1394 double A2 = SQ(A);
1395 double d = gt->d;
1396 double d2 = SQ(d);
1397
1398 // this expression calculated from wxMaxima using subsitutions for
1399 // D=delta-epsilon and A=alpha*delta.
1400 double f1 =
1401 - (8*A*A2) * D*D2
1402 + (12*d * A2) * D2
1403 + (12 * delta * A2 + (6*d - 6*d2)*A) * D
1404 - (6 * d * delta * A + d*d2 - 3*d2 + 2*d);
1405
1406 res += gt->n * pow(tau,gt->t) * pow(delta, d - 3.)
1407 * f1
1408 * exp(-(gt->alpha * D2 + gt->beta * T2));
1409 ++gt;
1410 }
1411 #endif
1412
1413 #if 0
1414 /* critical terms */
1415 n = data->nc;
1416 ct = &(data->ct[0]);
1417 for(i=0; i<n; ++i){
1418 #ifdef RESID_DEBUG
1419 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);
1420 #endif
1421
1422 DEFINE_DELTA;
1423 DEFINE_DELB;
1424 DEFINE_DPSIDDELTA;
1425 DEFINE_DDELDDELTA;
1426 DEFINE_DDELBDDELTA;
1427
1428 DEFINE_D2PSIDDELTA2;
1429 DEFINE_D2DELDDELTA2;
1430 DEFINE_D2DELBDDELTA2;
1431
1432 DEFINE_D3PSIDDELTA3;
1433 DEFINE_D3DELDDELTA3;
1434 DEFINE_D3DELBDDELTA3;
1435
1436 sum = ct->n * (
1437 delta * (DELB*d3PSIddelta3 + 3 * dDELbddelta * d2PSIddelta2 + 3 * d2DELbddelta2 * dPSIddelta + PSI * d3DELbddelta3)
1438 + 3 * (DELB*d2PSIddelta2 + 2 * dDELbddelta * dPSIddelta + d2DELbddelta2 * PSI)
1439 );
1440
1441 res += sum;
1442 ++ct;
1443 }
1444 #endif
1445
1446 return res;
1447 }
1448
1449
1450
1451

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