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

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