/[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 1989 - (show annotations) (download) (as text)
Tue Feb 3 07:29:29 2009 UTC (13 years, 5 months ago) by jpye
File MIME type: text/x-csrc
File size: 24958 byte(s)
Added helmholtz_w and helmholtz_cp.
Added test for helmholtz_w with water.
This exposes the fact that helm_resid_deldel and helm_resid_deltau need critical terms added.
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 Function to calculate isochoric heat capacity from the Helmholtz free energy
193 EOS given temperature and mass density.
194
195 @param T temperature in K
196 @param rho mass density in kg/m続
197 @return Isochoric specific heat capacity in J/kg/K.
198 */
199 double helmholtz_cv(double T, double rho, const HelmholtzData *data){
200 double tau = data->T_star / T;
201 double delta = rho / data->rho_star;
202
203 return - data->R * tau*tau * (helm_ideal_tautau(tau,data->ideal) + helm_resid_tautau(tau,delta,data));
204 }
205
206 /**
207 Function to calculate isobaric heat capacity from the Helmholtz free energy
208 EOS given temperature and mass density.
209
210 @param T temperature in K
211 @param rho mass density in kg/m続
212 @return Isobaric specific heat capacity in J/kg/K.
213 */
214 double helmholtz_cp(double T, double rho, const HelmholtzData *data){
215 double tau = data->T_star / T;
216 double delta = rho / data->rho_star;
217
218 double phir_d = helm_resid_del(tau,delta,data);
219 double phir_dd = helm_resid_deldel(tau,delta,data);
220 double phir_dt = helm_resid_deltau(tau,delta,data);
221
222 /* note similarities with helmholtz_w */
223 double temp1 = 1 + 2*delta*phir_d + SQ(delta)*phir_dd;
224 double temp2 = 1 + delta*phir_d - delta*tau*phir_dt;
225 double temp3 = SQ(tau)*(helm_ideal_tautau(tau,data->ideal) + helm_resid_tautau(tau,delta,data));
226
227 return data->R * (-temp3 + SQ(temp2)/temp1);
228 }
229
230
231 /**
232 Function to calculate the speed of sound in a fluid from the Helmholtz free
233 energy EOS, given temperature and mass density.
234
235 @param T temperature in K
236 @param rho mass density in kg/m続
237 @return Speed of sound in m/s.
238 */
239 double helmholtz_w(double T, double rho, const HelmholtzData *data){
240
241 double tau = data->T_star / T;
242 double delta = rho / data->rho_star;
243
244 double phir_d = helm_resid_del(tau,delta,data);
245 double phir_dd = helm_resid_deldel(tau,delta,data);
246 double phir_dt = helm_resid_deltau(tau,delta,data);
247
248 /* note similarities with helmholtz_cp */
249 double temp1 = 1 + 2*delta*phir_d + SQ(delta)*phir_dd;
250 double temp2 = 1 + delta*phir_d - delta*tau*phir_dt;
251 double temp3 = SQ(tau)*(helm_ideal_tautau(tau,data->ideal) + helm_resid_tautau(tau,delta,data));
252
253 return sqrt(data->R * T *(temp1 - SQ(temp2)/temp3));
254
255 }
256
257 /**
258 Calculation zero-pressure specific heat capacity
259 */
260 double helmholtz_cp0(double T, const HelmholtzData *data){
261 double val = helm_cp0(T,data->ideal);
262 #if 0
263 fprintf(stderr,"val = %f\n",val);
264 #endif
265 return val;
266 }
267
268 /**
269 Calculate partial derivative of p with respect to T, with rho constant
270 */
271 double helmholtz_dpdT_rho(double T, double rho, const HelmholtzData *data){
272 double tau = data->T_star / T;
273 double delta = rho / data->rho_star;
274
275 double phir_del = helm_resid_del(tau,delta,data);
276 double phir_deltau = helm_resid_deltau(tau,delta,data);
277 #ifdef TEST
278 assert(!isinf(phir_del));
279 assert(!isinf(phir_deltau));
280 assert(!isnan(phir_del));
281 assert(!isnan(phir_deltau));
282 assert(!isnan(data->R));
283 assert(!isnan(rho));
284 assert(!isnan(tau));
285 #endif
286
287 double res = data->R * rho * (1 + delta*phir_del - delta*tau*phir_deltau);
288
289 #ifdef TEST
290 assert(!isnan(res));
291 assert(!isinf(res));
292 #endif
293 return res;
294 }
295
296 /**
297 Calculate partial derivative of p with respect to rho, with T constant
298 */
299 double helmholtz_dpdrho_T(double T, double rho, const HelmholtzData *data){
300 double tau = data->T_star / T;
301 double delta = rho / data->rho_star;
302
303 double phir_del = helm_resid_del(tau,delta,data);
304 double phir_deldel = helm_resid_deldel(tau,delta,data);
305 #ifdef TEST
306 assert(!isinf(phir_del));
307 assert(!isinf(phir_deldel));
308 #endif
309 return data->R * T * (1 + 2*delta*phir_del + delta*delta* phir_deldel);
310 }
311
312 /**
313 Calculate partial derivative of h with respect to T, with rho constant
314 */
315 double helmholtz_dhdT_rho(double T, double rho, const HelmholtzData *data){
316 double tau = data->T_star / T;
317 double delta = rho / data->rho_star;
318
319 double phir_del = helm_resid_del(tau,delta,data);
320 double phir_deltau = helm_resid_deltau(tau,delta,data);
321 double phir_tautau = helm_resid_tautau(tau,delta,data);
322 double phi0_tautau = helm_ideal_tautau(tau,data->ideal);
323
324 //fprintf(stderr,"phir_del = %f, phir_deltau = %f, phir_tautau = %f, phi0_tautau = %f\n",phir_del,phir_deltau,phir_tautau,phi0_tautau);
325
326 //return (helmholtz_h(T+0.01,rho,data) - helmholtz_h(T,rho,data)) / 0.01;
327 return data->R * (1. + delta*phir_del - tau*tau*(phi0_tautau + phir_tautau) - delta*tau*phir_deltau);
328 }
329
330 /**
331 Calculate partial derivative of h with respect to rho, with T constant
332 */
333 double helmholtz_dhdrho_T(double T, double rho, const HelmholtzData *data){
334 double tau = data->T_star / T;
335 double delta = rho / data->rho_star;
336
337 double phir_del = helm_resid_del(tau,delta,data);
338 double phir_deltau = helm_resid_deltau(tau,delta,data);
339 double phir_deldel = helm_resid_deldel(tau,delta,data);
340
341 return data->R * T / rho * (tau*delta*(0 + phir_deltau) + delta * phir_del + SQ(delta)*phir_deldel);
342 }
343
344
345 /**
346 Calculate partial derivative of u with respect to T, with rho constant
347 */
348 double helmholtz_dudT_rho(double T, double rho, const HelmholtzData *data){
349 double tau = data->T_star / T;
350 double delta = rho / data->rho_star;
351
352 double phir_tautau = helm_resid_tautau(tau,delta,data);
353 double phi0_tautau = helm_ideal_tautau(tau,data->ideal);
354
355 return -data->R * SQ(tau) * (phi0_tautau + phir_tautau);
356 }
357
358 /**
359 Calculate partial derivative of u with respect to rho, with T constant
360 */
361 double helmholtz_dudrho_T(double T, double rho, const HelmholtzData *data){
362 double tau = data->T_star / T;
363 double delta = rho / data->rho_star;
364
365 double phir_deltau = helm_resid_deltau(tau,delta,data);
366
367 return data->R * T / rho * (tau * delta * phir_deltau);
368 }
369
370 /*---------------------------------------------
371 UTILITY FUNCTION(S)
372 */
373
374 /* ipow: public domain by Mark Stephen with suggestions by Keiichi Nakasato */
375 static double ipow(double x, int n){
376 double t = 1.0;
377
378 if(!n)return 1.0; /* At the top. x^0 = 1 */
379
380 if (n < 0){
381 n = -n;
382 x = 1.0/x; /* error if x == 0. Good */
383 } /* ZTC/SC returns inf, which is even better */
384
385 if (x == 0.0)return 0.0;
386
387 do{
388 if(n & 1)t *= x;
389 n /= 2; /* KN prefers if (n/=2) x*=x; This avoids an */
390 x *= x; /* unnecessary but benign multiplication on */
391 }while(n); /* the last pass, but the comparison is always
392 true _except_ on the last pass. */
393
394 return t;
395 }
396
397 //#define RESID_DEBUG
398
399 /**
400 Residual part of helmholtz function.
401 */
402 double helm_resid(double tau, double delta, const HelmholtzData *data){
403 double dell,ldell, term, sum, res = 0;
404 unsigned n, i;
405 const HelmholtzPowTerm *pt;
406 const HelmholtzGausTerm *gt;
407 const HelmholtzCritTerm *ct;
408
409 n = data->np;
410 pt = &(data->pt[0]);
411
412 #ifdef RESID_DEBUG
413 fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
414 #endif
415
416 /* power terms */
417 sum = 0;
418 dell = ipow(delta,pt->l);
419 ldell = pt->l * dell;
420 unsigned oldl;
421 for(i=0; i<n; ++i){
422 term = pt->a * pow(tau, pt->t) * ipow(delta, pt->d);
423 sum += term;
424 #ifdef RESID_DEBUG
425 fprintf(stderr,"i = %d, a=%e, t=%f, d=%d, term = %f, sum = %f",i,pt->a,pt->t,pt->d,term,sum);
426 if(pt->l==0){
427 fprintf(stderr,",row=%e\n",term);
428 }else{
429 fprintf(stderr,",row=%e\n,",term*exp(-dell));
430 }
431 #endif
432 oldl = pt->l;
433 ++pt;
434 if(i+1==n || oldl != pt->l){
435 if(oldl == 0){
436 #ifdef RESID_DEBUG
437 fprintf(stderr,"linear ");
438 #endif
439 res += sum;
440 }else{
441 #ifdef RESID_DEBUG
442 fprintf(stderr,"exp dell=%f, exp(-dell)=%f sum=%f: ",dell,exp(-dell),sum);
443 #endif
444 res += sum * exp(-dell);
445 }
446 #ifdef RESID_DEBUG
447 fprintf(stderr,"i = %d, res = %f\n",i,res);
448 #endif
449 sum = 0;
450 dell = ipow(delta,pt->l);
451 ldell = pt->l*dell;
452 }
453 }
454
455 /* gaussian terms */
456 n = data->ng;
457 //fprintf(stderr,"THERE ARE %d GAUSSIAN TERMS\n",n);
458 gt = &(data->gt[0]);
459 for(i=0; i<n; ++i){
460 #ifdef RESID_DEBUG
461 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);
462 #endif
463 double d1 = delta - gt->epsilon;
464 double t1 = tau - gt->gamma;
465 double e1 = -gt->alpha*d1*d1 - gt->beta*t1*t1;
466 sum = gt->n * pow(tau,gt->t) * pow(delta,gt->d) * exp(e1);
467 //fprintf(stderr,"sum = %f\n",sum);
468 res += sum;
469 ++gt;
470 }
471
472 /* critical terms */
473 n = data->nc;
474 ct = &(data->ct[0]);
475 for(i=0; i<n; ++i){
476 #ifdef RESID_DEBUG
477 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);
478 #endif
479 double d1 = delta - 1.;
480 double t1 = tau - 1.;
481 double theta = (1. - tau) + ct->A * pow(d1*d1, 0.5/ct->beta);
482 double psi = exp(-ct->C*d1*d1 - ct->D*t1*t1);
483 double DELTA = theta*theta + ct->B* pow(d1*d1, ct->a);
484 sum = ct->n * pow(DELTA, ct->b) * delta * psi;
485 res += sum;
486 ++ct;
487 }
488
489 #ifdef RESID_DEBUG
490 fprintf(stderr,"CALCULATED RESULT FOR phir = %f\n",res);
491 #endif
492 return res;
493 }
494
495 /*=================== FIRST DERIVATIVES =======================*/
496
497 /**
498 Derivative of the helmholtz residual function with respect to
499 delta.
500 */
501 double helm_resid_del(double tau,double delta, const HelmholtzData *data){
502 double sum = 0, res = 0;
503 double dell, ldell;
504 unsigned n, i;
505 const HelmholtzPowTerm *pt;
506 const HelmholtzGausTerm *gt;
507 const HelmholtzCritTerm *ct;
508
509 #ifdef RESID_DEBUG
510 fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
511 #endif
512
513 /* power terms */
514 n = data->np;
515 pt = &(data->pt[0]);
516 dell = ipow(delta,pt->l);
517 ldell = pt->l * dell;
518 unsigned oldl;
519 for(i=0; i<n; ++i){
520 sum += pt->a * pow(tau, pt->t) * ipow(delta, pt->d - 1) * (pt->d - ldell);
521 oldl = pt->l;
522 ++pt;
523 if(i+1==n || oldl != pt->l){
524 if(oldl == 0){
525 res += sum;
526 }else{
527 res += sum * exp(-dell);
528 }
529 sum = 0;
530 dell = ipow(delta,pt->l);
531 ldell = pt->l*dell;
532 }
533 }
534
535 /* gaussian terms */
536 n = data->ng;
537 gt = &(data->gt[0]);
538 for(i=0; i<n; ++i){
539 #ifdef RESID_DEBUG
540 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);
541 #endif
542 sum = - gt->n * pow(tau,gt->t) * pow(delta, -1. + gt->d)
543 * (2. * gt->alpha * delta * (delta - gt->epsilon) - gt->d)
544 * exp(-(gt->alpha * SQ(delta-gt->epsilon) + gt->beta*SQ(tau-gt->gamma)));
545 res += sum;
546 #ifdef RESID_DEBUG
547 fprintf(stderr,"sum = %f --> res = %f\n",sum,res);
548 #endif
549 ++gt;
550 }
551
552 /* critical terms */
553 n = data->nc;
554 ct = &(data->ct[0]);
555 for(i=0; i<n; ++i){
556 #ifdef RESID_DEBUG
557 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);
558 #endif
559 double d1 = delta - 1.;
560 double t1 = tau - 1.;
561 double theta = (1. - tau) + ct->A * pow(d1*d1, 0.5/ct->beta);
562 double psi = exp(-ct->C*d1*d1 - ct->D*t1*t1);
563 double DELTA = theta*theta + ct->B* pow(d1*d1, ct->a);
564
565 double dpsiddelta = -2. * ct->C * d1 * psi;
566
567 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));
568
569 double dDELbddelta = ct->b * pow(DELTA,ct->b - 1.) * dDELddelta;
570
571 sum = ct->n * (pow(DELTA, ct->b) * (psi + delta * dpsiddelta) + dDELbddelta * delta * psi);
572 res += sum;
573 ++ct;
574 }
575
576
577 return res;
578 }
579
580 /**
581 Derivative of the helmholtz residual function with respect to
582 tau.
583 */
584 double helm_resid_tau(double tau,double delta,const HelmholtzData *data){
585
586 double sum;
587 double res = 0;
588 double delX;
589 unsigned l;
590 unsigned n, i;
591 const HelmholtzPowTerm *pt;
592 const HelmholtzGausTerm *gt;
593 const HelmholtzCritTerm *ct;
594
595 n = data->np;
596 pt = &(data->pt[0]);
597
598 delX = 1;
599
600 l = 0;
601 sum = 0;
602 for(i=0; i<n; ++i){
603 if(pt->t){
604 //fprintf(stderr,"i = %d, a = %e, t = %f, d = %d, l = %d\n",i+1, pt->a, pt->t, pt->d, pt->l);
605 sum += pt->a * pow(tau, pt->t - 1) * ipow(delta, pt->d) * pt->t;
606 }
607 ++pt;
608 //fprintf(stderr,"l = %d\n",l);
609 if(i+1==n || l != pt->l){
610 if(l==0){
611 //fprintf(stderr,"Adding non-exp term\n");
612 res += sum;
613 }else{
614 //fprintf(stderr,"Adding exp term with l = %d, delX = %e\n",l,delX);
615 res += sum * exp(-delX);
616 }
617 /* set l to new value */
618 if(i+1!=n){
619 l = pt->l;
620 //fprintf(stderr,"New l = %d\n",l);
621 delX = ipow(delta,l);
622 sum = 0;
623 }
624 }
625 }
626
627 //#define RESID_DEBUG
628 /* gaussian terms */
629 n = data->ng;
630 gt = &(data->gt[0]);
631 for(i=0; i<n; ++i){
632 #ifdef RESID_DEBUG
633 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);
634 #endif
635
636 double val2;
637 val2 = -gt->n * pow(tau,gt->t - 1.) * pow(delta, gt->d)
638 * (2. * gt->beta * tau * (tau - gt->gamma) - gt->t)
639 * exp(-(gt->alpha * SQ(delta-gt->epsilon) + gt->beta*SQ(tau-gt->gamma)));
640 res += val2;
641 #ifdef RESID_DEBUG
642 fprintf(stderr,"res = %f\n",res);
643 #endif
644
645 ++gt;
646 }
647
648 /* critical terms */
649 n = data->nc;
650 ct = &(data->ct[0]);
651 for(i=0; i<n; ++i){
652 #ifdef RESID_DEBUG
653 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);
654 #endif
655 double d1 = delta - 1.;
656 double t1 = tau - 1.;
657 double theta = (1. - tau) + ct->A * pow(d1*d1, 0.5/ct->beta);
658 double psi = exp(-ct->C*d1*d1 - ct->D*t1*t1);
659 double DELTA = theta*theta + ct->B* pow(d1*d1, ct->a);
660
661 double dDELbdtau = -2. * theta * ct->b * pow(DELTA, ct->b - 1);
662
663 double dpsidtau = -2. * ct->D * t1 * psi;
664
665 sum = ct->n * delta * (dDELbdtau * psi + pow(DELTA, ct->b) * dpsidtau);
666 res += sum;
667 ++ct;
668 }
669
670 return res;
671 }
672
673
674 /*=================== SECOND DERIVATIVES =======================*/
675
676 /**
677 Mixed derivative of the helmholtz residual function with respect to
678 delta and tau.
679 */
680 double helm_resid_deltau(double tau,double delta,const HelmholtzData *data){
681 double dell,ldell, term, sum = 0, res = 0;
682 unsigned n, i;
683 const HelmholtzPowTerm *pt;
684 const HelmholtzGausTerm *gt;
685
686 /* power terms */
687 n = data->np;
688 pt = &(data->pt[0]);
689 dell = ipow(delta,pt->l);
690 ldell = pt->l * dell;
691 unsigned oldl;
692 for(i=0; i<n; ++i){
693 sum += pt->a * pt->t * pow(tau, pt->t - 1) * ipow(delta, pt->d - 1) * (pt->d - ldell);
694 oldl = pt->l;
695 ++pt;
696 if(i+1==n || oldl != pt->l){
697 if(oldl == 0){
698 res += sum;
699 }else{
700 res += sum * exp(-dell);
701 }
702 sum = 0;
703 dell = ipow(delta,pt->l);
704 ldell = pt->l*dell;
705 }
706 }
707
708 #ifdef TEST
709 assert(!isinf(res));
710 #endif
711
712 /* gaussian terms */
713 n = data->ng;
714 gt = &(data->gt[0]);
715 for(i=0; i<n; ++i){
716 #ifdef RESID_DEBUG
717 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);
718 #endif
719 double d1 = delta - gt->epsilon;
720 double t1 = tau - gt->gamma;
721 double e1 = -gt->alpha*SQ(d1) - gt->beta*SQ(t1);
722
723 double f1 = gt->t - 2*gt->beta*tau*(tau - gt->gamma);
724 double g1 = gt->d - 2*gt->alpha*delta*(delta - gt->epsilon);
725
726 sum = gt->n * f1 * pow(tau,gt->t-1) * g1 * pow(delta,gt->d-1) * exp(e1);
727
728 //fprintf(stderr,"sum = %f\n",sum);
729 res += sum;
730 #ifdef TEST
731 assert(!isinf(res));
732 #endif
733 ++gt;
734 }
735
736 /* FIXME add critical terms calculation */
737
738 #ifdef RESID_DEBUG
739 fprintf(stderr,"phir = %f\n",res);
740 #endif
741
742 #ifdef TEST
743 assert(!isnan(res));
744 assert(!isinf(res));
745 #endif
746 return res;
747 }
748
749 /**
750 Second derivative of helmholtz residual function with respect to
751 delta (twice).
752
753 FIXME this function is WRONG. (UPDATE? is this still true? Think not)
754 */
755 double helm_resid_deldel(double tau,double delta,const HelmholtzData *data){
756 double sum = 0, res = 0;
757 double dell, ldell;
758 unsigned n, i;
759 const HelmholtzPowTerm *pt;
760 const HelmholtzGausTerm *gt;
761
762
763 #ifdef RESID_DEBUG
764 fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
765 #endif
766
767 /* power terms */
768 n = data->np;
769 pt = &(data->pt[0]);
770 dell = ipow(delta,pt->l);
771 ldell = pt->l * dell;
772 unsigned oldl;
773 for(i=0; i<n; ++i){
774 double lpart = pt->l ? SQ(ldell) + ldell*(1. - 2*pt->d - pt->l) : 0;
775 sum += pt->a * pow(tau, pt->t) * ipow(delta, pt->d - 2) * (pt->d*(pt->d - 1) + lpart);
776 oldl = pt->l;
777 ++pt;
778 if(i+1==n || oldl != pt->l){
779 if(oldl == 0){
780 res += sum;
781 }else{
782 res += sum * exp(-dell);
783 }
784 sum = 0;
785 dell = ipow(delta,pt->l);
786 ldell = pt->l*dell;
787 }
788 }
789
790 /* gaussian terms */
791 n = data->ng;
792 //fprintf(stderr,"THERE ARE %d GAUSSIAN TERMS\n",n);
793 gt = &(data->gt[0]);
794 for(i=0; i<n; ++i){
795 double s1 = SQ(delta - gt->epsilon);
796 double f1 = gt->d*(gt->d - 1)
797 + 2.*gt->alpha*delta * (
798 delta * (2. * gt->alpha * s1 - 1)
799 - 2. * gt->d * (delta - gt->epsilon)
800 );
801 res += gt->n * pow(tau,gt->t) * pow(delta, gt->d - 2.)
802 * f1
803 * exp(-(gt->alpha * s1 + gt->beta*SQ(tau-gt->gamma)));
804 ++gt;
805 }
806
807 /* FIXME add critical terms calculation */
808
809 return res;
810 }
811
812
813
814 /**
815 Residual part of helmholtz function.
816 */
817 double helm_resid_tautau(double tau, double delta, const HelmholtzData *data){
818 double dell,ldell, term, sum, res = 0;
819 unsigned n, i;
820 const HelmholtzPowTerm *pt;
821 const HelmholtzGausTerm *gt;
822 const HelmholtzCritTerm *ct;
823
824 n = data->np;
825 pt = &(data->pt[0]);
826
827 #ifdef RESID_DEBUG
828 fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
829 #endif
830
831 /* power terms */
832 sum = 0;
833 dell = ipow(delta,pt->l);
834 ldell = pt->l * dell;
835 unsigned oldl;
836 for(i=0; i<n; ++i){
837 term = pt->a * pt->t * (pt->t - 1) * pow(tau, pt->t - 2) * ipow(delta, pt->d);
838 sum += term;
839 #ifdef RESID_DEBUG
840 fprintf(stderr,"i = %d, a=%e, t=%f, d=%d, term = %f, sum = %f",i,pt->a,pt->t,pt->d,term,sum);
841 if(pt->l==0){
842 fprintf(stderr,",row=%e\n",term);
843 }else{
844 fprintf(stderr,",row=%e\n,",term*exp(-dell));
845 }
846 #endif
847 oldl = pt->l;
848 ++pt;
849 if(i+1==n || oldl != pt->l){
850 if(oldl == 0){
851 #ifdef RESID_DEBUG
852 fprintf(stderr,"linear ");
853 #endif
854 res += sum;
855 }else{
856 #ifdef RESID_DEBUG
857 fprintf(stderr,"exp dell=%f, exp(-dell)=%f sum=%f: ",dell,exp(-dell),sum);
858 #endif
859 res += sum * exp(-dell);
860 }
861 #ifdef RESID_DEBUG
862 fprintf(stderr,"i = %d, res = %f\n",i,res);
863 #endif
864 sum = 0;
865 dell = ipow(delta,pt->l);
866 ldell = pt->l*dell;
867 }
868 }
869
870 /* gaussian terms */
871 n = data->ng;
872 //fprintf(stderr,"THERE ARE %d GAUSSIAN TERMS\n",n);
873 gt = &(data->gt[0]);
874 for(i=0; i<n; ++i){
875 #ifdef RESID_DEBUG
876 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);
877 #endif
878 double d1 = delta - gt->epsilon;
879 double t1 = tau - gt->gamma;
880 double f1 = gt->t*(gt->t - 1) + 4. * gt->beta * tau * (tau * (gt->beta*SQ(t1) - 0.5) - t1*gt->t);
881 double e1 = -gt->alpha*SQ(d1) - gt->beta*SQ(t1);
882 sum = gt->n * f1 * pow(tau,gt->t - 2) * pow(delta,gt->d) * exp(e1);
883 //fprintf(stderr,"sum = %f\n",sum);
884 res += sum;
885 ++gt;
886 }
887
888 /* critical terms */
889 n = data->nc;
890 ct = &(data->ct[0]);
891 for(i=0; i<n; ++i){
892 #ifdef RESID_DEBUG
893 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);
894 #endif
895 double d1 = delta - 1.;
896 double t1 = tau - 1.;
897 double theta = (1. - tau) + ct->A * pow(d1*d1, 0.5/ct->beta);
898 double psi = exp(-ct->C*d1*d1 - ct->D*t1*t1);
899 double DELTA = theta*theta + ct->B* pow(d1*d1, ct->a);
900
901 double d2DELbdtau2 = 2. * ct->b * pow(DELTA, ct->b - 1) + 4. * theta*theta * ct->b * (ct->b - 1) * pow(DELTA, ct->b - 2);
902 double dDELbdtau = -2. * theta * ct->b * pow(DELTA, ct->b - 1);
903 double dpsidtau = -2. * ct->D * t1 * psi;
904 double d2psidtau2 = 2. * ct->D * psi * (2. * ct->D * t1*t1 -1.);
905
906 sum = ct->n * delta * (d2DELbdtau2 * psi + 2 * dDELbdtau*dpsidtau + pow(DELTA, ct->b) * d2psidtau2);
907 res += sum;
908 ++ct;
909 }
910
911 #ifdef RESID_DEBUG
912 fprintf(stderr,"phir_tautau = %f\n",res);
913 #endif
914 return res;
915 }
916

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