/[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 2114 - (show annotations) (download) (as text)
Tue Dec 8 08:08:00 2009 UTC (10 years, 7 months ago) by jpye
File MIME type: text/x-csrc
File size: 30911 byte(s)
Starting some more work on saturation curves.
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 Calculation zero-pressure specific heat capacity
250 */
251 double helmholtz_cp0(double T, const HelmholtzData *data){
252 double val = helm_cp0(T,data->ideal);
253 #if 0
254 fprintf(stderr,"val = %f\n",val);
255 #endif
256 return val;
257 }
258
259 /**
260 Solve density given temperature and pressure and bounds for density,
261 used in solving the saturation curve. We use a single-variable Newton
262 method to solve for the desired pressure by iteratively adjusting density.
263
264 @param rho_l lower bound on density
265 @param rho_u upper bound on density
266 @param rho (returned) solved value of density
267 @return 0 on success (1 = didn't converge within permitted iterations)
268 */
269 static int helm_find_rho_tp(double T, double p, double *rho
270 , const double rho_l, const double rho_u, const HelmholtzData *data
271 ){
272 double rho_1 = 0.5 *(rho_l + rho_u);
273 double p_1, dpdrho;
274 double niter = 0, maxiter = 100;
275 int res = 1; /* 1 = exceeded iterations, 0 = converged */
276
277 #ifdef TEST
278 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);
279 #endif
280
281 #if 0
282 double r;
283 for(r = rho_l; r <= rho_u; r += (rho_u-rho_l)/20.){
284 p_1 = helmholtz_p(T, r, data);
285 fprintf(stderr,"T = %f, rho = %f --> p = %f MPa\n",T, r, p_1/1e6);
286 }
287 #endif
288
289 while(res){
290 if(++niter>maxiter)break;
291 p_1 = helmholtz_p(T, rho_1, data);
292 if(p_1 < p){
293 p_1 = 0.5*(p_1 + p);
294 }
295 //fprintf(stderr,"T = %f, rho = %f --> p = %f MPa, err = %f %%\n",T, rho_1, p_1/1e6, (p - p_1)/p *100);
296 dpdrho = helmholtz_dpdrho_T(T, rho_1, data);
297 if(dpdrho < 0){
298 if(rho_l < data->rho_star){
299 /* looking for gas solution */
300 rho_1 = 0.5*(*rho + rho_l);
301 }
302 }
303 rho_1 += (p - p_1)/dpdrho;
304 if(rho_1 > rho_u){
305 rho_1 = 0.5*(*rho + rho_u);
306 //fprintf(stderr,"UPPERBOUND ");
307 }
308 if(rho_1 < rho_l){
309 rho_1 = 0.5*(*rho + rho_l);
310 //fprintf(stderr,"LOWERBOUND ");
311 }
312 *rho = rho_1;
313
314 if(fabs((p - p_1)/p) < 1e-7){
315 #ifdef TEST
316 fprintf(stderr,"Converged to p = %f MPa with rho = %f\n", p/1e6, rho_1);
317 #endif
318 res = 0;
319 }
320 }
321
322 return res;
323 }
324
325 /**
326 Calculation of saturation conditions given temperature
327
328 @param T temperature [K]
329 @param rho_f (returned) saturated liquid density [kg/m続]
330 @param rho_g (returned) saturated gas density [kg/m続]
331 @param p pressure p_sat(T) [Pa]
332 @return 0 on success
333 */
334 int helmholtz_sat_t(double T, double *p, double *rho_f, double *rho_g, const HelmholtzData *data){
335 double tau = data->T_star / T;
336
337 if(T >= data->T_star){
338 #ifdef TEST
339 fprintf(stderr,"ERROR: temperature exceeds critical temperature in helmholtz_sat_t.\n");
340 #endif
341 /* return some reasonable values */
342 *rho_f = data->rho_star;
343 *rho_g = data->rho_star;
344 *p = helmholtz_p(data->T_star, data->rho_star, data);
345 return 1; /* error status */
346 }
347
348 /* get a first estimate of saturation pressure using acentric factor */
349
350 /* critical pressure */
351 #ifdef TEST
352 fprintf(stderr,"T_c = %f, rho_c = %f\n",data->T_star, data->rho_star);
353 #endif
354 double p_c = helmholtz_p(data->T_star, data->rho_star, data);
355
356 #ifdef TEST
357 fprintf(stderr,"Critical pressure = %f MPa\n",p_c/1.e6);
358 fprintf(stderr,"Acentric factor = %f\n",data->omega);
359 #endif
360
361
362 /* maybe use http://dx.doi.org/10.1016/S0009-2509(02)00017-9 for this part? */
363 /* FIXME need to cite this formula */
364 *p = p_c * pow(10.,(data->omega + 1.)*-7./3.*(tau - 1.));
365
366 #ifdef TEST
367 fprintf(stderr,"Estimated p_sat(T=%f) = %f MPa\n",T,(*p)/1e6);
368 #endif
369
370 if(tau < 1.01){
371 /* close to critical point: need a different approach */
372 #ifdef TEST
373 fprintf(stderr,"ERROR: not implemented\n");
374 #endif
375 return 1;
376 }else{
377 double niter = 0;
378 (void)niter;
379
380 int res = helm_find_rho_tp(T, *p, rho_f, data->rho_star,data->rho_star*10, data);
381 if(res){
382 #ifdef TEST
383 fprintf(stderr,"ERROR: failed to solve rho_f\n");
384 #endif
385 }
386
387 res = helm_find_rho_tp(T, *p, rho_g, 0.001*data->rho_star, data->rho_star, data);
388 if(res){
389 #ifdef TEST
390 fprintf(stderr,"ERROR: failed to solve rho_g\n");
391 #endif
392 }
393
394 #ifdef TEST
395 fprintf(stderr,"p = %f MPa: rho_f = %f, rho_g = %f\n", *p, *rho_f, *rho_g);
396 #endif
397
398 double LHS = *p/data->R/T*(1./(*rho_g) - 1./(*rho_f)) - log((*rho_f)/(*rho_g));
399 double delta_f = (*rho_f) / data->rho_star;
400 double delta_g = (*rho_g) / data->rho_star;
401 double RHS = helm_resid(delta_f,tau,data) - helm_resid(delta_g,tau,data);
402
403 double err = LHS - RHS;
404 (void)err;
405
406 #ifdef TEST
407 fprintf(stderr,"LHS = %f, RHS = %f, err = %f\n",LHS, RHS, err);
408 #endif
409 /* away from critical point... */
410 *rho_f = data->rho_star;
411 *rho_g = data->rho_star;
412 #ifdef TEST
413 fprintf(stderr,"ERROR: not implemented\n");
414 exit(1);
415 #endif
416 return 1;
417 }
418 }
419
420 /*----------------------------------------------------------------------------
421 PARTIAL DERIVATIVES
422 */
423
424 /**
425 Calculate partial derivative of p with respect to T, with rho constant
426 */
427 double helmholtz_dpdT_rho(double T, double rho, const HelmholtzData *data){
428 DEFINE_TD;
429
430 double phir_del = helm_resid_del(tau,delta,data);
431 double phir_deltau = helm_resid_deltau(tau,delta,data);
432 #ifdef TEST
433 assert(!isinf(phir_del));
434 assert(!isinf(phir_deltau));
435 assert(!isnan(phir_del));
436 assert(!isnan(phir_deltau));
437 assert(!isnan(data->R));
438 assert(!isnan(rho));
439 assert(!isnan(tau));
440 #endif
441
442 double res = data->R * rho * (1 + delta*phir_del - delta*tau*phir_deltau);
443
444 #ifdef TEST
445 assert(!isnan(res));
446 assert(!isinf(res));
447 #endif
448 return res;
449 }
450
451 /**
452 Calculate partial derivative of p with respect to rho, with T constant
453 */
454 double helmholtz_dpdrho_T(double T, double rho, const HelmholtzData *data){
455 DEFINE_TD;
456
457 double phir_del = helm_resid_del(tau,delta,data);
458 double phir_deldel = helm_resid_deldel(tau,delta,data);
459 #ifdef TEST
460 assert(!isinf(phir_del));
461 assert(!isinf(phir_deldel));
462 #endif
463 return data->R * T * (1 + 2*delta*phir_del + SQ(delta)*phir_deldel);
464 }
465
466 /**
467 Calculate partial derivative of h with respect to T, with rho constant
468 */
469 double helmholtz_dhdT_rho(double T, double rho, const HelmholtzData *data){
470 DEFINE_TD;
471
472 double phir_del = helm_resid_del(tau,delta,data);
473 double phir_deltau = helm_resid_deltau(tau,delta,data);
474 double phir_tautau = helm_resid_tautau(tau,delta,data);
475 double phi0_tautau = helm_ideal_tautau(tau,data->ideal);
476
477 //fprintf(stderr,"phir_del = %f, phir_deltau = %f, phir_tautau = %f, phi0_tautau = %f\n",phir_del,phir_deltau,phir_tautau,phi0_tautau);
478
479 //return (helmholtz_h(T+0.01,rho,data) - helmholtz_h(T,rho,data)) / 0.01;
480 return data->R * (1. + delta*phir_del - SQ(tau)*(phi0_tautau + phir_tautau) - delta*tau*phir_deltau);
481 }
482
483 /**
484 Calculate partial derivative of h with respect to rho, with T constant
485 */
486 double helmholtz_dhdrho_T(double T, double rho, const HelmholtzData *data){
487 DEFINE_TD;
488
489 double phir_del = helm_resid_del(tau,delta,data);
490 double phir_deltau = helm_resid_deltau(tau,delta,data);
491 double phir_deldel = helm_resid_deldel(tau,delta,data);
492
493 return data->R * T / rho * (tau*delta*(0 + phir_deltau) + delta * phir_del + SQ(delta)*phir_deldel);
494 }
495
496
497 /**
498 Calculate partial derivative of u with respect to T, with rho constant
499 */
500 double helmholtz_dudT_rho(double T, double rho, const HelmholtzData *data){
501 DEFINE_TD;
502
503 double phir_tautau = helm_resid_tautau(tau,delta,data);
504 double phi0_tautau = helm_ideal_tautau(tau,data->ideal);
505
506 return -data->R * SQ(tau) * (phi0_tautau + phir_tautau);
507 }
508
509 /**
510 Calculate partial derivative of u with respect to rho, with T constant
511 */
512 double helmholtz_dudrho_T(double T, double rho, const HelmholtzData *data){
513 DEFINE_TD;
514
515 double phir_deltau = helm_resid_deltau(tau,delta,data);
516
517 return data->R * T / rho * (tau * delta * phir_deltau);
518 }
519
520 /*---------------------------------------------
521 UTILITY FUNCTION(S)
522 */
523
524 /* ipow: public domain by Mark Stephen with suggestions by Keiichi Nakasato */
525 static double ipow(double x, int n){
526 double t = 1.0;
527
528 if(!n)return 1.0; /* At the top. x^0 = 1 */
529
530 if (n < 0){
531 n = -n;
532 x = 1.0/x; /* error if x == 0. Good */
533 } /* ZTC/SC returns inf, which is even better */
534
535 if (x == 0.0)return 0.0;
536
537 do{
538 if(n & 1)t *= x;
539 n /= 2; /* KN prefers if (n/=2) x*=x; This avoids an */
540 x *= x; /* unnecessary but benign multiplication on */
541 }while(n); /* the last pass, but the comparison is always
542 true _except_ on the last pass. */
543
544 return t;
545 }
546
547 //#define RESID_DEBUG
548
549 /*
550 We avoid duplication by using the following #defines for common code in
551 calculation of critical terms.
552 */
553 #define DEFINE_DELTA \
554 double d1 = delta - 1.; \
555 double t1 = tau - 1.; \
556 double d12 = SQ(d1); \
557 double theta = (1. - tau) + ct->A * pow(d12, 0.5/ct->beta); \
558 double psi = exp(-ct->C*d12 - ct->D*SQ(t1)); \
559 double DELTA = SQ(theta) + ct->B* pow(d12, ct->a)
560
561 #define DEFINE_DPSIDDELTA \
562 double dpsiddelta = -2. * ct->C * d1 * psi
563
564 #define DEFINE_DDELDDELTA \
565 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))
566
567 #define DEFINE_DDELBDTAU \
568 double dDELbdtau = -2. * theta * ct->b * pow(DELTA, ct->b - 1)
569
570 #define DEFINE_DPSIDTAU \
571 double dpsidtau = -2. * ct->D * t1 * psi
572
573 #define DEFINE_DDELBDDELTA \
574 double dDELbddelta = (DELTA==0?0:ct->b * pow(DELTA,ct->b - 1.) * dDELddelta)
575
576 /**
577 Residual part of helmholtz function.
578 */
579 double helm_resid(double tau, double delta, const HelmholtzData *data){
580 double dell,ldell, term, sum, res = 0;
581 unsigned n, i;
582 const HelmholtzPowTerm *pt;
583 const HelmholtzGausTerm *gt;
584 const HelmholtzCritTerm *ct;
585
586 n = data->np;
587 pt = &(data->pt[0]);
588
589 #ifdef RESID_DEBUG
590 fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
591 #endif
592
593 /* power terms */
594 sum = 0;
595 dell = ipow(delta,pt->l);
596 ldell = pt->l * dell;
597 unsigned oldl;
598 for(i=0; i<n; ++i){
599 term = pt->a * pow(tau, pt->t) * ipow(delta, pt->d);
600 sum += term;
601 #ifdef RESID_DEBUG
602 fprintf(stderr,"i = %d, a=%e, t=%f, d=%d, term = %f, sum = %f",i,pt->a,pt->t,pt->d,term,sum);
603 if(pt->l==0){
604 fprintf(stderr,",row=%e\n",term);
605 }else{
606 fprintf(stderr,",row=%e\n,",term*exp(-dell));
607 }
608 #endif
609 oldl = pt->l;
610 ++pt;
611 if(i+1==n || oldl != pt->l){
612 if(oldl == 0){
613 #ifdef RESID_DEBUG
614 fprintf(stderr,"linear ");
615 #endif
616 res += sum;
617 }else{
618 #ifdef RESID_DEBUG
619 fprintf(stderr,"exp dell=%f, exp(-dell)=%f sum=%f: ",dell,exp(-dell),sum);
620 #endif
621 res += sum * exp(-dell);
622 }
623 #ifdef RESID_DEBUG
624 fprintf(stderr,"i = %d, res = %f\n",i,res);
625 #endif
626 sum = 0;
627 dell = ipow(delta,pt->l);
628 ldell = pt->l*dell;
629 }
630 }
631
632 /* gaussian terms */
633 n = data->ng;
634 //fprintf(stderr,"THERE ARE %d GAUSSIAN TERMS\n",n);
635 gt = &(data->gt[0]);
636 for(i=0; i<n; ++i){
637 #ifdef RESID_DEBUG
638 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);
639 #endif
640 double d1 = delta - gt->epsilon;
641 double t1 = tau - gt->gamma;
642 double e1 = -gt->alpha*SQ(d1) - gt->beta*SQ(t1);
643 sum = gt->n * pow(tau,gt->t) * pow(delta,gt->d) * exp(e1);
644 //fprintf(stderr,"sum = %f\n",sum);
645 res += sum;
646 ++gt;
647 }
648
649 /* critical terms */
650 n = data->nc;
651 ct = &(data->ct[0]);
652 for(i=0; i<n; ++i){
653 #ifdef RESID_DEBUG
654 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);
655 #endif
656
657 DEFINE_DELTA;
658
659 sum = ct->n * pow(DELTA, ct->b) * delta * psi;
660 res += sum;
661 ++ct;
662 }
663
664 #ifdef RESID_DEBUG
665 fprintf(stderr,"CALCULATED RESULT FOR phir = %f\n",res);
666 #endif
667 return res;
668 }
669
670 /*=================== FIRST DERIVATIVES =======================*/
671
672 /**
673 Derivative of the helmholtz residual function with respect to
674 delta.
675 */
676 double helm_resid_del(double tau,double delta, const HelmholtzData *data){
677 double sum = 0, res = 0;
678 double dell, ldell;
679 unsigned n, i;
680 const HelmholtzPowTerm *pt;
681 const HelmholtzGausTerm *gt;
682 const HelmholtzCritTerm *ct;
683
684 #ifdef RESID_DEBUG
685 fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
686 #endif
687
688 /* power terms */
689 n = data->np;
690 pt = &(data->pt[0]);
691 dell = ipow(delta,pt->l);
692 ldell = pt->l * dell;
693 unsigned oldl;
694 for(i=0; i<n; ++i){
695 sum += pt->a * pow(tau, pt->t) * ipow(delta, pt->d - 1) * (pt->d - ldell);
696 oldl = pt->l;
697 ++pt;
698 if(i+1==n || oldl != pt->l){
699 if(oldl == 0){
700 res += sum;
701 }else{
702 res += sum * exp(-dell);
703 }
704 sum = 0;
705 dell = ipow(delta,pt->l);
706 ldell = pt->l*dell;
707 }
708 }
709
710 /* gaussian terms */
711 n = data->ng;
712 gt = &(data->gt[0]);
713 for(i=0; i<n; ++i){
714 #ifdef RESID_DEBUG
715 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);
716 #endif
717 sum = - gt->n * pow(tau,gt->t) * pow(delta, -1. + gt->d)
718 * (2. * gt->alpha * delta * (delta - gt->epsilon) - gt->d)
719 * exp(-(gt->alpha * SQ(delta-gt->epsilon) + gt->beta*SQ(tau-gt->gamma)));
720 res += sum;
721 #ifdef RESID_DEBUG
722 fprintf(stderr,"sum = %f --> res = %f\n",sum,res);
723 #endif
724 ++gt;
725 }
726
727 /* critical terms */
728 n = data->nc;
729 ct = &(data->ct[0]);
730 for(i=0; i<n; ++i){
731 #ifdef RESID_DEBUG
732 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);
733 #endif
734 DEFINE_DELTA;
735 DEFINE_DPSIDDELTA;
736 DEFINE_DDELDDELTA;
737 DEFINE_DDELBDDELTA;
738
739 #if 0
740 if(fabs(dpsiddelta) ==0)fprintf(stderr,"WARNING: dpsiddelta == 0\n");
741 if(fabs(dpsiddelta) ==0)fprintf(stderr,"WARNING: dpsiddelta == 0\n");
742 fprintf(stderr,"psi = %f\n",psi);
743 fprintf(stderr,"DELTA = %f\n",DELTA);
744
745 fprintf(stderr,"dDELddelta = %f\n",dDELddelta);
746 fprintf(stderr,"ct->b - 1. = %f\n",ct->b - 1.);
747 fprintf(stderr,"pow(DELTA,ct->b - 1.) = %f\n",pow(DELTA,ct->b - 1.));
748 assert(!isnan(pow(DELTA,ct->b - 1.)));
749 assert(!isnan(dDELddelta));
750 assert(!isnan(dDELbddelta));
751 //double dDELbddelta = ct->b * pow(DELTA,ct->b - 1.) * dDELddelta
752 fprintf(stderr,"sum = %f\n",sum);
753 if(isnan(sum))fprintf(stderr,"ERROR: sum isnan with i=%d at %d\n",i,__LINE__);
754 #endif
755 sum = ct->n * (pow(DELTA, ct->b) * (psi + delta * dpsiddelta) + dDELbddelta * delta * psi);
756 res += sum;
757
758 ++ct;
759 }
760
761 return res;
762 }
763
764 /**
765 Derivative of the helmholtz residual function with respect to
766 tau.
767 */
768 double helm_resid_tau(double tau,double delta,const HelmholtzData *data){
769
770 double sum;
771 double res = 0;
772 double delX;
773 unsigned l;
774 unsigned n, i;
775 const HelmholtzPowTerm *pt;
776 const HelmholtzGausTerm *gt;
777 const HelmholtzCritTerm *ct;
778
779 n = data->np;
780 pt = &(data->pt[0]);
781
782 delX = 1;
783
784 l = 0;
785 sum = 0;
786 for(i=0; i<n; ++i){
787 if(pt->t){
788 //fprintf(stderr,"i = %d, a = %e, t = %f, d = %d, l = %d\n",i+1, pt->a, pt->t, pt->d, pt->l);
789 sum += pt->a * pow(tau, pt->t - 1) * ipow(delta, pt->d) * pt->t;
790 }
791 ++pt;
792 //fprintf(stderr,"l = %d\n",l);
793 if(i+1==n || l != pt->l){
794 if(l==0){
795 //fprintf(stderr,"Adding non-exp term\n");
796 res += sum;
797 }else{
798 //fprintf(stderr,"Adding exp term with l = %d, delX = %e\n",l,delX);
799 res += sum * exp(-delX);
800 }
801 /* set l to new value */
802 if(i+1!=n){
803 l = pt->l;
804 //fprintf(stderr,"New l = %d\n",l);
805 delX = ipow(delta,l);
806 sum = 0;
807 }
808 }
809 }
810
811 //#define RESID_DEBUG
812 /* gaussian terms */
813 n = data->ng;
814 gt = &(data->gt[0]);
815 for(i=0; i<n; ++i){
816 #ifdef RESID_DEBUG
817 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);
818 #endif
819
820 double val2;
821 val2 = -gt->n * pow(tau,gt->t - 1.) * pow(delta, gt->d)
822 * (2. * gt->beta * tau * (tau - gt->gamma) - gt->t)
823 * exp(-(gt->alpha * SQ(delta-gt->epsilon) + gt->beta*SQ(tau-gt->gamma)));
824 res += val2;
825 #ifdef RESID_DEBUG
826 fprintf(stderr,"res = %f\n",res);
827 #endif
828
829 ++gt;
830 }
831
832 /* critical terms */
833 n = data->nc;
834 ct = &(data->ct[0]);
835 for(i=0; i<n; ++i){
836 #ifdef RESID_DEBUG
837 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);
838 #endif
839 DEFINE_DELTA;
840 DEFINE_DDELBDTAU;
841 DEFINE_DPSIDTAU;
842
843 sum = ct->n * delta * (dDELbdtau * psi + pow(DELTA, ct->b) * dpsidtau);
844 res += sum;
845 ++ct;
846 }
847
848 return res;
849 }
850
851
852 /*=================== SECOND DERIVATIVES =======================*/
853
854 /**
855 Mixed derivative of the helmholtz residual function with respect to
856 delta and tau.
857 */
858 double helm_resid_deltau(double tau,double delta,const HelmholtzData *data){
859 double dell,ldell, term, sum = 0, res = 0;
860 unsigned n, i;
861 const HelmholtzPowTerm *pt;
862 const HelmholtzGausTerm *gt;
863 const HelmholtzCritTerm *ct;
864
865 /* power terms */
866 n = data->np;
867 pt = &(data->pt[0]);
868 dell = ipow(delta,pt->l);
869 ldell = pt->l * dell;
870 unsigned oldl;
871 for(i=0; i<n; ++i){
872 sum += pt->a * pt->t * pow(tau, pt->t - 1) * ipow(delta, pt->d - 1) * (pt->d - ldell);
873 oldl = pt->l;
874 ++pt;
875 if(i+1==n || oldl != pt->l){
876 if(oldl == 0){
877 res += sum;
878 }else{
879 res += sum * exp(-dell);
880 }
881 sum = 0;
882 dell = ipow(delta,pt->l);
883 ldell = pt->l*dell;
884 }
885 }
886
887 #ifdef TEST
888 assert(!isinf(res));
889 #endif
890
891 /* gaussian terms */
892 n = data->ng;
893 gt = &(data->gt[0]);
894 for(i=0; i<n; ++i){
895 #ifdef RESID_DEBUG
896 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);
897 #endif
898 double d1 = delta - gt->epsilon;
899 double t1 = tau - gt->gamma;
900 double e1 = -gt->alpha*SQ(d1) - gt->beta*SQ(t1);
901
902 double f1 = gt->t - 2*gt->beta*tau*(tau - gt->gamma);
903 double g1 = gt->d - 2*gt->alpha*delta*(delta - gt->epsilon);
904
905 sum = gt->n * f1 * pow(tau,gt->t-1) * g1 * pow(delta,gt->d-1) * exp(e1);
906
907 //fprintf(stderr,"sum = %f\n",sum);
908 res += sum;
909 #ifdef TEST
910 assert(!isinf(res));
911 #endif
912 ++gt;
913 }
914
915 /* critical terms */
916 n = data->nc;
917 ct = &(data->ct[0]);
918 for(i=0; i<n; ++i){
919 #ifdef RESID_DEBUG
920 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);
921 #endif
922 DEFINE_DELTA;
923 DEFINE_DPSIDDELTA;
924 DEFINE_DDELBDTAU;
925 DEFINE_DDELDDELTA;
926
927 double d2DELbddeldtau = -ct->A * ct->b * 2./ct->beta * pow(DELTA,ct->b-1)*d1*pow(d12,0.5/ct->beta-1) \
928 - 2. * theta * ct->b * (ct->b - 1) * pow(DELTA,ct->b-2.) * dDELddelta;
929
930 double d2psiddeldtau = 4. * ct->C*ct->D*d1*t1*psi;
931
932 DEFINE_DPSIDTAU;
933
934 sum = ct->n * (pow(DELTA, ct->b) * (dpsidtau + delta * d2psiddeldtau) \
935 + delta *dDELbdtau*dpsidtau \
936 + dDELbdtau*(psi+delta*dpsiddelta) \
937 + d2DELbddeldtau*delta*psi
938 );
939 res += sum;
940 ++ct;
941 }
942
943 #ifdef RESID_DEBUG
944 fprintf(stderr,"phir = %f\n",res);
945 #endif
946
947 #ifdef TEST
948 assert(!isnan(res));
949 assert(!isinf(res));
950 #endif
951 return res;
952 }
953
954 /**
955 Second derivative of helmholtz residual function with respect to
956 delta (twice).
957 */
958 double helm_resid_deldel(double tau,double delta,const HelmholtzData *data){
959 double sum = 0, res = 0;
960 double dell, ldell;
961 unsigned n, i;
962 const HelmholtzPowTerm *pt;
963 const HelmholtzGausTerm *gt;
964 const HelmholtzCritTerm *ct;
965
966 #ifdef RESID_DEBUG
967 fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
968 #endif
969
970 /* power terms */
971 n = data->np;
972 pt = &(data->pt[0]);
973 dell = ipow(delta,pt->l);
974 ldell = pt->l * dell;
975 unsigned oldl;
976 for(i=0; i<n; ++i){
977 double lpart = pt->l ? SQ(ldell) + ldell*(1. - 2*pt->d - pt->l) : 0;
978 sum += pt->a * pow(tau, pt->t) * ipow(delta, pt->d - 2) * (pt->d*(pt->d - 1) + lpart);
979 oldl = pt->l;
980 ++pt;
981 if(i+1==n || oldl != pt->l){
982 if(oldl == 0){
983 res += sum;
984 }else{
985 res += sum * exp(-dell);
986 }
987 sum = 0;
988 dell = ipow(delta,pt->l);
989 ldell = pt->l*dell;
990 }
991 }
992
993 /* gaussian terms */
994 n = data->ng;
995 //fprintf(stderr,"THERE ARE %d GAUSSIAN TERMS\n",n);
996 gt = &(data->gt[0]);
997 for(i=0; i<n; ++i){
998 double s1 = SQ(delta - gt->epsilon);
999 double f1 = gt->d*(gt->d - 1)
1000 + 2.*gt->alpha*delta * (
1001 delta * (2. * gt->alpha * s1 - 1)
1002 - 2. * gt->d * (delta - gt->epsilon)
1003 );
1004 res += gt->n * pow(tau,gt->t) * pow(delta, gt->d - 2.)
1005 * f1
1006 * exp(-(gt->alpha * s1 + gt->beta*SQ(tau-gt->gamma)));
1007 ++gt;
1008 }
1009
1010 /* critical terms */
1011 n = data->nc;
1012 ct = &(data->ct[0]);
1013 for(i=0; i<n; ++i){
1014 #ifdef RESID_DEBUG
1015 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);
1016 #endif
1017
1018 DEFINE_DELTA;
1019 DEFINE_DPSIDDELTA;
1020 DEFINE_DDELDDELTA;
1021 DEFINE_DDELBDDELTA;
1022
1023 double powd12bm1 = pow(d12,0.5/ct->beta-1.);
1024
1025 double d2psiddel2 = (2.*ct->C*d12 - 1.)*2.*ct->C*psi;
1026
1027 double d2DELddel2 = 1./d1*dDELddelta + d12*(
1028 4.*ct->B*ct->a*(ct->a-1.)*pow(d12,ct->a-2.) \
1029 + 2.*SQ(ct->A)*SQ(1./ct->beta)*SQ(powd12bm1) \
1030 + ct->A*theta*4./ct->beta*(0.5/ct->beta-1.)*powd12bm1/d12
1031 );
1032
1033 double d2DELbddel2 = ct->b * (pow(DELTA,ct->b - 1.)*d2DELddel2 + (ct->b-1.)*pow(DELTA,ct->b-2.)*SQ(dDELddelta));
1034
1035 sum = ct->n * (pow(DELTA,ct->b)*(2.*dpsiddelta + delta*d2psiddel2) + 2.*dDELbddelta*(psi+delta*dpsiddelta) + d2DELbddel2*delta*psi);
1036
1037 res += sum;
1038 ++ct;
1039 }
1040
1041 return res;
1042 }
1043
1044
1045
1046 /**
1047 Residual part of helmholtz function.
1048 */
1049 double helm_resid_tautau(double tau, double delta, const HelmholtzData *data){
1050 double dell,ldell, term, sum, res = 0;
1051 unsigned n, i;
1052 const HelmholtzPowTerm *pt;
1053 const HelmholtzGausTerm *gt;
1054 const HelmholtzCritTerm *ct;
1055
1056 n = data->np;
1057 pt = &(data->pt[0]);
1058
1059 #ifdef RESID_DEBUG
1060 fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
1061 #endif
1062
1063 /* power terms */
1064 sum = 0;
1065 dell = ipow(delta,pt->l);
1066 ldell = pt->l * dell;
1067 unsigned oldl;
1068 for(i=0; i<n; ++i){
1069 term = pt->a * pt->t * (pt->t - 1) * pow(tau, pt->t - 2) * ipow(delta, pt->d);
1070 sum += term;
1071 #ifdef RESID_DEBUG
1072 fprintf(stderr,"i = %d, a=%e, t=%f, d=%d, term = %f, sum = %f",i,pt->a,pt->t,pt->d,term,sum);
1073 if(pt->l==0){
1074 fprintf(stderr,",row=%e\n",term);
1075 }else{
1076 fprintf(stderr,",row=%e\n,",term*exp(-dell));
1077 }
1078 #endif
1079 oldl = pt->l;
1080 ++pt;
1081 if(i+1==n || oldl != pt->l){
1082 if(oldl == 0){
1083 #ifdef RESID_DEBUG
1084 fprintf(stderr,"linear ");
1085 #endif
1086 res += sum;
1087 }else{
1088 #ifdef RESID_DEBUG
1089 fprintf(stderr,"exp dell=%f, exp(-dell)=%f sum=%f: ",dell,exp(-dell),sum);
1090 #endif
1091 res += sum * exp(-dell);
1092 }
1093 #ifdef RESID_DEBUG
1094 fprintf(stderr,"i = %d, res = %f\n",i,res);
1095 #endif
1096 sum = 0;
1097 dell = ipow(delta,pt->l);
1098 ldell = pt->l*dell;
1099 }
1100 }
1101
1102 /* gaussian terms */
1103 n = data->ng;
1104 //fprintf(stderr,"THERE ARE %d GAUSSIAN TERMS\n",n);
1105 gt = &(data->gt[0]);
1106 for(i=0; i<n; ++i){
1107 #ifdef RESID_DEBUG
1108 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);
1109 #endif
1110 double d1 = delta - gt->epsilon;
1111 double t1 = tau - gt->gamma;
1112 double f1 = gt->t*(gt->t - 1) + 4. * gt->beta * tau * (tau * (gt->beta*SQ(t1) - 0.5) - t1*gt->t);
1113 double e1 = -gt->alpha*SQ(d1) - gt->beta*SQ(t1);
1114 sum = gt->n * f1 * pow(tau,gt->t - 2) * pow(delta,gt->d) * exp(e1);
1115 //fprintf(stderr,"sum = %f\n",sum);
1116 res += sum;
1117 ++gt;
1118 }
1119
1120 /* critical terms */
1121 n = data->nc;
1122 ct = &(data->ct[0]);
1123 for(i=0; i<n; ++i){
1124 #ifdef RESID_DEBUG
1125 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);
1126 #endif
1127 DEFINE_DELTA;
1128 DEFINE_DDELBDTAU;
1129 DEFINE_DPSIDTAU;
1130
1131 double d2DELbdtau2 = 2. * ct->b * pow(DELTA, ct->b - 1) + 4. * SQ(theta) * ct->b * (ct->b - 1) * pow(DELTA, ct->b - 2);
1132
1133 double d2psidtau2 = 2. * ct->D * psi * (2. * ct->D * SQ(t1) -1.);
1134
1135 sum = ct->n * delta * (d2DELbdtau2 * psi + 2 * dDELbdtau*dpsidtau + pow(DELTA, ct->b) * d2psidtau2);
1136 res += sum;
1137 ++ct;
1138 }
1139
1140 #ifdef RESID_DEBUG
1141 fprintf(stderr,"phir_tautau = %f\n",res);
1142 #endif
1143 return res;
1144 }
1145

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