/[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 1986 - (show annotations) (download) (as text)
Tue Feb 3 05:19:15 2009 UTC (13 years, 5 months ago) by jpye
File MIME type: text/x-csrc
File size: 21149 byte(s)
Fixed helm_resid and helm_resid_del with critical terms.
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->C*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 const HelmholtzCritTerm *ct;
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 gt = &(data->gt[0]);
473 for(i=0; i<n; ++i){
474 #ifdef RESID_DEBUG
475 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);
476 #endif
477 sum = - gt->n * pow(tau,gt->t) * pow(delta, -1. + gt->d)
478 * (2. * gt->alpha * delta * (delta - gt->epsilon) - gt->d)
479 * exp(-(gt->alpha * SQ(delta-gt->epsilon) + gt->beta*SQ(tau-gt->gamma)));
480 res += sum;
481 #ifdef RESID_DEBUG
482 fprintf(stderr,"sum = %f --> res = %f\n",sum,res);
483 #endif
484 ++gt;
485 }
486
487 /* critical terms */
488 n = data->nc;
489 ct = &(data->ct[0]);
490 for(i=0; i<n; ++i){
491 #ifdef RESID_DEBUG
492 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);
493 #endif
494 double d1 = delta - 1.;
495 double t1 = tau - 1.;
496 double theta = (1. - tau) + ct->A * pow(d1*d1, 0.5/ct->beta);
497 double psi = exp(-ct->C*d1*d1 - ct->D*t1*t1);
498 double DELTA = theta*theta + ct->B* pow(d1*d1, ct->a);
499
500 double dpsiddelta = -2. * ct->C * d1 * psi;
501
502 double dDELddelta = d1 * (ct->A * theta * 2./ct->beta * pow(d1*d1, 0.5/ct->beta - 1) + 2* ct->B * ct->a * pow(d1*d1, ct->a - 1));
503
504 double dDELbddelta = ct->b * pow(DELTA,ct->b - 1.) * dDELddelta;
505
506 sum = ct->n * (pow(DELTA, ct->b) * (psi + delta * dpsiddelta) + dDELbddelta * delta * psi);
507 res += sum;
508 ++ct;
509 }
510
511
512 return res;
513 }
514
515 /**
516 Derivative of the helmholtz residual function with respect to
517 tau.
518 */
519 double helm_resid_tau(double tau,double delta,const HelmholtzData *data){
520
521 double sum;
522 double res = 0;
523 double delX;
524 unsigned l;
525 unsigned n, i;
526 const HelmholtzPowTerm *pt;
527 const HelmholtzGausTerm *gt;
528
529 n = data->np;
530 pt = &(data->pt[0]);
531
532 delX = 1;
533
534 l = 0;
535 sum = 0;
536 for(i=0; i<n; ++i){
537 if(pt->t){
538 //fprintf(stderr,"i = %d, a = %e, t = %f, d = %d, l = %d\n",i+1, pt->a, pt->t, pt->d, pt->l);
539 sum += pt->a * pow(tau, pt->t - 1) * ipow(delta, pt->d) * pt->t;
540 }
541 ++pt;
542 //fprintf(stderr,"l = %d\n",l);
543 if(i+1==n || l != pt->l){
544 if(l==0){
545 //fprintf(stderr,"Adding non-exp term\n");
546 res += sum;
547 }else{
548 //fprintf(stderr,"Adding exp term with l = %d, delX = %e\n",l,delX);
549 res += sum * exp(-delX);
550 }
551 /* set l to new value */
552 if(i+1!=n){
553 l = pt->l;
554 //fprintf(stderr,"New l = %d\n",l);
555 delX = ipow(delta,l);
556 sum = 0;
557 }
558 }
559 }
560
561 //#define RESID_DEBUG
562 /* gaussian terms */
563 n = data->ng;
564 gt = &(data->gt[0]);
565 for(i=0; i<n; ++i){
566 #ifdef RESID_DEBUG
567 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);
568 #endif
569
570 double val2;
571 val2 = -gt->n * pow(tau,gt->t - 1.) * pow(delta, gt->d)
572 * (2. * gt->beta * tau * (tau - gt->gamma) - gt->t)
573 * exp(-(gt->alpha * SQ(delta-gt->epsilon) + gt->beta*SQ(tau-gt->gamma)));
574 res += val2;
575 #ifdef RESID_DEBUG
576 fprintf(stderr,"res = %f\n",res);
577 #endif
578
579 ++gt;
580 }
581
582 /* FIXME add critical terms calculation */
583
584 return res;
585 }
586
587
588 /*=================== SECOND DERIVATIVES =======================*/
589
590 /**
591 Mixed derivative of the helmholtz residual function with respect to
592 delta and tau.
593 */
594 double helm_resid_deltau(double tau,double delta,const HelmholtzData *data){
595 double dell,ldell, term, sum = 0, res = 0;
596 unsigned n, i;
597 const HelmholtzPowTerm *pt;
598 const HelmholtzGausTerm *gt;
599
600 /* power terms */
601 n = data->np;
602 pt = &(data->pt[0]);
603 dell = ipow(delta,pt->l);
604 ldell = pt->l * dell;
605 unsigned oldl;
606 for(i=0; i<n; ++i){
607 sum += pt->a * pt->t * pow(tau, pt->t - 1) * ipow(delta, pt->d - 1) * (pt->d - ldell);
608 oldl = pt->l;
609 ++pt;
610 if(i+1==n || oldl != pt->l){
611 if(oldl == 0){
612 res += sum;
613 }else{
614 res += sum * exp(-dell);
615 }
616 sum = 0;
617 dell = ipow(delta,pt->l);
618 ldell = pt->l*dell;
619 }
620 }
621
622 #ifdef TEST
623 assert(!isinf(res));
624 #endif
625
626 /* gaussian terms */
627 n = data->ng;
628 gt = &(data->gt[0]);
629 for(i=0; i<n; ++i){
630 #ifdef RESID_DEBUG
631 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);
632 #endif
633 double d1 = delta - gt->epsilon;
634 double t1 = tau - gt->gamma;
635 double e1 = -gt->alpha*SQ(d1) - gt->beta*SQ(t1);
636
637 double f1 = gt->t - 2*gt->beta*tau*(tau - gt->gamma);
638 double g1 = gt->d - 2*gt->alpha*delta*(delta - gt->epsilon);
639
640 sum = gt->n * f1 * pow(tau,gt->t-1) * g1 * pow(delta,gt->d-1) * exp(e1);
641
642 //fprintf(stderr,"sum = %f\n",sum);
643 res += sum;
644 #ifdef TEST
645 assert(!isinf(res));
646 #endif
647 ++gt;
648 }
649
650 /* FIXME add critical terms calculation */
651
652 #ifdef RESID_DEBUG
653 fprintf(stderr,"phir = %f\n",res);
654 #endif
655
656 #ifdef TEST
657 assert(!isnan(res));
658 assert(!isinf(res));
659 #endif
660 return res;
661 }
662
663 /**
664 Second derivative of helmholtz residual function with respect to
665 delta (twice).
666
667 FIXME this function is WRONG.
668 */
669 double helm_resid_deldel(double tau,double delta,const HelmholtzData *data){
670 double sum = 0, res = 0;
671 double dell, ldell;
672 unsigned n, i;
673 const HelmholtzPowTerm *pt;
674 const HelmholtzGausTerm *gt;
675
676
677 #ifdef RESID_DEBUG
678 fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
679 #endif
680
681 /* power terms */
682 n = data->np;
683 pt = &(data->pt[0]);
684 dell = ipow(delta,pt->l);
685 ldell = pt->l * dell;
686 unsigned oldl;
687 for(i=0; i<n; ++i){
688 double lpart = pt->l ? SQ(ldell) + ldell*(1. - 2*pt->d - pt->l) : 0;
689 sum += pt->a * pow(tau, pt->t) * ipow(delta, pt->d - 2) * (pt->d*(pt->d - 1) + lpart);
690 oldl = pt->l;
691 ++pt;
692 if(i+1==n || oldl != pt->l){
693 if(oldl == 0){
694 res += sum;
695 }else{
696 res += sum * exp(-dell);
697 }
698 sum = 0;
699 dell = ipow(delta,pt->l);
700 ldell = pt->l*dell;
701 }
702 }
703
704 /* gaussian terms */
705 n = data->ng;
706 //fprintf(stderr,"THERE ARE %d GAUSSIAN TERMS\n",n);
707 gt = &(data->gt[0]);
708 for(i=0; i<n; ++i){
709 double s1 = SQ(delta - gt->epsilon);
710 double f1 = gt->d*(gt->d - 1)
711 + 2.*gt->alpha*delta * (
712 delta * (2. * gt->alpha * s1 - 1)
713 - 2. * gt->d * (delta - gt->epsilon)
714 );
715 res += gt->n * pow(tau,gt->t) * pow(delta, gt->d - 2.)
716 * f1
717 * exp(-(gt->alpha * s1 + gt->beta*SQ(tau-gt->gamma)));
718 ++gt;
719 }
720
721 /* FIXME add critical terms calculation */
722
723 return res;
724 }
725
726
727
728 /**
729 Residual part of helmholtz function.
730 */
731 double helm_resid_tautau(double tau, double delta, const HelmholtzData *data){
732 double dell,ldell, term, sum, res = 0;
733 unsigned n, i;
734 const HelmholtzPowTerm *pt;
735 const HelmholtzGausTerm *gt;
736
737 n = data->np;
738 pt = &(data->pt[0]);
739
740 #ifdef RESID_DEBUG
741 fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
742 #endif
743
744 /* power terms */
745 sum = 0;
746 dell = ipow(delta,pt->l);
747 ldell = pt->l * dell;
748 unsigned oldl;
749 for(i=0; i<n; ++i){
750 term = pt->a * pt->t * (pt->t - 1) * pow(tau, pt->t - 2) * ipow(delta, pt->d);
751 sum += term;
752 #ifdef RESID_DEBUG
753 fprintf(stderr,"i = %d, a=%e, t=%f, d=%d, term = %f, sum = %f",i,pt->a,pt->t,pt->d,term,sum);
754 if(pt->l==0){
755 fprintf(stderr,",row=%e\n",term);
756 }else{
757 fprintf(stderr,",row=%e\n,",term*exp(-dell));
758 }
759 #endif
760 oldl = pt->l;
761 ++pt;
762 if(i+1==n || oldl != pt->l){
763 if(oldl == 0){
764 #ifdef RESID_DEBUG
765 fprintf(stderr,"linear ");
766 #endif
767 res += sum;
768 }else{
769 #ifdef RESID_DEBUG
770 fprintf(stderr,"exp dell=%f, exp(-dell)=%f sum=%f: ",dell,exp(-dell),sum);
771 #endif
772 res += sum * exp(-dell);
773 }
774 #ifdef RESID_DEBUG
775 fprintf(stderr,"i = %d, res = %f\n",i,res);
776 #endif
777 sum = 0;
778 dell = ipow(delta,pt->l);
779 ldell = pt->l*dell;
780 }
781 }
782
783 /* gaussian terms */
784 n = data->ng;
785 //fprintf(stderr,"THERE ARE %d GAUSSIAN TERMS\n",n);
786 gt = &(data->gt[0]);
787 for(i=0; i<n; ++i){
788 #ifdef RESID_DEBUG
789 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);
790 #endif
791 double d1 = delta - gt->epsilon;
792 double t1 = tau - gt->gamma;
793 double f1 = gt->t*(gt->t - 1) + 4. * gt->beta * tau * (tau * (gt->beta*SQ(t1) - 0.5) - t1*gt->t);
794 double e1 = -gt->alpha*SQ(d1) - gt->beta*SQ(t1);
795 sum = gt->n * f1 * pow(tau,gt->t - 2) * pow(delta,gt->d) * exp(e1);
796 //fprintf(stderr,"sum = %f\n",sum);
797 res += sum;
798 ++gt;
799 }
800
801 /* FIXME add critical terms calculation */
802
803 #ifdef RESID_DEBUG
804 fprintf(stderr,"phir_tautau = %f\n",res);
805 #endif
806 return res;
807 }
808

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