/[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 1994 - (show annotations) (download) (as text)
Wed Feb 4 03:52:36 2009 UTC (11 years, 6 months ago) by jpye
File MIME type: text/x-csrc
File size: 26360 byte(s)
Still some error near critical point for calculation of speed of sound.
However this error is not reconcilable with both IAPWS95 and REFPROP8, because those
two sources give different values. The value returned by FPROPS is midway between
the two.
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 * SQ(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 + SQ(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 - SQ(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 We avoid duplication by using the following #defines for common code in
401 calculation of critical terms.
402 */
403 #define DEFINE_DELTA \
404 double d1 = delta - 1.; \
405 double t1 = tau - 1.; \
406 double d12 = SQ(d1); \
407 double theta = (1. - tau) + ct->A * pow(d12, 0.5/ct->beta); \
408 double psi = exp(-ct->C*d12 - ct->D*SQ(t1)); \
409 double DELTA = SQ(theta) + ct->B* pow(d12, ct->a)
410
411 #define DEFINE_DPSIDDELTA \
412 double dpsiddelta = -2. * ct->C * d1 * psi
413
414 #define DEFINE_DDELDDELTA \
415 double dDELddelta = d1 * (ct->A * theta * 2./ct->beta * pow(d12, 0.5/ct->beta - 1) + 2* ct->B * ct->a * pow(d12, ct->a - 1))
416
417 #define DEFINE_DDELBDTAU \
418 double dDELbdtau = -2. * theta * ct->b * pow(DELTA, ct->b - 1)
419
420 #define DEFINE_DPSIDTAU \
421 double dpsidtau = -2. * ct->D * t1 * psi
422
423 #define DEFINE_DDELBDDELTA \
424 double dDELbddelta = ct->b * pow(DELTA,ct->b - 1.) * dDELddelta
425
426 /**
427 Residual part of helmholtz function.
428 */
429 double helm_resid(double tau, double delta, const HelmholtzData *data){
430 double dell,ldell, term, sum, res = 0;
431 unsigned n, i;
432 const HelmholtzPowTerm *pt;
433 const HelmholtzGausTerm *gt;
434 const HelmholtzCritTerm *ct;
435
436 n = data->np;
437 pt = &(data->pt[0]);
438
439 #ifdef RESID_DEBUG
440 fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
441 #endif
442
443 /* power terms */
444 sum = 0;
445 dell = ipow(delta,pt->l);
446 ldell = pt->l * dell;
447 unsigned oldl;
448 for(i=0; i<n; ++i){
449 term = pt->a * pow(tau, pt->t) * ipow(delta, pt->d);
450 sum += term;
451 #ifdef RESID_DEBUG
452 fprintf(stderr,"i = %d, a=%e, t=%f, d=%d, term = %f, sum = %f",i,pt->a,pt->t,pt->d,term,sum);
453 if(pt->l==0){
454 fprintf(stderr,",row=%e\n",term);
455 }else{
456 fprintf(stderr,",row=%e\n,",term*exp(-dell));
457 }
458 #endif
459 oldl = pt->l;
460 ++pt;
461 if(i+1==n || oldl != pt->l){
462 if(oldl == 0){
463 #ifdef RESID_DEBUG
464 fprintf(stderr,"linear ");
465 #endif
466 res += sum;
467 }else{
468 #ifdef RESID_DEBUG
469 fprintf(stderr,"exp dell=%f, exp(-dell)=%f sum=%f: ",dell,exp(-dell),sum);
470 #endif
471 res += sum * exp(-dell);
472 }
473 #ifdef RESID_DEBUG
474 fprintf(stderr,"i = %d, res = %f\n",i,res);
475 #endif
476 sum = 0;
477 dell = ipow(delta,pt->l);
478 ldell = pt->l*dell;
479 }
480 }
481
482 /* gaussian terms */
483 n = data->ng;
484 //fprintf(stderr,"THERE ARE %d GAUSSIAN TERMS\n",n);
485 gt = &(data->gt[0]);
486 for(i=0; i<n; ++i){
487 #ifdef RESID_DEBUG
488 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);
489 #endif
490 double d1 = delta - gt->epsilon;
491 double t1 = tau - gt->gamma;
492 double e1 = -gt->alpha*SQ(d1) - gt->beta*SQ(t1);
493 sum = gt->n * pow(tau,gt->t) * pow(delta,gt->d) * exp(e1);
494 //fprintf(stderr,"sum = %f\n",sum);
495 res += sum;
496 ++gt;
497 }
498
499 /* critical terms */
500 n = data->nc;
501 ct = &(data->ct[0]);
502 for(i=0; i<n; ++i){
503 #ifdef RESID_DEBUG
504 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);
505 #endif
506
507 DEFINE_DELTA;
508
509 sum = ct->n * pow(DELTA, ct->b) * delta * psi;
510 res += sum;
511 ++ct;
512 }
513
514 #ifdef RESID_DEBUG
515 fprintf(stderr,"CALCULATED RESULT FOR phir = %f\n",res);
516 #endif
517 return res;
518 }
519
520 /*=================== FIRST DERIVATIVES =======================*/
521
522 /**
523 Derivative of the helmholtz residual function with respect to
524 delta.
525 */
526 double helm_resid_del(double tau,double delta, const HelmholtzData *data){
527 double sum = 0, res = 0;
528 double dell, ldell;
529 unsigned n, i;
530 const HelmholtzPowTerm *pt;
531 const HelmholtzGausTerm *gt;
532 const HelmholtzCritTerm *ct;
533
534 #ifdef RESID_DEBUG
535 fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
536 #endif
537
538 /* power terms */
539 n = data->np;
540 pt = &(data->pt[0]);
541 dell = ipow(delta,pt->l);
542 ldell = pt->l * dell;
543 unsigned oldl;
544 for(i=0; i<n; ++i){
545 sum += pt->a * pow(tau, pt->t) * ipow(delta, pt->d - 1) * (pt->d - ldell);
546 oldl = pt->l;
547 ++pt;
548 if(i+1==n || oldl != pt->l){
549 if(oldl == 0){
550 res += sum;
551 }else{
552 res += sum * exp(-dell);
553 }
554 sum = 0;
555 dell = ipow(delta,pt->l);
556 ldell = pt->l*dell;
557 }
558 }
559
560 /* gaussian terms */
561 n = data->ng;
562 gt = &(data->gt[0]);
563 for(i=0; i<n; ++i){
564 #ifdef RESID_DEBUG
565 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);
566 #endif
567 sum = - gt->n * pow(tau,gt->t) * pow(delta, -1. + gt->d)
568 * (2. * gt->alpha * delta * (delta - gt->epsilon) - gt->d)
569 * exp(-(gt->alpha * SQ(delta-gt->epsilon) + gt->beta*SQ(tau-gt->gamma)));
570 res += sum;
571 #ifdef RESID_DEBUG
572 fprintf(stderr,"sum = %f --> res = %f\n",sum,res);
573 #endif
574 ++gt;
575 }
576
577 /* critical terms */
578 n = data->nc;
579 ct = &(data->ct[0]);
580 for(i=0; i<n; ++i){
581 #ifdef RESID_DEBUG
582 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);
583 #endif
584 DEFINE_DELTA;
585 DEFINE_DPSIDDELTA;
586 DEFINE_DDELDDELTA;
587 DEFINE_DDELBDDELTA;
588
589 sum = ct->n * (pow(DELTA, ct->b) * (psi + delta * dpsiddelta) + dDELbddelta * delta * psi);
590 res += sum;
591 ++ct;
592 }
593
594
595 return res;
596 }
597
598 /**
599 Derivative of the helmholtz residual function with respect to
600 tau.
601 */
602 double helm_resid_tau(double tau,double delta,const HelmholtzData *data){
603
604 double sum;
605 double res = 0;
606 double delX;
607 unsigned l;
608 unsigned n, i;
609 const HelmholtzPowTerm *pt;
610 const HelmholtzGausTerm *gt;
611 const HelmholtzCritTerm *ct;
612
613 n = data->np;
614 pt = &(data->pt[0]);
615
616 delX = 1;
617
618 l = 0;
619 sum = 0;
620 for(i=0; i<n; ++i){
621 if(pt->t){
622 //fprintf(stderr,"i = %d, a = %e, t = %f, d = %d, l = %d\n",i+1, pt->a, pt->t, pt->d, pt->l);
623 sum += pt->a * pow(tau, pt->t - 1) * ipow(delta, pt->d) * pt->t;
624 }
625 ++pt;
626 //fprintf(stderr,"l = %d\n",l);
627 if(i+1==n || l != pt->l){
628 if(l==0){
629 //fprintf(stderr,"Adding non-exp term\n");
630 res += sum;
631 }else{
632 //fprintf(stderr,"Adding exp term with l = %d, delX = %e\n",l,delX);
633 res += sum * exp(-delX);
634 }
635 /* set l to new value */
636 if(i+1!=n){
637 l = pt->l;
638 //fprintf(stderr,"New l = %d\n",l);
639 delX = ipow(delta,l);
640 sum = 0;
641 }
642 }
643 }
644
645 //#define RESID_DEBUG
646 /* gaussian terms */
647 n = data->ng;
648 gt = &(data->gt[0]);
649 for(i=0; i<n; ++i){
650 #ifdef RESID_DEBUG
651 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);
652 #endif
653
654 double val2;
655 val2 = -gt->n * pow(tau,gt->t - 1.) * pow(delta, gt->d)
656 * (2. * gt->beta * tau * (tau - gt->gamma) - gt->t)
657 * exp(-(gt->alpha * SQ(delta-gt->epsilon) + gt->beta*SQ(tau-gt->gamma)));
658 res += val2;
659 #ifdef RESID_DEBUG
660 fprintf(stderr,"res = %f\n",res);
661 #endif
662
663 ++gt;
664 }
665
666 /* critical terms */
667 n = data->nc;
668 ct = &(data->ct[0]);
669 for(i=0; i<n; ++i){
670 #ifdef RESID_DEBUG
671 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);
672 #endif
673 DEFINE_DELTA;
674 DEFINE_DDELBDTAU;
675 DEFINE_DPSIDTAU;
676
677 sum = ct->n * delta * (dDELbdtau * psi + pow(DELTA, ct->b) * dpsidtau);
678 res += sum;
679 ++ct;
680 }
681
682 return res;
683 }
684
685
686 /*=================== SECOND DERIVATIVES =======================*/
687
688 /**
689 Mixed derivative of the helmholtz residual function with respect to
690 delta and tau.
691 */
692 double helm_resid_deltau(double tau,double delta,const HelmholtzData *data){
693 double dell,ldell, term, sum = 0, res = 0;
694 unsigned n, i;
695 const HelmholtzPowTerm *pt;
696 const HelmholtzGausTerm *gt;
697 const HelmholtzCritTerm *ct;
698
699 /* power terms */
700 n = data->np;
701 pt = &(data->pt[0]);
702 dell = ipow(delta,pt->l);
703 ldell = pt->l * dell;
704 unsigned oldl;
705 for(i=0; i<n; ++i){
706 sum += pt->a * pt->t * pow(tau, pt->t - 1) * ipow(delta, pt->d - 1) * (pt->d - ldell);
707 oldl = pt->l;
708 ++pt;
709 if(i+1==n || oldl != pt->l){
710 if(oldl == 0){
711 res += sum;
712 }else{
713 res += sum * exp(-dell);
714 }
715 sum = 0;
716 dell = ipow(delta,pt->l);
717 ldell = pt->l*dell;
718 }
719 }
720
721 #ifdef TEST
722 assert(!isinf(res));
723 #endif
724
725 /* gaussian terms */
726 n = data->ng;
727 gt = &(data->gt[0]);
728 for(i=0; i<n; ++i){
729 #ifdef RESID_DEBUG
730 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);
731 #endif
732 double d1 = delta - gt->epsilon;
733 double t1 = tau - gt->gamma;
734 double e1 = -gt->alpha*SQ(d1) - gt->beta*SQ(t1);
735
736 double f1 = gt->t - 2*gt->beta*tau*(tau - gt->gamma);
737 double g1 = gt->d - 2*gt->alpha*delta*(delta - gt->epsilon);
738
739 sum = gt->n * f1 * pow(tau,gt->t-1) * g1 * pow(delta,gt->d-1) * exp(e1);
740
741 //fprintf(stderr,"sum = %f\n",sum);
742 res += sum;
743 #ifdef TEST
744 assert(!isinf(res));
745 #endif
746 ++gt;
747 }
748
749 /* critical terms */
750 n = data->nc;
751 ct = &(data->ct[0]);
752 for(i=0; i<n; ++i){
753 #ifdef RESID_DEBUG
754 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);
755 #endif
756 DEFINE_DELTA;
757 DEFINE_DPSIDDELTA;
758 DEFINE_DDELBDTAU;
759 DEFINE_DDELDDELTA;
760
761 double d2DELbddeldtau = -ct->A * ct->b * 2./ct->beta * pow(DELTA,ct->b-1)*d1*pow(d12,0.5/ct->beta-1) \
762 - 2. * theta * ct->b * (ct->b - 1) * pow(DELTA,ct->b-2.) * dDELddelta;
763
764 double d2psiddeldtau = 4. * ct->C*ct->D*d1*t1*psi;
765
766 DEFINE_DPSIDTAU;
767
768 sum = ct->n * (pow(DELTA, ct->b) * (dpsidtau + delta * d2psiddeldtau) \
769 + delta *dDELbdtau*dpsidtau \
770 + dDELbdtau*(psi+delta*dpsiddelta) \
771 + d2DELbddeldtau*delta*psi
772 );
773 res += sum;
774 ++ct;
775 }
776
777 #ifdef RESID_DEBUG
778 fprintf(stderr,"phir = %f\n",res);
779 #endif
780
781 #ifdef TEST
782 assert(!isnan(res));
783 assert(!isinf(res));
784 #endif
785 return res;
786 }
787
788 /**
789 Second derivative of helmholtz residual function with respect to
790 delta (twice).
791 */
792 double helm_resid_deldel(double tau,double delta,const HelmholtzData *data){
793 double sum = 0, res = 0;
794 double dell, ldell;
795 unsigned n, i;
796 const HelmholtzPowTerm *pt;
797 const HelmholtzGausTerm *gt;
798 const HelmholtzCritTerm *ct;
799
800 #ifdef RESID_DEBUG
801 fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
802 #endif
803
804 /* power terms */
805 n = data->np;
806 pt = &(data->pt[0]);
807 dell = ipow(delta,pt->l);
808 ldell = pt->l * dell;
809 unsigned oldl;
810 for(i=0; i<n; ++i){
811 double lpart = pt->l ? SQ(ldell) + ldell*(1. - 2*pt->d - pt->l) : 0;
812 sum += pt->a * pow(tau, pt->t) * ipow(delta, pt->d - 2) * (pt->d*(pt->d - 1) + lpart);
813 oldl = pt->l;
814 ++pt;
815 if(i+1==n || oldl != pt->l){
816 if(oldl == 0){
817 res += sum;
818 }else{
819 res += sum * exp(-dell);
820 }
821 sum = 0;
822 dell = ipow(delta,pt->l);
823 ldell = pt->l*dell;
824 }
825 }
826
827 /* gaussian terms */
828 n = data->ng;
829 //fprintf(stderr,"THERE ARE %d GAUSSIAN TERMS\n",n);
830 gt = &(data->gt[0]);
831 for(i=0; i<n; ++i){
832 double s1 = SQ(delta - gt->epsilon);
833 double f1 = gt->d*(gt->d - 1)
834 + 2.*gt->alpha*delta * (
835 delta * (2. * gt->alpha * s1 - 1)
836 - 2. * gt->d * (delta - gt->epsilon)
837 );
838 res += gt->n * pow(tau,gt->t) * pow(delta, gt->d - 2.)
839 * f1
840 * exp(-(gt->alpha * s1 + gt->beta*SQ(tau-gt->gamma)));
841 ++gt;
842 }
843
844 /* critical terms */
845 n = data->nc;
846 ct = &(data->ct[0]);
847 for(i=0; i<n; ++i){
848 #ifdef RESID_DEBUG
849 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);
850 #endif
851
852 DEFINE_DELTA;
853 DEFINE_DPSIDDELTA;
854 DEFINE_DDELDDELTA;
855 DEFINE_DDELBDDELTA;
856
857 double powd12bm1 = pow(d12,0.5/ct->beta-1.);
858
859 double d2psiddel2 = (2.*ct->C*d12 - 1.)*2.*ct->C*psi;
860
861 double d2DELddel2 = 1./d1*dDELddelta + d12*(
862 4.*ct->B*ct->a*(ct->a-1.)*pow(d12,ct->a-2.) \
863 + 2.*SQ(ct->A)*SQ(1./ct->beta)*SQ(powd12bm1) \
864 + ct->A*theta*4./ct->beta*(0.5/ct->beta-1.)*powd12bm1/d12
865 );
866
867 double d2DELbddel2 = ct->b * (pow(DELTA,ct->b - 1.)*d2DELddel2 + (ct->b-1.)*pow(DELTA,ct->b-2.)*SQ(dDELddelta));
868
869 sum = ct->n * (pow(DELTA,ct->b)*(2.*dpsiddelta + delta*d2psiddel2) + 2.*dDELbddelta*(psi+delta*dpsiddelta) + d2DELbddel2*delta*psi);
870
871 res += sum;
872 ++ct;
873 }
874
875 return res;
876 }
877
878
879
880 /**
881 Residual part of helmholtz function.
882 */
883 double helm_resid_tautau(double tau, double delta, const HelmholtzData *data){
884 double dell,ldell, term, sum, res = 0;
885 unsigned n, i;
886 const HelmholtzPowTerm *pt;
887 const HelmholtzGausTerm *gt;
888 const HelmholtzCritTerm *ct;
889
890 n = data->np;
891 pt = &(data->pt[0]);
892
893 #ifdef RESID_DEBUG
894 fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
895 #endif
896
897 /* power terms */
898 sum = 0;
899 dell = ipow(delta,pt->l);
900 ldell = pt->l * dell;
901 unsigned oldl;
902 for(i=0; i<n; ++i){
903 term = pt->a * pt->t * (pt->t - 1) * pow(tau, pt->t - 2) * ipow(delta, pt->d);
904 sum += term;
905 #ifdef RESID_DEBUG
906 fprintf(stderr,"i = %d, a=%e, t=%f, d=%d, term = %f, sum = %f",i,pt->a,pt->t,pt->d,term,sum);
907 if(pt->l==0){
908 fprintf(stderr,",row=%e\n",term);
909 }else{
910 fprintf(stderr,",row=%e\n,",term*exp(-dell));
911 }
912 #endif
913 oldl = pt->l;
914 ++pt;
915 if(i+1==n || oldl != pt->l){
916 if(oldl == 0){
917 #ifdef RESID_DEBUG
918 fprintf(stderr,"linear ");
919 #endif
920 res += sum;
921 }else{
922 #ifdef RESID_DEBUG
923 fprintf(stderr,"exp dell=%f, exp(-dell)=%f sum=%f: ",dell,exp(-dell),sum);
924 #endif
925 res += sum * exp(-dell);
926 }
927 #ifdef RESID_DEBUG
928 fprintf(stderr,"i = %d, res = %f\n",i,res);
929 #endif
930 sum = 0;
931 dell = ipow(delta,pt->l);
932 ldell = pt->l*dell;
933 }
934 }
935
936 /* gaussian terms */
937 n = data->ng;
938 //fprintf(stderr,"THERE ARE %d GAUSSIAN TERMS\n",n);
939 gt = &(data->gt[0]);
940 for(i=0; i<n; ++i){
941 #ifdef RESID_DEBUG
942 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);
943 #endif
944 double d1 = delta - gt->epsilon;
945 double t1 = tau - gt->gamma;
946 double f1 = gt->t*(gt->t - 1) + 4. * gt->beta * tau * (tau * (gt->beta*SQ(t1) - 0.5) - t1*gt->t);
947 double e1 = -gt->alpha*SQ(d1) - gt->beta*SQ(t1);
948 sum = gt->n * f1 * pow(tau,gt->t - 2) * pow(delta,gt->d) * exp(e1);
949 //fprintf(stderr,"sum = %f\n",sum);
950 res += sum;
951 ++gt;
952 }
953
954 /* critical terms */
955 n = data->nc;
956 ct = &(data->ct[0]);
957 for(i=0; i<n; ++i){
958 #ifdef RESID_DEBUG
959 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);
960 #endif
961 DEFINE_DELTA;
962 DEFINE_DDELBDTAU;
963 DEFINE_DPSIDTAU;
964
965 double d2DELbdtau2 = 2. * ct->b * pow(DELTA, ct->b - 1) + 4. * SQ(theta) * ct->b * (ct->b - 1) * pow(DELTA, ct->b - 2);
966
967 double d2psidtau2 = 2. * ct->D * psi * (2. * ct->D * SQ(t1) -1.);
968
969 sum = ct->n * delta * (d2DELbdtau2 * psi + 2 * dDELbdtau*dpsidtau + pow(DELTA, ct->b) * d2psidtau2);
970 res += sum;
971 ++ct;
972 }
973
974 #ifdef RESID_DEBUG
975 fprintf(stderr,"phir_tautau = %f\n",res);
976 #endif
977 return res;
978 }
979

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