/[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 1985 - (show annotations) (download) (as text)
Tue Feb 3 01:15:17 2009 UTC (13 years, 5 months ago) by jpye
File MIME type: text/x-csrc
File size: 20373 byte(s)
Added critical term calculation in helm_resid.
Added gaussian and critical term values for water from IAPWS95.
Added some test data from IAPWS95 to water.c tests.
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 #define SQ(X) ((X)*(X))
39
40 /* forward decls */
41
42 #include "helmholtz_impl.h"
43
44 /**
45 Function to calculate pressure from Helmholtz free energy EOS, given temperature
46 and mass density.
47
48 @param T temperature in K
49 @param rho mass density in kg/m続
50 @return pressure in Pa???
51 */
52 double helmholtz_p(double T, double rho, const HelmholtzData *data){
53
54 double tau = data->T_star / T;
55 double delta = rho / data->rho_star;
56
57 #ifdef TEST
58 assert(data->rho_star!=0);
59 assert(T!=0);
60 assert(!isnan(tau));
61 assert(!isnan(delta));
62 assert(!isnan(data->R));
63
64 //fprintf(stderr,"p calc: T = %f\n",T);
65 //fprintf(stderr,"p calc: tau = %f\n",tau);
66 //fprintf(stderr,"p calc: rho = %f\n",rho);
67 //fprintf(stderr,"p calc: delta = %f\n",delta);
68 //fprintf(stderr,"p calc: R*T*rho = %f\n",data->R * T * rho);
69
70 //fprintf(stderr,"T = %f\n", T);
71 //fprintf(stderr,"rhob = %f, rhob* = %f, delta = %f\n", rho/data->M, data->rho_star/data->M, delta);
72 #endif
73
74 return data->R * T * rho * (1 + delta * helm_resid_del(tau,delta,data));
75 }
76
77 /**
78 Function to calculate internal energy from Helmholtz free energy EOS, given
79 temperature and mass density.
80
81 @param T temperature in K
82 @param rho mass density in kg/m続
83 @return internal energy in ???
84 */
85 double helmholtz_u(double T, double rho, const HelmholtzData *data){
86
87 double tau = data->T_star / T;
88 double delta = rho / data->rho_star;
89
90 #ifdef TEST
91 assert(data->rho_star!=0);
92 assert(T!=0);
93 assert(!isnan(tau));
94 assert(!isnan(delta));
95 assert(!isnan(data->R));
96 #endif
97
98 #if 0
99 fprintf(stderr,"ideal_tau = %f\n",helm_ideal_tau(tau,delta,data->ideal));
100 fprintf(stderr,"resid_tau = %f\n",helm_resid_tau(tau,delta,data));
101 fprintf(stderr,"R T = %f\n",data->R * data->T_star);
102 #endif
103
104 return data->R * data->T_star * (helm_ideal_tau(tau,delta,data->ideal) + helm_resid_tau(tau,delta,data));
105 }
106
107 /**
108 Function to calculate enthalpy from Helmholtz free energy EOS, given
109 temperature and mass density.
110
111 @param T temperature in K
112 @param rho mass density in kg/m続
113 @return enthalpy in J/kg
114 */
115 double helmholtz_h(double T, double rho, const HelmholtzData *data){
116
117 double tau = data->T_star / T;
118 double delta = rho / data->rho_star;
119
120 #ifdef TEST
121 assert(data->rho_star!=0);
122 assert(T!=0);
123 assert(!isnan(tau));
124 assert(!isnan(delta));
125 assert(!isnan(data->R));
126 #endif
127
128 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));
129 }
130
131 /**
132 Function to calculate entropy from Helmholtz free energy EOS, given
133 temperature and mass density.
134
135 @param T temperature in K
136 @param rho mass density in kg/m続
137 @return entropy in J/kgK
138 */
139 double helmholtz_s(double T, double rho, const HelmholtzData *data){
140
141 double tau = data->T_star / T;
142 double delta = rho / data->rho_star;
143
144 #ifdef ENTROPY_DEBUG
145 assert(data->rho_star!=0);
146 assert(T!=0);
147 assert(!isnan(tau));
148 assert(!isnan(delta));
149 assert(!isnan(data->R));
150
151 fprintf(stderr,"helm_ideal_tau = %f\n",helm_ideal_tau(tau,delta,data->ideal));
152 fprintf(stderr,"helm_resid_tau = %f\n",helm_resid_tau(tau,delta,data));
153 fprintf(stderr,"helm_ideal = %f\n",helm_ideal(tau,delta,data->ideal));
154 fprintf(stderr,"helm_resid = %f\n",helm_resid(tau,delta,data));
155 #endif
156 return data->R * (
157 tau * (helm_ideal_tau(tau,delta,data->ideal) + helm_resid_tau(tau,delta,data))
158 - (helm_ideal(tau,delta,data->ideal) + helm_resid(tau,delta,data))
159 );
160 }
161
162 /**
163 Function to calculate Helmholtz energy from the Helmholtz free energy EOS,
164 given temperature and mass density.
165
166 @param T temperature in K
167 @param rho mass density in kg/m続
168 @return Helmholtz energy 'a', in J/kg
169 */
170 double helmholtz_a(double T, double rho, const HelmholtzData *data){
171
172 double tau = data->T_star / T;
173 double delta = rho / data->rho_star;
174
175 #ifdef TEST
176 assert(data->rho_star!=0);
177 assert(T!=0);
178 assert(!isnan(tau));
179 assert(!isnan(delta));
180 assert(!isnan(data->R));
181 #endif
182
183 #ifdef HELMHOLTZ_DEBUG
184 fprintf(stderr,"helmholtz_a: T = %f, rho = %f\n",T,rho);
185 fprintf(stderr,"multiplying by RT = %f\n",data->R*T);
186 #endif
187
188 return data->R * T * (helm_ideal(tau,delta,data->ideal) + helm_resid(tau,delta,data));
189 }
190
191
192 /**
193 Calculation zero-pressure specific heat capacity
194 */
195 double helmholtz_cp0(double T, const HelmholtzData *data){
196 double val = helm_cp0(T,data->ideal);
197 #if 0
198 fprintf(stderr,"val = %f\n",val);
199 #endif
200 return val;
201 }
202
203 /**
204 Calculate partial derivative of p with respect to T, with rho constant
205 */
206 double helmholtz_dpdT_rho(double T, double rho, const HelmholtzData *data){
207 double tau = data->T_star / T;
208 double delta = rho / data->rho_star;
209
210 double phir_del = helm_resid_del(tau,delta,data);
211 double phir_deltau = helm_resid_deltau(tau,delta,data);
212 #ifdef TEST
213 assert(!isinf(phir_del));
214 assert(!isinf(phir_deltau));
215 assert(!isnan(phir_del));
216 assert(!isnan(phir_deltau));
217 assert(!isnan(data->R));
218 assert(!isnan(rho));
219 assert(!isnan(tau));
220 #endif
221
222 double res = data->R * rho * (1 + delta*phir_del - delta*tau*phir_deltau);
223
224 #ifdef TEST
225 assert(!isnan(res));
226 assert(!isinf(res));
227 #endif
228 return res;
229 }
230
231 /**
232 Calculate partial derivative of p with respect to rho, with T constant
233 */
234 double helmholtz_dpdrho_T(double T, double rho, const HelmholtzData *data){
235 double tau = data->T_star / T;
236 double delta = rho / data->rho_star;
237
238 double phir_del = helm_resid_del(tau,delta,data);
239 double phir_deldel = helm_resid_deldel(tau,delta,data);
240 #ifdef TEST
241 assert(!isinf(phir_del));
242 assert(!isinf(phir_deldel));
243 #endif
244 return data->R * T * (1 + 2*delta*phir_del + delta*delta* phir_deldel);
245 }
246
247 /**
248 Calculate partial derivative of h with respect to T, with rho constant
249 */
250 double helmholtz_dhdT_rho(double T, double rho, const HelmholtzData *data){
251 double tau = data->T_star / T;
252 double delta = rho / data->rho_star;
253
254 double phir_del = helm_resid_del(tau,delta,data);
255 double phir_deltau = helm_resid_deltau(tau,delta,data);
256 double phir_tautau = helm_resid_tautau(tau,delta,data);
257 double phi0_tautau = helm_ideal_tautau(tau,data->ideal);
258
259 //fprintf(stderr,"phir_del = %f, phir_deltau = %f, phir_tautau = %f, phi0_tautau = %f\n",phir_del,phir_deltau,phir_tautau,phi0_tautau);
260
261 //return (helmholtz_h(T+0.01,rho,data) - helmholtz_h(T,rho,data)) / 0.01;
262 return data->R * (1. + delta*phir_del - tau*tau*(phi0_tautau + phir_tautau) - delta*tau*phir_deltau);
263 }
264
265 /**
266 Calculate partial derivative of h with respect to rho, with T constant
267 */
268 double helmholtz_dhdrho_T(double T, double rho, const HelmholtzData *data){
269 double tau = data->T_star / T;
270 double delta = rho / data->rho_star;
271
272 double phir_del = helm_resid_del(tau,delta,data);
273 double phir_deltau = helm_resid_deltau(tau,delta,data);
274 double phir_deldel = helm_resid_deldel(tau,delta,data);
275
276 return data->R * T / rho * (tau*delta*(0 + phir_deltau) + delta * phir_del + SQ(delta)*phir_deldel);
277 }
278
279
280 /**
281 Calculate partial derivative of u with respect to T, with rho constant
282 */
283 double helmholtz_dudT_rho(double T, double rho, const HelmholtzData *data){
284 double tau = data->T_star / T;
285 double delta = rho / data->rho_star;
286
287 double phir_tautau = helm_resid_tautau(tau,delta,data);
288 double phi0_tautau = helm_ideal_tautau(tau,data->ideal);
289
290 return -data->R * SQ(tau) * (phi0_tautau + phir_tautau);
291 }
292
293 /**
294 Calculate partial derivative of u with respect to rho, with T constant
295 */
296 double helmholtz_dudrho_T(double T, double rho, const HelmholtzData *data){
297 double tau = data->T_star / T;
298 double delta = rho / data->rho_star;
299
300 double phir_deltau = helm_resid_deltau(tau,delta,data);
301
302 return data->R * T / rho * (tau * delta * phir_deltau);
303 }
304
305 /*---------------------------------------------
306 UTILITY FUNCTION(S)
307 */
308
309 /* ipow: public domain by Mark Stephen with suggestions by Keiichi Nakasato */
310 static double ipow(double x, int n){
311 double t = 1.0;
312
313 if(!n)return 1.0; /* At the top. x^0 = 1 */
314
315 if (n < 0){
316 n = -n;
317 x = 1.0/x; /* error if x == 0. Good */
318 } /* ZTC/SC returns inf, which is even better */
319
320 if (x == 0.0)return 0.0;
321
322 do{
323 if(n & 1)t *= x;
324 n /= 2; /* KN prefers if (n/=2) x*=x; This avoids an */
325 x *= x; /* unnecessary but benign multiplication on */
326 }while(n); /* the last pass, but the comparison is always
327 true _except_ on the last pass. */
328
329 return t;
330 }
331
332 //#define RESID_DEBUG
333
334 /**
335 Residual part of helmholtz function.
336 */
337 double helm_resid(double tau, double delta, const HelmholtzData *data){
338 double dell,ldell, term, sum, res = 0;
339 unsigned n, i;
340 const HelmholtzPowTerm *pt;
341 const HelmholtzGausTerm *gt;
342 const HelmholtzCritTerm *ct;
343
344 n = data->np;
345 pt = &(data->pt[0]);
346
347 #ifdef RESID_DEBUG
348 fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
349 #endif
350
351 /* power terms */
352 sum = 0;
353 dell = ipow(delta,pt->l);
354 ldell = pt->l * dell;
355 unsigned oldl;
356 for(i=0; i<n; ++i){
357 term = pt->a * pow(tau, pt->t) * ipow(delta, pt->d);
358 sum += term;
359 #ifdef RESID_DEBUG
360 fprintf(stderr,"i = %d, a=%e, t=%f, d=%d, term = %f, sum = %f",i,pt->a,pt->t,pt->d,term,sum);
361 if(pt->l==0){
362 fprintf(stderr,",row=%e\n",term);
363 }else{
364 fprintf(stderr,",row=%e\n,",term*exp(-dell));
365 }
366 #endif
367 oldl = pt->l;
368 ++pt;
369 if(i+1==n || oldl != pt->l){
370 if(oldl == 0){
371 #ifdef RESID_DEBUG
372 fprintf(stderr,"linear ");
373 #endif
374 res += sum;
375 }else{
376 #ifdef RESID_DEBUG
377 fprintf(stderr,"exp dell=%f, exp(-dell)=%f sum=%f: ",dell,exp(-dell),sum);
378 #endif
379 res += sum * exp(-dell);
380 }
381 #ifdef RESID_DEBUG
382 fprintf(stderr,"i = %d, res = %f\n",i,res);
383 #endif
384 sum = 0;
385 dell = ipow(delta,pt->l);
386 ldell = pt->l*dell;
387 }
388 }
389
390 /* gaussian terms */
391 n = data->ng;
392 //fprintf(stderr,"THERE ARE %d GAUSSIAN TERMS\n",n);
393 gt = &(data->gt[0]);
394 for(i=0; i<n; ++i){
395 #ifdef RESID_DEBUG
396 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);
397 #endif
398 double d1 = delta - gt->epsilon;
399 double t1 = tau - gt->gamma;
400 double e1 = -gt->alpha*d1*d1 - gt->beta*t1*t1;
401 sum = gt->n * pow(tau,gt->t) * pow(delta,gt->d) * exp(e1);
402 //fprintf(stderr,"sum = %f\n",sum);
403 res += sum;
404 ++gt;
405 }
406
407 /* critical terms */
408 n = data->nc;
409 ct = &(data->ct[0]);
410 for(i=0; i<n; ++i){
411 #ifdef RESID_DEBUG
412 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);
413 #endif
414 double d1 = delta - 1.;
415 double t1 = tau - 1.;
416 double theta = (1. - tau) + ct->A * pow(d1*d1, 0.5/ct->beta);
417 double psi = exp(-ct->A*d1*d1 - ct->D*t1*t1);
418 double DELTA = theta*theta + ct->B* pow(d1*d1, ct->a);
419 sum = ct->n * pow(DELTA, ct->b) * delta * psi;
420 res += sum;
421 ++ct;
422 }
423
424 #ifdef RESID_DEBUG
425 fprintf(stderr,"CALCULATED RESULT FOR phir = %f\n",res);
426 #endif
427 return res;
428 }
429
430 /*=================== FIRST DERIVATIVES =======================*/
431
432 /**
433 Derivative of the helmholtz residual function with respect to
434 delta.
435 */
436 double helm_resid_del(double tau,double delta, const HelmholtzData *data){
437 double sum = 0, res = 0;
438 double dell, ldell;
439 unsigned n, i;
440 const HelmholtzPowTerm *pt;
441 const HelmholtzGausTerm *gt;
442
443
444 #ifdef RESID_DEBUG
445 fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
446 #endif
447
448 /* power terms */
449 n = data->np;
450 pt = &(data->pt[0]);
451 dell = ipow(delta,pt->l);
452 ldell = pt->l * dell;
453 unsigned oldl;
454 for(i=0; i<n; ++i){
455 sum += pt->a * pow(tau, pt->t) * ipow(delta, pt->d - 1) * (pt->d - ldell);
456 oldl = pt->l;
457 ++pt;
458 if(i+1==n || oldl != pt->l){
459 if(oldl == 0){
460 res += sum;
461 }else{
462 res += sum * exp(-dell);
463 }
464 sum = 0;
465 dell = ipow(delta,pt->l);
466 ldell = pt->l*dell;
467 }
468 }
469
470 /* gaussian terms */
471 n = data->ng;
472 //fprintf(stderr,"THERE ARE %d GAUSSIAN TERMS\n",n);
473 gt = &(data->gt[0]);
474 for(i=0; i<n; ++i){
475 #ifdef RESID_DEBUG
476 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);
477 #endif
478 double val2;
479 val2 = - gt->n * pow(tau,gt->t) * pow(delta, -1. + gt->d)
480 * (2. * gt->alpha * delta * (delta - gt->epsilon) - gt->d)
481 * exp(-(gt->alpha * SQ(delta-gt->epsilon) + gt->beta*SQ(tau-gt->gamma)));
482 res += val2;
483 #ifdef RESID_DEBUG
484 fprintf(stderr,"val2 = %f --> res = %f\n",val2,res);
485 #endif
486 ++gt;
487 }
488
489 /* FIXME add critical terms calculation */
490
491 return res;
492 }
493
494 /**
495 Derivative of the helmholtz residual function with respect to
496 tau.
497 */
498 double helm_resid_tau(double tau,double delta,const HelmholtzData *data){
499
500 double sum;
501 double res = 0;
502 double delX;
503 unsigned l;
504 unsigned n, i;
505 const HelmholtzPowTerm *pt;
506 const HelmholtzGausTerm *gt;
507
508 n = data->np;
509 pt = &(data->pt[0]);
510
511 delX = 1;
512
513 l = 0;
514 sum = 0;
515 for(i=0; i<n; ++i){
516 if(pt->t){
517 //fprintf(stderr,"i = %d, a = %e, t = %f, d = %d, l = %d\n",i+1, pt->a, pt->t, pt->d, pt->l);
518 sum += pt->a * pow(tau, pt->t - 1) * ipow(delta, pt->d) * pt->t;
519 }
520 ++pt;
521 //fprintf(stderr,"l = %d\n",l);
522 if(i+1==n || l != pt->l){
523 if(l==0){
524 //fprintf(stderr,"Adding non-exp term\n");
525 res += sum;
526 }else{
527 //fprintf(stderr,"Adding exp term with l = %d, delX = %e\n",l,delX);
528 res += sum * exp(-delX);
529 }
530 /* set l to new value */
531 if(i+1!=n){
532 l = pt->l;
533 //fprintf(stderr,"New l = %d\n",l);
534 delX = ipow(delta,l);
535 sum = 0;
536 }
537 }
538 }
539
540 //#define RESID_DEBUG
541 /* gaussian terms */
542 n = data->ng;
543 gt = &(data->gt[0]);
544 for(i=0; i<n; ++i){
545 #ifdef RESID_DEBUG
546 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);
547 #endif
548
549 double val2;
550 val2 = -gt->n * pow(tau,gt->t - 1.) * pow(delta, gt->d)
551 * (2. * gt->beta * tau * (tau - gt->gamma) - gt->t)
552 * exp(-(gt->alpha * SQ(delta-gt->epsilon) + gt->beta*SQ(tau-gt->gamma)));
553 res += val2;
554 #ifdef RESID_DEBUG
555 fprintf(stderr,"res = %f\n",res);
556 #endif
557
558 ++gt;
559 }
560
561 /* FIXME add critical terms calculation */
562
563 return res;
564 }
565
566
567 /*=================== SECOND DERIVATIVES =======================*/
568
569 /**
570 Mixed derivative of the helmholtz residual function with respect to
571 delta and tau.
572 */
573 double helm_resid_deltau(double tau,double delta,const HelmholtzData *data){
574 double dell,ldell, term, sum = 0, res = 0;
575 unsigned n, i;
576 const HelmholtzPowTerm *pt;
577 const HelmholtzGausTerm *gt;
578
579 /* power terms */
580 n = data->np;
581 pt = &(data->pt[0]);
582 dell = ipow(delta,pt->l);
583 ldell = pt->l * dell;
584 unsigned oldl;
585 for(i=0; i<n; ++i){
586 sum += pt->a * pt->t * pow(tau, pt->t - 1) * ipow(delta, pt->d - 1) * (pt->d - ldell);
587 oldl = pt->l;
588 ++pt;
589 if(i+1==n || oldl != pt->l){
590 if(oldl == 0){
591 res += sum;
592 }else{
593 res += sum * exp(-dell);
594 }
595 sum = 0;
596 dell = ipow(delta,pt->l);
597 ldell = pt->l*dell;
598 }
599 }
600
601 #ifdef TEST
602 assert(!isinf(res));
603 #endif
604
605 /* gaussian terms */
606 n = data->ng;
607 gt = &(data->gt[0]);
608 for(i=0; i<n; ++i){
609 #ifdef RESID_DEBUG
610 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);
611 #endif
612 double d1 = delta - gt->epsilon;
613 double t1 = tau - gt->gamma;
614 double e1 = -gt->alpha*SQ(d1) - gt->beta*SQ(t1);
615
616 double f1 = gt->t - 2*gt->beta*tau*(tau - gt->gamma);
617 double g1 = gt->d - 2*gt->alpha*delta*(delta - gt->epsilon);
618
619 sum = gt->n * f1 * pow(tau,gt->t-1) * g1 * pow(delta,gt->d-1) * exp(e1);
620
621 //fprintf(stderr,"sum = %f\n",sum);
622 res += sum;
623 #ifdef TEST
624 assert(!isinf(res));
625 #endif
626 ++gt;
627 }
628
629 /* FIXME add critical terms calculation */
630
631 #ifdef RESID_DEBUG
632 fprintf(stderr,"phir = %f\n",res);
633 #endif
634
635 #ifdef TEST
636 assert(!isnan(res));
637 assert(!isinf(res));
638 #endif
639 return res;
640 }
641
642 /**
643 Second derivative of helmholtz residual function with respect to
644 delta (twice).
645
646 FIXME this function is WRONG.
647 */
648 double helm_resid_deldel(double tau,double delta,const HelmholtzData *data){
649 double sum = 0, res = 0;
650 double dell, ldell;
651 unsigned n, i;
652 const HelmholtzPowTerm *pt;
653 const HelmholtzGausTerm *gt;
654
655
656 #ifdef RESID_DEBUG
657 fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
658 #endif
659
660 /* power terms */
661 n = data->np;
662 pt = &(data->pt[0]);
663 dell = ipow(delta,pt->l);
664 ldell = pt->l * dell;
665 unsigned oldl;
666 for(i=0; i<n; ++i){
667 double lpart = pt->l ? SQ(ldell) + ldell*(1. - 2*pt->d - pt->l) : 0;
668 sum += pt->a * pow(tau, pt->t) * ipow(delta, pt->d - 2) * (pt->d*(pt->d - 1) + lpart);
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 //fprintf(stderr,"THERE ARE %d GAUSSIAN TERMS\n",n);
686 gt = &(data->gt[0]);
687 for(i=0; i<n; ++i){
688 double s1 = SQ(delta - gt->epsilon);
689 double f1 = gt->d*(gt->d - 1)
690 + 2.*gt->alpha*delta * (
691 delta * (2. * gt->alpha * s1 - 1)
692 - 2. * gt->d * (delta - gt->epsilon)
693 );
694 res += gt->n * pow(tau,gt->t) * pow(delta, gt->d - 2.)
695 * f1
696 * exp(-(gt->alpha * s1 + gt->beta*SQ(tau-gt->gamma)));
697 ++gt;
698 }
699
700 /* FIXME add critical terms calculation */
701
702 return res;
703 }
704
705
706
707 /**
708 Residual part of helmholtz function.
709 */
710 double helm_resid_tautau(double tau, double delta, const HelmholtzData *data){
711 double dell,ldell, term, sum, res = 0;
712 unsigned n, i;
713 const HelmholtzPowTerm *pt;
714 const HelmholtzGausTerm *gt;
715
716 n = data->np;
717 pt = &(data->pt[0]);
718
719 #ifdef RESID_DEBUG
720 fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
721 #endif
722
723 /* power terms */
724 sum = 0;
725 dell = ipow(delta,pt->l);
726 ldell = pt->l * dell;
727 unsigned oldl;
728 for(i=0; i<n; ++i){
729 term = pt->a * pt->t * (pt->t - 1) * pow(tau, pt->t - 2) * ipow(delta, pt->d);
730 sum += term;
731 #ifdef RESID_DEBUG
732 fprintf(stderr,"i = %d, a=%e, t=%f, d=%d, term = %f, sum = %f",i,pt->a,pt->t,pt->d,term,sum);
733 if(pt->l==0){
734 fprintf(stderr,",row=%e\n",term);
735 }else{
736 fprintf(stderr,",row=%e\n,",term*exp(-dell));
737 }
738 #endif
739 oldl = pt->l;
740 ++pt;
741 if(i+1==n || oldl != pt->l){
742 if(oldl == 0){
743 #ifdef RESID_DEBUG
744 fprintf(stderr,"linear ");
745 #endif
746 res += sum;
747 }else{
748 #ifdef RESID_DEBUG
749 fprintf(stderr,"exp dell=%f, exp(-dell)=%f sum=%f: ",dell,exp(-dell),sum);
750 #endif
751 res += sum * exp(-dell);
752 }
753 #ifdef RESID_DEBUG
754 fprintf(stderr,"i = %d, res = %f\n",i,res);
755 #endif
756 sum = 0;
757 dell = ipow(delta,pt->l);
758 ldell = pt->l*dell;
759 }
760 }
761
762 /* gaussian terms */
763 n = data->ng;
764 //fprintf(stderr,"THERE ARE %d GAUSSIAN TERMS\n",n);
765 gt = &(data->gt[0]);
766 for(i=0; i<n; ++i){
767 #ifdef RESID_DEBUG
768 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);
769 #endif
770 double d1 = delta - gt->epsilon;
771 double t1 = tau - gt->gamma;
772 double f1 = gt->t*(gt->t - 1) + 4. * gt->beta * tau * (tau * (gt->beta*SQ(t1) - 0.5) - t1*gt->t);
773 double e1 = -gt->alpha*SQ(d1) - gt->beta*SQ(t1);
774 sum = gt->n * f1 * pow(tau,gt->t - 2) * pow(delta,gt->d) * exp(e1);
775 //fprintf(stderr,"sum = %f\n",sum);
776 res += sum;
777 ++gt;
778 }
779
780 /* FIXME add critical terms calculation */
781
782 #ifdef RESID_DEBUG
783 fprintf(stderr,"phir_tautau = %f\n",res);
784 #endif
785 return res;
786 }
787

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