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

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