/[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 2670 - (show annotations) (download) (as text)
Wed Jan 23 05:58:37 2013 UTC (5 years, 5 months ago) by jpye
File MIME type: text/x-csrc
File size: 40283 byte(s)
Add test of saturation curve convergence for all helmholtz fluids.
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, see <http://www.gnu.org/licenses/>.
16 *//** @file
17 Implementation of the reduced molar Helmholtz free energy equation of state.
18
19 For nomenclature see Tillner-Roth, Harms-Watzenberg and Baehr, Eine neue
20 Fundamentalgleichung f端r Ammoniak.
21
22 John Pye, 29 Jul 2008.
23 */
24
25 #include <math.h>
26 #include <stdio.h>
27 #include <stdlib.h>
28
29 #include "helmholtz.h"
30 #include "ideal_impl.h"
31 #include "sat.h"
32 #include "cp0.h"
33 #include "refstate.h"
34
35
36
37 /* these are the 'raw' functions, they don't do phase equilibrium. */
38 PropEvalFn helmholtz_p;
39 PropEvalFn helmholtz_u;
40 PropEvalFn helmholtz_h;
41 PropEvalFn helmholtz_s;
42 PropEvalFn helmholtz_a;
43 PropEvalFn helmholtz_g;
44 PropEvalFn helmholtz_cp;
45 PropEvalFn helmholtz_cv;
46 PropEvalFn helmholtz_w;
47 PropEvalFn helmholtz_dpdrho_T;
48 PropEvalFn helmholtz_alphap;
49 PropEvalFn helmholtz_betap;
50 SatEvalFn helmholtz_sat;
51
52 double helmholtz_dpdT_rho(double T, double rho, const FluidData *data, FpropsError *err);
53 double helmholtz_d2pdrho2_T(double T, double rho, const FluidData *data, FpropsError *err);
54
55 double helmholtz_dhdT_rho(double T, double rho, const FluidData *data, FpropsError *err);
56 double helmholtz_dhdrho_T(double T, double rho, const FluidData *data, FpropsError *err);
57
58 double helmholtz_dudT_rho(double T, double rho, const FluidData *data, FpropsError *err);
59 double helmholtz_dudrho_T(double T, double rho, const FluidData *data, FpropsError *err);
60
61
62 //#define HELM_DEBUG
63 #define HELM_ERRORS
64 //#define SAT_DEBUG
65
66 #ifdef HELM_DEBUG
67 # include "color.h"
68 # define MSG(FMT, ...) \
69 color_on(stderr,ASC_FG_BRIGHTRED);\
70 fprintf(stderr,"%s:%d: ",__FILE__,__LINE__);\
71 color_on(stderr,ASC_FG_BRIGHTBLUE);\
72 fprintf(stderr,"%s: ",__func__);\
73 color_off(stderr);\
74 fprintf(stderr,FMT "\n",##__VA_ARGS__)
75 #else
76 # define MSG(ARGS...) ((void)0)
77 #endif
78
79 /* TODO centralise declaration of our error-reporting function somehow...? */
80 #ifdef HELM_ERRORS
81 # include "color.h"
82 # define ERRMSG(STR,...) \
83 color_on(stderr,ASC_FG_BRIGHTRED);\
84 fprintf(stderr,"ERROR:");\
85 color_off(stderr);\
86 fprintf(stderr," %s:%d:" STR "\n", __func__, __LINE__ ,##__VA_ARGS__)
87 #else
88 # define ERRMSG(ARGS...) ((void)0)
89 #endif
90
91
92 //#define ASSERT_DEBUG
93 #ifdef ASSERT_DEBUG
94 # include <assert.h>
95 #else
96 # define assert(ARGS...)
97 #endif
98
99 //#define RESID_DEBUG
100
101 #define INCLUDE_THIRD_DERIV_CODE
102
103 /* macros and forward decls */
104
105 #define SQ(X) ((X)*(X))
106
107 #include "helmholtz_impl.h"
108
109 /* shortcut take us straight into the correct data structure for Helmholtz correlation calculations */
110 #define HD data->corr.helm
111 #define HD_R data->R
112 #define HD_CP0 data->cp0
113
114 /* calculate tau and delta using a macro -- is used in most functions */
115 #define DEFINE_TD \
116 double tau = data->corr.helm->T_star / T; \
117 double delta = rho / data->corr.helm->rho_star
118
119
120 PureFluid *helmholtz_prepare(const EosData *E, const ReferenceState *ref){
121 PureFluid *P = FPROPS_NEW(PureFluid);
122
123 if(E->type != FPROPS_HELMHOLTZ){
124 fprintf(stderr,"%s: Error: invalid EOS data, wrong type\n",__func__);
125 return NULL;
126 }
127
128 MSG("Fluid '%s' with T_t = %f", E->name, E->data.helm->T_t);
129
130 P->data = FPROPS_NEW(FluidData);
131 P->data->corr.helm = FPROPS_NEW(HelmholtzRunData);
132
133 /* metadata */
134 /* FIXME strings should be copied, not just referenced */
135 P->name = E->name;
136 P->source = E->source;
137 P->type = E->type;
138 MSG("name = %s",P->name);
139
140 /* common data across all correlation types */
141 #define I E->data.helm
142 P->data->M = I->M;
143 if(I->R == 0){
144 P->data->R = R_UNIVERSAL / I->M;
145 }else{
146 P->data->R = I->R;
147 }
148 P->data->T_t = I->T_t;
149 P->data->T_c = I->T_c;
150 P->data->p_c = 0; // we calculate this later...
151 P->data->rho_c = I->rho_c;
152 P->data->omega = I->omega;
153 P->data->cp0 = cp0_prepare(E->data.helm->ideal, P->data->R, P->data->T_c);
154
155 /* data specific to helmholtz correlations */
156 #define H P->data->corr.helm
157 H->rho_star = I->rho_star;
158 H->T_star = I->T_star;
159 H->np = I->np;
160 // FIXME copy et, ct, pt to runtime struct?
161 H->pt = I->pt;
162 H->ng = I->ng;
163 H->gt = I->gt;
164 H->nc = I->nc;
165 H->ct = I->ct;
166 MSG("np = %d, ng = %d, nc = %d, T_t = %f",H->np,H->ng,H->nc,I->T_t);
167 #undef H
168
169 /* function pointers... more to come still? */
170 #define FN(VAR) P->VAR##_fn = &helmholtz_##VAR
171 FN(p); FN(u); FN(h); FN(s); FN(a); FN(g); FN(cp); FN(cv); FN(w);
172 FN(alphap); FN(betap); FN(dpdrho_T);
173 FN(sat);
174 #undef FN
175
176 FpropsError err = 0;
177
178 /* calculate critical pressure (doesn't require h0, s0) */
179 MSG("Calculating critical pressure at T_c = %f K, rho_c = %f kg/m3",P->data->T_c, P->data->rho_c);
180 P->data->p_c = helmholtz_p(P->data->T_c, P->data->rho_c, P->data, &err);
181 if(err){
182 fprintf(stderr,"Failed to calculate critical pressure\n");
183 FPROPS_FREE(P->data);
184 FPROPS_FREE(P->data->corr.helm);
185 return NULL;
186 }
187 if(P->data->p_c <= 0){
188 fprintf(stderr,"Calculated a critical pressure <= 0! (value = %f)\n",P->data->p_c);
189 //return NULL;
190 }
191
192 // fix up the reference point now...
193 if(ref == NULL){
194 // use the provided ReferenceState, or the default one otherwise.
195 ref = &(I->ref);
196 }
197 int res = fprops_set_reference_state(P,ref);
198 if(res){
199 fprintf(stderr,"Unable to apply reference state (type %d, err %d)\n",ref->type,res);
200 return NULL;
201 }
202
203 #undef I
204 return P;
205 }
206
207 /**
208 Function to calculate pressure from Helmholtz free energy EOS, given temperature
209 and mass density.
210
211 @param T temperature in K
212 @param rho mass density in kg/m続
213 @return pressure in Pa
214 */
215 double helmholtz_p(double T, double rho, const FluidData *data, FpropsError *err){
216 DEFINE_TD;
217
218 assert(HD->rho_star!=0);
219 assert(T!=0);
220 assert(!isnan(T));
221 assert(!isnan(rho));
222 assert(!isnan(HD_R));
223
224 //fprintf(stderr,"p calc: T = %f\n",T);
225 //fprintf(stderr,"p calc: tau = %f\n",tau);
226 //fprintf(stderr,"p calc: rho = %f\n",rho);
227 //fprintf(stderr,"p calc: delta = %f\n",delta);
228 //fprintf(stderr,"p calc: R*T*rho = %f\n",HD_R * T * rho);
229
230 //fprintf(stderr,"T = %f\n", T);
231 //fprintf(stderr,"rhob = %f, rhob* = %f, delta = %f\n", rho/HD->M, HD->rho_star/HD->M, delta);
232
233 double p = HD_R * T * rho * (1 + delta * helm_resid_del(tau,delta,HD));
234 #if 0
235 if(isnan(p)){
236 fprintf(stderr,"T = %.12e, rho = %.12e\n",T,rho);
237 }
238 #endif
239 //abort();
240 if(isnan(p))*err = FPROPS_NUMERIC_ERROR;
241 return p;
242 }
243
244 /**
245 Function to calculate internal energy from Helmholtz free energy EOS, given
246 temperature and mass density.
247
248 @param T temperature in K
249 @param rho mass density in kg/m続
250 @return internal energy in ???
251 */
252 double helmholtz_u(double T, double rho, const FluidData *data, FpropsError *err){
253 DEFINE_TD;
254
255 #ifdef TEST
256 assert(HD->rho_star!=0);
257 assert(T!=0);
258 assert(!isnan(tau));
259 assert(!isnan(delta));
260 assert(!isnan(HD_R));
261 #endif
262
263 #if 0
264 fprintf(stderr,"ideal_tau = %f\n",ideal_phi_tau(tau,delta,HD_CP0));
265 fprintf(stderr,"resid_tau = %f\n",helm_resid_tau(tau,delta,HD));
266 fprintf(stderr,"R T = %f\n",HD_R * HD->T_star);
267 #endif
268
269 return HD_R * HD->T_star * (ideal_phi_tau(tau,delta,HD_CP0) + helm_resid_tau(tau,delta,HD));
270 }
271
272 /**
273 Function to calculate enthalpy from Helmholtz free energy EOS, given
274 temperature and mass density.
275
276 @param T temperature in K
277 @param rho mass density in kg/m続
278 @return enthalpy in J/kg
279 */
280 double helmholtz_h(double T, double rho, const FluidData *data, FpropsError *err){
281 DEFINE_TD;
282
283 //#ifdef TEST
284 assert(HD->rho_star!=0);
285 assert(T!=0);
286 assert(!isnan(tau));
287 assert(!isnan(delta));
288 assert(!isnan(HD_R));
289 //#endif
290 double h = HD_R * T * (1 + tau * (ideal_phi_tau(tau,delta,HD_CP0) + helm_resid_tau(tau,delta,HD)) \
291 + delta*helm_resid_del(tau,delta,HD));
292 assert(!isnan(h));
293 return h;
294 }
295
296 /**
297 Function to calculate entropy from Helmholtz free energy EOS, given
298 temperature and mass density.
299
300 @param T temperature in K
301 @param rho mass density in kg/m続
302 @return entropy in J/kgK
303 */
304 double helmholtz_s(double T, double rho, const FluidData *data, FpropsError *err){
305 DEFINE_TD;
306
307 #ifdef ENTROPY_DEBUG
308 assert(HD->rho_star!=0);
309 assert(T!=0);
310 assert(!isnan(tau));
311 assert(!isnan(delta));
312 assert(!isnan(HD_R));
313
314 fprintf(stderr,"ideal_phi_tau = %f\n",ideal_phi_tau(tau,delta,HD_CP0));
315 fprintf(stderr,"helm_resid_tau = %f\n",helm_resid_tau(tau,delta,HD));
316 fprintf(stderr,"ideal_phi = %f\n",ideal_phi(tau,delta,HD_CP0));
317 fprintf(stderr,"helm_resid = %f\n",helm_resid(tau,delta,HD));
318 #endif
319 return HD_R * (
320 tau * (ideal_phi_tau(tau,delta,HD_CP0) + helm_resid_tau(tau,delta,HD))
321 - (ideal_phi(tau,delta,HD_CP0) + helm_resid(tau,delta,HD))
322 );
323 }
324
325 /**
326 Function to calculate Helmholtz energy from the Helmholtz free energy EOS,
327 given temperature and mass density.
328
329 @param T temperature in K
330 @param rho mass density in kg/m続
331 @return Helmholtz energy 'a', in J/kg
332 */
333 double helmholtz_a(double T, double rho, const FluidData *data, FpropsError *err){
334 DEFINE_TD;
335
336 #ifdef TEST
337 assert(HD->rho_star!=0);
338 assert(T!=0);
339 assert(!isnan(tau));
340 assert(!isnan(delta));
341 assert(!isnan(HD_R));
342 #endif
343
344 #ifdef HELMHOLTZ_DEBUG
345 fprintf(stderr,"helmholtz_a: T = %f, rho = %f\n",T,rho);
346 fprintf(stderr,"multiplying by RT = %f\n",HD_R*T);
347 #endif
348
349 return HD_R * T * (ideal_phi(tau,delta,HD_CP0) + helm_resid(tau,delta,HD));
350 }
351
352 /**
353 Function to calculate isochoric heat capacity from the Helmholtz free energy
354 EOS given temperature and mass density.
355
356 @param T temperature in K
357 @param rho mass density in kg/m続
358 @return Isochoric specific heat capacity in J/kg/K.
359 */
360 double helmholtz_cv(double T, double rho, const FluidData *data, FpropsError *err){
361 DEFINE_TD;
362
363 return - HD_R * SQ(tau) * (ideal_phi_tautau(tau,HD_CP0) + helm_resid_tautau(tau,delta,HD));
364 }
365
366 /**
367 Function to calculate isobaric heat capacity from the Helmholtz free energy
368 EOS given temperature and mass density.
369
370 @param T temperature in K
371 @param rho mass density in kg/m続
372 @return Isobaric specific heat capacity in J/kg/K.
373 */
374 double helmholtz_cp(double T, double rho, const FluidData *data, FpropsError *err){
375 DEFINE_TD;
376
377 double phir_d = helm_resid_del(tau,delta,HD);
378 double phir_dd = helm_resid_deldel(tau,delta,HD);
379 double phir_dt = helm_resid_deltau(tau,delta,HD);
380
381 /* note similarities with helmholtz_w */
382 double temp1 = 1 + 2*delta*phir_d + SQ(delta)*phir_dd;
383 double temp2 = 1 + delta*phir_d - delta*tau*phir_dt;
384 double temp3 = -SQ(tau)*(ideal_phi_tautau(tau,HD_CP0) + helm_resid_tautau(tau,delta,HD));
385
386 return HD_R * (temp3 + SQ(temp2)/temp1);
387 }
388
389
390 /**
391 Function to calculate the speed of sound in a fluid from the Helmholtz free
392 energy EOS, given temperature and mass density.
393
394 @param T temperature in K
395 @param rho mass density in kg/m続
396 @return Speed of sound in m/s.
397 */
398 double helmholtz_w(double T, double rho, const FluidData *data, FpropsError *err){
399 DEFINE_TD;
400
401 double phir_d = helm_resid_del(tau,delta,HD);
402 double phir_dd = helm_resid_deldel(tau,delta,HD);
403 double phir_dt = helm_resid_deltau(tau,delta,HD);
404
405 /* note similarities with helmholtz_cp */
406 double temp1 = 1. + 2.*delta*phir_d + SQ(delta)*phir_dd;
407 double temp2 = 1. + delta*phir_d - delta*tau*phir_dt;
408 double temp3 = -SQ(tau)*(ideal_phi_tautau(tau,HD_CP0) + helm_resid_tautau(tau,delta,HD));
409
410 return sqrt(HD_R * T * (temp1 + SQ(temp2)/temp3));
411
412 }
413
414 /**
415 Function to calculate the Gibbs energy fluid from the Helmholtz free
416 energy EOS, given temperature and mass density.
417
418 @param T temperature in K
419 @param rho mass density in kg/m続
420 @return Gibbs energy, in J/kg.
421 */
422 double helmholtz_g(double T, double rho, const FluidData *data, FpropsError *err){
423 DEFINE_TD;
424
425 double phir_d = helm_resid_del(tau,delta,HD);
426 double phir = helm_resid(tau,delta,HD);
427 double phi0 = ideal_phi(tau,delta,HD_CP0);
428
429 return HD_R * T * (phi0 + phir + 1. + delta * phir_d);
430 }
431
432 /**
433 alpha_p function from IAPWS Advisory Note 3, used in calculation of
434 partial property derivatives.
435 */
436 double helmholtz_alphap(double T, double rho, const FluidData *data, FpropsError *err){
437 DEFINE_TD;
438 double phir_d = helm_resid_del(tau,delta,HD);
439 double phir_dt = helm_resid_deltau(tau,delta,HD);
440 return 1./T * (1. - delta*tau*phir_dt/(1 + delta*phir_d));
441 }
442
443 /**
444 beta_p function from IAPWS Advisory Note 3 , used in calculation of partial
445 property derivatives.
446 */
447 double helmholtz_betap(double T, double rho, const FluidData *data, FpropsError *err){
448 DEFINE_TD;
449 double phir_d = helm_resid_del(tau,delta,HD);
450 double phir_dd = helm_resid_deldel(tau,delta,HD);
451 return rho*(1. + (delta*phir_d + SQ(delta)*phir_dd)/(1+delta*phir_d));
452 }
453
454 /*----------------------------------------------------------------------------
455 PARTIAL DERIVATIVES
456 */
457
458 /**
459 Calculate partial derivative of p with respect to T, with rho constant
460 */
461 double helmholtz_dpdT_rho(double T, double rho, const FluidData *data, FpropsError *err){
462 DEFINE_TD;
463
464 double phir_del = helm_resid_del(tau,delta,HD);
465 double phir_deltau = helm_resid_deltau(tau,delta,HD);
466 #ifdef TEST
467 assert(!isinf(phir_del));
468 assert(!isinf(phir_deltau));
469 assert(!isnan(phir_del));
470 assert(!isnan(phir_deltau));
471 assert(!isnan(HD_R));
472 assert(!isnan(rho));
473 assert(!isnan(tau));
474 #endif
475
476 double res = HD_R * rho * (1 + delta*phir_del - delta*tau*phir_deltau);
477
478 #ifdef TEST
479 assert(!isnan(res));
480 assert(!isinf(res));
481 #endif
482 return res;
483 }
484
485 /**
486 Calculate partial derivative of p with respect to rho, with T constant
487 */
488 double helmholtz_dpdrho_T(double T, double rho, const FluidData *data, FpropsError *err){
489 DEFINE_TD;
490 //MSG("...");
491 double phir_del = helm_resid_del(tau,delta,HD);
492 double phir_deldel = helm_resid_deldel(tau,delta,HD);
493 #ifdef TEST
494 assert(!isinf(phir_del));
495 assert(!isinf(phir_deldel));
496 #endif
497 return HD_R * T * (1 + 2*delta*phir_del + SQ(delta)*phir_deldel);
498 }
499
500
501 double helmholtz_d2pdrho2_T(double T, double rho, const FluidData *data, FpropsError *err){
502 DEFINE_TD;
503
504 double phir_del = helm_resid_del(tau,delta,HD);
505 double phir_deldel = helm_resid_deldel(tau,delta,HD);
506 double phir_deldeldel = helm_resid_deldeldel(tau,delta,HD);
507 #ifdef TEST
508 assert(!isinf(phir_del));
509 assert(!isinf(phir_deldel));
510 assert(!isinf(phir_deldeldel));
511 #endif
512
513 return HD_R * T / rho * delta * (2*phir_del + delta*(4*phir_deldel + delta*phir_deldeldel));
514 }
515
516 /**
517 Calculate partial derivative of h with respect to T, with rho constant
518 */
519 double helmholtz_dhdT_rho(double T, double rho, const FluidData *data, FpropsError *err){
520 DEFINE_TD;
521
522 double phir_del = helm_resid_del(tau,delta,HD);
523 double phir_deltau = helm_resid_deltau(tau,delta,HD);
524 double phir_tautau = helm_resid_tautau(tau,delta,HD);
525 double phi0_tautau = ideal_phi_tautau(tau,HD_CP0);
526
527 //fprintf(stderr,"phir_del = %f, phir_deltau = %f, phir_tautau = %f, phi0_tautau = %f\n",phir_del,phir_deltau,phir_tautau,phi0_tautau);
528
529 //return (helmholtz_h(T+0.01,rho,data) - helmholtz_h(T,rho,data)) / 0.01;
530 return HD_R * (1. + delta*phir_del - SQ(tau)*(phi0_tautau + phir_tautau) - delta*tau*phir_deltau);
531 }
532
533 /**
534 Calculate partial derivative of h with respect to rho, with T constant
535 */
536 double helmholtz_dhdrho_T(double T, double rho, const FluidData *data, FpropsError *err){
537 DEFINE_TD;
538
539 double phir_del = helm_resid_del(tau,delta,HD);
540 double phir_deltau = helm_resid_deltau(tau,delta,HD);
541 double phir_deldel = helm_resid_deldel(tau,delta,HD);
542
543 return HD_R * T / rho * (tau*delta*(0 + phir_deltau) + delta * phir_del + SQ(delta)*phir_deldel);
544 }
545
546
547 /**
548 Calculate partial derivative of u with respect to T, with rho constant
549 */
550 double helmholtz_dudT_rho(double T, double rho, const FluidData *data, FpropsError *err){
551 DEFINE_TD;
552
553 double phir_tautau = helm_resid_tautau(tau,delta,HD);
554 double phi0_tautau = ideal_phi_tautau(tau,HD_CP0);
555
556 return -HD_R * SQ(tau) * (phi0_tautau + phir_tautau);
557 }
558
559
560 /**
561 Calculate partial derivative of u with respect to rho, with T constant
562 */
563 double helmholtz_dudrho_T(double T, double rho, const FluidData *data, FpropsError *err){
564 DEFINE_TD;
565
566 double phir_deltau = helm_resid_deltau(tau,delta,HD);
567
568 return HD_R * T / rho * (tau * delta * phir_deltau);
569 }
570
571
572 /**
573 Solve saturation condition for a specified temperature using approach of
574 Akasaka, but adapted for general use to non-helmholtz property correlations.
575 @param T temperature [K]
576 @param psat_out output, saturation pressure [Pa]
577 @param rhof_out output, saturated liquid density [kg/m^3]
578 @param rhog_out output, saturated vapour density [kg/m^3]
579 @param d helmholtz data object for the fluid in question.
580 @return 0 on success, non-zero on error (eg algorithm failed to converge, T out of range, etc.)
581 */
582 double helmholtz_sat(double T, double *rhof_out, double * rhog_out, const FluidData *data, FpropsError *err){
583 if(T < data->T_t - 1e-8){
584 ERRMSG("Input temperature %f K is below triple-point temperature %f K",T,data->T_t);
585 return FPROPS_RANGE_ERROR;
586 }
587
588 if(T > data->T_c + 1e-8){
589 ERRMSG("Input temperature is above critical point temperature");
590 *err = FPROPS_RANGE_ERROR;
591 }
592
593 // we're at the critical point
594 if(fabs(T - data->T_c) < 1e-9){
595 *rhof_out = data->rho_c;
596 *rhog_out = data->rho_c;
597 return data->p_c;
598 }
599
600 // FIXME at present step-length multiplier is set to 0.4 just because of
601 // ONE FLUID, ethanol. Probably our initial guess data isn't good enough,
602 // or maybe there's a problem with the acentric factor or something like
603 // that. This factor 0.4 will be slowing down the whole system, so it's not
604 // good. TODO XXX.
605
606 // initial guesses for liquid and vapour density
607 double rhof = 1.1 * fprops_rhof_T_rackett(T,data);
608 double rhog= 0.9 * fprops_rhog_T_chouaieb(T,data);
609 double R = data->R;
610 double pc = data->p_c;
611
612 #ifdef SAT_DEBUG
613 MSG("initial guess rho_f = %f, rho_g = %f",rhof,rhog);
614 MSG("calculating at T = %.12e",T);
615 #endif
616
617 int i = 0;
618 while(i++ < 200){
619 assert(!isnan(rhog));
620 assert(!isnan(rhof));
621 #ifdef SAT_DEBUG
622 MSG("iter %d: T = %f, rhof = %f, rhog = %f",i,T, rhof, rhog);
623 #endif
624
625 double pf = helmholtz_p(T,rhof,data,err);
626 double pg = helmholtz_p(T,rhog,data,err);
627 double gf = helmholtz_a(T,rhof,data,err) + pf/rhof;
628 double gg = helmholtz_a(T,rhog,data,err) + pg/rhog;
629 double dpdrf = helmholtz_dpdrho_T(T,rhof,data,err);
630 double dpdrg = helmholtz_dpdrho_T(T,rhog,data,err);
631
632 // jacobian for [F;G](rhof, rhog) --- derivatives wrt rhof and rhog
633 double F = (pf - pg)/pc;
634 double G = (gf - gg)/R/T;
635
636 if(fabs(F) + fabs(G) < 1e-12){
637 //fprintf(stderr,"%s: CONVERGED\n",__func__);
638 *rhof_out = rhof;
639 *rhog_out = rhog;
640 return helmholtz_p(T, *rhog_out, data, err);
641 /* SUCCESS */
642 }
643
644 double Ff = dpdrf/pc;
645 double Fg = -dpdrg/pc;
646 //MSG("Ff = %e, Fg = %e",Ff,Fg);
647
648 double Gf = dpdrf/rhof/R/T;
649 double Gg = -dpdrg/rhog/R/T;
650 //MSG("Gf = %e, Gg = %e",Gf,Gg);
651
652 double DET = Ff*Gg - Fg*Gf;
653 //MSG("DET = %f",DET);
654
655 // 'gamma' needs to be increased to 0.5 for water to solve correctly (see 'test/sat.c')
656 #define gamma 1
657 rhof += gamma/DET * (Fg*G - Gg*F);
658 rhog += gamma/DET * ( Gf*F - Ff*G);
659 #undef gamma
660
661 assert(!isnan(rhof));
662 assert(!isnan(rhog));
663
664 if(rhog < 0)rhog = -0.5*rhog;
665 if(rhof < 0)rhof = -0.5*rhof;
666 }
667 *rhof_out = rhof;
668 *rhog_out = rhog;
669 *err = FPROPS_SAT_CVGC_ERROR;
670 ERRMSG("Not converged: with T = %e (rhof=%f, rhog=%f).",T,*rhof_out,*rhog_out);
671 return helmholtz_p(T, rhog, data, err);
672 }
673
674
675
676
677 /*---------------------------------------------
678 UTILITY FUNCTION(S)
679 */
680
681 /* ipow: public domain by Mark Stephen with suggestions by Keiichi Nakasato */
682 static double ipow(double x, int n){
683 double t = 1.0;
684
685 if(!n)return 1.0; /* At the top. x^0 = 1 */
686
687 if (n < 0){
688 n = -n;
689 x = 1.0/x; /* error if x == 0. Good */
690 } /* ZTC/SC returns inf, which is even better */
691
692 if (x == 0.0)return 0.0;
693
694 do{
695 if(n & 1)t *= x;
696 n /= 2; /* KN prefers if (n/=2) x*=x; This avoids an */
697 x *= x; /* unnecessary but benign multiplication on */
698 }while(n); /* the last pass, but the comparison is always
699 true _except_ on the last pass. */
700
701 return t;
702 }
703
704 /* maxima expressions:
705 Psi(delta) := exp(-C*(delta-1)^2 -D*(tau-1)^2);
706 theta(delta) := (1-tau) + A*((delta-1)^2)^(1/(2*beta));
707 Delta(delta):= theta(delta)^2 + B*((delta-1)^2)^a;
708 n*Delta(delta)^b*delta*Psi(delta);
709 diff(%,delta,3);
710 yikes, that's scary! break down into steps.
711 */
712
713 #undef HD
714
715 /*
716 We avoid duplication by using the following #defines for common code in
717 calculation of critical terms.
718 */
719 #define DEFINE_DELTA \
720 double d1 = delta - 1.; \
721 double t1 = tau - 1.; \
722 double d12 = SQ(d1); \
723 double theta = (1. - tau) + ct->A * pow(d12, 0.5/ct->beta); \
724 double PSI = exp(-ct->C*d12 - ct->D*SQ(t1)); \
725 double DELTA = SQ(theta) + ct->B* pow(d12, ct->a)
726
727 #define DEFINE_DELB \
728 double DELB = pow(DELTA,ct->b)
729
730 #define DEFINE_DPSIDDELTA \
731 double dPSIddelta = -2. * ct->C * d1 * PSI
732
733 #define DEFINE_DDELDDELTA \
734 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))
735
736 #define DEFINE_DDELBDTAU \
737 double dDELbdtau = (DELTA == 0) ? 0 : -2. * theta * ct->b * (DELB/DELTA);\
738 assert(!__isnan(dDELbdtau))
739
740 #define DEFINE_DPSIDTAU \
741 double dPSIdtau = -2. * ct->D * t1 * PSI
742
743 #define DEFINE_DDELBDDELTA \
744 double dDELbddelta = (DELTA==0?0:ct->b * (DELB/DELTA) * dDELddelta)
745
746 #define DEFINE_D2DELDDELTA2 \
747 double powd12bm1 = pow(d12,0.5/ct->beta-1.); \
748 double d2DELddelta2 = 1./d1*dDELddelta + d12*( \
749 4.*ct->B*ct->a*(ct->a-1.)*pow(d12,ct->a-2.) \
750 + 2.*SQ(ct->A)*SQ(1./ct->beta)*SQ(powd12bm1) \
751 + ct->A*theta*4./ct->beta*(0.5/ct->beta-1.)*powd12bm1/d12 \
752 )
753
754 #define DEFINE_D2DELBDDELTA2 \
755 double d2DELbddelta2 = ct->b * ( (DELB/DELTA)*d2DELddelta2 + (ct->b-1.)*(DELB/SQ(DELTA)*SQ(dDELddelta)))
756
757 #define DEFINE_D2PSIDDELTA2 \
758 double d2PSIddelta2 = (2.*ct->C*d12 - 1.)*2.*ct->C * PSI
759
760 #define DEFINE_D3PSIDDELTA3 \
761 double d3PSIddelta3 = -4. * d1 * SQ(ct->C) * (2.*d12*ct->C - 3.) * PSI
762
763 #define DEFINE_D3DELDDELTA3 \
764 double d3DELddelta3 = 1./(d1*d12*ct->beta*SQ(ct->beta)) * (\
765 4*ct->B*ct->a*(1.+ct->a*(2*ct->a-3))*SQ(ct->beta)*pow(d12,ct->a)\
766 + ct->A * (1.+ct->beta*(2*ct->beta-3))*pow(d12,0.5/ct->beta)\
767 )
768
769 #define DEFINE_D3DELBDDELTA3 \
770 double d3DELbddelta3 = ct->b / (DELTA*SQ(DELTA)) * ( \
771 (2+ct->b*(ct->b-3))*dDELddelta*SQ(dDELddelta)*DELB \
772 + DELB*SQ(DELTA)*d3DELddelta3 \
773 + 3*(ct->b-1) * DELB * DELTA * dDELddelta * d2DELddelta2 \
774 )
775
776 /**
777 Residual part of helmholtz function.
778 */
779 double helm_resid(double tau, double delta, const HelmholtzRunData *HD){
780 double dell, term, sum, res = 0;
781 unsigned n, i;
782 const HelmholtzPowTerm *pt;
783 const HelmholtzGausTerm *gt;
784 const HelmholtzCritTerm *ct;
785
786 n = HD->np;
787 pt = &(HD->pt[0]);
788
789 //MSG("tau=%f, del=%f",tau,delta);
790 //if(isinf(tau))abort();
791
792 /* power terms */
793 sum = 0;
794 dell = ipow(delta,pt->l);
795 //ldell = pt->l * dell;
796 unsigned oldl;
797 for(i=0; i<n; ++i){
798 term = pt->a * pow(tau, pt->t) * ipow(delta, pt->d);
799 sum += term;
800 #ifdef RESID_DEBUG
801 MSG("i = %d, a=%e, t=%f, d=%d, l=%u, term = %f, sum = %f",i,pt->a,pt->t,pt->d,pt->l,term,sum);
802 if(pt->l==0){
803 MSG(",row=%e\n",term);
804 }else{
805 MSG(",row=%e\n",term*exp(-dell));
806 }
807 #endif
808 oldl = pt->l;
809 ++pt;
810 if(i+1==n || oldl != pt->l){
811 if(oldl == 0){
812 #ifdef RESID_DEBUG
813 MSG(" linear ");
814 #endif
815 res += sum;
816 }else{
817 #ifdef RESID_DEBUG
818 MSG(" %sEXP dell=%f, exp(-dell)=%f sum=%f: ",(i+1==n?"LAST ":""),dell,exp(-dell),sum);
819 #endif
820 res += sum * exp(-dell);
821 }
822 #ifdef RESID_DEBUG
823 MSG("i = %d, res = %f\n",i,res);
824 #endif
825 sum = 0;
826 if(i+1<n){
827 #ifdef RESID_DEBUG
828 MSG(" next delta = %.12e, l = %u\n",delta, pt->l);
829 #endif
830 dell = (delta==0 ? 0 : ipow(delta,pt->l));
831 //ldell = pt->l*dell;
832 }
833 }
834 }
835 assert(!__isnan(res));
836
837 /* gaussian terms */
838 n = HD->ng;
839 //fprintf(stderr,"THERE ARE %d GAUSSIAN TERMS\n",n);
840 gt = &(HD->gt[0]);
841 for(i=0; i<n; ++i){
842 #ifdef RESID_DEBUG
843 MSG("i = %d, GAUSSIAN, n = %e, t = %f, d = %f, alpha = %f, beta = %f, gamma = %f, epsilon = %f",i+1, gt->n, gt->t, gt->d, gt->alpha, gt->beta, gt->gamma, gt->epsilon);
844 #endif
845 double d1 = delta - gt->epsilon;
846 double t1 = tau - gt->gamma;
847 double e1 = -gt->alpha*SQ(d1) - gt->beta*SQ(t1);
848 sum = gt->n * pow(tau,gt->t) * pow(delta,gt->d) * exp(e1);
849 //fprintf(stderr,"sum = %f\n",sum);
850 res += sum;
851 ++gt;
852 }
853 assert(!__isnan(res));
854
855 /* critical terms */
856 n = HD->nc;
857 ct = &(HD->ct[0]);
858 for(i=0; i<n; ++i){
859 #ifdef RESID_DEBUG
860 MSG("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);
861 #endif
862
863 DEFINE_DELTA;
864 DEFINE_DELB;
865
866 sum = ct->n * DELB * delta * PSI;
867 res += sum;
868 ++ct;
869 }
870 assert(!__isnan(res));
871
872 #ifdef RESID_DEBUG
873 fprintf(stderr,"CALCULATED RESULT FOR phir = %f\n",res);
874 #endif
875 return res;
876 }
877
878 /*=================== FIRST DERIVATIVES =======================*/
879
880 /**
881 Derivative of the helmholtz residual function with respect to
882 delta.
883 */
884 double helm_resid_del(double tau,double delta, const HelmholtzRunData *HD){
885 double sum = 0, res = 0;
886 double dell, ldell;
887 unsigned n, i;
888 const HelmholtzPowTerm *pt;
889 const HelmholtzGausTerm *gt;
890 const HelmholtzCritTerm *ct;
891
892 #ifdef RESID_DEBUG
893 fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
894 #endif
895
896 /* power terms */
897 n = HD->np;
898 pt = &(HD->pt[0]);
899 dell = ipow(delta,pt->l);
900 ldell = pt->l * dell;
901 unsigned oldl;
902 for(i=0; i<n; ++i){
903 sum += pt->a * pow(tau, pt->t) * ipow(delta, pt->d - 1) * (pt->d - ldell);
904 oldl = pt->l;
905 ++pt;
906 if(i+1==n || oldl != pt->l){
907 if(oldl == 0){
908 res += sum;
909 }else{
910 res += sum * exp(-dell);
911 }
912 sum = 0;
913 if(i+1<n){
914 dell = (delta==0 ? 0 : ipow(delta,pt->l));
915 ldell = pt->l*dell;
916 }
917 }
918 }
919
920 /* gaussian terms */
921 n = HD->ng;
922 gt = &(HD->gt[0]);
923 for(i=0; i<n; ++i){
924 #ifdef RESID_DEBUG
925 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);
926 #endif
927 sum = - gt->n * pow(tau,gt->t) * pow(delta, -1. + gt->d)
928 * (2. * gt->alpha * delta * (delta - gt->epsilon) - gt->d)
929 * exp(-(gt->alpha * SQ(delta-gt->epsilon) + gt->beta*SQ(tau-gt->gamma)));
930 res += sum;
931 #ifdef RESID_DEBUG
932 fprintf(stderr,"sum = %f --> res = %f\n",sum,res);
933 #endif
934 ++gt;
935 }
936
937 /* critical terms */
938 n = HD->nc;
939 ct = &(HD->ct[0]);
940 for(i=0; i<n; ++i){
941 #ifdef RESID_DEBUG
942 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);
943 #endif
944 DEFINE_DELTA;
945 DEFINE_DELB;
946 DEFINE_DPSIDDELTA;
947 DEFINE_DDELDDELTA;
948 DEFINE_DDELBDDELTA;
949
950 #if 0
951 if(fabs(dpsiddelta) ==0)fprintf(stderr,"WARNING: dpsiddelta == 0\n");
952 if(fabs(dpsiddelta) ==0)fprintf(stderr,"WARNING: dpsiddelta == 0\n");
953 fprintf(stderr,"psi = %f\n",psi);
954 fprintf(stderr,"DELTA = %f\n",DELTA);
955
956 fprintf(stderr,"dDELddelta = %f\n",dDELddelta);
957 fprintf(stderr,"ct->b - 1. = %f\n",ct->b - 1.);
958 fprintf(stderr,"pow(DELTA,ct->b - 1.) = %f\n",pow(DELTA,ct->b - 1.));
959 assert(!isnan(pow(DELTA,ct->b - 1.)));
960 assert(!isnan(dDELddelta));
961 assert(!isnan(dDELbddelta));
962 //double dDELbddelta = ct->b * pow(DELTA,ct->b - 1.) * dDELddelta
963 fprintf(stderr,"sum = %f\n",sum);
964 if(isnan(sum))fprintf(stderr,"ERROR: sum isnan with i=%d at %d\n",i,__LINE__);
965 #endif
966 sum = ct->n * (DELB * (PSI + delta * dPSIddelta) + dDELbddelta * delta * PSI);
967 res += sum;
968
969 ++ct;
970 }
971
972 return res;
973 }
974
975 /**
976 Derivative of the helmholtz residual function with respect to
977 tau.
978 */
979 double helm_resid_tau(double tau,double delta,const HelmholtzRunData *HD){
980
981 double sum;
982 double res = 0;
983 double delX;
984 unsigned l;
985 unsigned n, i;
986 const HelmholtzPowTerm *pt;
987 const HelmholtzGausTerm *gt;
988 const HelmholtzCritTerm *ct;
989
990 n = HD->np;
991 pt = &(HD->pt[0]);
992
993 delX = 1;
994
995 l = 0;
996 sum = 0;
997 for(i=0; i<n; ++i){
998 if(pt->t){
999 //fprintf(stderr,"i = %d, a = %e, t = %f, d = %d, l = %d\n",i+1, pt->a, pt->t, pt->d, pt->l);
1000 sum += pt->a * pow(tau, pt->t - 1) * ipow(delta, pt->d) * pt->t;
1001 }
1002 ++pt;
1003 //fprintf(stderr,"l = %d\n",l);
1004 if(i+1==n || l != pt->l){
1005 if(l==0){
1006 //fprintf(stderr,"Adding non-exp term\n");
1007 res += sum;
1008 }else{
1009 //fprintf(stderr,"Adding exp term with l = %d, delX = %e\n",l,delX);
1010 res += sum * exp(-delX);
1011 }
1012 /* set l to new value */
1013 if(i+1!=n){
1014 l = pt->l;
1015 //fprintf(stderr,"New l = %d\n",l);
1016 delX = ipow(delta,l);
1017 sum = 0;
1018 }
1019 }
1020 }
1021 assert(!__isnan(res));
1022
1023 //#define RESID_DEBUG
1024 /* gaussian terms */
1025 n = HD->ng;
1026 gt = &(HD->gt[0]);
1027 for(i=0; i<n; ++i){
1028 #ifdef RESID_DEBUG
1029 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);
1030 #endif
1031
1032 double val2;
1033 val2 = -gt->n * pow(tau,gt->t - 1.) * pow(delta, gt->d)
1034 * (2. * gt->beta * tau * (tau - gt->gamma) - gt->t)
1035 * exp(-(gt->alpha * SQ(delta-gt->epsilon) + gt->beta*SQ(tau-gt->gamma)));
1036 res += val2;
1037 #ifdef RESID_DEBUG
1038 fprintf(stderr,"res = %f\n",res);
1039 #endif
1040
1041 ++gt;
1042 }
1043 assert(!__isnan(res));
1044
1045 /* critical terms */
1046 n = HD->nc;
1047 ct = &(HD->ct[0]);
1048 for(i=0; i<n; ++i){
1049 #ifdef RESID_DEBUG
1050 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);
1051 #endif
1052 DEFINE_DELTA;
1053 DEFINE_DELB;
1054 DEFINE_DDELBDTAU;
1055 DEFINE_DPSIDTAU;
1056
1057 sum = ct->n * delta * (dDELbdtau * PSI + DELB * dPSIdtau);
1058 res += sum;
1059 ++ct;
1060 }
1061 assert(!__isnan(res));
1062
1063 return res;
1064 }
1065
1066
1067 /*=================== SECOND DERIVATIVES =======================*/
1068
1069 /**
1070 Mixed derivative of the helmholtz residual function with respect to
1071 delta and tau.
1072 */
1073 double helm_resid_deltau(double tau,double delta,const HelmholtzRunData *HD){
1074 double dell,ldell, sum = 0, res = 0;
1075 unsigned n, i;
1076 const HelmholtzPowTerm *pt;
1077 const HelmholtzGausTerm *gt;
1078 const HelmholtzCritTerm *ct;
1079
1080 /* power terms */
1081 n = HD->np;
1082 pt = &(HD->pt[0]);
1083 dell = ipow(delta,pt->l);
1084 ldell = pt->l * dell;
1085 unsigned oldl;
1086 for(i=0; i<n; ++i){
1087 sum += pt->a * pt->t * pow(tau, pt->t - 1) * ipow(delta, pt->d - 1) * (pt->d - ldell);
1088 oldl = pt->l;
1089 ++pt;
1090 if(i+1==n || oldl != pt->l){
1091 if(oldl == 0){
1092 res += sum;
1093 }else{
1094 res += sum * exp(-dell);
1095 }
1096 sum = 0;
1097 if(i+1<n){
1098 dell = ipow(delta,pt->l);
1099 ldell = pt->l*dell;
1100 }
1101 }
1102 }
1103
1104 #ifdef TEST
1105 assert(!isinf(res));
1106 #endif
1107
1108 /* gaussian terms */
1109 n = HD->ng;
1110 gt = &(HD->gt[0]);
1111 for(i=0; i<n; ++i){
1112 #ifdef RESID_DEBUG
1113 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);
1114 #endif
1115 double d1 = delta - gt->epsilon;
1116 double t1 = tau - gt->gamma;
1117 double e1 = -gt->alpha*SQ(d1) - gt->beta*SQ(t1);
1118
1119 double f1 = gt->t - 2*gt->beta*tau*(tau - gt->gamma);
1120 double g1 = gt->d - 2*gt->alpha*delta*(delta - gt->epsilon);
1121
1122 sum = gt->n * f1 * pow(tau,gt->t-1) * g1 * pow(delta,gt->d-1) * exp(e1);
1123
1124 //fprintf(stderr,"sum = %f\n",sum);
1125 res += sum;
1126 #ifdef TEST
1127 assert(!isinf(res));
1128 #endif
1129 ++gt;
1130 }
1131
1132 /* critical terms */
1133 n = HD->nc;
1134 ct = &(HD->ct[0]);
1135 for(i=0; i<n; ++i){
1136 #ifdef RESID_DEBUG
1137 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);
1138 #endif
1139 DEFINE_DELTA;
1140 DEFINE_DELB;
1141 DEFINE_DPSIDDELTA;
1142 DEFINE_DDELBDTAU;
1143 DEFINE_DDELDDELTA;
1144
1145 double d2DELbddeldtau = -ct->A * ct->b * 2./ct->beta * (DELB/DELTA)*d1*pow(d12,0.5/ct->beta-1) \
1146 - 2. * theta * ct->b * (ct->b - 1) * (DELB/SQ(DELTA)) * dDELddelta;
1147
1148 double d2PSIddeldtau = 4. * ct->C*ct->D*d1*t1*PSI;
1149
1150 DEFINE_DPSIDTAU;
1151
1152 sum = ct->n * (DELB * (dPSIdtau + delta * d2PSIddeldtau) \
1153 + delta *dDELbdtau*dPSIdtau \
1154 + dDELbdtau*(PSI+delta*dPSIddelta) \
1155 + d2DELbddeldtau*delta*PSI
1156 );
1157 res += sum;
1158 ++ct;
1159 }
1160
1161 #ifdef RESID_DEBUG
1162 fprintf(stderr,"phir = %f\n",res);
1163 #endif
1164
1165 #ifdef TEST
1166 assert(!isnan(res));
1167 assert(!isinf(res));
1168 #endif
1169 return res;
1170 }
1171
1172 /**
1173 Second derivative of helmholtz residual function with respect to
1174 delta (twice).
1175 */
1176 double helm_resid_deldel(double tau,double delta,const HelmholtzRunData *HD){
1177 double sum = 0, res = 0;
1178 double dell, ldell;
1179 unsigned n, i;
1180 const HelmholtzPowTerm *pt;
1181 const HelmholtzGausTerm *gt;
1182 const HelmholtzCritTerm *ct;
1183
1184 #ifdef RESID_DEBUG
1185 fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
1186 #endif
1187
1188 /* power terms */
1189 n = HD->np;
1190 pt = &(HD->pt[0]);
1191 dell = ipow(delta,pt->l);
1192 ldell = pt->l * dell;
1193 unsigned oldl;
1194 for(i=0; i<n; ++i){
1195 double lpart = pt->l ? SQ(ldell) + ldell*(1. - 2*pt->d - pt->l) : 0;
1196 sum += pt->a * pow(tau, pt->t) * ipow(delta, pt->d - 2) * (pt->d*(pt->d - 1) + lpart);
1197 oldl = pt->l;
1198 ++pt;
1199 if(i+1==n || oldl != pt->l){
1200 if(oldl == 0){
1201 res += sum;
1202 }else{
1203 res += sum * exp(-dell);
1204 }
1205 sum = 0;
1206 if(i+1<n){
1207 dell = ipow(delta,pt->l);
1208 ldell = pt->l*dell;
1209 }
1210 }
1211 }
1212 if(isnan(res)){
1213 fprintf(stderr,"tau = %.12e, del = %.12e\n",tau,delta);
1214 }
1215 assert(!__isnan(res));
1216
1217 /* gaussian terms */
1218 n = HD->ng;
1219 //fprintf(stderr,"THERE ARE %d GAUSSIAN TERMS\n",n);
1220 gt = &(HD->gt[0]);
1221 for(i=0; i<n; ++i){
1222 double s1 = SQ(delta - gt->epsilon);
1223 double f1 = gt->d*(gt->d - 1)
1224 + 2.*gt->alpha*delta * (
1225 delta * (2. * gt->alpha * s1 - 1)
1226 - 2. * gt->d * (delta - gt->epsilon)
1227 );
1228 res += gt->n * pow(tau,gt->t) * pow(delta, gt->d - 2.)
1229 * f1
1230 * exp(-(gt->alpha * s1 + gt->beta*SQ(tau-gt->gamma)));
1231 ++gt;
1232 }
1233 assert(!__isnan(res));
1234
1235 /* critical terms */
1236 n = HD->nc;
1237 ct = &(HD->ct[0]);
1238 for(i=0; i<n; ++i){
1239 #ifdef RESID_DEBUG
1240 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);
1241 #endif
1242
1243 DEFINE_DELTA;
1244 DEFINE_DELB;
1245 DEFINE_DPSIDDELTA;
1246 DEFINE_DDELDDELTA;
1247 DEFINE_DDELBDDELTA;
1248
1249 DEFINE_D2DELDDELTA2;
1250 DEFINE_D2DELBDDELTA2;
1251
1252 DEFINE_D2PSIDDELTA2;
1253
1254 sum = ct->n * (DELB*(2.*dPSIddelta + delta*d2PSIddelta2) + 2.*dDELbddelta*(PSI+delta*dPSIddelta) + d2DELbddelta2*delta*PSI);
1255
1256 res += sum;
1257 ++ct;
1258 }
1259 assert(!__isnan(res));
1260
1261 return res;
1262 }
1263
1264
1265
1266 /**
1267 Residual part of helmholtz function.
1268 */
1269 double helm_resid_tautau(double tau, double delta, const HelmholtzRunData *HD){
1270 double dell, term, sum, res = 0;
1271 unsigned n, i;
1272 const HelmholtzPowTerm *pt;
1273 const HelmholtzGausTerm *gt;
1274 const HelmholtzCritTerm *ct;
1275
1276 n = HD->np;
1277 pt = &(HD->pt[0]);
1278
1279 #ifdef RESID_DEBUG
1280 fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
1281 #endif
1282
1283 /* power terms */
1284 sum = 0;
1285 dell = ipow(delta,pt->l);
1286 //ldell = pt->l * dell;
1287 unsigned oldl;
1288 for(i=0; i<n; ++i){
1289 term = pt->a * pt->t * (pt->t - 1) * pow(tau, pt->t - 2) * ipow(delta, pt->d);
1290 sum += term;
1291 #ifdef RESID_DEBUG
1292 fprintf(stderr,"i = %d, a=%e, t=%f, d=%d, term = %f, sum = %f",i,pt->a,pt->t,pt->d,term,sum);
1293 if(pt->l==0){
1294 fprintf(stderr,",row=%e\n",term);
1295 }else{
1296 fprintf(stderr,",row=%e\n,",term*exp(-dell));
1297 }
1298 #endif
1299 oldl = pt->l;
1300 ++pt;
1301 if(i+1==n || oldl != pt->l){
1302 if(oldl == 0){
1303 #ifdef RESID_DEBUG
1304 fprintf(stderr,"linear ");
1305 #endif
1306 res += sum;
1307 }else{
1308 #ifdef RESID_DEBUG
1309 fprintf(stderr,"exp dell=%f, exp(-dell)=%f sum=%f: ",dell,exp(-dell),sum);
1310 #endif
1311 res += sum * exp(-dell);
1312 }
1313 #ifdef RESID_DEBUG
1314 fprintf(stderr,"i = %d, res = %f\n",i,res);
1315 #endif
1316 sum = 0;
1317 if(i+1<n){
1318 dell = ipow(delta,pt->l);
1319 //ldell = pt->l*dell;
1320 }
1321 }
1322 }
1323
1324 /* gaussian terms */
1325 n = HD->ng;
1326 //fprintf(stderr,"THERE ARE %d GAUSSIAN TERMS\n",n);
1327 gt = &(HD->gt[0]);
1328 for(i=0; i<n; ++i){
1329 #ifdef RESID_DEBUG
1330 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);
1331 #endif
1332 double d1 = delta - gt->epsilon;
1333 double t1 = tau - gt->gamma;
1334 double f1 = gt->t*(gt->t - 1) + 4. * gt->beta * tau * (tau * (gt->beta*SQ(t1) - 0.5) - t1*gt->t);
1335 double e1 = -gt->alpha*SQ(d1) - gt->beta*SQ(t1);
1336 sum = gt->n * f1 * pow(tau,gt->t - 2) * pow(delta,gt->d) * exp(e1);
1337 //fprintf(stderr,"sum = %f\n",sum);
1338 res += sum;
1339 ++gt;
1340 }
1341
1342 /* critical terms */
1343 n = HD->nc;
1344 ct = &(HD->ct[0]);
1345 for(i=0; i<n; ++i){
1346 #ifdef RESID_DEBUG
1347 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);
1348 #endif
1349 DEFINE_DELTA;
1350 DEFINE_DELB;
1351 DEFINE_DDELBDTAU;
1352 DEFINE_DPSIDTAU;
1353
1354 double d2DELbdtau2 = 2. * ct->b * (DELB/DELTA) + 4. * SQ(theta) * ct->b * (ct->b - 1) * (DELB/SQ(DELTA));
1355
1356 double d2PSIdtau2 = 2. * ct->D * PSI * (2. * ct->D * SQ(t1) -1.);
1357
1358 sum = ct->n * delta * (d2DELbdtau2 * PSI + 2 * dDELbdtau*dPSIdtau + DELB * d2PSIdtau2);
1359 res += sum;
1360 ++ct;
1361 }
1362
1363 #ifdef RESID_DEBUG
1364 fprintf(stderr,"phir_tautau = %f\n",res);
1365 #endif
1366 return res;
1367 }
1368
1369 /* === THIRD DERIVATIVES (this is getting boring now) === */
1370
1371 #ifdef INCLUDE_THIRD_DERIV_CODE
1372 /**
1373 Third derivative of helmholtz residual function, with respect to
1374 delta (thrice).
1375 */
1376 double helm_resid_deldeldel(double tau,double delta,const HelmholtzRunData *HD){
1377 double sum = 0, res = 0;
1378 double D;
1379 unsigned n, i;
1380 const HelmholtzPowTerm *pt;
1381 const HelmholtzGausTerm *gt;
1382 const HelmholtzCritTerm *ct;
1383
1384 #ifdef RESID_DEBUG
1385 fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
1386 #endif
1387
1388 #if 1
1389 /* major shortcut, but not efficient */
1390 double ddel = 0.0000000001;
1391 return (helm_resid_deldel(tau,delta+ddel,HD) - helm_resid_deldel(tau,delta,HD))/ddel;
1392 #endif
1393
1394 /* seem to be errors in the following, still haven't tracked them all down. */
1395
1396 #if 1
1397 /* wxmaxima code:
1398 a*delta^d*%e^(-delta^l)*tau^t
1399 diff(%,delta,3);
1400 */
1401 /* power terms */
1402 n = HD->np;
1403 pt = &(HD->pt[0]);
1404 D = ipow(delta,pt->l);
1405 unsigned oldl;
1406 for(i=0; i<n; ++i){
1407 double d = pt->d;
1408 double l = pt->l;
1409 double lpart = pt->l
1410 ? D*((D-1)*(D-2)-1) * l*SQ(l)
1411 + 3*D*(1-d)*(D-1) * SQ(l)
1412 + D*(3*SQ(d-1)-1) * l
1413 : 0;
1414 sum += pt->a * pow(tau,pt->t) * ipow(delta, d-3) * (d*(d-1)*(d-2) + lpart);
1415 oldl = pt->l;
1416 ++pt;
1417 if(i+1==n || oldl != pt->l){
1418 if(oldl == 0){
1419 res += sum; // note special meaning of l==0 case: no exponential
1420 }else{
1421 res += sum * exp(-D);
1422 }
1423 sum = 0;
1424 D = ipow(delta,pt->l);
1425 }
1426 }
1427
1428 //fprintf(stderr,"DELDELDEL fiff = %f, sum = %f ",fdiff, res);
1429 #endif
1430
1431 #if 1
1432 /* gaussian terms */
1433 n = HD->ng;
1434 //fprintf(stderr,"THERE ARE %d GAUSSIAN TERMS\n",n);
1435 gt = &(HD->gt[0]);
1436 for(i=0; i<n; ++i){
1437 double D = delta - gt->epsilon;
1438 double D2 = SQ(D);
1439 double T2 = SQ(tau - gt->gamma);
1440 double A = gt->alpha * delta;
1441 double A2 = SQ(A);
1442 double d = gt->d;
1443 double d2 = SQ(d);
1444
1445 // this expression calculated from wxMaxima using subsitutions for
1446 // D=delta-epsilon and A=alpha*delta.
1447 double f1 =
1448 - (8*A*A2) * D*D2
1449 + (12*d * A2) * D2
1450 + (12 * delta * A2 + (6*d - 6*d2)*A) * D
1451 - (6 * d * delta * A + d*d2 - 3*d2 + 2*d);
1452
1453 res += gt->n * pow(tau,gt->t) * pow(delta, d - 3.)
1454 * f1
1455 * exp(-(gt->alpha * D2 + gt->beta * T2));
1456 ++gt;
1457 }
1458 #endif
1459
1460 #if 1
1461 /* critical terms */
1462 n = HD->nc;
1463 ct = &(HD->ct[0]);
1464 for(i=0; i<n; ++i){
1465 #ifdef RESID_DEBUG
1466 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);
1467 #endif
1468
1469 DEFINE_DELTA;
1470 DEFINE_DELB;
1471 DEFINE_DPSIDDELTA;
1472 DEFINE_DDELDDELTA;
1473 DEFINE_DDELBDDELTA;
1474
1475 DEFINE_D2PSIDDELTA2;
1476 DEFINE_D2DELDDELTA2;
1477 DEFINE_D2DELBDDELTA2;
1478
1479 DEFINE_D3PSIDDELTA3;
1480 DEFINE_D3DELDDELTA3;
1481 DEFINE_D3DELBDDELTA3;
1482
1483 sum = ct->n * (
1484 delta * (DELB*d3PSIddelta3 + 3 * dDELbddelta * d2PSIddelta2 + 3 * d2DELbddelta2 * dPSIddelta + PSI * d3DELbddelta3)
1485 + 3 * (DELB*d2PSIddelta2 + 2 * dDELbddelta * dPSIddelta + d2DELbddelta2 * PSI)
1486 );
1487
1488 res += sum;
1489 ++ct;
1490 }
1491 #endif
1492
1493 return res;
1494 }
1495
1496 #endif
1497
1498

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