/[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 2215 - (show annotations) (download) (as text)
Mon Jul 12 23:11:29 2010 UTC (13 years, 11 months ago) by jpye
File MIME type: text/x-csrc
File size: 31420 byte(s)
Copying HongKe's methane correlation to trunk (pending testing?) and fixing helmholtz_g as suggested by Shrikantch (also needs double-checking).
Fixed problem with GSL in fprops/SConstruct.
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
32 #ifdef TEST
33 #include <assert.h>
34 #include <stdlib.h>
35 #include <stdio.h>
36 #endif
37
38 /* macros and forward decls */
39
40 #define SQ(X) ((X)*(X))
41
42 #include "helmholtz_impl.h"
43
44 /* calculate tau and delta using a macro -- is used in most functions */
45 #define DEFINE_TD \
46 double tau = data->T_star / T; \
47 double delta = rho / data->rho_star
48
49 /**
50 Function to calculate pressure from Helmholtz free energy EOS, given temperature
51 and mass density.
52
53 @param T temperature in K
54 @param rho mass density in kg/m続
55 @return pressure in Pa???
56 */
57 double helmholtz_p(double T, double rho, const HelmholtzData *data){
58 DEFINE_TD;
59
60 #ifdef TEST
61 assert(data->rho_star!=0);
62 assert(T!=0);
63 assert(!isnan(tau));
64 assert(!isnan(delta));
65 assert(!isnan(data->R));
66
67 //fprintf(stderr,"p calc: T = %f\n",T);
68 //fprintf(stderr,"p calc: tau = %f\n",tau);
69 //fprintf(stderr,"p calc: rho = %f\n",rho);
70 //fprintf(stderr,"p calc: delta = %f\n",delta);
71 //fprintf(stderr,"p calc: R*T*rho = %f\n",data->R * T * rho);
72
73 //fprintf(stderr,"T = %f\n", T);
74 //fprintf(stderr,"rhob = %f, rhob* = %f, delta = %f\n", rho/data->M, data->rho_star/data->M, delta);
75 #endif
76
77 return data->R * T * rho * (1 + delta * helm_resid_del(tau,delta,data));
78 }
79
80 /**
81 Function to calculate internal energy from Helmholtz free energy EOS, given
82 temperature and mass density.
83
84 @param T temperature in K
85 @param rho mass density in kg/m続
86 @return internal energy in ???
87 */
88 double helmholtz_u(double T, double rho, const HelmholtzData *data){
89 DEFINE_TD;
90
91 #ifdef TEST
92 assert(data->rho_star!=0);
93 assert(T!=0);
94 assert(!isnan(tau));
95 assert(!isnan(delta));
96 assert(!isnan(data->R));
97 #endif
98
99 #if 0
100 fprintf(stderr,"ideal_tau = %f\n",helm_ideal_tau(tau,delta,data->ideal));
101 fprintf(stderr,"resid_tau = %f\n",helm_resid_tau(tau,delta,data));
102 fprintf(stderr,"R T = %f\n",data->R * data->T_star);
103 #endif
104
105 return data->R * data->T_star * (helm_ideal_tau(tau,delta,data->ideal) + helm_resid_tau(tau,delta,data));
106 }
107
108 /**
109 Function to calculate enthalpy from Helmholtz free energy EOS, given
110 temperature and mass density.
111
112 @param T temperature in K
113 @param rho mass density in kg/m続
114 @return enthalpy in J/kg
115 */
116 double helmholtz_h(double T, double rho, const HelmholtzData *data){
117 DEFINE_TD;
118
119 #ifdef TEST
120 assert(data->rho_star!=0);
121 assert(T!=0);
122 assert(!isnan(tau));
123 assert(!isnan(delta));
124 assert(!isnan(data->R));
125 #endif
126
127 return data->R * T * (1 + tau * (helm_ideal_tau(tau,delta,data->ideal) + helm_resid_tau(tau,delta,data)) + delta*helm_resid_del(tau,delta,data));
128 }
129
130 /**
131 Function to calculate entropy from Helmholtz free energy EOS, given
132 temperature and mass density.
133
134 @param T temperature in K
135 @param rho mass density in kg/m続
136 @return entropy in J/kgK
137 */
138 double helmholtz_s(double T, double rho, const HelmholtzData *data){
139 DEFINE_TD;
140
141 #ifdef ENTROPY_DEBUG
142 assert(data->rho_star!=0);
143 assert(T!=0);
144 assert(!isnan(tau));
145 assert(!isnan(delta));
146 assert(!isnan(data->R));
147
148 fprintf(stderr,"helm_ideal_tau = %f\n",helm_ideal_tau(tau,delta,data->ideal));
149 fprintf(stderr,"helm_resid_tau = %f\n",helm_resid_tau(tau,delta,data));
150 fprintf(stderr,"helm_ideal = %f\n",helm_ideal(tau,delta,data->ideal));
151 fprintf(stderr,"helm_resid = %f\n",helm_resid(tau,delta,data));
152 #endif
153 return data->R * (
154 tau * (helm_ideal_tau(tau,delta,data->ideal) + helm_resid_tau(tau,delta,data))
155 - (helm_ideal(tau,delta,data->ideal) + helm_resid(tau,delta,data))
156 );
157 }
158
159 /**
160 Function to calculate Helmholtz energy from the Helmholtz free energy EOS,
161 given temperature and mass density.
162
163 @param T temperature in K
164 @param rho mass density in kg/m続
165 @return Helmholtz energy 'a', in J/kg
166 */
167 double helmholtz_a(double T, double rho, const HelmholtzData *data){
168 DEFINE_TD;
169
170 #ifdef TEST
171 assert(data->rho_star!=0);
172 assert(T!=0);
173 assert(!isnan(tau));
174 assert(!isnan(delta));
175 assert(!isnan(data->R));
176 #endif
177
178 #ifdef HELMHOLTZ_DEBUG
179 fprintf(stderr,"helmholtz_a: T = %f, rho = %f\n",T,rho);
180 fprintf(stderr,"multiplying by RT = %f\n",data->R*T);
181 #endif
182
183 return data->R * T * (helm_ideal(tau,delta,data->ideal) + helm_resid(tau,delta,data));
184 }
185
186 /**
187 Function to calculate isochoric heat capacity from the Helmholtz free energy
188 EOS given temperature and mass density.
189
190 @param T temperature in K
191 @param rho mass density in kg/m続
192 @return Isochoric specific heat capacity in J/kg/K.
193 */
194 double helmholtz_cv(double T, double rho, const HelmholtzData *data){
195 DEFINE_TD;
196
197 return - data->R * SQ(tau) * (helm_ideal_tautau(tau,data->ideal) + helm_resid_tautau(tau,delta,data));
198 }
199
200 /**
201 Function to calculate isobaric heat capacity from the Helmholtz free energy
202 EOS given temperature and mass density.
203
204 @param T temperature in K
205 @param rho mass density in kg/m続
206 @return Isobaric specific heat capacity in J/kg/K.
207 */
208 double helmholtz_cp(double T, double rho, const HelmholtzData *data){
209 DEFINE_TD;
210
211 double phir_d = helm_resid_del(tau,delta,data);
212 double phir_dd = helm_resid_deldel(tau,delta,data);
213 double phir_dt = helm_resid_deltau(tau,delta,data);
214
215 /* note similarities with helmholtz_w */
216 double temp1 = 1 + 2*delta*phir_d + SQ(delta)*phir_dd;
217 double temp2 = 1 + delta*phir_d - delta*tau*phir_dt;
218 double temp3 = -SQ(tau)*(helm_ideal_tautau(tau,data->ideal) + helm_resid_tautau(tau,delta,data));
219
220 return data->R * (temp3 + SQ(temp2)/temp1);
221 }
222
223
224 /**
225 Function to calculate the speed of sound in a fluid from the Helmholtz free
226 energy EOS, given temperature and mass density.
227
228 @param T temperature in K
229 @param rho mass density in kg/m続
230 @return Speed of sound in m/s.
231 */
232 double helmholtz_w(double T, double rho, const HelmholtzData *data){
233 DEFINE_TD;
234
235 double phir_d = helm_resid_del(tau,delta,data);
236 double phir_dd = helm_resid_deldel(tau,delta,data);
237 double phir_dt = helm_resid_deltau(tau,delta,data);
238
239 /* note similarities with helmholtz_cp */
240 double temp1 = 1. + 2.*delta*phir_d + SQ(delta)*phir_dd;
241 double temp2 = 1. + delta*phir_d - delta*tau*phir_dt;
242 double temp3 = -SQ(tau)*(helm_ideal_tautau(tau,data->ideal) + helm_resid_tautau(tau,delta,data));
243
244 return sqrt(data->R * T * (temp1 + SQ(temp2)/temp3));
245
246 }
247
248 /**
249 Function to calculate the Gibbs energy fluid from the Helmholtz free
250 energy EOS, given temperature and mass density.
251
252 @param T temperature in K
253 @param rho mass density in kg/m続
254 @return Gibbs energy, in J/kg.
255 */
256 double helmholtz_g(double T, double rho, const HelmholtzData *data){
257 DEFINE_TD;
258
259 double phir_d = helm_resid_del(tau,delta,data);
260 double phir = helm_resid(tau,delta,data);
261 double phi0 = helm_ideal(tau,delta,data->ideal);
262
263 return data->R * T * (phi0 + phir + 1. - delta * phir_d);
264 }
265
266 /**
267 Calculation zero-pressure specific heat capacity
268 */
269 double helmholtz_cp0(double T, const HelmholtzData *data){
270 double val = helm_cp0(T,data->ideal);
271 #if 0
272 fprintf(stderr,"val = %f\n",val);
273 #endif
274 return val;
275 }
276
277 /**
278 Solve density given temperature and pressure and bounds for density,
279 used in solving the saturation curve. We use a single-variable Newton
280 method to solve for the desired pressure by iteratively adjusting density.
281
282 @param rho_l lower bound on density
283 @param rho_u upper bound on density
284 @param rho (returned) solved value of density
285 @return 0 on success (1 = didn't converge within permitted iterations)
286 */
287 static int helm_find_rho_tp(double T, double p, double *rho
288 , const double rho_l, const double rho_u, const HelmholtzData *data
289 ){
290 double rho_1 = 0.5 *(rho_l + rho_u);
291 double p_1, dpdrho;
292 double niter = 0, maxiter = 100;
293 int res = 1; /* 1 = exceeded iterations, 0 = converged */
294
295 #ifdef TEST
296 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);
297 #endif
298
299 #if 0
300 double r;
301 for(r = rho_l; r <= rho_u; r += (rho_u-rho_l)/20.){
302 p_1 = helmholtz_p(T, r, data);
303 fprintf(stderr,"T = %f, rho = %f --> p = %f MPa\n",T, r, p_1/1e6);
304 }
305 #endif
306
307 while(res){
308 if(++niter>maxiter)break;
309 p_1 = helmholtz_p(T, rho_1, data);
310 if(p_1 < p){
311 p_1 = 0.5*(p_1 + p);
312 }
313 //fprintf(stderr,"T = %f, rho = %f --> p = %f MPa, err = %f %%\n",T, rho_1, p_1/1e6, (p - p_1)/p *100);
314 dpdrho = helmholtz_dpdrho_T(T, rho_1, data);
315 if(dpdrho < 0){
316 if(rho_l < data->rho_star){
317 /* looking for gas solution */
318 rho_1 = 0.5*(*rho + rho_l);
319 }
320 }
321 rho_1 += (p - p_1)/dpdrho;
322 if(rho_1 > rho_u){
323 rho_1 = 0.5*(*rho + rho_u);
324 //fprintf(stderr,"UPPERBOUND ");
325 }
326 if(rho_1 < rho_l){
327 rho_1 = 0.5*(*rho + rho_l);
328 //fprintf(stderr,"LOWERBOUND ");
329 }
330 *rho = rho_1;
331
332 if(fabs((p - p_1)/p) < 1e-7){
333 #ifdef TEST
334 fprintf(stderr,"Converged to p = %f MPa with rho = %f\n", p/1e6, rho_1);
335 #endif
336 res = 0;
337 }
338 }
339
340 return res;
341 }
342
343 /**
344 Calculation of saturation conditions given temperature
345
346 @param T temperature [K]
347 @param rho_f (returned) saturated liquid density [kg/m続]
348 @param rho_g (returned) saturated gas density [kg/m続]
349 @param p pressure p_sat(T) [Pa]
350 @return 0 on success
351 */
352 int helmholtz_sat_t(double T, double *p, double *rho_f, double *rho_g, const HelmholtzData *data){
353 double tau = data->T_star / T;
354
355 if(T >= data->T_star){
356 #ifdef TEST
357 fprintf(stderr,"ERROR: temperature exceeds critical temperature in helmholtz_sat_t.\n");
358 #endif
359 /* return some reasonable values */
360 *rho_f = data->rho_star;
361 *rho_g = data->rho_star;
362 *p = helmholtz_p(data->T_star, data->rho_star, data);
363 return 1; /* error status */
364 }
365
366 /* get a first estimate of saturation pressure using acentric factor */
367
368 /* critical pressure */
369 #ifdef TEST
370 fprintf(stderr,"T_c = %f, rho_c = %f\n",data->T_star, data->rho_star);
371 #endif
372 double p_c = helmholtz_p(data->T_star, data->rho_star, data);
373
374 #ifdef TEST
375 fprintf(stderr,"Critical pressure = %f MPa\n",p_c/1.e6);
376 fprintf(stderr,"Acentric factor = %f\n",data->omega);
377 #endif
378
379
380 /* maybe use http://dx.doi.org/10.1016/S0009-2509(02)00017-9 for this part? */
381 /* FIXME need to cite this formula */
382 *p = p_c * pow(10.,(data->omega + 1.)*-7./3.*(tau - 1.));
383
384 #ifdef TEST
385 fprintf(stderr,"Estimated p_sat(T=%f) = %f MPa\n",T,(*p)/1e6);
386 #endif
387
388 if(tau < 1.01){
389 /* close to critical point: need a different approach */
390 #ifdef TEST
391 fprintf(stderr,"ERROR: not implemented\n");
392 #endif
393 return 1;
394 }else{
395 double niter = 0;
396 (void)niter;
397
398 int res = helm_find_rho_tp(T, *p, rho_f, data->rho_star,data->rho_star*10, data);
399 if(res){
400 #ifdef TEST
401 fprintf(stderr,"ERROR: failed to solve rho_f\n");
402 #endif
403 }
404
405 res = helm_find_rho_tp(T, *p, rho_g, 0.001*data->rho_star, data->rho_star, data);
406 if(res){
407 #ifdef TEST
408 fprintf(stderr,"ERROR: failed to solve rho_g\n");
409 #endif
410 }
411
412 #ifdef TEST
413 fprintf(stderr,"p = %f MPa: rho_f = %f, rho_g = %f\n", *p, *rho_f, *rho_g);
414 #endif
415
416 double LHS = *p/data->R/T*(1./(*rho_g) - 1./(*rho_f)) - log((*rho_f)/(*rho_g));
417 double delta_f = (*rho_f) / data->rho_star;
418 double delta_g = (*rho_g) / data->rho_star;
419 double RHS = helm_resid(delta_f,tau,data) - helm_resid(delta_g,tau,data);
420
421 double err = LHS - RHS;
422 (void)err;
423
424 #ifdef TEST
425 fprintf(stderr,"LHS = %f, RHS = %f, err = %f\n",LHS, RHS, err);
426 #endif
427 /* away from critical point... */
428 *rho_f = data->rho_star;
429 *rho_g = data->rho_star;
430 #ifdef TEST
431 fprintf(stderr,"ERROR: not implemented\n");
432 exit(1);
433 #endif
434 return 1;
435 }
436 }
437
438 /*----------------------------------------------------------------------------
439 PARTIAL DERIVATIVES
440 */
441
442 /**
443 Calculate partial derivative of p with respect to T, with rho constant
444 */
445 double helmholtz_dpdT_rho(double T, double rho, const HelmholtzData *data){
446 DEFINE_TD;
447
448 double phir_del = helm_resid_del(tau,delta,data);
449 double phir_deltau = helm_resid_deltau(tau,delta,data);
450 #ifdef TEST
451 assert(!isinf(phir_del));
452 assert(!isinf(phir_deltau));
453 assert(!isnan(phir_del));
454 assert(!isnan(phir_deltau));
455 assert(!isnan(data->R));
456 assert(!isnan(rho));
457 assert(!isnan(tau));
458 #endif
459
460 double res = data->R * rho * (1 + delta*phir_del - delta*tau*phir_deltau);
461
462 #ifdef TEST
463 assert(!isnan(res));
464 assert(!isinf(res));
465 #endif
466 return res;
467 }
468
469 /**
470 Calculate partial derivative of p with respect to rho, with T constant
471 */
472 double helmholtz_dpdrho_T(double T, double rho, const HelmholtzData *data){
473 DEFINE_TD;
474
475 double phir_del = helm_resid_del(tau,delta,data);
476 double phir_deldel = helm_resid_deldel(tau,delta,data);
477 #ifdef TEST
478 assert(!isinf(phir_del));
479 assert(!isinf(phir_deldel));
480 #endif
481 return data->R * T * (1 + 2*delta*phir_del + SQ(delta)*phir_deldel);
482 }
483
484 /**
485 Calculate partial derivative of h with respect to T, with rho constant
486 */
487 double helmholtz_dhdT_rho(double T, double rho, const HelmholtzData *data){
488 DEFINE_TD;
489
490 double phir_del = helm_resid_del(tau,delta,data);
491 double phir_deltau = helm_resid_deltau(tau,delta,data);
492 double phir_tautau = helm_resid_tautau(tau,delta,data);
493 double phi0_tautau = helm_ideal_tautau(tau,data->ideal);
494
495 //fprintf(stderr,"phir_del = %f, phir_deltau = %f, phir_tautau = %f, phi0_tautau = %f\n",phir_del,phir_deltau,phir_tautau,phi0_tautau);
496
497 //return (helmholtz_h(T+0.01,rho,data) - helmholtz_h(T,rho,data)) / 0.01;
498 return data->R * (1. + delta*phir_del - SQ(tau)*(phi0_tautau + phir_tautau) - delta*tau*phir_deltau);
499 }
500
501 /**
502 Calculate partial derivative of h with respect to rho, with T constant
503 */
504 double helmholtz_dhdrho_T(double T, double rho, const HelmholtzData *data){
505 DEFINE_TD;
506
507 double phir_del = helm_resid_del(tau,delta,data);
508 double phir_deltau = helm_resid_deltau(tau,delta,data);
509 double phir_deldel = helm_resid_deldel(tau,delta,data);
510
511 return data->R * T / rho * (tau*delta*(0 + phir_deltau) + delta * phir_del + SQ(delta)*phir_deldel);
512 }
513
514
515 /**
516 Calculate partial derivative of u with respect to T, with rho constant
517 */
518 double helmholtz_dudT_rho(double T, double rho, const HelmholtzData *data){
519 DEFINE_TD;
520
521 double phir_tautau = helm_resid_tautau(tau,delta,data);
522 double phi0_tautau = helm_ideal_tautau(tau,data->ideal);
523
524 return -data->R * SQ(tau) * (phi0_tautau + phir_tautau);
525 }
526
527 /**
528 Calculate partial derivative of u with respect to rho, with T constant
529 */
530 double helmholtz_dudrho_T(double T, double rho, const HelmholtzData *data){
531 DEFINE_TD;
532
533 double phir_deltau = helm_resid_deltau(tau,delta,data);
534
535 return data->R * T / rho * (tau * delta * phir_deltau);
536 }
537
538 /*---------------------------------------------
539 UTILITY FUNCTION(S)
540 */
541
542 /* ipow: public domain by Mark Stephen with suggestions by Keiichi Nakasato */
543 static double ipow(double x, int n){
544 double t = 1.0;
545
546 if(!n)return 1.0; /* At the top. x^0 = 1 */
547
548 if (n < 0){
549 n = -n;
550 x = 1.0/x; /* error if x == 0. Good */
551 } /* ZTC/SC returns inf, which is even better */
552
553 if (x == 0.0)return 0.0;
554
555 do{
556 if(n & 1)t *= x;
557 n /= 2; /* KN prefers if (n/=2) x*=x; This avoids an */
558 x *= x; /* unnecessary but benign multiplication on */
559 }while(n); /* the last pass, but the comparison is always
560 true _except_ on the last pass. */
561
562 return t;
563 }
564
565 //#define RESID_DEBUG
566
567 /*
568 We avoid duplication by using the following #defines for common code in
569 calculation of critical terms.
570 */
571 #define DEFINE_DELTA \
572 double d1 = delta - 1.; \
573 double t1 = tau - 1.; \
574 double d12 = SQ(d1); \
575 double theta = (1. - tau) + ct->A * pow(d12, 0.5/ct->beta); \
576 double psi = exp(-ct->C*d12 - ct->D*SQ(t1)); \
577 double DELTA = SQ(theta) + ct->B* pow(d12, ct->a)
578
579 #define DEFINE_DPSIDDELTA \
580 double dpsiddelta = -2. * ct->C * d1 * psi
581
582 #define DEFINE_DDELDDELTA \
583 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))
584
585 #define DEFINE_DDELBDTAU \
586 double dDELbdtau = -2. * theta * ct->b * pow(DELTA, ct->b - 1)
587
588 #define DEFINE_DPSIDTAU \
589 double dpsidtau = -2. * ct->D * t1 * psi
590
591 #define DEFINE_DDELBDDELTA \
592 double dDELbddelta = (DELTA==0?0:ct->b * pow(DELTA,ct->b - 1.) * dDELddelta)
593
594 /**
595 Residual part of helmholtz function.
596 */
597 double helm_resid(double tau, double delta, const HelmholtzData *data){
598 double dell,ldell, term, sum, res = 0;
599 unsigned n, i;
600 const HelmholtzPowTerm *pt;
601 const HelmholtzGausTerm *gt;
602 const HelmholtzCritTerm *ct;
603
604 n = data->np;
605 pt = &(data->pt[0]);
606
607 #ifdef RESID_DEBUG
608 fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
609 #endif
610
611 /* power terms */
612 sum = 0;
613 dell = ipow(delta,pt->l);
614 ldell = pt->l * dell;
615 unsigned oldl;
616 for(i=0; i<n; ++i){
617 term = pt->a * pow(tau, pt->t) * ipow(delta, pt->d);
618 sum += term;
619 #ifdef RESID_DEBUG
620 fprintf(stderr,"i = %d, a=%e, t=%f, d=%d, term = %f, sum = %f",i,pt->a,pt->t,pt->d,term,sum);
621 if(pt->l==0){
622 fprintf(stderr,",row=%e\n",term);
623 }else{
624 fprintf(stderr,",row=%e\n,",term*exp(-dell));
625 }
626 #endif
627 oldl = pt->l;
628 ++pt;
629 if(i+1==n || oldl != pt->l){
630 if(oldl == 0){
631 #ifdef RESID_DEBUG
632 fprintf(stderr,"linear ");
633 #endif
634 res += sum;
635 }else{
636 #ifdef RESID_DEBUG
637 fprintf(stderr,"exp dell=%f, exp(-dell)=%f sum=%f: ",dell,exp(-dell),sum);
638 #endif
639 res += sum * exp(-dell);
640 }
641 #ifdef RESID_DEBUG
642 fprintf(stderr,"i = %d, res = %f\n",i,res);
643 #endif
644 sum = 0;
645 dell = ipow(delta,pt->l);
646 ldell = pt->l*dell;
647 }
648 }
649
650 /* gaussian terms */
651 n = data->ng;
652 //fprintf(stderr,"THERE ARE %d GAUSSIAN TERMS\n",n);
653 gt = &(data->gt[0]);
654 for(i=0; i<n; ++i){
655 #ifdef RESID_DEBUG
656 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);
657 #endif
658 double d1 = delta - gt->epsilon;
659 double t1 = tau - gt->gamma;
660 double e1 = -gt->alpha*SQ(d1) - gt->beta*SQ(t1);
661 sum = gt->n * pow(tau,gt->t) * pow(delta,gt->d) * exp(e1);
662 //fprintf(stderr,"sum = %f\n",sum);
663 res += sum;
664 ++gt;
665 }
666
667 /* critical terms */
668 n = data->nc;
669 ct = &(data->ct[0]);
670 for(i=0; i<n; ++i){
671 #ifdef RESID_DEBUG
672 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);
673 #endif
674
675 DEFINE_DELTA;
676
677 sum = ct->n * pow(DELTA, ct->b) * delta * psi;
678 res += sum;
679 ++ct;
680 }
681
682 #ifdef RESID_DEBUG
683 fprintf(stderr,"CALCULATED RESULT FOR phir = %f\n",res);
684 #endif
685 return res;
686 }
687
688 /*=================== FIRST DERIVATIVES =======================*/
689
690 /**
691 Derivative of the helmholtz residual function with respect to
692 delta.
693 */
694 double helm_resid_del(double tau,double delta, const HelmholtzData *data){
695 double sum = 0, res = 0;
696 double dell, ldell;
697 unsigned n, i;
698 const HelmholtzPowTerm *pt;
699 const HelmholtzGausTerm *gt;
700 const HelmholtzCritTerm *ct;
701
702 #ifdef RESID_DEBUG
703 fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
704 #endif
705
706 /* power terms */
707 n = data->np;
708 pt = &(data->pt[0]);
709 dell = ipow(delta,pt->l);
710 ldell = pt->l * dell;
711 unsigned oldl;
712 for(i=0; i<n; ++i){
713 sum += pt->a * pow(tau, pt->t) * ipow(delta, pt->d - 1) * (pt->d - ldell);
714 oldl = pt->l;
715 ++pt;
716 if(i+1==n || oldl != pt->l){
717 if(oldl == 0){
718 res += sum;
719 }else{
720 res += sum * exp(-dell);
721 }
722 sum = 0;
723 dell = ipow(delta,pt->l);
724 ldell = pt->l*dell;
725 }
726 }
727
728 /* gaussian terms */
729 n = data->ng;
730 gt = &(data->gt[0]);
731 for(i=0; i<n; ++i){
732 #ifdef RESID_DEBUG
733 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);
734 #endif
735 sum = - gt->n * pow(tau,gt->t) * pow(delta, -1. + gt->d)
736 * (2. * gt->alpha * delta * (delta - gt->epsilon) - gt->d)
737 * exp(-(gt->alpha * SQ(delta-gt->epsilon) + gt->beta*SQ(tau-gt->gamma)));
738 res += sum;
739 #ifdef RESID_DEBUG
740 fprintf(stderr,"sum = %f --> res = %f\n",sum,res);
741 #endif
742 ++gt;
743 }
744
745 /* critical terms */
746 n = data->nc;
747 ct = &(data->ct[0]);
748 for(i=0; i<n; ++i){
749 #ifdef RESID_DEBUG
750 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);
751 #endif
752 DEFINE_DELTA;
753 DEFINE_DPSIDDELTA;
754 DEFINE_DDELDDELTA;
755 DEFINE_DDELBDDELTA;
756
757 #if 0
758 if(fabs(dpsiddelta) ==0)fprintf(stderr,"WARNING: dpsiddelta == 0\n");
759 if(fabs(dpsiddelta) ==0)fprintf(stderr,"WARNING: dpsiddelta == 0\n");
760 fprintf(stderr,"psi = %f\n",psi);
761 fprintf(stderr,"DELTA = %f\n",DELTA);
762
763 fprintf(stderr,"dDELddelta = %f\n",dDELddelta);
764 fprintf(stderr,"ct->b - 1. = %f\n",ct->b - 1.);
765 fprintf(stderr,"pow(DELTA,ct->b - 1.) = %f\n",pow(DELTA,ct->b - 1.));
766 assert(!isnan(pow(DELTA,ct->b - 1.)));
767 assert(!isnan(dDELddelta));
768 assert(!isnan(dDELbddelta));
769 //double dDELbddelta = ct->b * pow(DELTA,ct->b - 1.) * dDELddelta
770 fprintf(stderr,"sum = %f\n",sum);
771 if(isnan(sum))fprintf(stderr,"ERROR: sum isnan with i=%d at %d\n",i,__LINE__);
772 #endif
773 sum = ct->n * (pow(DELTA, ct->b) * (psi + delta * dpsiddelta) + dDELbddelta * delta * psi);
774 res += sum;
775
776 ++ct;
777 }
778
779 return res;
780 }
781
782 /**
783 Derivative of the helmholtz residual function with respect to
784 tau.
785 */
786 double helm_resid_tau(double tau,double delta,const HelmholtzData *data){
787
788 double sum;
789 double res = 0;
790 double delX;
791 unsigned l;
792 unsigned n, i;
793 const HelmholtzPowTerm *pt;
794 const HelmholtzGausTerm *gt;
795 const HelmholtzCritTerm *ct;
796
797 n = data->np;
798 pt = &(data->pt[0]);
799
800 delX = 1;
801
802 l = 0;
803 sum = 0;
804 for(i=0; i<n; ++i){
805 if(pt->t){
806 //fprintf(stderr,"i = %d, a = %e, t = %f, d = %d, l = %d\n",i+1, pt->a, pt->t, pt->d, pt->l);
807 sum += pt->a * pow(tau, pt->t - 1) * ipow(delta, pt->d) * pt->t;
808 }
809 ++pt;
810 //fprintf(stderr,"l = %d\n",l);
811 if(i+1==n || l != pt->l){
812 if(l==0){
813 //fprintf(stderr,"Adding non-exp term\n");
814 res += sum;
815 }else{
816 //fprintf(stderr,"Adding exp term with l = %d, delX = %e\n",l,delX);
817 res += sum * exp(-delX);
818 }
819 /* set l to new value */
820 if(i+1!=n){
821 l = pt->l;
822 //fprintf(stderr,"New l = %d\n",l);
823 delX = ipow(delta,l);
824 sum = 0;
825 }
826 }
827 }
828
829 //#define RESID_DEBUG
830 /* gaussian terms */
831 n = data->ng;
832 gt = &(data->gt[0]);
833 for(i=0; i<n; ++i){
834 #ifdef RESID_DEBUG
835 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);
836 #endif
837
838 double val2;
839 val2 = -gt->n * pow(tau,gt->t - 1.) * pow(delta, gt->d)
840 * (2. * gt->beta * tau * (tau - gt->gamma) - gt->t)
841 * exp(-(gt->alpha * SQ(delta-gt->epsilon) + gt->beta*SQ(tau-gt->gamma)));
842 res += val2;
843 #ifdef RESID_DEBUG
844 fprintf(stderr,"res = %f\n",res);
845 #endif
846
847 ++gt;
848 }
849
850 /* critical terms */
851 n = data->nc;
852 ct = &(data->ct[0]);
853 for(i=0; i<n; ++i){
854 #ifdef RESID_DEBUG
855 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);
856 #endif
857 DEFINE_DELTA;
858 DEFINE_DDELBDTAU;
859 DEFINE_DPSIDTAU;
860
861 sum = ct->n * delta * (dDELbdtau * psi + pow(DELTA, ct->b) * dpsidtau);
862 res += sum;
863 ++ct;
864 }
865
866 return res;
867 }
868
869
870 /*=================== SECOND DERIVATIVES =======================*/
871
872 /**
873 Mixed derivative of the helmholtz residual function with respect to
874 delta and tau.
875 */
876 double helm_resid_deltau(double tau,double delta,const HelmholtzData *data){
877 double dell,ldell, term, sum = 0, res = 0;
878 unsigned n, i;
879 const HelmholtzPowTerm *pt;
880 const HelmholtzGausTerm *gt;
881 const HelmholtzCritTerm *ct;
882
883 /* power terms */
884 n = data->np;
885 pt = &(data->pt[0]);
886 dell = ipow(delta,pt->l);
887 ldell = pt->l * dell;
888 unsigned oldl;
889 for(i=0; i<n; ++i){
890 sum += pt->a * pt->t * pow(tau, pt->t - 1) * ipow(delta, pt->d - 1) * (pt->d - ldell);
891 oldl = pt->l;
892 ++pt;
893 if(i+1==n || oldl != pt->l){
894 if(oldl == 0){
895 res += sum;
896 }else{
897 res += sum * exp(-dell);
898 }
899 sum = 0;
900 dell = ipow(delta,pt->l);
901 ldell = pt->l*dell;
902 }
903 }
904
905 #ifdef TEST
906 assert(!isinf(res));
907 #endif
908
909 /* gaussian terms */
910 n = data->ng;
911 gt = &(data->gt[0]);
912 for(i=0; i<n; ++i){
913 #ifdef RESID_DEBUG
914 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);
915 #endif
916 double d1 = delta - gt->epsilon;
917 double t1 = tau - gt->gamma;
918 double e1 = -gt->alpha*SQ(d1) - gt->beta*SQ(t1);
919
920 double f1 = gt->t - 2*gt->beta*tau*(tau - gt->gamma);
921 double g1 = gt->d - 2*gt->alpha*delta*(delta - gt->epsilon);
922
923 sum = gt->n * f1 * pow(tau,gt->t-1) * g1 * pow(delta,gt->d-1) * exp(e1);
924
925 //fprintf(stderr,"sum = %f\n",sum);
926 res += sum;
927 #ifdef TEST
928 assert(!isinf(res));
929 #endif
930 ++gt;
931 }
932
933 /* critical terms */
934 n = data->nc;
935 ct = &(data->ct[0]);
936 for(i=0; i<n; ++i){
937 #ifdef RESID_DEBUG
938 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);
939 #endif
940 DEFINE_DELTA;
941 DEFINE_DPSIDDELTA;
942 DEFINE_DDELBDTAU;
943 DEFINE_DDELDDELTA;
944
945 double d2DELbddeldtau = -ct->A * ct->b * 2./ct->beta * pow(DELTA,ct->b-1)*d1*pow(d12,0.5/ct->beta-1) \
946 - 2. * theta * ct->b * (ct->b - 1) * pow(DELTA,ct->b-2.) * dDELddelta;
947
948 double d2psiddeldtau = 4. * ct->C*ct->D*d1*t1*psi;
949
950 DEFINE_DPSIDTAU;
951
952 sum = ct->n * (pow(DELTA, ct->b) * (dpsidtau + delta * d2psiddeldtau) \
953 + delta *dDELbdtau*dpsidtau \
954 + dDELbdtau*(psi+delta*dpsiddelta) \
955 + d2DELbddeldtau*delta*psi
956 );
957 res += sum;
958 ++ct;
959 }
960
961 #ifdef RESID_DEBUG
962 fprintf(stderr,"phir = %f\n",res);
963 #endif
964
965 #ifdef TEST
966 assert(!isnan(res));
967 assert(!isinf(res));
968 #endif
969 return res;
970 }
971
972 /**
973 Second derivative of helmholtz residual function with respect to
974 delta (twice).
975 */
976 double helm_resid_deldel(double tau,double delta,const HelmholtzData *data){
977 double sum = 0, res = 0;
978 double dell, ldell;
979 unsigned n, i;
980 const HelmholtzPowTerm *pt;
981 const HelmholtzGausTerm *gt;
982 const HelmholtzCritTerm *ct;
983
984 #ifdef RESID_DEBUG
985 fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
986 #endif
987
988 /* power terms */
989 n = data->np;
990 pt = &(data->pt[0]);
991 dell = ipow(delta,pt->l);
992 ldell = pt->l * dell;
993 unsigned oldl;
994 for(i=0; i<n; ++i){
995 double lpart = pt->l ? SQ(ldell) + ldell*(1. - 2*pt->d - pt->l) : 0;
996 sum += pt->a * pow(tau, pt->t) * ipow(delta, pt->d - 2) * (pt->d*(pt->d - 1) + lpart);
997 oldl = pt->l;
998 ++pt;
999 if(i+1==n || oldl != pt->l){
1000 if(oldl == 0){
1001 res += sum;
1002 }else{
1003 res += sum * exp(-dell);
1004 }
1005 sum = 0;
1006 dell = ipow(delta,pt->l);
1007 ldell = pt->l*dell;
1008 }
1009 }
1010
1011 /* gaussian terms */
1012 n = data->ng;
1013 //fprintf(stderr,"THERE ARE %d GAUSSIAN TERMS\n",n);
1014 gt = &(data->gt[0]);
1015 for(i=0; i<n; ++i){
1016 double s1 = SQ(delta - gt->epsilon);
1017 double f1 = gt->d*(gt->d - 1)
1018 + 2.*gt->alpha*delta * (
1019 delta * (2. * gt->alpha * s1 - 1)
1020 - 2. * gt->d * (delta - gt->epsilon)
1021 );
1022 res += gt->n * pow(tau,gt->t) * pow(delta, gt->d - 2.)
1023 * f1
1024 * exp(-(gt->alpha * s1 + gt->beta*SQ(tau-gt->gamma)));
1025 ++gt;
1026 }
1027
1028 /* critical terms */
1029 n = data->nc;
1030 ct = &(data->ct[0]);
1031 for(i=0; i<n; ++i){
1032 #ifdef RESID_DEBUG
1033 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);
1034 #endif
1035
1036 DEFINE_DELTA;
1037 DEFINE_DPSIDDELTA;
1038 DEFINE_DDELDDELTA;
1039 DEFINE_DDELBDDELTA;
1040
1041 double powd12bm1 = pow(d12,0.5/ct->beta-1.);
1042
1043 double d2psiddel2 = (2.*ct->C*d12 - 1.)*2.*ct->C*psi;
1044
1045 double d2DELddel2 = 1./d1*dDELddelta + d12*(
1046 4.*ct->B*ct->a*(ct->a-1.)*pow(d12,ct->a-2.) \
1047 + 2.*SQ(ct->A)*SQ(1./ct->beta)*SQ(powd12bm1) \
1048 + ct->A*theta*4./ct->beta*(0.5/ct->beta-1.)*powd12bm1/d12
1049 );
1050
1051 double d2DELbddel2 = ct->b * (pow(DELTA,ct->b - 1.)*d2DELddel2 + (ct->b-1.)*pow(DELTA,ct->b-2.)*SQ(dDELddelta));
1052
1053 sum = ct->n * (pow(DELTA,ct->b)*(2.*dpsiddelta + delta*d2psiddel2) + 2.*dDELbddelta*(psi+delta*dpsiddelta) + d2DELbddel2*delta*psi);
1054
1055 res += sum;
1056 ++ct;
1057 }
1058
1059 return res;
1060 }
1061
1062
1063
1064 /**
1065 Residual part of helmholtz function.
1066 */
1067 double helm_resid_tautau(double tau, double delta, const HelmholtzData *data){
1068 double dell,ldell, term, sum, res = 0;
1069 unsigned n, i;
1070 const HelmholtzPowTerm *pt;
1071 const HelmholtzGausTerm *gt;
1072 const HelmholtzCritTerm *ct;
1073
1074 n = data->np;
1075 pt = &(data->pt[0]);
1076
1077 #ifdef RESID_DEBUG
1078 fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
1079 #endif
1080
1081 /* power terms */
1082 sum = 0;
1083 dell = ipow(delta,pt->l);
1084 ldell = pt->l * dell;
1085 unsigned oldl;
1086 for(i=0; i<n; ++i){
1087 term = pt->a * pt->t * (pt->t - 1) * pow(tau, pt->t - 2) * ipow(delta, pt->d);
1088 sum += term;
1089 #ifdef RESID_DEBUG
1090 fprintf(stderr,"i = %d, a=%e, t=%f, d=%d, term = %f, sum = %f",i,pt->a,pt->t,pt->d,term,sum);
1091 if(pt->l==0){
1092 fprintf(stderr,",row=%e\n",term);
1093 }else{
1094 fprintf(stderr,",row=%e\n,",term*exp(-dell));
1095 }
1096 #endif
1097 oldl = pt->l;
1098 ++pt;
1099 if(i+1==n || oldl != pt->l){
1100 if(oldl == 0){
1101 #ifdef RESID_DEBUG
1102 fprintf(stderr,"linear ");
1103 #endif
1104 res += sum;
1105 }else{
1106 #ifdef RESID_DEBUG
1107 fprintf(stderr,"exp dell=%f, exp(-dell)=%f sum=%f: ",dell,exp(-dell),sum);
1108 #endif
1109 res += sum * exp(-dell);
1110 }
1111 #ifdef RESID_DEBUG
1112 fprintf(stderr,"i = %d, res = %f\n",i,res);
1113 #endif
1114 sum = 0;
1115 dell = ipow(delta,pt->l);
1116 ldell = pt->l*dell;
1117 }
1118 }
1119
1120 /* gaussian terms */
1121 n = data->ng;
1122 //fprintf(stderr,"THERE ARE %d GAUSSIAN TERMS\n",n);
1123 gt = &(data->gt[0]);
1124 for(i=0; i<n; ++i){
1125 #ifdef RESID_DEBUG
1126 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);
1127 #endif
1128 double d1 = delta - gt->epsilon;
1129 double t1 = tau - gt->gamma;
1130 double f1 = gt->t*(gt->t - 1) + 4. * gt->beta * tau * (tau * (gt->beta*SQ(t1) - 0.5) - t1*gt->t);
1131 double e1 = -gt->alpha*SQ(d1) - gt->beta*SQ(t1);
1132 sum = gt->n * f1 * pow(tau,gt->t - 2) * pow(delta,gt->d) * exp(e1);
1133 //fprintf(stderr,"sum = %f\n",sum);
1134 res += sum;
1135 ++gt;
1136 }
1137
1138 /* critical terms */
1139 n = data->nc;
1140 ct = &(data->ct[0]);
1141 for(i=0; i<n; ++i){
1142 #ifdef RESID_DEBUG
1143 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);
1144 #endif
1145 DEFINE_DELTA;
1146 DEFINE_DDELBDTAU;
1147 DEFINE_DPSIDTAU;
1148
1149 double d2DELbdtau2 = 2. * ct->b * pow(DELTA, ct->b - 1) + 4. * SQ(theta) * ct->b * (ct->b - 1) * pow(DELTA, ct->b - 2);
1150
1151 double d2psidtau2 = 2. * ct->D * psi * (2. * ct->D * SQ(t1) -1.);
1152
1153 sum = ct->n * delta * (d2DELbdtau2 * psi + 2 * dDELbdtau*dpsidtau + pow(DELTA, ct->b) * d2psidtau2);
1154 res += sum;
1155 ++ct;
1156 }
1157
1158 #ifdef RESID_DEBUG
1159 fprintf(stderr,"phir_tautau = %f\n",res);
1160 #endif
1161 return res;
1162 }
1163

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