/[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 2233 - (show annotations) (download) (as text)
Wed Jul 28 09:09:11 2010 UTC (13 years, 10 months ago) by jpye
File MIME type: text/x-csrc
File size: 35084 byte(s)
helm_resid_deldeldel implementation complete... now for testing.
1 /* ASCEND modelling environment
2 Copyright (C) 2008-2009 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 /* macros and forward decls */
39
40 #define SQ(X) ((X)*(X))
41
42 #include "helmholtz_impl.h"
43
44 /* calculate tau and delta using a macro -- is used in most functions */
45 #define DEFINE_TD \
46 double tau = data->T_star / T; \
47 double delta = rho / data->rho_star
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 DEFINE_TD;
59
60 #ifdef TEST
61 assert(data->rho_star!=0);
62 assert(T!=0);
63 assert(!isnan(tau));
64 assert(!isnan(delta));
65 assert(!isnan(data->R));
66
67 //fprintf(stderr,"p calc: T = %f\n",T);
68 //fprintf(stderr,"p calc: tau = %f\n",tau);
69 //fprintf(stderr,"p calc: rho = %f\n",rho);
70 //fprintf(stderr,"p calc: delta = %f\n",delta);
71 //fprintf(stderr,"p calc: R*T*rho = %f\n",data->R * T * rho);
72
73 //fprintf(stderr,"T = %f\n", T);
74 //fprintf(stderr,"rhob = %f, rhob* = %f, delta = %f\n", rho/data->M, data->rho_star/data->M, delta);
75 #endif
76
77 return data->R * T * rho * (1 + delta * helm_resid_del(tau,delta,data));
78 }
79
80 /**
81 Function to calculate internal energy from Helmholtz free energy EOS, given
82 temperature and mass density.
83
84 @param T temperature in K
85 @param rho mass density in kg/m続
86 @return internal energy in ???
87 */
88 double helmholtz_u(double T, double rho, const HelmholtzData *data){
89 DEFINE_TD;
90
91 #ifdef TEST
92 assert(data->rho_star!=0);
93 assert(T!=0);
94 assert(!isnan(tau));
95 assert(!isnan(delta));
96 assert(!isnan(data->R));
97 #endif
98
99 #if 0
100 fprintf(stderr,"ideal_tau = %f\n",helm_ideal_tau(tau,delta,data->ideal));
101 fprintf(stderr,"resid_tau = %f\n",helm_resid_tau(tau,delta,data));
102 fprintf(stderr,"R T = %f\n",data->R * data->T_star);
103 #endif
104
105 return data->R * data->T_star * (helm_ideal_tau(tau,delta,data->ideal) + helm_resid_tau(tau,delta,data));
106 }
107
108 /**
109 Function to calculate enthalpy from Helmholtz free energy EOS, given
110 temperature and mass density.
111
112 @param T temperature in K
113 @param rho mass density in kg/m続
114 @return enthalpy in J/kg
115 */
116 double helmholtz_h(double T, double rho, const HelmholtzData *data){
117 DEFINE_TD;
118
119 #ifdef TEST
120 assert(data->rho_star!=0);
121 assert(T!=0);
122 assert(!isnan(tau));
123 assert(!isnan(delta));
124 assert(!isnan(data->R));
125 #endif
126
127 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));
128 }
129
130 /**
131 Function to calculate entropy from Helmholtz free energy EOS, given
132 temperature and mass density.
133
134 @param T temperature in K
135 @param rho mass density in kg/m続
136 @return entropy in J/kgK
137 */
138 double helmholtz_s(double T, double rho, const HelmholtzData *data){
139 DEFINE_TD;
140
141 #ifdef ENTROPY_DEBUG
142 assert(data->rho_star!=0);
143 assert(T!=0);
144 assert(!isnan(tau));
145 assert(!isnan(delta));
146 assert(!isnan(data->R));
147
148 fprintf(stderr,"helm_ideal_tau = %f\n",helm_ideal_tau(tau,delta,data->ideal));
149 fprintf(stderr,"helm_resid_tau = %f\n",helm_resid_tau(tau,delta,data));
150 fprintf(stderr,"helm_ideal = %f\n",helm_ideal(tau,delta,data->ideal));
151 fprintf(stderr,"helm_resid = %f\n",helm_resid(tau,delta,data));
152 #endif
153 return data->R * (
154 tau * (helm_ideal_tau(tau,delta,data->ideal) + helm_resid_tau(tau,delta,data))
155 - (helm_ideal(tau,delta,data->ideal) + helm_resid(tau,delta,data))
156 );
157 }
158
159 /**
160 Function to calculate Helmholtz energy from the Helmholtz free energy EOS,
161 given temperature and mass density.
162
163 @param T temperature in K
164 @param rho mass density in kg/m続
165 @return Helmholtz energy 'a', in J/kg
166 */
167 double helmholtz_a(double T, double rho, const HelmholtzData *data){
168 DEFINE_TD;
169
170 #ifdef TEST
171 assert(data->rho_star!=0);
172 assert(T!=0);
173 assert(!isnan(tau));
174 assert(!isnan(delta));
175 assert(!isnan(data->R));
176 #endif
177
178 #ifdef HELMHOLTZ_DEBUG
179 fprintf(stderr,"helmholtz_a: T = %f, rho = %f\n",T,rho);
180 fprintf(stderr,"multiplying by RT = %f\n",data->R*T);
181 #endif
182
183 return data->R * T * (helm_ideal(tau,delta,data->ideal) + helm_resid(tau,delta,data));
184 }
185
186 /**
187 Function to calculate isochoric heat capacity from the Helmholtz free energy
188 EOS given temperature and mass density.
189
190 @param T temperature in K
191 @param rho mass density in kg/m続
192 @return Isochoric specific heat capacity in J/kg/K.
193 */
194 double helmholtz_cv(double T, double rho, const HelmholtzData *data){
195 DEFINE_TD;
196
197 return - data->R * SQ(tau) * (helm_ideal_tautau(tau,data->ideal) + helm_resid_tautau(tau,delta,data));
198 }
199
200 /**
201 Function to calculate isobaric heat capacity from the Helmholtz free energy
202 EOS given temperature and mass density.
203
204 @param T temperature in K
205 @param rho mass density in kg/m続
206 @return Isobaric specific heat capacity in J/kg/K.
207 */
208 double helmholtz_cp(double T, double rho, const HelmholtzData *data){
209 DEFINE_TD;
210
211 double phir_d = helm_resid_del(tau,delta,data);
212 double phir_dd = helm_resid_deldel(tau,delta,data);
213 double phir_dt = helm_resid_deltau(tau,delta,data);
214
215 /* note similarities with helmholtz_w */
216 double temp1 = 1 + 2*delta*phir_d + SQ(delta)*phir_dd;
217 double temp2 = 1 + delta*phir_d - delta*tau*phir_dt;
218 double temp3 = -SQ(tau)*(helm_ideal_tautau(tau,data->ideal) + helm_resid_tautau(tau,delta,data));
219
220 return data->R * (temp3 + SQ(temp2)/temp1);
221 }
222
223
224 /**
225 Function to calculate the speed of sound in a fluid from the Helmholtz free
226 energy EOS, given temperature and mass density.
227
228 @param T temperature in K
229 @param rho mass density in kg/m続
230 @return Speed of sound in m/s.
231 */
232 double helmholtz_w(double T, double rho, const HelmholtzData *data){
233 DEFINE_TD;
234
235 double phir_d = helm_resid_del(tau,delta,data);
236 double phir_dd = helm_resid_deldel(tau,delta,data);
237 double phir_dt = helm_resid_deltau(tau,delta,data);
238
239 /* note similarities with helmholtz_cp */
240 double temp1 = 1. + 2.*delta*phir_d + SQ(delta)*phir_dd;
241 double temp2 = 1. + delta*phir_d - delta*tau*phir_dt;
242 double temp3 = -SQ(tau)*(helm_ideal_tautau(tau,data->ideal) + helm_resid_tautau(tau,delta,data));
243
244 return sqrt(data->R * T * (temp1 + SQ(temp2)/temp3));
245
246 }
247
248 /**
249 Function to calculate the Gibbs energy fluid from the Helmholtz free
250 energy EOS, given temperature and mass density.
251
252 @param T temperature in K
253 @param rho mass density in kg/m続
254 @return Gibbs energy, in J/kg.
255 */
256 double helmholtz_g(double T, double rho, const HelmholtzData *data){
257 DEFINE_TD;
258
259 double phir_d = helm_resid_del(tau,delta,data);
260 double phir = helm_resid(tau,delta,data);
261 double phi0 = helm_ideal(tau,delta,data->ideal);
262
263 return data->R * T * (phi0 + phir + 1. - delta * phir_d);
264 }
265
266 /**
267 Calculation zero-pressure specific heat capacity
268 */
269 double helmholtz_cp0(double T, const HelmholtzData *data){
270 double val = helm_cp0(T,data->ideal);
271 #if 0
272 fprintf(stderr,"val = %f\n",val);
273 #endif
274 return val;
275 }
276
277 /**
278 Solve density given temperature and pressure and bounds for density,
279 used in solving the saturation curve. We use a single-variable Newton
280 method to solve for the desired pressure by iteratively adjusting density.
281
282 @param rho_l lower bound on density
283 @param rho_u upper bound on density
284 @param rho (returned) solved value of density
285 @return 0 on success (1 = didn't converge within permitted iterations)
286 */
287 static int helm_find_rho_tp(double T, double p, double *rho
288 , const double rho_l, const double rho_u, const HelmholtzData *data
289 ){
290 double rho_1 = 0.5 *(rho_l + rho_u);
291 double p_1, dpdrho;
292 double niter = 0, maxiter = 100;
293 int res = 1; /* 1 = exceeded iterations, 0 = converged */
294
295 #ifdef TEST
296 fprintf(stderr,"\nHELM_FIND_RHO_TP: rho_l = %f, rho_u = %f (v_u = %f, v_l = %f)\n", rho_l, rho_u, 1./rho_l, 1./rho_u);
297 #endif
298
299 #if 0
300 double r;
301 for(r = rho_l; r <= rho_u; r += (rho_u-rho_l)/20.){
302 p_1 = helmholtz_p(T, r, data);
303 fprintf(stderr,"T = %f, rho = %f --> p = %f MPa\n",T, r, p_1/1e6);
304 }
305 #endif
306
307 while(res){
308 if(++niter>maxiter)break;
309 p_1 = helmholtz_p(T, rho_1, data);
310 if(p_1 < p){
311 p_1 = 0.5*(p_1 + p);
312 }
313 //fprintf(stderr,"T = %f, rho = %f --> p = %f MPa, err = %f %%\n",T, rho_1, p_1/1e6, (p - p_1)/p *100);
314 dpdrho = helmholtz_dpdrho_T(T, rho_1, data);
315 if(dpdrho < 0){
316 if(rho_l < data->rho_star){
317 /* looking for gas solution */
318 rho_1 = 0.5*(*rho + rho_l);
319 }
320 }
321 rho_1 += (p - p_1)/dpdrho;
322 if(rho_1 > rho_u){
323 rho_1 = 0.5*(*rho + rho_u);
324 //fprintf(stderr,"UPPERBOUND ");
325 }
326 if(rho_1 < rho_l){
327 rho_1 = 0.5*(*rho + rho_l);
328 //fprintf(stderr,"LOWERBOUND ");
329 }
330 *rho = rho_1;
331
332 if(fabs((p - p_1)/p) < 1e-7){
333 #ifdef TEST
334 fprintf(stderr,"Converged to p = %f MPa with rho = %f\n", p/1e6, rho_1);
335 #endif
336 res = 0;
337 }
338 }
339
340 return res;
341 }
342
343 /**
344 Calculation of saturation conditions given temperature
345
346 @param T temperature [K]
347 @param rho_f (returned) saturated liquid density [kg/m続]
348 @param rho_g (returned) saturated gas density [kg/m続]
349 @param p pressure p_sat(T) [Pa]
350 @return 0 on success
351 */
352 int helmholtz_sat_t(double T, double *p, double *rho_f, double *rho_g, const HelmholtzData *data){
353 double tau = data->T_star / T;
354
355 if(T >= data->T_star){
356 #ifdef TEST
357 fprintf(stderr,"ERROR: temperature exceeds critical temperature in helmholtz_sat_t.\n");
358 #endif
359 /* return some reasonable values */
360 *rho_f = data->rho_star;
361 *rho_g = data->rho_star;
362 *p = helmholtz_p(data->T_star, data->rho_star, data);
363 return 1; /* error status */
364 }
365
366 /* get a first estimate of saturation pressure using acentric factor */
367
368 /* critical pressure */
369 #ifdef TEST
370 fprintf(stderr,"T_c = %f, rho_c = %f\n",data->T_star, data->rho_star);
371 #endif
372 double p_c = helmholtz_p(data->T_star, data->rho_star, data);
373
374 #ifdef TEST
375 fprintf(stderr,"Critical pressure = %f MPa\n",p_c/1.e6);
376 fprintf(stderr,"Acentric factor = %f\n",data->omega);
377 #endif
378
379
380 /* maybe use http://dx.doi.org/10.1016/S0009-2509(02)00017-9 for this part? */
381 /* FIXME need to cite this formula */
382 *p = p_c * pow(10.,(data->omega + 1.)*-7./3.*(tau - 1.));
383
384 #ifdef TEST
385 fprintf(stderr,"Estimated p_sat(T=%f) = %f MPa\n",T,(*p)/1e6);
386 #endif
387
388 if(tau < 1.01){
389 /* close to critical point: need a different approach */
390 #ifdef TEST
391 fprintf(stderr,"ERROR: not implemented\n");
392 #endif
393 return 1;
394 }else{
395 double niter = 0;
396 (void)niter;
397
398 int res = helm_find_rho_tp(T, *p, rho_f, data->rho_star,data->rho_star*10, data);
399 if(res){
400 #ifdef TEST
401 fprintf(stderr,"ERROR: failed to solve rho_f\n");
402 #endif
403 }
404
405 res = helm_find_rho_tp(T, *p, rho_g, 0.001*data->rho_star, data->rho_star, data);
406 if(res){
407 #ifdef TEST
408 fprintf(stderr,"ERROR: failed to solve rho_g\n");
409 #endif
410 }
411
412 #ifdef TEST
413 fprintf(stderr,"p = %f MPa: rho_f = %f, rho_g = %f\n", *p, *rho_f, *rho_g);
414 #endif
415
416 double LHS = *p/data->R/T*(1./(*rho_g) - 1./(*rho_f)) - log((*rho_f)/(*rho_g));
417 double delta_f = (*rho_f) / data->rho_star;
418 double delta_g = (*rho_g) / data->rho_star;
419 double RHS = helm_resid(delta_f,tau,data) - helm_resid(delta_g,tau,data);
420
421 double err = LHS - RHS;
422 (void)err;
423
424 #ifdef TEST
425 fprintf(stderr,"LHS = %f, RHS = %f, err = %f\n",LHS, RHS, err);
426 #endif
427 /* away from critical point... */
428 *rho_f = data->rho_star;
429 *rho_g = data->rho_star;
430 #ifdef TEST
431 fprintf(stderr,"ERROR: not implemented\n");
432 exit(1);
433 #endif
434 return 1;
435 }
436 }
437
438 /*----------------------------------------------------------------------------
439 PARTIAL DERIVATIVES
440 */
441
442 /**
443 Calculate partial derivative of p with respect to T, with rho constant
444 */
445 double helmholtz_dpdT_rho(double T, double rho, const HelmholtzData *data){
446 DEFINE_TD;
447
448 double phir_del = helm_resid_del(tau,delta,data);
449 double phir_deltau = helm_resid_deltau(tau,delta,data);
450 #ifdef TEST
451 assert(!isinf(phir_del));
452 assert(!isinf(phir_deltau));
453 assert(!isnan(phir_del));
454 assert(!isnan(phir_deltau));
455 assert(!isnan(data->R));
456 assert(!isnan(rho));
457 assert(!isnan(tau));
458 #endif
459
460 double res = data->R * rho * (1 + delta*phir_del - delta*tau*phir_deltau);
461
462 #ifdef TEST
463 assert(!isnan(res));
464 assert(!isinf(res));
465 #endif
466 return res;
467 }
468
469 /**
470 Calculate partial derivative of p with respect to rho, with T constant
471 */
472 double helmholtz_dpdrho_T(double T, double rho, const HelmholtzData *data){
473 DEFINE_TD;
474
475 double phir_del = helm_resid_del(tau,delta,data);
476 double phir_deldel = helm_resid_deldel(tau,delta,data);
477 #ifdef TEST
478 assert(!isinf(phir_del));
479 assert(!isinf(phir_deldel));
480 #endif
481 return data->R * T * (1 + 2*delta*phir_del + SQ(delta)*phir_deldel);
482 }
483
484 /**
485 Calculate partial derivative of h with respect to T, with rho constant
486 */
487 double helmholtz_dhdT_rho(double T, double rho, const HelmholtzData *data){
488 DEFINE_TD;
489
490 double phir_del = helm_resid_del(tau,delta,data);
491 double phir_deltau = helm_resid_deltau(tau,delta,data);
492 double phir_tautau = helm_resid_tautau(tau,delta,data);
493 double phi0_tautau = helm_ideal_tautau(tau,data->ideal);
494
495 //fprintf(stderr,"phir_del = %f, phir_deltau = %f, phir_tautau = %f, phi0_tautau = %f\n",phir_del,phir_deltau,phir_tautau,phi0_tautau);
496
497 //return (helmholtz_h(T+0.01,rho,data) - helmholtz_h(T,rho,data)) / 0.01;
498 return data->R * (1. + delta*phir_del - SQ(tau)*(phi0_tautau + phir_tautau) - delta*tau*phir_deltau);
499 }
500
501 /**
502 Calculate partial derivative of h with respect to rho, with T constant
503 */
504 double helmholtz_dhdrho_T(double T, double rho, const HelmholtzData *data){
505 DEFINE_TD;
506
507 double phir_del = helm_resid_del(tau,delta,data);
508 double phir_deltau = helm_resid_deltau(tau,delta,data);
509 double phir_deldel = helm_resid_deldel(tau,delta,data);
510
511 return data->R * T / rho * (tau*delta*(0 + phir_deltau) + delta * phir_del + SQ(delta)*phir_deldel);
512 }
513
514
515 /**
516 Calculate partial derivative of u with respect to T, with rho constant
517 */
518 double helmholtz_dudT_rho(double T, double rho, const HelmholtzData *data){
519 DEFINE_TD;
520
521 double phir_tautau = helm_resid_tautau(tau,delta,data);
522 double phi0_tautau = helm_ideal_tautau(tau,data->ideal);
523
524 return -data->R * SQ(tau) * (phi0_tautau + phir_tautau);
525 }
526
527 /**
528 Calculate partial derivative of u with respect to rho, with T constant
529 */
530 double helmholtz_dudrho_T(double T, double rho, const HelmholtzData *data){
531 DEFINE_TD;
532
533 double phir_deltau = helm_resid_deltau(tau,delta,data);
534
535 return data->R * T / rho * (tau * delta * phir_deltau);
536 }
537
538 /*---------------------------------------------
539 UTILITY FUNCTION(S)
540 */
541
542 /* ipow: public domain by Mark Stephen with suggestions by Keiichi Nakasato */
543 static double ipow(double x, int n){
544 double t = 1.0;
545
546 if(!n)return 1.0; /* At the top. x^0 = 1 */
547
548 if (n < 0){
549 n = -n;
550 x = 1.0/x; /* error if x == 0. Good */
551 } /* ZTC/SC returns inf, which is even better */
552
553 if (x == 0.0)return 0.0;
554
555 do{
556 if(n & 1)t *= x;
557 n /= 2; /* KN prefers if (n/=2) x*=x; This avoids an */
558 x *= x; /* unnecessary but benign multiplication on */
559 }while(n); /* the last pass, but the comparison is always
560 true _except_ on the last pass. */
561
562 return t;
563 }
564
565 //#define RESID_DEBUG
566
567 /* maxima expressions:
568 Psi(delta) := exp(-C*(delta-1)^2 -D*(tau-1)^2);
569 theta(delta) := (1-tau) + A*((delta-1)^2)^(1/(2*beta));
570 Delta(delta):= theta(delta)^2 + B*((delta-1)^2)^a;
571 n*Delta(delta)^b*delta*Psi(delta);
572 diff(%,delta,3);
573 yikes, that's scary! break down into steps.
574 */
575
576 /*
577 We avoid duplication by using the following #defines for common code in
578 calculation of critical terms.
579 */
580 #define DEFINE_DELTA \
581 double d1 = delta - 1.; \
582 double t1 = tau - 1.; \
583 double d12 = SQ(d1); \
584 double theta = (1. - tau) + ct->A * pow(d12, 0.5/ct->beta); \
585 double PSI = exp(-ct->C*d12 - ct->D*SQ(t1)); \
586 double DELTA = SQ(theta) + ct->B* pow(d12, ct->a)
587
588 #define DEFINE_DELB \
589 double DELB = pow(DELTA,ct->b)
590
591 #define DEFINE_DPSIDDELTA \
592 double dPSIddelta = -2. * ct->C * d1 * PSI
593
594 #define DEFINE_DDELDDELTA \
595 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))
596
597 #define DEFINE_DDELBDTAU \
598 double dDELbdtau = -2. * theta * ct->b * (DELB/DELTA)
599
600 #define DEFINE_DPSIDTAU \
601 double dPSIdtau = -2. * ct->D * t1 * PSI
602
603 #define DEFINE_DDELBDDELTA \
604 double dDELbddelta = (DELTA==0?0:ct->b * (DELB/DELTA) * dDELddelta)
605
606 #define DEFINE_D2DELDDELTA2 \
607 double powd12bm1 = pow(d12,0.5/ct->beta-1.); \
608 double d2DELddelta2 = 1./d1*dDELddelta + d12*( \
609 4.*ct->B*ct->a*(ct->a-1.)*pow(d12,ct->a-2.) \
610 + 2.*SQ(ct->A)*SQ(1./ct->beta)*SQ(powd12bm1) \
611 + ct->A*theta*4./ct->beta*(0.5/ct->beta-1.)*powd12bm1/d12 \
612 )
613
614 #define DEFINE_D2DELBDDELTA2 \
615 double d2DELbddelta2 = ct->b * ( (DELB/DELTA)*d2DELddelta2 + (ct->b-1.)*(DELB/SQ(DELTA)*SQ(dDELddelta)))
616
617 #define DEFINE_D2PSIDDELTA2 \
618 double d2PSIddelta2 = (2.*ct->C*d12 - 1.)*2.*ct->C * PSI
619
620 #define DEFINE_D3PSIDDELTA3 \
621 double d3PSIddelta3 = -4. * d1 * SQ(ct->C) * (2.*d12*ct->C - 3.) * PSI
622
623 #define DEFINE_D3DELDDELTA3 \
624 double d3DELddelta3 = 1./(d1*d12*ct->beta*SQ(ct->beta)) * (\
625 4*ct->B*ct->a*(1.+ct->a*(2*ct->a-3))*SQ(ct->beta)*pow(d12,ct->a)\
626 + ct->A * (1.+ct->beta*(2*ct->beta-3))*pow(d12,0.5/ct->beta)\
627 )
628
629 #define DEFINE_D3DELBDDELTA3 \
630 double d3DELbddelta3 = ct->b / (DELTA*SQ(DELTA)) * ( \
631 (2+ct->b*(ct->b-3))*dDELddelta*SQ(dDELddelta)*DELB \
632 + DELB*SQ(DELTA)*d3DELddelta3 \
633 + 3*(ct->b-1) * DELB * DELTA * dDELddelta * d2DELddelta2 \
634 )
635
636 /**
637 Residual part of helmholtz function.
638 */
639 double helm_resid(double tau, double delta, const HelmholtzData *data){
640 double dell,ldell, term, sum, res = 0;
641 unsigned n, i;
642 const HelmholtzPowTerm *pt;
643 const HelmholtzGausTerm *gt;
644 const HelmholtzCritTerm *ct;
645
646 n = data->np;
647 pt = &(data->pt[0]);
648
649 #ifdef RESID_DEBUG
650 fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
651 #endif
652
653 /* power terms */
654 sum = 0;
655 dell = ipow(delta,pt->l);
656 ldell = pt->l * dell;
657 unsigned oldl;
658 for(i=0; i<n; ++i){
659 term = pt->a * pow(tau, pt->t) * ipow(delta, pt->d);
660 sum += term;
661 #ifdef RESID_DEBUG
662 fprintf(stderr,"i = %d, a=%e, t=%f, d=%d, term = %f, sum = %f",i,pt->a,pt->t,pt->d,term,sum);
663 if(pt->l==0){
664 fprintf(stderr,",row=%e\n",term);
665 }else{
666 fprintf(stderr,",row=%e\n,",term*exp(-dell));
667 }
668 #endif
669 oldl = pt->l;
670 ++pt;
671 if(i+1==n || oldl != pt->l){
672 if(oldl == 0){
673 #ifdef RESID_DEBUG
674 fprintf(stderr,"linear ");
675 #endif
676 res += sum;
677 }else{
678 #ifdef RESID_DEBUG
679 fprintf(stderr,"exp dell=%f, exp(-dell)=%f sum=%f: ",dell,exp(-dell),sum);
680 #endif
681 res += sum * exp(-dell);
682 }
683 #ifdef RESID_DEBUG
684 fprintf(stderr,"i = %d, res = %f\n",i,res);
685 #endif
686 sum = 0;
687 dell = ipow(delta,pt->l);
688 ldell = pt->l*dell;
689 }
690 }
691
692 /* gaussian terms */
693 n = data->ng;
694 //fprintf(stderr,"THERE ARE %d GAUSSIAN TERMS\n",n);
695 gt = &(data->gt[0]);
696 for(i=0; i<n; ++i){
697 #ifdef RESID_DEBUG
698 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);
699 #endif
700 double d1 = delta - gt->epsilon;
701 double t1 = tau - gt->gamma;
702 double e1 = -gt->alpha*SQ(d1) - gt->beta*SQ(t1);
703 sum = gt->n * pow(tau,gt->t) * pow(delta,gt->d) * exp(e1);
704 //fprintf(stderr,"sum = %f\n",sum);
705 res += sum;
706 ++gt;
707 }
708
709 /* critical terms */
710 n = data->nc;
711 ct = &(data->ct[0]);
712 for(i=0; i<n; ++i){
713 #ifdef RESID_DEBUG
714 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);
715 #endif
716
717 DEFINE_DELTA;
718 DEFINE_DELB;
719
720 sum = ct->n * DELB * delta * PSI;
721 res += sum;
722 ++ct;
723 }
724
725 #ifdef RESID_DEBUG
726 fprintf(stderr,"CALCULATED RESULT FOR phir = %f\n",res);
727 #endif
728 return res;
729 }
730
731 /*=================== FIRST DERIVATIVES =======================*/
732
733 /**
734 Derivative of the helmholtz residual function with respect to
735 delta.
736 */
737 double helm_resid_del(double tau,double delta, const HelmholtzData *data){
738 double sum = 0, res = 0;
739 double dell, ldell;
740 unsigned n, i;
741 const HelmholtzPowTerm *pt;
742 const HelmholtzGausTerm *gt;
743 const HelmholtzCritTerm *ct;
744
745 #ifdef RESID_DEBUG
746 fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
747 #endif
748
749 /* power terms */
750 n = data->np;
751 pt = &(data->pt[0]);
752 dell = ipow(delta,pt->l);
753 ldell = pt->l * dell;
754 unsigned oldl;
755 for(i=0; i<n; ++i){
756 sum += pt->a * pow(tau, pt->t) * ipow(delta, pt->d - 1) * (pt->d - ldell);
757 oldl = pt->l;
758 ++pt;
759 if(i+1==n || oldl != pt->l){
760 if(oldl == 0){
761 res += sum;
762 }else{
763 res += sum * exp(-dell);
764 }
765 sum = 0;
766 dell = ipow(delta,pt->l);
767 ldell = pt->l*dell;
768 }
769 }
770
771 /* gaussian terms */
772 n = data->ng;
773 gt = &(data->gt[0]);
774 for(i=0; i<n; ++i){
775 #ifdef RESID_DEBUG
776 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);
777 #endif
778 sum = - gt->n * pow(tau,gt->t) * pow(delta, -1. + gt->d)
779 * (2. * gt->alpha * delta * (delta - gt->epsilon) - gt->d)
780 * exp(-(gt->alpha * SQ(delta-gt->epsilon) + gt->beta*SQ(tau-gt->gamma)));
781 res += sum;
782 #ifdef RESID_DEBUG
783 fprintf(stderr,"sum = %f --> res = %f\n",sum,res);
784 #endif
785 ++gt;
786 }
787
788 /* critical terms */
789 n = data->nc;
790 ct = &(data->ct[0]);
791 for(i=0; i<n; ++i){
792 #ifdef RESID_DEBUG
793 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);
794 #endif
795 DEFINE_DELTA;
796 DEFINE_DELB;
797 DEFINE_DPSIDDELTA;
798 DEFINE_DDELDDELTA;
799 DEFINE_DDELBDDELTA;
800
801 #if 0
802 if(fabs(dpsiddelta) ==0)fprintf(stderr,"WARNING: dpsiddelta == 0\n");
803 if(fabs(dpsiddelta) ==0)fprintf(stderr,"WARNING: dpsiddelta == 0\n");
804 fprintf(stderr,"psi = %f\n",psi);
805 fprintf(stderr,"DELTA = %f\n",DELTA);
806
807 fprintf(stderr,"dDELddelta = %f\n",dDELddelta);
808 fprintf(stderr,"ct->b - 1. = %f\n",ct->b - 1.);
809 fprintf(stderr,"pow(DELTA,ct->b - 1.) = %f\n",pow(DELTA,ct->b - 1.));
810 assert(!isnan(pow(DELTA,ct->b - 1.)));
811 assert(!isnan(dDELddelta));
812 assert(!isnan(dDELbddelta));
813 //double dDELbddelta = ct->b * pow(DELTA,ct->b - 1.) * dDELddelta
814 fprintf(stderr,"sum = %f\n",sum);
815 if(isnan(sum))fprintf(stderr,"ERROR: sum isnan with i=%d at %d\n",i,__LINE__);
816 #endif
817 sum = ct->n * (DELB * (PSI + delta * dPSIddelta) + dDELbddelta * delta * PSI);
818 res += sum;
819
820 ++ct;
821 }
822
823 return res;
824 }
825
826 /**
827 Derivative of the helmholtz residual function with respect to
828 tau.
829 */
830 double helm_resid_tau(double tau,double delta,const HelmholtzData *data){
831
832 double sum;
833 double res = 0;
834 double delX;
835 unsigned l;
836 unsigned n, i;
837 const HelmholtzPowTerm *pt;
838 const HelmholtzGausTerm *gt;
839 const HelmholtzCritTerm *ct;
840
841 n = data->np;
842 pt = &(data->pt[0]);
843
844 delX = 1;
845
846 l = 0;
847 sum = 0;
848 for(i=0; i<n; ++i){
849 if(pt->t){
850 //fprintf(stderr,"i = %d, a = %e, t = %f, d = %d, l = %d\n",i+1, pt->a, pt->t, pt->d, pt->l);
851 sum += pt->a * pow(tau, pt->t - 1) * ipow(delta, pt->d) * pt->t;
852 }
853 ++pt;
854 //fprintf(stderr,"l = %d\n",l);
855 if(i+1==n || l != pt->l){
856 if(l==0){
857 //fprintf(stderr,"Adding non-exp term\n");
858 res += sum;
859 }else{
860 //fprintf(stderr,"Adding exp term with l = %d, delX = %e\n",l,delX);
861 res += sum * exp(-delX);
862 }
863 /* set l to new value */
864 if(i+1!=n){
865 l = pt->l;
866 //fprintf(stderr,"New l = %d\n",l);
867 delX = ipow(delta,l);
868 sum = 0;
869 }
870 }
871 }
872
873 //#define RESID_DEBUG
874 /* gaussian terms */
875 n = data->ng;
876 gt = &(data->gt[0]);
877 for(i=0; i<n; ++i){
878 #ifdef RESID_DEBUG
879 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);
880 #endif
881
882 double val2;
883 val2 = -gt->n * pow(tau,gt->t - 1.) * pow(delta, gt->d)
884 * (2. * gt->beta * tau * (tau - gt->gamma) - gt->t)
885 * exp(-(gt->alpha * SQ(delta-gt->epsilon) + gt->beta*SQ(tau-gt->gamma)));
886 res += val2;
887 #ifdef RESID_DEBUG
888 fprintf(stderr,"res = %f\n",res);
889 #endif
890
891 ++gt;
892 }
893
894 /* critical terms */
895 n = data->nc;
896 ct = &(data->ct[0]);
897 for(i=0; i<n; ++i){
898 #ifdef RESID_DEBUG
899 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);
900 #endif
901 DEFINE_DELTA;
902 DEFINE_DELB;
903 DEFINE_DDELBDTAU;
904 DEFINE_DPSIDTAU;
905
906 sum = ct->n * delta * (dDELbdtau * PSI + DELB * dPSIdtau);
907 res += sum;
908 ++ct;
909 }
910
911 return res;
912 }
913
914
915 /*=================== SECOND DERIVATIVES =======================*/
916
917 /**
918 Mixed derivative of the helmholtz residual function with respect to
919 delta and tau.
920 */
921 double helm_resid_deltau(double tau,double delta,const HelmholtzData *data){
922 double dell,ldell, term, sum = 0, res = 0;
923 unsigned n, i;
924 const HelmholtzPowTerm *pt;
925 const HelmholtzGausTerm *gt;
926 const HelmholtzCritTerm *ct;
927
928 /* power terms */
929 n = data->np;
930 pt = &(data->pt[0]);
931 dell = ipow(delta,pt->l);
932 ldell = pt->l * dell;
933 unsigned oldl;
934 for(i=0; i<n; ++i){
935 sum += pt->a * pt->t * pow(tau, pt->t - 1) * ipow(delta, pt->d - 1) * (pt->d - ldell);
936 oldl = pt->l;
937 ++pt;
938 if(i+1==n || oldl != pt->l){
939 if(oldl == 0){
940 res += sum;
941 }else{
942 res += sum * exp(-dell);
943 }
944 sum = 0;
945 dell = ipow(delta,pt->l);
946 ldell = pt->l*dell;
947 }
948 }
949
950 #ifdef TEST
951 assert(!isinf(res));
952 #endif
953
954 /* gaussian terms */
955 n = data->ng;
956 gt = &(data->gt[0]);
957 for(i=0; i<n; ++i){
958 #ifdef RESID_DEBUG
959 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);
960 #endif
961 double d1 = delta - gt->epsilon;
962 double t1 = tau - gt->gamma;
963 double e1 = -gt->alpha*SQ(d1) - gt->beta*SQ(t1);
964
965 double f1 = gt->t - 2*gt->beta*tau*(tau - gt->gamma);
966 double g1 = gt->d - 2*gt->alpha*delta*(delta - gt->epsilon);
967
968 sum = gt->n * f1 * pow(tau,gt->t-1) * g1 * pow(delta,gt->d-1) * exp(e1);
969
970 //fprintf(stderr,"sum = %f\n",sum);
971 res += sum;
972 #ifdef TEST
973 assert(!isinf(res));
974 #endif
975 ++gt;
976 }
977
978 /* critical terms */
979 n = data->nc;
980 ct = &(data->ct[0]);
981 for(i=0; i<n; ++i){
982 #ifdef RESID_DEBUG
983 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);
984 #endif
985 DEFINE_DELTA;
986 DEFINE_DELB;
987 DEFINE_DPSIDDELTA;
988 DEFINE_DDELBDTAU;
989 DEFINE_DDELDDELTA;
990
991 double d2DELbddeldtau = -ct->A * ct->b * 2./ct->beta * (DELB/DELTA)*d1*pow(d12,0.5/ct->beta-1) \
992 - 2. * theta * ct->b * (ct->b - 1) * (DELB/SQ(DELTA)) * dDELddelta;
993
994 double d2PSIddeldtau = 4. * ct->C*ct->D*d1*t1*PSI;
995
996 DEFINE_DPSIDTAU;
997
998 sum = ct->n * (DELB * (dPSIdtau + delta * d2PSIddeldtau) \
999 + delta *dDELbdtau*dPSIdtau \
1000 + dDELbdtau*(PSI+delta*dPSIddelta) \
1001 + d2DELbddeldtau*delta*PSI
1002 );
1003 res += sum;
1004 ++ct;
1005 }
1006
1007 #ifdef RESID_DEBUG
1008 fprintf(stderr,"phir = %f\n",res);
1009 #endif
1010
1011 #ifdef TEST
1012 assert(!isnan(res));
1013 assert(!isinf(res));
1014 #endif
1015 return res;
1016 }
1017
1018 /**
1019 Second derivative of helmholtz residual function with respect to
1020 delta (twice).
1021 */
1022 double helm_resid_deldel(double tau,double delta,const HelmholtzData *data){
1023 double sum = 0, res = 0;
1024 double dell, ldell;
1025 unsigned n, i;
1026 const HelmholtzPowTerm *pt;
1027 const HelmholtzGausTerm *gt;
1028 const HelmholtzCritTerm *ct;
1029
1030 #ifdef RESID_DEBUG
1031 fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
1032 #endif
1033
1034 /* power terms */
1035 n = data->np;
1036 pt = &(data->pt[0]);
1037 dell = ipow(delta,pt->l);
1038 ldell = pt->l * dell;
1039 unsigned oldl;
1040 for(i=0; i<n; ++i){
1041 double lpart = pt->l ? SQ(ldell) + ldell*(1. - 2*pt->d - pt->l) : 0;
1042 sum += pt->a * pow(tau, pt->t) * ipow(delta, pt->d - 2) * (pt->d*(pt->d - 1) + lpart);
1043 oldl = pt->l;
1044 ++pt;
1045 if(i+1==n || oldl != pt->l){
1046 if(oldl == 0){
1047 res += sum;
1048 }else{
1049 res += sum * exp(-dell);
1050 }
1051 sum = 0;
1052 dell = ipow(delta,pt->l);
1053 ldell = pt->l*dell;
1054 }
1055 }
1056
1057 /* gaussian terms */
1058 n = data->ng;
1059 //fprintf(stderr,"THERE ARE %d GAUSSIAN TERMS\n",n);
1060 gt = &(data->gt[0]);
1061 for(i=0; i<n; ++i){
1062 double s1 = SQ(delta - gt->epsilon);
1063 double f1 = gt->d*(gt->d - 1)
1064 + 2.*gt->alpha*delta * (
1065 delta * (2. * gt->alpha * s1 - 1)
1066 - 2. * gt->d * (delta - gt->epsilon)
1067 );
1068 res += gt->n * pow(tau,gt->t) * pow(delta, gt->d - 2.)
1069 * f1
1070 * exp(-(gt->alpha * s1 + gt->beta*SQ(tau-gt->gamma)));
1071 ++gt;
1072 }
1073
1074 /* critical terms */
1075 n = data->nc;
1076 ct = &(data->ct[0]);
1077 for(i=0; i<n; ++i){
1078 #ifdef RESID_DEBUG
1079 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);
1080 #endif
1081
1082 DEFINE_DELTA;
1083 DEFINE_DELB;
1084 DEFINE_DPSIDDELTA;
1085 DEFINE_DDELDDELTA;
1086 DEFINE_DDELBDDELTA;
1087
1088 DEFINE_D2DELDDELTA2;
1089 DEFINE_D2DELBDDELTA2;
1090
1091 DEFINE_D2PSIDDELTA2;
1092
1093 sum = ct->n * (DELB*(2.*dPSIddelta + delta*d2PSIddelta2) + 2.*dDELbddelta*(PSI+delta*dPSIddelta) + d2DELbddelta2*delta*PSI);
1094
1095 res += sum;
1096 ++ct;
1097 }
1098
1099 return res;
1100 }
1101
1102
1103
1104 /**
1105 Residual part of helmholtz function.
1106 */
1107 double helm_resid_tautau(double tau, double delta, const HelmholtzData *data){
1108 double dell,ldell, term, sum, res = 0;
1109 unsigned n, i;
1110 const HelmholtzPowTerm *pt;
1111 const HelmholtzGausTerm *gt;
1112 const HelmholtzCritTerm *ct;
1113
1114 n = data->np;
1115 pt = &(data->pt[0]);
1116
1117 #ifdef RESID_DEBUG
1118 fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
1119 #endif
1120
1121 /* power terms */
1122 sum = 0;
1123 dell = ipow(delta,pt->l);
1124 ldell = pt->l * dell;
1125 unsigned oldl;
1126 for(i=0; i<n; ++i){
1127 term = pt->a * pt->t * (pt->t - 1) * pow(tau, pt->t - 2) * ipow(delta, pt->d);
1128 sum += term;
1129 #ifdef RESID_DEBUG
1130 fprintf(stderr,"i = %d, a=%e, t=%f, d=%d, term = %f, sum = %f",i,pt->a,pt->t,pt->d,term,sum);
1131 if(pt->l==0){
1132 fprintf(stderr,",row=%e\n",term);
1133 }else{
1134 fprintf(stderr,",row=%e\n,",term*exp(-dell));
1135 }
1136 #endif
1137 oldl = pt->l;
1138 ++pt;
1139 if(i+1==n || oldl != pt->l){
1140 if(oldl == 0){
1141 #ifdef RESID_DEBUG
1142 fprintf(stderr,"linear ");
1143 #endif
1144 res += sum;
1145 }else{
1146 #ifdef RESID_DEBUG
1147 fprintf(stderr,"exp dell=%f, exp(-dell)=%f sum=%f: ",dell,exp(-dell),sum);
1148 #endif
1149 res += sum * exp(-dell);
1150 }
1151 #ifdef RESID_DEBUG
1152 fprintf(stderr,"i = %d, res = %f\n",i,res);
1153 #endif
1154 sum = 0;
1155 dell = ipow(delta,pt->l);
1156 ldell = pt->l*dell;
1157 }
1158 }
1159
1160 /* gaussian terms */
1161 n = data->ng;
1162 //fprintf(stderr,"THERE ARE %d GAUSSIAN TERMS\n",n);
1163 gt = &(data->gt[0]);
1164 for(i=0; i<n; ++i){
1165 #ifdef RESID_DEBUG
1166 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);
1167 #endif
1168 double d1 = delta - gt->epsilon;
1169 double t1 = tau - gt->gamma;
1170 double f1 = gt->t*(gt->t - 1) + 4. * gt->beta * tau * (tau * (gt->beta*SQ(t1) - 0.5) - t1*gt->t);
1171 double e1 = -gt->alpha*SQ(d1) - gt->beta*SQ(t1);
1172 sum = gt->n * f1 * pow(tau,gt->t - 2) * pow(delta,gt->d) * exp(e1);
1173 //fprintf(stderr,"sum = %f\n",sum);
1174 res += sum;
1175 ++gt;
1176 }
1177
1178 /* critical terms */
1179 n = data->nc;
1180 ct = &(data->ct[0]);
1181 for(i=0; i<n; ++i){
1182 #ifdef RESID_DEBUG
1183 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);
1184 #endif
1185 DEFINE_DELTA;
1186 DEFINE_DELB;
1187 DEFINE_DDELBDTAU;
1188 DEFINE_DPSIDTAU;
1189
1190 double d2DELbdtau2 = 2. * ct->b * (DELB/DELTA) + 4. * SQ(theta) * ct->b * (ct->b - 1) * (DELB/SQ(DELTA));
1191
1192 double d2PSIdtau2 = 2. * ct->D * PSI * (2. * ct->D * SQ(t1) -1.);
1193
1194 sum = ct->n * delta * (d2DELbdtau2 * PSI + 2 * dDELbdtau*dPSIdtau + DELB * d2PSIdtau2);
1195 res += sum;
1196 ++ct;
1197 }
1198
1199 #ifdef RESID_DEBUG
1200 fprintf(stderr,"phir_tautau = %f\n",res);
1201 #endif
1202 return res;
1203 }
1204
1205 /* === THIRD DERIVATIVES (this is getting boring now) === */
1206
1207 /**
1208 Third derivative of helmholtz residual function, with respect to
1209 delta (thrice).
1210 */
1211 double helm_resid_deldeldel(double tau,double delta,const HelmholtzData *data){
1212 double sum = 0, res = 0;
1213 double dell, L, L2;
1214 unsigned n, i;
1215 const HelmholtzPowTerm *pt;
1216 const HelmholtzGausTerm *gt;
1217 const HelmholtzCritTerm *ct;
1218
1219 #ifdef RESID_DEBUG
1220 fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
1221 #endif
1222
1223 /* power terms */
1224 n = data->np;
1225 pt = &(data->pt[0]);
1226 dell = ipow(delta,pt->l);
1227 L = pt->l * dell;
1228 L2 = SQ(L);
1229 unsigned oldl;
1230 for(i=0; i<n; ++i){
1231 double lpart = pt->l
1232 ? L2*L + 3*L2*(1 - pt->d - pt->l) + L*(SQ(pt->l) + 3*(pt->d - 1)*pt->l + 3*SQ(pt->d) - 6*pt->d + 1)
1233 : 0;
1234 sum += pt->a * pow(tau,pt->t) * ipow(delta, pt->d - 3) * (pt->d*(pt->d - 1)*(pt->d - 2) + lpart);
1235 oldl = pt->l;
1236 ++pt;
1237 if(i+1==n || oldl != pt->l){
1238 if(oldl == 0){
1239 res += sum; // note special meaning of l==0 case: no exponential
1240 }else{
1241 res += sum * exp(-dell);
1242 }
1243 sum = 0;
1244 dell = ipow(delta,pt->l);
1245 L = pt->l*dell;
1246 L2 = SQ(L);
1247 }
1248 }
1249
1250 /* gaussian terms */
1251 n = data->ng;
1252 //fprintf(stderr,"THERE ARE %d GAUSSIAN TERMS\n",n);
1253 gt = &(data->gt[0]);
1254 for(i=0; i<n; ++i){
1255 double D = delta - gt->epsilon;
1256 double D2 = SQ(D);
1257 double T2 = SQ(tau - gt->gamma);
1258 double A = gt->alpha * delta;
1259 double A2 = SQ(A);
1260 double d = gt->d;
1261 double d2 = SQ(d);
1262
1263 // this expression calculated from wxMaxima using subsitutions for
1264 // D=delta-epsilon and A=alpha*delta.
1265 double f1 =
1266 - (8*A*A2) * D*D2
1267 + (12*d * A2) * D2
1268 + (12 * delta * A2 + (6*d - 6*d2)*A) * D
1269 - (6 * d * delta * A + d*d2 - 3*d2 + 2*d);
1270
1271 res += gt->n * pow(tau,gt->t) * pow(delta, d - 3.)
1272 * f1
1273 * exp(-(gt->alpha * D2 + gt->beta * T2));
1274 ++gt;
1275 }
1276
1277
1278 #if 1
1279 /* critical terms */
1280 n = data->nc;
1281 ct = &(data->ct[0]);
1282 for(i=0; i<n; ++i){
1283 #ifdef RESID_DEBUG
1284 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);
1285 #endif
1286
1287 DEFINE_DELTA;
1288 DEFINE_DELB;
1289 DEFINE_DPSIDDELTA;
1290 DEFINE_DDELDDELTA;
1291 DEFINE_DDELBDDELTA;
1292
1293 DEFINE_D2PSIDDELTA2;
1294 DEFINE_D2DELDDELTA2;
1295 DEFINE_D2DELBDDELTA2;
1296
1297 DEFINE_D3PSIDDELTA3;
1298 DEFINE_D3DELDDELTA3;
1299 DEFINE_D3DELBDDELTA3;
1300
1301 sum = ct->n * (
1302 delta * (DELB*d3PSIddelta3 + 3 * dDELbddelta * d2PSIddelta2 + 3 * d2DELbddelta2 * dPSIddelta + PSI * d3DELbddelta3)
1303 + 3 * (DELB*d2PSIddelta2 + 2 * dDELbddelta * dPSIddelta + d2DELbddelta2 * PSI)
1304 );
1305
1306 res += sum;
1307 ++ct;
1308 }
1309 #endif
1310
1311 return res;
1312 }
1313
1314
1315
1316

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