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

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