/[ascend]/trunk/models/johnpye/fprops/helmholtz.c
ViewVC logotype

Annotation of /trunk/models/johnpye/fprops/helmholtz.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1907 - (hide annotations) (download) (as text)
Sat Sep 27 09:45:06 2008 UTC (13 years, 10 months ago) by jpye
File MIME type: text/x-csrc
File size: 15576 byte(s)
Working on dpdrho_T, still looks like problems with Gaussian terms.
1 jpye 1810 /* ASCEND modelling environment
2 jpye 1831 Copyright (C) 2008 Carnegie Mellon University
3 jpye 1810
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 jpye 1822 For nomenclature see Tillner-Roth, Harms-Watzenberg and Baehr, Eine neue
22     Fundamentalgleichung f端r Ammoniak.
23    
24 jpye 1810 John Pye, 29 Jul 2008.
25     */
26    
27 jpye 1822 #include <math.h>
28 jpye 1810
29 jpye 1822 #include "helmholtz.h"
30 jpye 1847 #include "ideal_impl.h"
31 jpye 1822
32 jpye 1826 #ifdef TEST
33     #include <assert.h>
34     #include <stdlib.h>
35     #include <stdio.h>
36     #endif
37    
38 jpye 1905 #define SQ(X) ((X)*(X))
39    
40 jpye 1822 /* forward decls */
41    
42 jpye 1826 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 jpye 1822
48 jpye 1810 /**
49     Function to calculate pressure from Helmholtz free energy EOS, given temperature
50 jpye 1824 and mass density.
51 jpye 1810
52     @param T temperature in K
53 jpye 1824 @param rho mass density in kg/m続
54     @return pressure in Pa???
55 jpye 1810 */
56 jpye 1826 double helmholtz_p(double T, double rho, const HelmholtzData *data){
57 jpye 1822
58     double tau = data->T_star / T;
59     double delta = rho / data->rho_star;
60 jpye 1810
61 jpye 1838 #ifdef TEST
62     assert(data->rho_star!=0);
63     assert(T!=0);
64     assert(!isnan(tau));
65     assert(!isnan(delta));
66     assert(!isnan(data->R));
67 jpye 1839
68 jpye 1840 //fprintf(stderr,"p calc: T = %f\n",T);
69     //fprintf(stderr,"p calc: tau = %f\n",tau);
70     //fprintf(stderr,"p calc: rho = %f\n",rho);
71     //fprintf(stderr,"p calc: delta = %f\n",delta);
72     //fprintf(stderr,"p calc: R*T*rho = %f\n",data->R * T * rho);
73 jpye 1841
74     //fprintf(stderr,"T = %f\n", T);
75     //fprintf(stderr,"rhob = %f, rhob* = %f, delta = %f\n", rho/data->M, data->rho_star/data->M, delta);
76 jpye 1838 #endif
77 jpye 1839
78 jpye 1841 return data->R * T * rho * (1 + delta * helm_resid_del(tau,delta,data));
79 jpye 1822 }
80 jpye 1810
81 jpye 1824 /**
82     Function to calculate internal energy from Helmholtz free energy EOS, given
83     temperature and mass density.
84    
85     @param T temperature in K
86     @param rho mass density in kg/m続
87     @return internal energy in ???
88     */
89 jpye 1826 double helmholtz_u(double T, double rho, const HelmholtzData *data){
90 jpye 1824
91     double tau = data->T_star / T;
92     double delta = rho / data->rho_star;
93    
94 jpye 1832 #ifdef TEST
95 jpye 1838 assert(data->rho_star!=0);
96     assert(T!=0);
97     assert(!isnan(tau));
98     assert(!isnan(delta));
99     assert(!isnan(data->R));
100     #endif
101    
102 jpye 1857 #if 0
103 jpye 1835 fprintf(stderr,"ideal_tau = %f\n",helm_ideal_tau(tau,delta,data->ideal));
104 jpye 1827 fprintf(stderr,"resid_tau = %f\n",helm_resid_tau(tau,delta,data));
105     fprintf(stderr,"R T = %f\n",data->R * data->T_star);
106 jpye 1832 #endif
107 jpye 1827
108 jpye 1835 return data->R * data->T_star * (helm_ideal_tau(tau,delta,data->ideal) + helm_resid_tau(tau,delta,data));
109 jpye 1824 }
110    
111     /**
112     Function to calculate enthalpy from Helmholtz free energy EOS, given
113     temperature and mass density.
114    
115     @param T temperature in K
116     @param rho mass density in kg/m続
117 jpye 1829 @return enthalpy in J/kg
118 jpye 1824 */
119 jpye 1826 double helmholtz_h(double T, double rho, const HelmholtzData *data){
120 jpye 1824
121     double tau = data->T_star / T;
122     double delta = rho / data->rho_star;
123    
124 jpye 1838 #ifdef TEST
125     assert(data->rho_star!=0);
126     assert(T!=0);
127     assert(!isnan(tau));
128     assert(!isnan(delta));
129     assert(!isnan(data->R));
130     #endif
131    
132 jpye 1835 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));
133 jpye 1824 }
134    
135 jpye 1829 /**
136     Function to calculate entropy from Helmholtz free energy EOS, given
137     temperature and mass density.
138    
139     @param T temperature in K
140     @param rho mass density in kg/m続
141     @return entropy in J/kgK
142     */
143     double helmholtz_s(double T, double rho, const HelmholtzData *data){
144    
145     double tau = data->T_star / T;
146     double delta = rho / data->rho_star;
147    
148 jpye 1870 #ifdef ENTROPY_DEBUG
149 jpye 1838 assert(data->rho_star!=0);
150     assert(T!=0);
151     assert(!isnan(tau));
152     assert(!isnan(delta));
153     assert(!isnan(data->R));
154    
155 jpye 1845 fprintf(stderr,"helm_ideal_tau = %f\n",helm_ideal_tau(tau,delta,data->ideal));
156     fprintf(stderr,"helm_resid_tau = %f\n",helm_resid_tau(tau,delta,data));
157     fprintf(stderr,"helm_ideal = %f\n",helm_ideal(tau,delta,data->ideal));
158     fprintf(stderr,"helm_resid = %f\n",helm_resid(tau,delta,data));
159 jpye 1847 #endif
160 jpye 1829 return data->R * (
161 jpye 1835 tau * (helm_ideal_tau(tau,delta,data->ideal) + helm_resid_tau(tau,delta,data))
162 jpye 1871 - (helm_ideal(tau,delta,data->ideal) + helm_resid(tau,delta,data))
163 jpye 1829 );
164     }
165    
166 jpye 1850 /**
167 jpye 1863 Function to calculate Helmholtz energy from the Helmholtz free energy EOS,
168     given temperature and mass density.
169    
170     @param T temperature in K
171     @param rho mass density in kg/m続
172     @return Helmholtz energy 'a', in J/kg
173     */
174     double helmholtz_a(double T, double rho, const HelmholtzData *data){
175    
176     double tau = data->T_star / T;
177     double delta = rho / data->rho_star;
178    
179     #ifdef TEST
180     assert(data->rho_star!=0);
181     assert(T!=0);
182     assert(!isnan(tau));
183     assert(!isnan(delta));
184     assert(!isnan(data->R));
185     #endif
186    
187 jpye 1868 #ifdef HELMHOLTZ_DEBUG
188 jpye 1863 fprintf(stderr,"helmholtz_a: T = %f, rho = %f\n",T,rho);
189 jpye 1865 fprintf(stderr,"multiplying by RT = %f\n",data->R*T);
190 jpye 1863 #endif
191    
192 jpye 1870 return data->R * T * (helm_ideal(tau,delta,data->ideal) + helm_resid(tau,delta,data));
193 jpye 1863 }
194    
195    
196     /**
197 jpye 1850 Calculation zero-pressure specific heat capacity
198     */
199     double helmholtz_cp0(double T, const HelmholtzData *data){
200     double val = helm_cp0(T,data->ideal);
201 jpye 1855 #if 0
202 jpye 1850 fprintf(stderr,"val = %f\n",val);
203 jpye 1851 #endif
204 jpye 1850 return val;
205     }
206    
207 jpye 1905 /**
208     Calculate partial derivative of p with respect to T, with rho constant
209     */
210     double helmholtz_dpdT_rho(double T, double rho, const HelmholtzData *data){
211     double tau = data->T_star / T;
212     double delta = rho / data->rho_star;
213 jpye 1906
214     double phir_del = helm_resid_del(tau,delta,data);
215     double phir_deltau = helm_resid_deltau(tau,delta,data);
216     assert(!isinf(phir_del));
217     assert(!isinf(phir_deltau));
218     assert(!isnan(phir_del));
219     assert(!isnan(phir_deltau));
220     assert(!isnan(data->R));
221     assert(!isnan(rho));
222     assert(!isnan(tau));
223    
224     double res = data->R * rho * (1 + delta*phir_del - delta*tau*phir_deltau);
225    
226     assert(!isnan(res));
227     assert(!isinf(res));
228     return res;
229 jpye 1905 }
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 jpye 1906
238     double phir_del = helm_resid_del(tau,delta,data);
239     double phir_deldel = helm_resid_deldel(tau,delta,data);
240     assert(!isinf(phir_del));
241     assert(!isinf(phir_deldel));
242 jpye 1905
243 jpye 1906 return data->R * T * (1 + 2*delta*phir_del + delta*delta* phir_deldel);
244 jpye 1905 }
245    
246 jpye 1824 /*---------------------------------------------
247 jpye 1826 UTILITY FUNCTION(S)
248     */
249    
250     /* ipow: public domain by Mark Stephen with suggestions by Keiichi Nakasato */
251     static double ipow(double x, int n){
252     double t = 1.0;
253    
254 jpye 1829 if(!n)return 1.0; /* At the top. x^0 = 1 */
255 jpye 1826
256     if (n < 0){
257     n = -n;
258     x = 1.0/x; /* error if x == 0. Good */
259     } /* ZTC/SC returns inf, which is even better */
260    
261     if (x == 0.0)return 0.0;
262    
263     do{
264     if(n & 1)t *= x;
265     n /= 2; /* KN prefers if (n/=2) x*=x; This avoids an */
266     x *= x; /* unnecessary but benign multiplication on */
267     }while(n); /* the last pass, but the comparison is always
268     true _except_ on the last pass. */
269    
270     return t;
271     }
272    
273 jpye 1889 //#define RESID_DEBUG
274 jpye 1872
275 jpye 1824 /**
276 jpye 1845 Residual part of helmholtz function.
277 jpye 1822 */
278 jpye 1826 double helm_resid(double tau, double delta, const HelmholtzData *data){
279 jpye 1867 double dell,ldell, term, sum, res = 0;
280 jpye 1845 unsigned n, i;
281     const HelmholtzPowTerm *pt;
282 jpye 1887 const HelmholtzGausTerm *gt;
283 jpye 1845
284     n = data->np;
285     pt = &(data->pt[0]);
286    
287 jpye 1868 #ifdef RESID_DEBUG
288 jpye 1867 fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
289     #endif
290    
291 jpye 1845 /* power terms */
292     sum = 0;
293 jpye 1865 dell = ipow(delta,pt->l);
294     ldell = pt->l * dell;
295 jpye 1845 unsigned oldl;
296     for(i=0; i<n; ++i){
297 jpye 1867 term = pt->a * pow(tau, pt->t) * ipow(delta, pt->d);
298     sum += term;
299 jpye 1868 #ifdef RESID_DEBUG
300 jpye 1867 fprintf(stderr,"i = %d, a=%e, t=%f, d=%d, term = %f, sum = %f",i,pt->a,pt->t,pt->d,term,sum);
301     if(pt->l==0){
302     fprintf(stderr,",row=%e\n",term);
303     }else{
304     fprintf(stderr,",row=%e\n,",term*exp(-dell));
305     }
306     #endif
307 jpye 1845 oldl = pt->l;
308     ++pt;
309     if(i+1==n || oldl != pt->l){
310     if(oldl == 0){
311 jpye 1868 #ifdef RESID_DEBUG
312 jpye 1865 fprintf(stderr,"linear ");
313 jpye 1867 #endif
314 jpye 1845 res += sum;
315     }else{
316 jpye 1868 #ifdef RESID_DEBUG
317 jpye 1865 fprintf(stderr,"exp dell=%f, exp(-dell)=%f sum=%f: ",dell,exp(-dell),sum);
318 jpye 1867 #endif
319 jpye 1865 res += sum * exp(-dell);
320 jpye 1845 }
321 jpye 1868 #ifdef RESID_DEBUG
322 jpye 1865 fprintf(stderr,"i = %d, res = %f\n",i,res);
323 jpye 1867 #endif
324 jpye 1845 sum = 0;
325 jpye 1865 dell = ipow(delta,pt->l);
326     ldell = pt->l*dell;
327 jpye 1845 }
328     }
329    
330 jpye 1887 /* gaussian terms */
331     n = data->ng;
332 jpye 1889 //fprintf(stderr,"THERE ARE %d GAUSSIAN TERMS\n",n);
333 jpye 1887 gt = &(data->gt[0]);
334     for(i=0; i<n; ++i){
335 jpye 1868 #ifdef RESID_DEBUG
336 jpye 1887 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);
337     #endif
338     double d1 = delta - gt->epsilon;
339     double t1 = tau - gt->gamma;
340     double e1 = -gt->alpha*d1*d1 - gt->beta*t1*t1;
341     sum = gt->n * pow(tau,gt->t) * pow(delta,gt->d) * exp(e1);
342     //fprintf(stderr,"sum = %f\n",sum);
343     res += sum;
344     ++gt;
345     }
346    
347     #ifdef RESID_DEBUG
348 jpye 1864 fprintf(stderr,"phir = %f\n",res);
349     #endif
350 jpye 1845 return res;
351     }
352    
353 jpye 1822 /**
354     Derivative of the helmholtz residual function with respect to
355     delta.
356     */
357 jpye 1826 double helm_resid_del(double tau,double delta, const HelmholtzData *data){
358 jpye 1905 double sum = 0, res = 0;
359 jpye 1844 double dell, ldell;
360 jpye 1838 unsigned n, i;
361     const HelmholtzPowTerm *pt;
362 jpye 1887 const HelmholtzGausTerm *gt;
363 jpye 1810
364 jpye 1822
365 jpye 1898 #ifdef RESID_DEBUG
366     fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
367     #endif
368    
369 jpye 1905 /* power terms */
370     n = data->np;
371     pt = &(data->pt[0]);
372 jpye 1844 dell = ipow(delta,pt->l);
373     ldell = pt->l * dell;
374     unsigned oldl;
375 jpye 1841 for(i=0; i<n; ++i){
376 jpye 1844 sum += pt->a * pow(tau, pt->t) * ipow(delta, pt->d - 1) * (pt->d - ldell);
377     oldl = pt->l;
378 jpye 1841 ++pt;
379 jpye 1844 if(i+1==n || oldl != pt->l){
380     if(oldl == 0){
381 jpye 1836 res += sum;
382     }else{
383 jpye 1844 res += sum * exp(-dell);
384 jpye 1836 }
385 jpye 1844 sum = 0;
386     dell = ipow(delta,pt->l);
387     ldell = pt->l*dell;
388 jpye 1836 }
389 jpye 1822 }
390    
391 jpye 1887 /* gaussian terms */
392     n = data->ng;
393 jpye 1889 //fprintf(stderr,"THERE ARE %d GAUSSIAN TERMS\n",n);
394 jpye 1887 gt = &(data->gt[0]);
395     for(i=0; i<n; ++i){
396     #ifdef RESID_DEBUG
397 jpye 1891 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);
398 jpye 1840 #endif
399 jpye 1898 double val2;
400     val2 = - gt->n * pow(tau,gt->t) * pow(delta, -1. + gt->d)
401     * (2. * gt->alpha * delta * (delta - gt->epsilon) - gt->d)
402     * exp(-(gt->alpha * SQ(delta-gt->epsilon) + gt->beta*SQ(tau-gt->gamma)));
403     res += val2;
404 jpye 1891 #ifdef RESID_DEBUG
405 jpye 1898 fprintf(stderr,"val2 = %f --> res = %f\n",val2,res);
406 jpye 1891 #endif
407 jpye 1887 ++gt;
408     }
409 jpye 1838
410 jpye 1836 return res;
411 jpye 1822 }
412    
413 jpye 1824 /**
414     Derivative of the helmholtz residual function with respect to
415     tau.
416     */
417 jpye 1826 double helm_resid_tau(double tau,double delta,const HelmholtzData *data){
418 jpye 1824
419     double sum;
420 jpye 1827 double res = 0;
421 jpye 1832 double delX;
422     unsigned l;
423 jpye 1845 unsigned n, i;
424 jpye 1838 const HelmholtzPowTerm *pt;
425 jpye 1891 const HelmholtzGausTerm *gt;
426 jpye 1822
427 jpye 1845 n = data->np;
428 jpye 1838 pt = &(data->pt[0]);
429 jpye 1832
430     delX = 1;
431    
432     l = 0;
433     sum = 0;
434 jpye 1845 for(i=0; i<n; ++i){
435 jpye 1838 if(pt->t){
436     //fprintf(stderr,"i = %d, a = %e, t = %f, d = %d, l = %d\n",i+1, pt->a, pt->t, pt->d, pt->l);
437     sum += pt->a * pow(tau, pt->t - 1) * ipow(delta, pt->d) * pt->t;
438 jpye 1832 }
439 jpye 1838 ++pt;
440 jpye 1832 //fprintf(stderr,"l = %d\n",l);
441 jpye 1845 if(i+1==n || l != pt->l){
442 jpye 1832 if(l==0){
443     //fprintf(stderr,"Adding non-exp term\n");
444     res += sum;
445     }else{
446     //fprintf(stderr,"Adding exp term with l = %d, delX = %e\n",l,delX);
447     res += sum * exp(-delX);
448     }
449     /* set l to new value */
450 jpye 1845 if(i+1!=n){
451 jpye 1838 l = pt->l;
452 jpye 1832 //fprintf(stderr,"New l = %d\n",l);
453     delX = ipow(delta,l);
454     sum = 0;
455     }
456     }
457     }
458    
459 jpye 1897 //#define RESID_DEBUG
460 jpye 1891 /* gaussian terms */
461     n = data->ng;
462     gt = &(data->gt[0]);
463     for(i=0; i<n; ++i){
464     #ifdef RESID_DEBUG
465     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);
466     #endif
467 jpye 1897
468     double val2;
469 jpye 1900 val2 = -gt->n * pow(tau,gt->t - 1.) * pow(delta, gt->d)
470     * (2. * gt->beta * tau * (tau - gt->gamma) - gt->t)
471 jpye 1898 * exp(-(gt->alpha * SQ(delta-gt->epsilon) + gt->beta*SQ(tau-gt->gamma)));
472     res += val2;
473 jpye 1891 #ifdef RESID_DEBUG
474 jpye 1897 fprintf(stderr,"res = %f\n",res);
475 jpye 1891 #endif
476 jpye 1897
477 jpye 1891 ++gt;
478     }
479    
480 jpye 1827 return res;
481 jpye 1824 }
482    
483    
484    
485 jpye 1826 /**
486     Mixed derivative of the helmholtz residual function with respect to
487 jpye 1907 delta and tau.
488 jpye 1868 */
489 jpye 1826 double helm_resid_deltau(double tau,double delta,const HelmholtzData *data){
490 jpye 1906 double dell,ldell, term, sum = 0, res = 0;
491 jpye 1905 unsigned n, i;
492     const HelmholtzPowTerm *pt;
493     const HelmholtzGausTerm *gt;
494 jpye 1824
495 jpye 1905 /* power terms */
496     n = data->np;
497     pt = &(data->pt[0]);
498     dell = ipow(delta,pt->l);
499     ldell = pt->l * dell;
500     unsigned oldl;
501     for(i=0; i<n; ++i){
502     sum += pt->a * pt->t * pow(tau, pt->t - 1) * ipow(delta, pt->d - 1) * (pt->d - ldell);
503     oldl = pt->l;
504 jpye 1838 ++pt;
505 jpye 1905 if(i+1==n || oldl != pt->l){
506     if(oldl == 0){
507     res += sum;
508     }else{
509     res += sum * exp(-dell);
510     }
511     sum = 0;
512     dell = ipow(delta,pt->l);
513     ldell = pt->l*dell;
514     }
515 jpye 1826 }
516    
517 jpye 1906 assert(!isinf(res));
518    
519 jpye 1905 /* gaussian terms */
520     n = data->ng;
521     gt = &(data->gt[0]);
522     for(i=0; i<n; ++i){
523     #ifdef RESID_DEBUG
524     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);
525     #endif
526     double d1 = delta - gt->epsilon;
527     double t1 = tau - gt->gamma;
528     double e1 = -gt->alpha*SQ(d1) - gt->beta*SQ(t1);
529 jpye 1826
530 jpye 1905 double f1 = gt->t - 2*gt->beta*tau*(tau - gt->gamma);
531     double g1 = gt->d - 2*gt->alpha*delta*(delta - gt->epsilon);
532 jpye 1826
533 jpye 1905 sum = gt->n * f1 * pow(tau,gt->t-1) * g1 * pow(delta,gt->d-1) * exp(e1);
534 jpye 1906
535 jpye 1905 //fprintf(stderr,"sum = %f\n",sum);
536     res += sum;
537 jpye 1906 assert(!isinf(res));
538    
539 jpye 1905 ++gt;
540 jpye 1826 }
541    
542 jpye 1905 #ifdef RESID_DEBUG
543     fprintf(stderr,"phir = %f\n",res);
544     #endif
545 jpye 1906
546     assert(!isnan(res));
547     assert(!isinf(res));
548 jpye 1905 return res;
549 jpye 1826 }
550    
551     /**
552     Second derivative of helmholtz residual function with respect to
553     delta (twice).
554 jpye 1890
555     FIXME this function is WRONG.
556 jpye 1826 */
557     double helm_resid_deldel(double tau,double delta,const HelmholtzData *data){
558 jpye 1907 double sum = 0, res = 0;
559     double dell, ldell;
560     unsigned n, i;
561     const HelmholtzPowTerm *pt;
562     const HelmholtzGausTerm *gt;
563 jpye 1826
564    
565 jpye 1907 #ifdef RESID_DEBUG
566     fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
567     #endif
568 jpye 1826
569 jpye 1907 /* power terms */
570     n = data->np;
571     pt = &(data->pt[0]);
572     dell = ipow(delta,pt->l);
573     ldell = pt->l * dell;
574     unsigned oldl;
575     for(i=0; i<n; ++i){
576     double lpart = pt->l ? SQ(ldell) + ldell*(1. - 2*pt->d - pt->l) : 0;
577     sum += pt->a * pow(tau, pt->t) * ipow(delta, pt->d - 2) * (pt->d*(pt->d - 1) + lpart);
578     oldl = pt->l;
579 jpye 1838 ++pt;
580 jpye 1907 if(i+1==n || oldl != pt->l){
581     if(oldl == 0){
582     res += sum;
583     }else{
584     res += sum * exp(-dell);
585     }
586     sum = 0;
587     dell = ipow(delta,pt->l);
588     ldell = pt->l*dell;
589     }
590 jpye 1826 }
591    
592 jpye 1907 /* gaussian terms */
593     n = data->ng;
594     //fprintf(stderr,"THERE ARE %d GAUSSIAN TERMS\n",n);
595     gt = &(data->gt[0]);
596     for(i=0; i<n; ++i){
597     double s1 = SQ(delta - gt->epsilon);
598     double f1 = 2*delta*gt->alpha *(2*gt->d*gt->epsilon
599     - delta * (2*gt->d + 1 - 2 * gt->alpha * s1));
600     res += - gt->n * pow(tau,gt->t) * pow(delta, -1. + gt->d)
601     * f1
602     * exp(-(gt->alpha * s1 + gt->beta*SQ(tau-gt->gamma)));
603     ++gt;
604 jpye 1826 }
605    
606 jpye 1907 return res;
607 jpye 1826 }
608    

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