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

Contents of /branches/sid/models/johnpye/fprops/helmholtz.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3024 - (show annotations) (download) (as text)
Thu Jul 23 04:16:35 2015 UTC (4 years, 1 month ago) by sid
File MIME type: text/x-csrc
File size: 56901 byte(s)
first attempt at two phase

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_d2pdT2_rho(double T, double rho, const FluidData *data, FpropsError *err);
54 double helmholtz_dpdrho_T(double T, double rho, const FluidData *data, FpropsError *err);
55 double helmholtz_d2pdrho2_T(double T, double rho, const FluidData *data, FpropsError *err);
56 double helmholtz_d2pdTdrho(double T, double rho, const FluidData *data, FpropsError *err);
57
58 double helmholtz_dhdT_rho(double T, double rho, const FluidData *data, FpropsError *err);
59 double helmholtz_d2hdT2_rho(double T, double rho, const FluidData *data, FpropsError *err);
60 double helmholtz_dhdrho_T(double T, double rho, const FluidData *data, FpropsError *err);
61 double helmholtz_d2hdrho2_T(double T, double rho, const FluidData *data, FpropsError *err);
62 double helmholtz_d2hdTdrho(double T, double rho, const FluidData *data, FpropsError *err);
63
64 double helmholtz_dsdT_rho(double T, double rho, const FluidData *data, FpropsError *err);
65 double helmholtz_d2sdT2_rho(double T, double rho, const FluidData *data, FpropsError *err);
66 double helmholtz_dsdrho_T(double T, double rho, const FluidData *data, FpropsError *err);
67 double helmholtz_d2sdrho2_T(double T, double rho, const FluidData *data, FpropsError *err);
68 double helmholtz_d2sdTdrho(double T, double rho, const FluidData *data, FpropsError *err);
69
70
71 double helmholtz_dudT_rho(double T, double rho, const FluidData *data, FpropsError *err);
72 double helmholtz_d2udT2_rho(double T, double rho, const FluidData *data, FpropsError *err);
73 double helmholtz_dudrho_T(double T, double rho, const FluidData *data, FpropsError *err);
74 double helmholtz_d2udrho2_T(double T, double rho, const FluidData *data, FpropsError *err);
75 double helmholtz_d2udTdrho(double T, double rho, const FluidData *data, FpropsError *err);
76
77
78 double helmholtz_dgdT_rho(double T, double rho, const FluidData *data, FpropsError *err);
79 double helmholtz_d2gdT2_rho(double T, double rho, const FluidData *data, FpropsError *err);
80 double helmholtz_dgdrho_T(double T, double rho, const FluidData *data, FpropsError *err);
81 double helmholtz_d2gdrho2_T(double T, double rho, const FluidData *data, FpropsError *err);
82 double helmholtz_d2gdTdrho(double T, double rho, const FluidData *data, FpropsError *err);
83
84
85
86 //#define HELM_DEBUG
87 #define HELM_ERRORS
88 //#define SAT_DEBUG
89
90 #ifdef HELM_DEBUG
91 # include "color.h"
92 # define MSG FPROPS_MSG
93 #else
94 # define MSG(ARGS...) ((void)0)
95 #endif
96
97 /* TODO centralise declaration of our error-reporting function somehow...? */
98 #ifdef HELM_ERRORS
99 # include "color.h"
100 # define ERRMSG FPROPS_ERRMSG
101 #else
102 # define ERRMSG(ARGS...) ((void)0)
103 #endif
104
105
106 //#define ASSERT_DEBUG
107 #ifdef ASSERT_DEBUG
108 # include <assert.h>
109 #else
110 # define assert(ARGS...)
111 #endif
112
113 //#define RESID_DEBUG
114
115 #define INCLUDE_THIRD_DERIV_CODE
116
117 /* macros and forward decls */
118
119 #define SQ(X) ((X)*(X))
120 #define CU(X) ((X)*(X)*(X))
121
122 #include "helmholtz_impl.h"
123
124 /* shortcut take us straight into the correct data structure for Helmholtz correlation calculations */
125 #define HD data->corr.helm
126 #define HD_R data->R
127 #define HD_CP0 data->cp0
128
129 /* calculate tau and delta using a macro -- is used in most functions */
130 #define DEFINE_TD \
131 double tau = data->corr.helm->T_star / T; \
132 double delta = rho / data->corr.helm->rho_star
133
134
135
136
137 PureFluid *helmholtz_prepare(const EosData *E, const ReferenceState *ref){
138 PureFluid *P = FPROPS_NEW(PureFluid);
139
140 if(E->type != FPROPS_HELMHOLTZ){
141 ERRMSG("invalid EOS data, wrong type");
142 return NULL;
143 }
144
145 MSG("Fluid '%s' with T_t = %f", E->name, E->data.helm->T_t);
146
147 P->data = FPROPS_NEW(FluidData);
148 P->data->corr.helm = FPROPS_NEW(HelmholtzRunData);
149
150 /* metadata */
151 /* FIXME strings should be copied, not just referenced */
152 P->name = E->name;
153 P->source = E->source;
154 P->type = E->type;
155 MSG("name = %s",P->name);
156
157 /* common data across all correlation types */
158 #define I E->data.helm
159 P->data->M = I->M;
160 if(I->R == 0){
161 P->data->R = R_UNIVERSAL / I->M;
162 }else{
163 P->data->R = I->R;
164 }
165 P->data->T_t = I->T_t;
166 P->data->T_c = I->T_c;
167 P->data->p_c = 0; // we calculate this later...
168 P->data->rho_c = I->rho_c;
169 P->data->omega = I->omega;
170 P->data->Tstar = I->T_c;
171 P->data->rhostar = I->rho_c;
172 P->data->cp0 = cp0_prepare(E->data.helm->ideal, P->data->R, P->data->T_c);
173
174 /* data specific to helmholtz correlations */
175 #define H P->data->corr.helm
176 H->rho_star = I->rho_star;
177 H->T_star = I->T_star;
178 H->np = I->np;
179 // FIXME copy et, ct, pt to runtime struct?, FIXME see helmholtz_destroy below.
180 H->pt = I->pt;
181 H->ng = I->ng;
182 H->gt = I->gt;
183 H->nc = I->nc;
184 H->ct = I->ct;
185 MSG("np = %d, ng = %d, nc = %d, T_t = %f",H->np,H->ng,H->nc,I->T_t);
186 #undef H
187
188 /* function pointers... more to come still? */
189 #define FN(VAR) P->VAR##_fn = &helmholtz_##VAR
190 FN(p); FN(u); FN(h); FN(s); FN(a); FN(g); FN(cp); FN(cv); FN(w);
191 FN(alphap); FN(betap);
192 FN(dpdrho_T);FN(d2pdrho2_T);FN(dpdT_rho);FN(d2pdT2_rho);FN(d2pdTdrho);
193 FN(dhdrho_T);FN(d2hdrho2_T);FN(dhdT_rho);FN(d2hdT2_rho);FN(d2hdTdrho);
194 FN(dsdrho_T);FN(d2sdrho2_T);FN(dsdT_rho);FN(d2sdT2_rho);FN(d2sdTdrho);
195 FN(dudrho_T);FN(d2udrho2_T);FN(dudT_rho);FN(d2udT2_rho);FN(d2udTdrho);
196 FN(dgdrho_T);FN(d2gdrho2_T);FN(dgdT_rho);FN(d2gdT2_rho);FN(d2gdTdrho);
197 FN(sat);
198 #undef FN
199
200 FpropsError err = 0;
201
202 /* calculate critical pressur
203 P->data->cp0 = cp0_prepare(E->data.helm->ideal, P->dae (doesn't require h0, s0) */
204 MSG("Calculating critical pressure at T_c = %f K, rho_c = %f kg/m3",P->data->T_c, P->data->rho_c);
205 P->data->p_c = helmholtz_p(P->data->T_c, P->data->rho_c, P->data, &err);
206 if(err){
207 ERRMSG("Failed to calculate critical pressure.");
208 FPROPS_FREE(P->data);
209 FPROPS_FREE(P->data->corr.helm);
210 return NULL;
211 }
212 if(P->data->p_c <= 0){
213 ERRMSG("Calculated a critical pressure <= 0! (value = %f)",P->data->p_c);
214 //return NULL;
215 }
216
217 // ref0 is not yet supported for this fluid type:
218 P->data->ref0 = (ReferenceState){FPROPS_REF_TPHG,{.tphg={298.15,0,NAN,NAN}}};
219
220 // fix up the reference point now...
221 if(ref == NULL){
222 // use the provided ReferenceState, or the default one otherwise.
223 ref = &(I->ref);
224 }
225 int res = fprops_set_reference_state(P,ref);
226 if(res){
227 ERRMSG("Unable to apply reference state (type %d, err %d)",ref->type,res);
228 return NULL;
229 }
230
231 P->table = FPROPS_NEW(Ttse); //sid
232
233 #undef I
234 return P;
235 }
236
237 void helmholtz_destroy(PureFluid *P){
238 assert(FPROPS_HELMHOLTZ == P->data);
239 cp0_destroy(P->data->cp0);
240 FPROPS_FREE(P->data->corr.helm);
241 FPROPS_FREE(P->data);
242 FPROPS_FREE(P);
243 }
244
245 /**
246 Function to calculate pressure from Helmholtz free energy EOS, given temperature
247 and mass density.
248
249 @param T temperature in K
250 @param rho mass density in kg/m続
251 @return pressure in Pa
252 */
253 double helmholtz_p(double T, double rho, const FluidData *data, FpropsError *err){
254 DEFINE_TD;
255
256 assert(HD->rho_star!=0);
257 assert(T!=0);
258 assert(!isnan(T));
259 assert(!isnan(rho));
260 assert(!isnan(HD_R));
261
262 //fprintf(stderr,"p calc: T = %f\n",T);
263 //fprintf(stderr,"p calc: tau = %f\n",tau);
264 //fprintf(stderr,"p calc: rho = %f\n",rho);
265 //fprintf(stderr,"p calc: delta = %f\n",delta);
266 //fprintf(stderr,"p calc: R*T*rho = %f\n",HD_R * T * rho);
267
268 //fprintf(stderr,"T = %f\n", T);
269 //fprintf(stderr,"rhob = %f, rhob* = %f, delta = %f\n", rho/HD->M, HD->rho_star/HD->M, delta);
270
271 double p = HD_R * T * rho * (1 + delta * helm_resid_del(tau,delta,HD));
272
273 //printf("%e\n",p);
274 #if 0
275 if(isnan(p)){
276 fprintf(stderr,"T = %.12e, rho = %.12e\n",T,rho);
277 }
278 #endif
279 //abort();
280 if(isnan(p))*err = FPROPS_NUMERIC_ERROR;
281 return p;
282 }
283
284 /**
285 Function to calculate internal energy from Helmholtz free energy EOS, given
286 temperature and mass density.
287
288 @param T temperature in K
289 @param rho mass density in kg/m続
290 @return internal energy in ???
291 */
292 double helmholtz_u(double T, double rho, const FluidData *data, FpropsError *err){
293 DEFINE_TD;
294
295 #ifdef TEST
296 assert(HD->rho_star!=0);
297 assert(T!=0);
298 assert(!isnan(tau));
299 assert(!isnan(delta));
300 assert(!isnan(HD_R));
301 #endif
302
303 #if 0
304 fprintf(stderr,"ideal_tau = %f\n",ideal_phi_tau(tau,delta,HD_CP0));
305 fprintf(stderr,"resid_tau = %f\n",helm_resid_tau(tau,delta,HD));
306 fprintf(stderr,"R T = %f\n",HD_R * HD->T_star);
307 #endif
308
309 return HD_R * HD->T_star * (ideal_phi_tau(tau,delta,HD_CP0) + helm_resid_tau(tau,delta,HD));
310 }
311
312 /**
313 Function to calculate enthalpy from Helmholtz free energy EOS, given
314 temperature and mass density.
315
316 @param T temperature in K
317 @param rho mass density in kg/m続
318 @return enthalpy in J/kg
319 */
320 double helmholtz_h(double T, double rho, const FluidData *data, FpropsError *err){
321 DEFINE_TD;
322
323 //#ifdef TEST
324 assert(HD->rho_star!=0);
325 assert(T!=0);
326 assert(!isnan(tau));
327 assert(!isnan(delta));
328 assert(!isnan(HD_R));
329 //#endif
330 double h = HD_R * T * (1 + tau * (ideal_phi_tau(tau,delta,HD_CP0) + helm_resid_tau(tau,delta,HD)) \
331 + delta*helm_resid_del(tau,delta,HD));
332 assert(!isnan(h));
333 return h;
334 }
335
336 /**
337 Function to calculate entropy from Helmholtz free energy EOS, given
338 temperature and mass density.
339
340 @param T temperature in K
341 @param rho mass density in kg/m続
342 @return entropy in J/kgK
343 */
344 double helmholtz_s(double T, double rho, const FluidData *data, FpropsError *err){
345 DEFINE_TD;
346
347 #ifdef ENTROPY_DEBUG
348 assert(HD->rho_star!=0);
349 assert(T!=0);
350 assert(!isnan(tau));
351 assert(!isnan(delta));
352 assert(!isnan(HD_R));
353
354 fprintf(stderr,"ideal_phi_tau = %f\n",ideal_phi_tau(tau,delta,HD_CP0));
355 fprintf(stderr,"helm_resid_tau = %f\n",helm_resid_tau(tau,delta,HD));
356 fprintf(stderr,"ideal_phi = %f\n",ideal_phi(tau,delta,HD_CP0));
357 fprintf(stderr,"helm_resid = %f\n",helm_resid(tau,delta,HD));
358 #endif
359 return HD_R * (
360 tau * (ideal_phi_tau(tau,delta,HD_CP0) + helm_resid_tau(tau,delta,HD))
361 - (ideal_phi(tau,delta,HD_CP0) + helm_resid(tau,delta,HD))
362 );
363 }
364
365 /**
366 Function to calculate Helmholtz energy from the Helmholtz free energy EOS,
367 given temperature and mass density.
368
369 @param T temperature in K
370 @param rho mass density in kg/m続
371 @return Helmholtz energy 'a', in J/kg
372 */
373 double helmholtz_a(double T, double rho, const FluidData *data, FpropsError *err){
374 DEFINE_TD;
375
376 #ifdef TEST
377 assert(HD->rho_star!=0);
378 assert(T!=0);
379 assert(!isnan(tau));
380 assert(!isnan(delta));
381 assert(!isnan(HD_R));
382 #endif
383
384 #ifdef HELMHOLTZ_DEBUG
385 fprintf(stderr,"helmholtz_a: T = %f, rho = %f\n",T,rho);
386 fprintf(stderr,"multiplying by RT = %f\n",HD_R*T);
387 #endif
388
389 return HD_R * T * (ideal_phi(tau,delta,HD_CP0) + helm_resid(tau,delta,HD));
390 }
391
392 /**
393 Function to calculate isochoric heat capacity from the Helmholtz free energy
394 EOS given temperature and mass density.
395
396 @param T temperature in K
397 @param rho mass density in kg/m続
398 @return Isochoric specific heat capacity in J/kg/K.
399 */
400 double helmholtz_cv(double T, double rho, const FluidData *data, FpropsError *err){
401 DEFINE_TD;
402
403 return - HD_R * SQ(tau) * (ideal_phi_tautau(tau,HD_CP0) + helm_resid_tautau(tau,delta,HD));
404 }
405
406 /**
407 Function to calculate isobaric heat capacity from the Helmholtz free energy
408 EOS given temperature and mass density.
409
410 @param T temperature in K
411 @param rho mass density in kg/m続
412 @return Isobaric specific heat capacity in J/kg/K.
413 */
414 double helmholtz_cp(double T, double rho, const FluidData *data, FpropsError *err){
415 DEFINE_TD;
416
417 double phir_d = helm_resid_del(tau,delta,HD);
418 double phir_dd = helm_resid_deldel(tau,delta,HD);
419 double phir_dt = helm_resid_deltau(tau,delta,HD);
420
421 /* note similarities with helmholtz_w */
422 double temp1 = 1 + 2*delta*phir_d + SQ(delta)*phir_dd;
423 double temp2 = 1 + delta*phir_d - delta*tau*phir_dt;
424 double temp3 = -SQ(tau)*(ideal_phi_tautau(tau,HD_CP0) + helm_resid_tautau(tau,delta,HD));
425
426 return HD_R * (temp3 + SQ(temp2)/temp1);
427 }
428
429
430 /**
431 Function to calculate the speed of sound in a fluid from the Helmholtz free
432 energy EOS, given temperature and mass density.
433
434 @param T temperature in K
435 @param rho mass density in kg/m続
436 @return Speed of sound in m/s.
437 */
438 double helmholtz_w(double T, double rho, const FluidData *data, FpropsError *err){
439 DEFINE_TD;
440
441 double phir_d = helm_resid_del(tau,delta,HD);
442 double phir_dd = helm_resid_deldel(tau,delta,HD);
443 double phir_dt = helm_resid_deltau(tau,delta,HD);
444
445 /* note similarities with helmholtz_cp */
446 double temp1 = 1. + 2.*delta*phir_d + SQ(delta)*phir_dd;
447 double temp2 = 1. + delta*phir_d - delta*tau*phir_dt;
448 double temp3 = -SQ(tau)*(ideal_phi_tautau(tau,HD_CP0) + helm_resid_tautau(tau,delta,HD));
449
450 return sqrt(HD_R * T * (temp1 + SQ(temp2)/temp3));
451
452 }
453
454 /**
455 Function to calculate the Gibbs energy fluid from the Helmholtz free
456 energy EOS, given temperature and mass density.
457
458 @param T temperature in K
459 @param rho mass density in kg/m続
460 @return Gibbs energy, in J/kg.
461 */
462 double helmholtz_g(double T, double rho, const FluidData *data, FpropsError *err){
463 DEFINE_TD;
464
465 double phir_d = helm_resid_del(tau,delta,HD);
466 double phir = helm_resid(tau,delta,HD);
467 double phi0 = ideal_phi(tau,delta,HD_CP0);
468
469 return HD_R * T * (phi0 + phir + 1. + delta * phir_d);
470 }
471
472 /**
473 alpha_p function from IAPWS Advisory Note 3, used in calculation of
474 partial property derivatives.
475 */
476 double helmholtz_alphap(double T, double rho, const FluidData *data, FpropsError *err){
477 DEFINE_TD;
478 double phir_d = helm_resid_del(tau,delta,HD);
479 double phir_dt = helm_resid_deltau(tau,delta,HD);
480 return 1./T * (1. - delta*tau*phir_dt/(1 + delta*phir_d));
481 }
482
483 /**
484 beta_p function from IAPWS Advisory Note 3 , used in calculation of partial
485 property derivatives.
486 */
487 double helmholtz_betap(double T, double rho, const FluidData *data, FpropsError *err){
488 DEFINE_TD;
489 double phir_d = helm_resid_del(tau,delta,HD);
490 double phir_dd = helm_resid_deldel(tau,delta,HD);
491 return rho*(1. + (delta*phir_d + SQ(delta)*phir_dd)/(1+delta*phir_d));
492 }
493
494 /*----------------------------------------------------------------------------
495 PARTIAL DERIVATIVES
496 */
497
498 /**
499 Calculate partial derivative of p with respect to T, with rho constant
500 */
501 double helmholtz_dpdT_rho(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_deltau = helm_resid_deltau(tau,delta,HD);
506 #ifdef TEST
507 assert(!isinf(phir_del));
508 assert(!isinf(phir_deltau));
509 assert(!isnan(phir_del));
510 assert(!isnan(phir_deltau));
511 assert(!isnan(HD_R));
512 assert(!isnan(rho));
513 assert(!isnan(tau));
514 #endif
515
516 double res = HD_R * rho * (1 + delta*phir_del - delta*tau*phir_deltau);
517
518 #ifdef TEST
519 assert(!isnan(res));
520 assert(!isinf(res));
521 #endif
522 return res;
523 }
524
525 /**
526 Calculate partial derivative of p with respect to rho, with T constant
527 */
528 double helmholtz_dpdrho_T(double T, double rho, const FluidData *data, FpropsError *err){
529 DEFINE_TD;
530 //MSG("...");
531 double phir_del = helm_resid_del(tau,delta,HD);
532 double phir_deldel = helm_resid_deldel(tau,delta,HD);
533 #ifdef TEST
534 assert(!isinf(phir_del));
535 assert(!isinf(phir_deldel));
536 #endif
537 double ret = HD_R * T * (1 + 2*delta*phir_del + SQ(delta)*phir_deldel);
538
539
540 //printf ("%e %e %e %e\n", phir_del, phir_deldel, HD_R, ret);
541 return ret;
542 }
543
544 /**
545 Calculate partial second order derivative of p with respect to rho, with T constant
546 */
547
548 double helmholtz_d2pdrho2_T(double T, double rho, const FluidData *data, FpropsError *err){
549 DEFINE_TD;
550
551 double phir_del = helm_resid_del(tau,delta,HD);
552 double phir_deldel = helm_resid_deldel(tau,delta,HD);
553 double phir_deldeldel = helm_resid_deldeldel(tau,delta,HD);
554 #ifdef TEST
555 assert(!isinf(phir_del));
556 assert(!isinf(phir_deldel));
557 assert(!isinf(phir_deldeldel));
558 #endif
559
560 return HD_R * T / rho * delta * (2*phir_del + delta*(4*phir_deldel + delta*phir_deldeldel));
561 }
562
563
564 /**
565 Calculate partial second order derivative of p with respect to T, with rho constant
566 */
567
568 double helmholtz_d2pdT2_rho(double T, double rho, const FluidData *data, FpropsError *err){
569 DEFINE_TD;
570
571 double phir_tautaudel = helm_resid_deltautau(tau,delta,HD);
572 #ifdef TEST
573 assert(!isnan(phir_tautaudel));
574 assert(!isinf(phir_tautaudel));
575 #endif
576 return HD_R * rho / T * SQ(tau) * delta * phir_tautaudel;
577
578 }
579
580 /**
581 Calculate partial second order mixed derivative of p with respect to rho, T
582 */
583
584 double helmholtz_d2pdTdrho(double T, double rho, const FluidData *data, FpropsError *err){
585 DEFINE_TD;
586 double phir_del = helm_resid_del(tau,delta,HD);
587 double phir_deldel = helm_resid_deldel(tau,delta,HD);
588 double phir_taudel = helm_resid_deltau(tau,delta,HD);
589 double phir_taudeldel = helm_resid_deldeltau(tau,delta,HD);
590 #ifdef TEST
591 assert(!isinf(phir_del));
592 assert(!isinf(phir_deldel));
593 assert(!isinf(phir_taudel));
594 assert(!isinf(phir_taudeldel));
595 assert(!isnan(phir_del));
596 assert(!isnan(phir_deldel));
597 assert(!isnan(phir_taudel));
598 assert(!isnan(phir_taudeldel));
599 #endif
600 double res = HD_R * ( 1 + 2*delta* phir_del + SQ(delta) * phir_deldel - 2*tau*delta*phir_taudel - tau*SQ(delta)*phir_taudeldel);
601
602 #ifdef TEST
603 assert(!isinf(res));
604 assert(!isnan(res));
605 #endif
606 return res;
607 }
608
609
610
611 /**
612 Calculate partial derivative of h with respect to T, with rho constant
613 */
614 double helmholtz_dhdT_rho(double T, double rho, const FluidData *data, FpropsError *err){
615 DEFINE_TD;
616
617 double phir_del = helm_resid_del(tau,delta,HD);
618 double phir_deltau = helm_resid_deltau(tau,delta,HD);
619 double phir_tautau = helm_resid_tautau(tau,delta,HD);
620 double phi0_tautau = ideal_phi_tautau(tau,HD_CP0);
621
622 //fprintf(stderr,"phir_del = %f, phir_deltau = %f, phir_tautau = %f, phi0_tautau = %f\n",phir_del,phir_deltau,phir_tautau,phi0_tautau);
623
624 //return (helmholtz_h(T+0.01,rho,data) - helmholtz_h(T,rho,data)) / 0.01;
625 return HD_R * (1. + delta*phir_del - SQ(tau)*(phi0_tautau + phir_tautau) - delta*tau*phir_deltau);
626 }
627
628 /**
629 Calculate partial derivative of h with respect to rho, with T constant
630 */
631 double helmholtz_dhdrho_T(double T, double rho, const FluidData *data, FpropsError *err){
632 DEFINE_TD;
633
634 double phir_del = helm_resid_del(tau,delta,HD);
635 double phir_deltau = helm_resid_deltau(tau,delta,HD);
636 double phir_deldel = helm_resid_deldel(tau,delta,HD);
637
638 return HD_R * T / rho * (tau*delta*(phir_deltau) + delta * phir_del + SQ(delta)*phir_deldel);
639 }
640
641
642
643
644 /**
645 Calculate partial second order derivative of h with respect to rho, with T constant
646 */
647
648 double helmholtz_d2hdrho2_T(double T, double rho, const FluidData *data, FpropsError *err){
649 DEFINE_TD;
650
651 double phir_deldel = helm_resid_deldel(tau,delta,HD);
652 double phir_deldeldel = helm_resid_deldeldel(tau,delta,HD);
653 double phir_taudeldel = helm_resid_deldeltau(tau,delta,HD);
654
655 #ifdef TEST
656 assert(!isinf(phir_deldel));
657 assert(!isinf(phir_deldeldel));
658 assert(!isinf(phir_taudeldel));
659 assert(!isnan(phir_deldel));
660 assert(!isnan(phir_deldeldel));
661 assert(!isnan(phir_taudeldel));
662 #endif
663
664 double res = HD_R * T / SQ(rho) * ( tau*SQ(delta)*phir_taudeldel + 2*SQ(delta)*phir_deldel + CU(delta)*phir_deldeldel );
665
666 #ifdef TEST
667 assert(!isinf(res));
668 assert(!isnan(res));
669 #endif
670 return res;
671 }
672
673 /**
674 Calculate partial second order derivative of h with respect to T, with rho constant
675 */
676
677 double helmholtz_d2hdT2_rho(double T, double rho, const FluidData *data, FpropsError *err){
678 DEFINE_TD;
679
680
681 double phi0_tautau = ideal_phi_tautau(tau,HD_CP0);
682 double phi0_tautautau = ideal_phi_tautautau(tau,HD_CP0);
683
684 double phir_tautau = helm_resid_tautau(tau,delta,HD);
685 double phir_tautaudel = helm_resid_deltautau(tau,delta,HD);
686 double phir_tautautau = helm_resid_tautautau(tau,delta,HD);
687
688 #ifdef TEST
689 assert(!isinf(phi0_tautau));
690 assert(!isinf(phi0_tautautau));
691 assert(!isinf(phir_tautau));
692 assert(!isinf(phir_tautaudel));
693 assert(!isinf(phir_tautautau));
694 assert(!isnan(phi0_tautau));
695 assert(!isnan(phi0_tautautau));
696 assert(!isnan(phir_tautau));
697 assert(!isnan(phir_tautaudel));
698 assert(!isnan(phir_tautautau));
699 #endif
700
701 double res = HD_R / T * ( CU(tau)*(phi0_tautautau+phir_tautautau) + 2*SQ(tau)*(phi0_tautau+phir_tautau) + SQ(tau)*delta*phir_tautaudel );
702
703 #ifdef TEST
704 assert(!isinf(res));
705 assert(!isnan(res));
706 #endif
707 return res;
708 }
709
710 /**
711 Calculate partial second order mixed derivative of h with respect to rho, T
712 */
713
714 double helmholtz_d2hdTdrho(double T, double rho, const FluidData *data, FpropsError *err){
715 DEFINE_TD;
716
717 double phir_deldel = helm_resid_deldel(tau,delta,HD);
718 double phir_tautaudel = helm_resid_deltautau(tau,delta,HD);
719 double phir_del = helm_resid_del(tau,delta,HD);
720 double phir_taudeldel = helm_resid_deldeltau(tau,delta,HD);
721 double phir_taudel = helm_resid_deltau(tau,delta,HD);
722
723 #ifdef TEST
724 assert(!isinf(phir_deldel));
725 assert(!isinf(phir_tautaudel));
726 assert(!isinf(phir_del));
727 assert(!isinf(phir_taudeldel));
728 assert(!isinf(phir_taudel));
729 assert(!isnan(phir_deldel));
730 assert(!isnan(phir_tautaudel));
731 assert(!isnan(phir_del));
732 assert(!isnan(phir_taudeldel));
733 assert(!isnan(phir_taudel));
734 #endif
735
736 double res = HD_R / rho * ( SQ(delta)*(phir_deldel) - SQ(tau)*delta*(phir_tautaudel) + delta*phir_del - SQ(delta)*tau*(phir_taudeldel) - tau*delta*(phir_taudel) );
737
738 #ifdef TEST
739 assert(!isinf(res));
740 assert(!isnan(res));
741 #endif
742 return res;
743 }
744
745
746
747
748
749
750
751
752
753
754 /**
755 Calculate partial derivative of u with respect to T, with rho constant
756 */
757 double helmholtz_dudT_rho(double T, double rho, const FluidData *data, FpropsError *err){
758 DEFINE_TD;
759
760 double phir_tautau = helm_resid_tautau(tau,delta,HD);
761 double phi0_tautau = ideal_phi_tautau(tau,HD_CP0);
762
763 return -HD_R * SQ(tau) * (phi0_tautau + phir_tautau);
764 }
765
766
767 /**
768 Calculate partial derivative of u with respect to rho, with T constant
769 */
770 double helmholtz_dudrho_T(double T, double rho, const FluidData *data, FpropsError *err){
771 DEFINE_TD;
772
773 double phir_deltau = helm_resid_deltau(tau,delta,HD);
774
775 return HD_R * T / rho * (tau * delta * phir_deltau);
776 }
777
778 /**
779 Calculate partial second order derivative of u with respect to rho, with T constant
780 */
781
782 double helmholtz_d2udrho2_T(double T, double rho, const FluidData *data, FpropsError *err){
783 DEFINE_TD;
784
785 double phir_taudeldel = helm_resid_deldeltau(tau,delta,HD);
786
787 #ifdef TEST
788 assert(!isinf(phir_taudeldel));
789 assert(!isnan(phir_taudeldel));
790 #endif
791
792 double res = HD_R * T / SQ(rho) * ( tau * SQ(delta) * phir_taudeldel);
793
794 #ifdef TEST
795 assert(!isinf(res));
796 assert(!isnan(res));
797 #endif
798 return res;
799 }
800
801 /**
802 Calculate partial second order derivative of u with respect to T, with rho constant
803 */
804
805 double helmholtz_d2udT2_rho(double T, double rho, const FluidData *data, FpropsError *err){
806 DEFINE_TD;
807
808
809 double phi0_tautau = ideal_phi_tautau(tau,HD_CP0);
810 double phi0_tautautau = ideal_phi_tautautau(tau,HD_CP0);
811
812 double phir_tautau = helm_resid_tautau(tau,delta,HD);
813 double phir_tautautau = helm_resid_tautautau(tau,delta,HD);
814
815 #ifdef TEST
816 assert(!isinf(phi0_tautau));
817 assert(!isinf(phi0_tautautau));
818 assert(!isinf(phir_tautau));
819 assert(!isinf(phir_tautautau));
820 assert(!isnan(phi0_tautau));
821 assert(!isnan(phi0_tautautau));
822 assert(!isnan(phir_tautau));
823 assert(!isnan(phir_tautautau));
824 #endif
825
826 double res = HD_R / T * ( CU(tau)*(phi0_tautautau+phir_tautautau) + 2*SQ(tau)*(phi0_tautau+phir_tautau) );
827
828 #ifdef TEST
829 assert(!isinf(res));
830 assert(!isnan(res));
831 #endif
832 return res;
833 }
834
835 /**
836 Calculate partial second order mixed derivative of u with respect to rho, T
837 */
838
839 double helmholtz_d2udTdrho(double T, double rho, const FluidData *data, FpropsError *err){
840 DEFINE_TD;
841
842 double phir_tautaudel = helm_resid_deltautau(tau,delta,HD);
843
844 #ifdef TEST
845 assert(!isinf(phir_tautaudel));
846 assert(!isnan(phir_tautaudel));
847 #endif
848
849 double res = HD_R / rho * ( -SQ(tau)*delta*phir_tautaudel);
850
851 #ifdef TEST
852 assert(!isinf(res));
853 assert(!isnan(res));
854 #endif
855 return res;
856 }
857
858
859
860
861
862
863
864
865
866 /**
867 Calculate partial derivative of s with respect to T, with rho constant
868 */
869 double helmholtz_dsdT_rho(double T, double rho, const FluidData *data, FpropsError *err){
870 DEFINE_TD;
871
872 double phi0_tautau = ideal_phi_tautau(tau,HD_CP0);
873 double phir_tautau = helm_resid_tautau(tau,delta,HD);
874
875 #ifdef TEST
876 assert(!isinf(phi0_tautau));
877 assert(!isinf(phir_tautau));
878 assert(!isnan(phi0_tautau));
879 assert(!isnan(phir_tautau));
880 #endif
881
882 double res = HD_R / T * ( -SQ(tau) * (phi0_tautau + phir_tautau) );
883
884 #ifdef TEST
885 assert(!isinf(res));
886 assert(!isnan(res));
887 #endif
888 return res;
889 }
890 /**
891 Calculate partial derivative of s with respect to rho, with T constant
892 */
893 double helmholtz_dsdrho_T(double T, double rho, const FluidData *data, FpropsError *err){
894 DEFINE_TD;
895
896 double phir_del = helm_resid_del(tau,delta,HD);
897 double phir_taudel = helm_resid_deltau(tau,delta,HD);
898
899 #ifdef TEST
900 assert(!isinf(phir_del));
901 assert(!isinf(phir_taudel));
902 assert(!isnan(phir_del));
903 assert(!isnan(phir_taudel));
904 #endif
905
906 double res = HD_R / rho * (-1 - delta*phir_del + tau*delta*phir_taudel);
907
908 #ifdef TEST
909 assert(!isinf(res));
910 assert(!isnan(res));
911 #endif
912 return res;
913 }
914
915 /**
916 Calculate partial second order derivative of s with respect to rho, with T constant
917 */
918 double helmholtz_d2sdrho2_T(double T, double rho, const FluidData *data, FpropsError *err){
919 DEFINE_TD;
920
921 double phir_taudeldel = helm_resid_deldeltau(tau,delta,HD);
922 double phir_deldel = helm_resid_deldel(tau,delta,HD);
923
924 #ifdef TEST
925 assert(!isinf(phir_taudeldel));
926 assert(!isinf(phir_deldel));
927 assert(!isnan(phir_taudeldel));
928 assert(!isnan(phir_deldel));
929 #endif
930
931 double res = (HD_R / SQ(rho)) * (1 - SQ(delta)*phir_deldel + tau*SQ(delta)*phir_taudeldel);
932
933 #ifdef TEST
934 assert(!isinf(res));
935 assert(!isnan(res));
936 #endif
937 return res;
938 }
939
940 /**
941 Calculate partial second order derivative of s with respect to T, with rho constant
942 */
943 double helmholtz_d2sdT2_rho(double T, double rho, const FluidData *data, FpropsError *err){
944 DEFINE_TD;
945
946 double phi0_tautau = ideal_phi_tautau(tau,HD_CP0);
947 double phi0_tautautau = ideal_phi_tautautau(tau,HD_CP0);
948
949 double phir_tautau = helm_resid_tautau(tau,delta,HD);
950 double phir_tautautau = helm_resid_tautautau(tau,delta,HD);
951
952 #ifdef TEST
953 assert(!isinf(phi0_tautau));
954 assert(!isinf(phi0_tautautau));
955 assert(!isinf(phir_tautau));
956 assert(!isinf(phir_tautautau));
957 assert(!isnan(phi0_tautau));
958 assert(!isnan(phi0_tautautau));
959 assert(!isnan(phir_tautau));
960 assert(!isnan(phir_tautautau));
961 #endif
962
963 double res = HD_R / SQ(T) * ( CU(tau)*(phi0_tautautau+phir_tautautau) + 3*SQ(tau)*(phi0_tautau+phir_tautau) );
964
965 #ifdef TEST
966 assert(!isinf(res));
967 assert(!isnan(res));
968 #endif
969 return res;
970 }
971
972
973 /**
974 Calculate partial second order mixed derivative of s with respect to rho, T
975 */
976 double helmholtz_d2sdTdrho(double T, double rho, const FluidData *data, FpropsError *err){
977 DEFINE_TD;
978
979 double phir_tautaudel = helm_resid_deltautau(tau,delta,HD);
980
981 #ifdef TEST
982 assert(!isinf(phir_tautaudel));
983 assert(!isnan(phir_tautaudel));
984 #endif
985
986 double res = HD_R / (T*rho) * ( -SQ(tau)*delta*phir_tautaudel);
987
988 #ifdef TEST
989 assert(!isinf(res));
990 assert(!isnan(res));
991 #endif
992 return res;
993 }
994
995
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
1007 /**
1008 Calculate partial derivative of g with respect to T, with rho constant
1009 */
1010 double helmholtz_dgdT_rho(double T, double rho, const FluidData *data, FpropsError *err){
1011 DEFINE_TD;
1012
1013 double phi0 = ideal_phi(tau,delta,HD_CP0);
1014 double phi0_tau = ideal_phi_tau(tau,delta,HD_CP0);
1015 double phir = helm_resid(tau,delta,HD);
1016 double phir_tau = helm_resid_tau(tau,delta,HD);
1017 double phir_del = helm_resid_del(tau,delta,HD);
1018 double phir_taudel = helm_resid_deltau(tau,delta,HD);
1019
1020 #ifdef TEST
1021 assert(!isinf(phi0));
1022 assert(!isinf(phi0_tau));
1023 assert(!isinf(phir));
1024 assert(!isinf(phir_tau));
1025 assert(!isinf(phir_del));
1026 assert(!isinf(phir_taudel));
1027 assert(!isnan(phi0));
1028 assert(!isnan(phi0_tau));
1029 assert(!isnan(phir));
1030 assert(!isnan(phir_tau));
1031 assert(!isnan(phir_del));
1032 assert(!isnan(phir_taudel));
1033 #endif
1034
1035 double res = HD_R * ( - tau*(phi0_tau + phir_tau) + (phi0 + phir) + (1 + delta*phir_del - tau*delta*phir_taudel) );
1036
1037 #ifdef TEST
1038 assert(!isinf(res));
1039 assert(!isnan(res));
1040 #endif
1041
1042 return res;
1043 }
1044
1045 /**
1046 Calculate partial derivative of g with respect to rho, with T constant
1047 */
1048 double helmholtz_dgdrho_T(double T, double rho, const FluidData *data, FpropsError *err){
1049 DEFINE_TD;
1050
1051 double phir_del = helm_resid_del(tau,delta,HD);
1052 double phir_deldel = helm_resid_deldel(tau,delta,HD);
1053
1054 #ifdef TEST
1055 assert(!isinf(phir_del));
1056 assert(!isinf(phir_deldel));
1057 assert(!isnan(phir_del));
1058 assert(!isnan(phir_deldel));
1059 #endif
1060
1061 double res = HD_R * T / rho * (1 + 2*delta*phir_del - SQ(delta)*phir_deldel) ;
1062
1063 #ifdef TEST
1064 assert(!isinf(res));
1065 assert(!isnan(res));
1066 #endif
1067
1068 return res;
1069 }
1070
1071 /**
1072 Calculate partial second order derivative of g with respect to rho, with T constant
1073 */
1074 double helmholtz_d2gdrho2_T(double T, double rho, const FluidData *data, FpropsError *err){
1075 DEFINE_TD;
1076
1077 double phir_deldel = helm_resid_deldel(tau,delta,HD);
1078 double phir_deldeldel = helm_resid_deldeldel(tau,delta,HD);
1079
1080 #ifdef TEST
1081 assert(!isinf(phir_deldel));
1082 assert(!isinf(phir_deldeldel));
1083 assert(!isnan(phir_deldel));
1084 assert(!isnan(phir_deldeldel));
1085 #endif
1086
1087 double res = HD_R * T / SQ(rho) * ( -1 + 3 * SQ(delta) * phir_deldel + CU(delta) * phir_deldeldel );
1088
1089 #ifdef TEST
1090 assert(!isinf(res));
1091 assert(!isnan(res));
1092 #endif
1093 return res;
1094 }
1095
1096 /**
1097 Calculate partial second order derivative of g with respect to T, with rho constant
1098 */
1099 double helmholtz_d2gdT2_rho(double T, double rho, const FluidData *data, FpropsError *err){
1100 DEFINE_TD;
1101
1102 double phi0_tautau = ideal_phi_tautau(tau,HD_CP0);
1103 double phir_tautau = helm_resid_tautau(tau,delta,HD);
1104 double phir_tautaudel = helm_resid_deltautau(tau,delta,HD);
1105
1106 #ifdef TEST
1107 assert(!isinf(phi0_tautau));
1108 assert(!isinf(phir_tautau));
1109 assert(!isinf(phir_tautaudel));
1110 assert(!isnan(phi0_tautau));
1111 assert(!isnan(phir_tautau));
1112 assert(!isnan(phir_tautaudel));
1113 #endif
1114
1115 double res = HD_R / T * ( SQ(tau)*(phi0_tautau+phir_tautau) + SQ(tau)*delta*phir_tautaudel );
1116
1117 #ifdef TEST
1118 assert(!isinf(res));
1119 assert(!isnan(res));
1120 #endif
1121 return res;
1122 }
1123
1124 /**
1125 Calculate partial second order mixed derivative of g with respect to rho, T
1126 */
1127 double helmholtz_d2gdTdrho(double T, double rho, const FluidData *data, FpropsError *err){
1128 DEFINE_TD;
1129
1130 double phir_deldel = helm_resid_deldel(tau,delta,HD);
1131 double phir_del = helm_resid_del(tau,delta,HD);
1132 double phir_taudeldel = helm_resid_deldeltau(tau,delta,HD);
1133 double phir_taudel = helm_resid_deltau(tau,delta,HD);
1134
1135 #ifdef TEST
1136 assert(!isinf(phir_deldel));
1137 assert(!isinf(phir_del));
1138 assert(!isinf(phir_taudeldel));
1139 assert(!isinf(phir_taudel));
1140 assert(!isnan(phir_deldel));
1141 assert(!isnan(phir_del));
1142 assert(!isnan(phir_taudeldel));
1143 assert(!isnan(phir_taudel));
1144 #endif
1145
1146 double res = HD_R / rho * ( 1 + 2*delta*phir_del - 2*tau*delta*phir_taudel + SQ(delta)*phir_deldel - tau*SQ(delta)*phir_taudeldel );
1147
1148 #ifdef TEST
1149 assert(!isinf(res));
1150 assert(!isnan(res));
1151 #endif
1152 return res;
1153 }
1154
1155
1156
1157
1158 /**
1159 Solve saturation condition for a specified temperature using approach of
1160 Akasaka, but adapted for general use to non-helmholtz property correlations.
1161 @param T temperature [K]
1162 @param psat_out output, saturation pressure [Pa]
1163 @param rhof_out output, saturated liquid density [kg/m^3]
1164 @param rhog_out output, saturated vapour density [kg/m^3]
1165 @param d helmholtz data object for the fluid in question.
1166 @return 0 on success, non-zero on error (eg algorithm failed to converge, T out of range, etc.)
1167 */
1168 double helmholtz_sat(double T, double *rhof_out, double * rhog_out, const FluidData *data, FpropsError *err){
1169 if(T < data->T_t - 1e-8){
1170 ERRMSG("Input temperature %f K is below triple-point temperature %f K",T,data->T_t);
1171 return FPROPS_RANGE_ERROR;
1172 }
1173 if(T > data->T_c + 1e-8){
1174 ERRMSG("Input temperature is above critical point temperature");
1175 *err = FPROPS_RANGE_ERROR;
1176 }
1177
1178 // we're at the critical point
1179 if(fabs(T - data->T_c) < 1e-9){
1180 *rhof_out = data->rho_c;
1181 *rhog_out = data->rho_c;
1182 return data->p_c;
1183 }
1184
1185 // FIXME at present step-length multiplier is set to 0.4 just because of
1186 // ONE FLUID, ethanol. Probably our initial guess data isn't good enough,
1187 // or maybe there's a problem with the acentric factor or something like
1188 // that. This factor 0.4 will be slowing down the whole system, so it's not
1189 // good. TODO XXX.
1190
1191 // initial guesses for liquid and vapour density
1192 double rhof = 1.1 * fprops_rhof_T_rackett(T,data);
1193 double rhog= 0.9 * fprops_rhog_T_chouaieb(T,data);
1194 double R = data->R;
1195 double pc = data->p_c;
1196
1197 #ifdef SAT_DEBUG
1198 MSG("initial guess rho_f = %f, rho_g = %f",rhof,rhog);
1199 MSG("calculating at T = %.12e",T);
1200 #endif
1201
1202 int i = 0;
1203 while(i++ < 200){
1204 assert(!isnan(rhog));
1205 assert(!isnan(rhof));
1206 #ifdef SAT_DEBUG
1207 MSG("iter %d: T = %f, rhof = %f, rhog = %f",i,T, rhof, rhog);
1208 #endif
1209
1210 double pf = helmholtz_p(T,rhof,data,err);
1211 double pg = helmholtz_p(T,rhog,data,err);
1212 double gf = helmholtz_a(T,rhof,data,err) + pf/rhof;
1213 double gg = helmholtz_a(T,rhog,data,err) + pg/rhog;
1214 double dpdrf = helmholtz_dpdrho_T(T,rhof,data,err);
1215 double dpdrg = helmholtz_dpdrho_T(T,rhog,data,err);
1216
1217 // jacobian for [F;G](rhof, rhog) --- derivatives wrt rhof and rhog
1218 double F = (pf - pg)/pc;
1219 double G = (gf - gg)/R/T;
1220
1221 if(fabs(F) + fabs(G) < 1e-12){
1222 //fprintf(stderr,"%s: CONVERGED\n",__func__);
1223 *rhof_out = rhof;
1224 *rhog_out = rhog;
1225 return helmholtz_p(T, *rhog_out, data, err);
1226 /* SUCCESS */
1227 }
1228
1229 double Ff = dpdrf/pc;
1230 double Fg = -dpdrg/pc;
1231 //MSG("Ff = %e, Fg = %e",Ff,Fg);
1232
1233 double Gf = dpdrf/rhof/R/T;
1234 double Gg = -dpdrg/rhog/R/T;
1235 //MSG("Gf = %e, Gg = %e",Gf,Gg);
1236
1237 double DET = Ff*Gg - Fg*Gf;
1238 //MSG("DET = %f",DET);
1239
1240 // 'gamma' needs to be increased to 0.5 for water to solve correctly (see 'test/sat.c')
1241 // 'gamma' needs to be not more than 0.4 for ethanol to solve correctly (see 'test/sat.c')
1242 #define gamma 0.40
1243 rhof += gamma/DET * (Fg*G - Gg*F);
1244 rhog += gamma/DET * ( Gf*F - Ff*G);
1245 #undef gamma
1246
1247 assert(!isnan(rhof));
1248 assert(!isnan(rhog));
1249
1250 if(rhog < 0)rhog = -0.5*rhog;
1251 if(rhof < 0)rhof = -0.5*rhof;
1252 }
1253 *rhof_out = rhof;
1254 *rhog_out = rhog;
1255 *err = FPROPS_SAT_CVGC_ERROR;
1256 ERRMSG("Not converged: with T = %e (rhof=%f, rhog=%f).",T,*rhof_out,*rhog_out);
1257 return helmholtz_p(T, rhog, data, err);
1258 }
1259
1260
1261
1262
1263 /*---------------------------------------------
1264 UTILITY FUNCTION(S)
1265 */
1266
1267 /* ipow: public domain by Mark Stephen with suggestions by Keiichi Nakasato */
1268 static double ipow(double x, int n){
1269 double t = 1.0;
1270
1271 if(!n)return 1.0; /* At the top. x^0 = 1 */
1272
1273 if (n < 0){
1274 n = -n;
1275 x = 1.0/x; /* error if x == 0. Good */
1276 } /* ZTC/SC returns inf, which is even better */
1277
1278 if (x == 0.0)return 0.0;
1279
1280 do{
1281 if(n & 1)t *= x;
1282 n /= 2; /* KN prefers if (n/=2) x*=x; This avoids an */
1283 x *= x; /* unnecessary but benign multiplication on */
1284 }while(n); /* the last pass, but the comparison is always
1285 true _except_ on the last pass. */
1286
1287 return t;
1288 }
1289
1290 /* maxima expressions:
1291 Psi(delta) := exp(-C*(delta-1)^2 -D*(tau-1)^2);
1292 theta(delta) := (1-tau) + A*((delta-1)^2)^(1/(2*beta));
1293 Delta(delta):= theta(delta)^2 + B*((delta-1)^2)^a;
1294 n*Delta(delta)^b*delta*Psi(delta);
1295 diff(%,delta,3);
1296 yikes, that's scary! break down into steps.
1297 */
1298
1299 #undef HD
1300
1301 /*
1302 We avoid duplication by using the following #defines for common code in
1303 calculation of critical terms.
1304 */
1305 #define DEFINE_DELTA \
1306 double d1 = delta - 1.; \
1307 double t1 = tau - 1.; \
1308 double d12 = SQ(d1); \
1309 double theta = (1. - tau) + ct->A * pow(d12, 0.5/ct->beta); \
1310 double PSI = exp(-ct->C*d12 - ct->D*SQ(t1)); \
1311 double DELTA = SQ(theta) + ct->B* pow(d12, ct->a)
1312
1313 #define DEFINE_DELB \
1314 double DELB = pow(DELTA,ct->b)
1315
1316 #define DEFINE_DPSIDDELTA \
1317 double dPSIddelta = -2. * ct->C * d1 * PSI
1318
1319 #define DEFINE_DDELDDELTA \
1320 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))
1321
1322 #define DEFINE_DDELBDTAU \
1323 double dDELbdtau = (DELTA == 0) ? 0 : -2. * theta * ct->b * (DELB/DELTA);\
1324 assert(!__isnan(dDELbdtau))
1325
1326 #define DEFINE_DPSIDTAU \
1327 double dPSIdtau = -2. * ct->D * t1 * PSI
1328
1329 #define DEFINE_DDELBDDELTA \
1330 double dDELbddelta = (DELTA==0?0:ct->b * (DELB/DELTA) * dDELddelta)
1331
1332 #define DEFINE_D2DELDDELTA2 \
1333 double powd12bm1 = pow(d12,0.5/ct->beta-1.); \
1334 double d2DELddelta2 = 1./d1*dDELddelta + d12*( \
1335 4.*ct->B*ct->a*(ct->a-1.)*pow(d12,ct->a-2.) \
1336 + 2.*SQ(ct->A)*SQ(1./ct->beta)*SQ(powd12bm1) \
1337 + ct->A*theta*4./ct->beta*(0.5/ct->beta-1.)*powd12bm1/d12 \
1338 )
1339
1340 #define DEFINE_D2DELBDDELTA2 \
1341 double d2DELbddelta2 = ct->b * ( (DELB/DELTA)*d2DELddelta2 + (ct->b-1.)*(DELB/SQ(DELTA)*SQ(dDELddelta)))
1342
1343 #define DEFINE_D2PSIDDELTA2 \
1344 double d2PSIddelta2 = (2.*ct->C*d12 - 1.)*2.*ct->C * PSI
1345
1346 #define DEFINE_D3PSIDDELTA3 \
1347 double d3PSIddelta3 = -4. * d1 * SQ(ct->C) * (2.*d12*ct->C - 3.) * PSI
1348
1349 #define DEFINE_D3DELDDELTA3 \
1350 double d3DELddelta3 = 1./(d1*d12*ct->beta*SQ(ct->beta)) * (\
1351 4*ct->B*ct->a*(1.+ct->a*(2*ct->a-3))*SQ(ct->beta)*pow(d12,ct->a)\
1352 + ct->A * (1.+ct->beta*(2*ct->beta-3))*pow(d12,0.5/ct->beta)\
1353 )
1354
1355 #define DEFINE_D3DELBDDELTA3 \
1356 double d3DELbddelta3 = ct->b / (DELTA*SQ(DELTA)) * ( \
1357 (2+ct->b*(ct->b-3))*dDELddelta*SQ(dDELddelta)*DELB \
1358 + DELB*SQ(DELTA)*d3DELddelta3 \
1359 + 3*(ct->b-1) * DELB * DELTA * dDELddelta * d2DELddelta2 \
1360 )
1361
1362 /**
1363 Residual part of helmholtz function.
1364 */
1365 double helm_resid(double tau, double delta, const HelmholtzRunData *HD){
1366 double dell, term, sum, res = 0;
1367 unsigned n, i;
1368 const HelmholtzPowTerm *pt;
1369 const HelmholtzGausTerm *gt;
1370 const HelmholtzCritTerm *ct;
1371
1372 n = HD->np;
1373 pt = &(HD->pt[0]);
1374
1375 //MSG("tau=%f, del=%f",tau,delta);
1376 //if(isinf(tau))abort();
1377
1378 /* power terms */
1379 sum = 0;
1380 dell = ipow(delta,pt->l);
1381 //ldell = pt->l * dell;
1382 unsigned oldl;
1383 for(i=0; i<n; ++i){
1384 term = pt->a * pow(tau, pt->t) * ipow(delta, pt->d);
1385 sum += term;
1386 #ifdef RESID_DEBUG
1387 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);
1388 if(pt->l==0){
1389 MSG(",row=%e\n",term);
1390 }else{
1391 MSG(",row=%e\n",term*exp(-dell));
1392 }
1393 #endif
1394 oldl = pt->l;
1395 ++pt;
1396 if(i+1==n || oldl != pt->l){
1397 if(oldl == 0){
1398 #ifdef RESID_DEBUG
1399 MSG(" linear ");
1400 #endif
1401 res += sum;
1402 }else{
1403 #ifdef RESID_DEBUG
1404 MSG(" %sEXP dell=%f, exp(-dell)=%f sum=%f: ",(i+1==n?"LAST ":""),dell,exp(-dell),sum);
1405 #endif
1406 res += sum * exp(-dell);
1407 }
1408 #ifdef RESID_DEBUG
1409 MSG("i = %d, res = %f\n",i,res);
1410 #endif
1411 sum = 0;
1412 if(i+1<n){
1413 #ifdef RESID_DEBUG
1414 MSG(" next delta = %.12e, l = %u\n",delta, pt->l);
1415 #endif
1416 dell = (delta==0 ? 0 : ipow(delta,pt->l));
1417 //ldell = pt->l*dell;
1418 }
1419 }
1420 }
1421 assert(!__isnan(res));
1422
1423 /* gaussian terms */
1424 n = HD->ng;
1425 //fprintf(stderr,"THERE ARE %d GAUSSIAN TERMS\n",n);
1426 gt = &(HD->gt[0]);
1427 for(i=0; i<n; ++i){
1428 #ifdef RESID_DEBUG
1429 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);
1430 #endif
1431 double d1 = delta - gt->epsilon;
1432 double t1 = tau - gt->gamma;
1433 double e1 = -gt->alpha*SQ(d1) - gt->beta*SQ(t1);
1434 sum = gt->n * pow(tau,gt->t) * pow(delta,gt->d) * exp(e1);
1435 //fprintf(stderr,"sum = %f\n",sum);
1436 res += sum;
1437 ++gt;
1438 }
1439 assert(!__isnan(res));
1440
1441 /* critical terms */
1442 n = HD->nc;
1443 ct = &(HD->ct[0]);
1444 for(i=0; i<n; ++i){
1445 #ifdef RESID_DEBUG
1446 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);
1447 #endif
1448 DEFINE_DELTA;
1449 DEFINE_DELB;
1450
1451 sum = ct->n * DELB * delta * PSI;
1452 res += sum;
1453 ++ct;
1454 }
1455 assert(!__isnan(res));
1456
1457 #ifdef RESID_DEBUG
1458 fprintf(stderr,"CALCULATED RESULT FOR phir = %f\n",res);
1459 #endif
1460 return res;
1461 }
1462
1463 /*=================== FIRST DERIVATIVES =======================*/
1464
1465 /**
1466 Derivative of the helmholtz residual function with respect to
1467 delta.
1468 */
1469 double helm_resid_del(double tau,double delta, const HelmholtzRunData *HD){
1470 double sum = 0, res = 0;
1471 double dell, ldell;
1472 unsigned n, i;
1473 const HelmholtzPowTerm *pt;
1474 const HelmholtzGausTerm *gt;
1475 const HelmholtzCritTerm *ct;
1476
1477 #ifdef RESID_DEBUG
1478 fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
1479 #endif
1480
1481 /* power terms */
1482 n = HD->np;
1483 pt = &(HD->pt[0]);
1484 dell = ipow(delta,pt->l);
1485 ldell = pt->l * dell;
1486 unsigned oldl;
1487 for(i=0; i<n; ++i){
1488 sum += pt->a * pow(tau, pt->t) * ipow(delta, pt->d - 1) * (pt->d - ldell);
1489 oldl = pt->l;
1490 ++pt;
1491 if(i+1==n || oldl != pt->l){
1492 if(oldl == 0){
1493 res += sum;
1494 }else{
1495 res += sum * exp(-dell);
1496 }
1497 sum = 0;
1498 if(i+1<n){
1499 dell = (delta==0 ? 0 : ipow(delta,pt->l));
1500 ldell = pt->l*dell;
1501 }
1502 }
1503 }
1504
1505 /* gaussian terms */
1506 n = HD->ng;
1507 gt = &(HD->gt[0]);
1508 for(i=0; i<n; ++i){
1509 #ifdef RESID_DEBUG
1510 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);
1511 #endif
1512 sum = - gt->n * pow(tau,gt->t) * pow(delta, -1. + gt->d)
1513 * (2. * gt->alpha * delta * (delta - gt->epsilon) - gt->d)
1514 * exp(-(gt->alpha * SQ(delta-gt->epsilon) + gt->beta*SQ(tau-gt->gamma)));
1515 res += sum;
1516 #ifdef RESID_DEBUG
1517 fprintf(stderr,"sum = %f --> res = %f\n",sum,res);
1518 #endif
1519 ++gt;
1520 }
1521
1522 /* critical terms */
1523 n = HD->nc;
1524 ct = &(HD->ct[0]);
1525 for(i=0; i<n; ++i){
1526 #ifdef RESID_DEBUG
1527 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);
1528 #endif
1529 DEFINE_DELTA;
1530 DEFINE_DELB;
1531 DEFINE_DPSIDDELTA;
1532 DEFINE_DDELDDELTA;
1533 DEFINE_DDELBDDELTA;
1534
1535 #if 0
1536 if(fabs(dpsiddelta) ==0)fprintf(stderr,"WARNING: dpsiddelta == 0\n");
1537 if(fabs(dpsiddelta) ==0)fprintf(stderr,"WARNING: dpsiddelta == 0\n");
1538 fprintf(stderr,"psi = %f\n",psi);
1539 fprintf(stderr,"DELTA = %f\n",DELTA);
1540
1541 fprintf(stderr,"dDELddelta = %f\n",dDELddelta);
1542 fprintf(stderr,"ct->b - 1. = %f\n",ct->b - 1.);
1543 fprintf(stderr,"pow(DELTA,ct->b - 1.) = %f\n",pow(DELTA,ct->b - 1.));
1544 assert(!isnan(pow(DELTA,ct->b - 1.)));
1545 assert(!isnan(dDELddelta));
1546 assert(!isnan(dDELbddelta));
1547 //double dDELbddelta = ct->b * pow(DELTA,ct->b - 1.) * dDELddelta
1548 fprintf(stderr,"sum = %f\n",sum);
1549 if(isnan(sum))fprintf(stderr,"ERROR: sum isnan with i=%d at %d\n",i,__LINE__);
1550 #endif
1551 sum = ct->n * (DELB * (PSI + delta * dPSIddelta) + dDELbddelta * delta * PSI);
1552 res += sum;
1553
1554 ++ct;
1555 }
1556
1557 return res;
1558 }
1559
1560 /**
1561 Derivative of the helmholtz residual function with respect to
1562 tau.
1563 */
1564 double helm_resid_tau(double tau,double delta,const HelmholtzRunData *HD){
1565
1566 double sum;
1567 double res = 0;
1568 double delX;
1569 unsigned l;
1570 unsigned n, i;
1571 const HelmholtzPowTerm *pt;
1572 const HelmholtzGausTerm *gt;
1573 const HelmholtzCritTerm *ct;
1574
1575 n = HD->np;
1576 pt = &(HD->pt[0]);
1577
1578 delX = 1;
1579
1580 l = 0;
1581 sum = 0;
1582 for(i=0; i<n; ++i){
1583 if(pt->t){
1584 //fprintf(stderr,"i = %d, a = %e, t = %f, d = %d, l = %d\n",i+1, pt->a, pt->t, pt->d, pt->l);
1585 sum += pt->a * pow(tau, pt->t - 1) * ipow(delta, pt->d) * pt->t;
1586 }
1587 ++pt;
1588 //fprintf(stderr,"l = %d\n",l);
1589 if(i+1==n || l != pt->l){
1590 if(l==0){
1591 //fprintf(stderr,"Adding non-exp term\n");
1592 res += sum;
1593 }else{
1594 //fprintf(stderr,"Adding exp term with l = %d, delX = %e\n",l,delX);
1595 res += sum * exp(-delX);
1596 }
1597 /* set l to new value */
1598 if(i+1!=n){
1599 l = pt->l;
1600 //fprintf(stderr,"New l = %d\n",l);
1601 delX = ipow(delta,l);
1602 sum = 0;
1603 }
1604 }
1605 }
1606 assert(!__isnan(res));
1607
1608 //#define RESID_DEBUG
1609 /* gaussian terms */
1610 n = HD->ng;
1611 gt = &(HD->gt[0]);
1612 for(i=0; i<n; ++i){
1613 #ifdef RESID_DEBUGrundata
1614 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);
1615 #endif
1616
1617 double val2;
1618 val2 = -gt->n * pow(tau,gt->t - 1.) * pow(delta, gt->d)
1619 * (2. * gt->beta * tau * (tau - gt->gamma) - gt->t)
1620 * exp(-(gt->alpha * SQ(delta-gt->epsilon) + gt->beta*SQ(tau-gt->gamma)));
1621 res += val2;
1622 #ifdef RESID_DEBUG
1623 fprintf(stderr,"res = %f\n",res);
1624 #endif
1625
1626 ++gt;
1627 }
1628 assert(!__isnan(res));
1629
1630 /* critical terms */
1631 n = HD->nc;
1632 ct = &(HD->ct[0]);
1633 for(i=0; i<n; ++i){
1634 #ifdef RESID_DEBUG
1635 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);
1636 #endif
1637 DEFINE_DELTA;
1638 DEFINE_DELB;
1639 DEFINE_DDELBDTAU;
1640 DEFINE_DPSIDTAU;
1641
1642 sum = ct->n * delta * (dDELbdtau * PSI + DELB * dPSIdtau);
1643 res += sum;
1644 ++ct;
1645 }
1646 assert(!__isnan(res));
1647
1648 return res;
1649 }
1650
1651
1652 /*=================== SECOND DERIVATIVES =======================*/
1653
1654 /**
1655 Mixed derivative of the helmholtz residual function with respect to
1656 delta and tau.
1657 */
1658 double helm_resid_deltau(double tau,double delta,const HelmholtzRunData *HD){
1659 double dell,ldell, sum = 0, res = 0;
1660 unsigned n, i;
1661 const HelmholtzPowTerm *pt;
1662 const HelmholtzGausTerm *gt;
1663 const HelmholtzCritTerm *ct;
1664
1665 /* power terms */
1666 n = HD->np;
1667 pt = &(HD->pt[0]);
1668 dell = ipow(delta,pt->l);
1669 ldell = pt->l * dell;
1670 unsigned oldl;
1671 for(i=0; i<n; ++i){
1672 sum += pt->a * pt->t * pow(tau, pt->t - 1) * ipow(delta, pt->d - 1) * (pt->d - ldell);
1673 oldl = pt->l;
1674 ++pt;
1675 if(i+1==n || oldl != pt->l){
1676 if(oldl == 0){
1677 res += sum;
1678 }else{
1679 res += sum * exp(-dell);
1680 }
1681 sum = 0;
1682 if(i+1<n){
1683 dell = ipow(delta,pt->l);
1684 ldell = pt->l*dell;
1685 }
1686 }
1687 }
1688
1689 #ifdef TEST
1690 assert(!isinf(res));
1691 #endif
1692
1693 /* gaussian terms */
1694 n = HD->ng;
1695 gt = &(HD->gt[0]);
1696 for(i=0; i<n; ++i){
1697 #ifdef RESID_DEBUG
1698 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);
1699 #endif
1700 double d1 = delta - gt->epsilon;
1701 double t1 = tau - gt->gamma;
1702 double e1 = -gt->alpha*SQ(d1) - gt->beta*SQ(t1);
1703
1704 double f1 = gt->t - 2*gt->beta*tau*(tau - gt->gamma);
1705 double g1 = gt->d - 2*gt->alpha*delta*(delta - gt->epsilon);
1706
1707 sum = gt->n * f1 * pow(tau,gt->t-1) * g1 * pow(delta,gt->d-1) * exp(e1);
1708
1709 //fprintf(stderr,"sum = %f\n",sum);
1710 res += sum;
1711 #ifdef TEST
1712 assert(!isinf(res));
1713 #endif
1714 ++gt;
1715 }
1716
1717 /* critical terms */
1718 n = HD->nc;
1719 ct = &(HD->ct[0]);
1720 for(i=0; i<n; ++i){
1721 #ifdef RESID_DEBUG
1722 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);
1723 #endif
1724 DEFINE_DELTA;
1725 DEFINE_DELB;
1726 DEFINE_DPSIDDELTA;
1727 DEFINE_DDELBDTAU;
1728 DEFINE_DDELDDELTA;
1729
1730 double d2DELbddeldtau = -ct->A * ct->b * 2./ct->beta * (DELB/DELTA)*d1*pow(d12,0.5/ct->beta-1) \
1731 - 2. * theta * ct->b * (ct->b - 1) * (DELB/SQ(DELTA)) * dDELddelta;
1732
1733 double d2PSIddeldtau = 4. * ct->C*ct->D*d1*t1*PSI;
1734
1735 DEFINE_DPSIDTAU;
1736
1737 sum = ct->n * (DELB * (dPSIdtau + delta * d2PSIddeldtau) \
1738 + delta *dDELbdtau*dPSIdtau \
1739 + dDELbdtau*(PSI+delta*dPSIddelta) \
1740 + d2DELbddeldtau*delta*PSI
1741 );
1742 res += sum;
1743 ++ct;
1744 }
1745
1746 #ifdef RESID_DEBUG
1747 fprintf(stderr,"phir = %f\n",res);
1748 #endif
1749
1750 #ifdef TEST
1751 assert(!isnan(res));
1752 assert(!isinf(res));
1753 #endif
1754 return res;
1755 }
1756
1757 /**
1758 Second derivative of helmholtz residual function with respect to
1759 delta (twice).
1760 */
1761 double helm_resid_deldel(double tau,double delta,const HelmholtzRunData *HD){
1762 double sum = 0, res = 0;
1763 double dell, ldell;
1764 unsigned n, i;
1765 const HelmholtzPowTerm *pt;
1766 const HelmholtzGausTerm *gt;
1767 const HelmholtzCritTerm *ct;
1768
1769 #ifdef RESID_DEBUG
1770 fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
1771 #endif
1772
1773 /* power terms */
1774 n = HD->np;
1775 pt = &(HD->pt[0]);
1776 dell = ipow(delta,pt->l);
1777 ldell = pt->l * dell;
1778 unsigned oldl;
1779 for(i=0; i<n; ++i){
1780 double lpart = pt->l ? SQ(ldell) + ldell*(1. - 2*pt->d - pt->l) : 0;
1781 sum += pt->a * pow(tau, pt->t) * ipow(delta, pt->d - 2) * (pt->d*(pt->d - 1) + lpart);
1782 oldl = pt->l;
1783 ++pt;
1784 if(i+1==n || oldl != pt->l){
1785 if(oldl == 0){
1786 res += sum;
1787 }else{
1788 res += sum * exp(-dell);
1789 }
1790 sum = 0;
1791 if(i+1<n){
1792 dell = ipow(delta,pt->l);
1793 ldell = pt->l*dell;
1794 }
1795 }
1796 }
1797 #if RESID_DEBUG
1798 if(isnan(res)){
1799 fprintf(stderr,"got NAN in %s: tau = %.12e, del = %.12e\n",__func__,tau,delta);
1800 }
1801 assert(!__isnan(res));
1802 #endif
1803
1804 /* gaussian terms */
1805 n = HD->ng;
1806 //fprintf(stderr,"THERE ARE %d GAUSSIAN TERMS\n",n);
1807 gt = &(HD->gt[0]);
1808 for(i=0; i<n; ++i){
1809 double s1 = SQ(delta - gt->epsilon);
1810 double f1 = gt->d*(gt->d - 1)
1811 + 2.*gt->alpha*delta * (
1812 delta * (2. * gt->alpha * s1 - 1)
1813 - 2. * gt->d * (delta - gt->epsilon)
1814 );
1815 res += gt->n * pow(tau,gt->t) * pow(delta, gt->d - 2.)
1816 * f1
1817 * exp(-(gt->alpha * s1 + gt->beta*SQ(tau-gt->gamma)));
1818 ++gt;
1819 }
1820 assert(!__isnan(res));
1821
1822 /* critical terms */
1823 n = HD->nc;
1824 ct = &(HD->ct[0]);
1825 for(i=0; i<n; ++i){
1826 #ifdef RESID_DEBUG
1827 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);
1828 #endif
1829
1830 DEFINE_DELTA;
1831 DEFINE_DELB;
1832 DEFINE_DPSIDDELTA;
1833 DEFINE_DDELDDELTA;
1834 DEFINE_DDELBDDELTA;
1835
1836 DEFINE_D2DELDDELTA2;
1837 DEFINE_D2DELBDDELTA2;
1838
1839 DEFINE_D2PSIDDELTA2;
1840
1841 sum = ct->n * (DELB*(2.*dPSIddelta + delta*d2PSIddelta2) + 2.*dDELbddelta*(PSI+delta*dPSIddelta) + d2DELbddelta2*delta*PSI);
1842
1843 res += sum;
1844 ++ct;
1845 }
1846 assert(!__isnan(res));
1847
1848 return res;
1849 }
1850
1851
1852
1853 /**
1854 Residual part of helmholtz function.
1855 */
1856 double helm_resid_tautau(double tau, double delta, const HelmholtzRunData *HD){
1857 double dell, term, sum, res = 0;
1858 unsigned n, i;
1859 const HelmholtzPowTerm *pt;
1860 const HelmholtzGausTerm *gt;
1861 const HelmholtzCritTerm *ct;
1862
1863 n = HD->np;
1864 pt = &(HD->pt[0]);
1865
1866 #ifdef RESID_DEBUG
1867 fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
1868 #endif
1869
1870 /* power terms */
1871 sum = 0;
1872 dell = ipow(delta,pt->l);
1873 //ldell = pt->l * dell;
1874 unsigned oldl;
1875 for(i=0; i<n; ++i){
1876 term = pt->a * pt->t * (pt->t - 1) * pow(tau, pt->t - 2) * ipow(delta, pt->d);
1877 sum += term;
1878 #ifdef RESID_DEBUG
1879 fprintf(stderr,"i = %d, a=%e, t=%f, d=%d, term = %f, sum = %f",i,pt->a,pt->t,pt->d,term,sum);
1880 if(pt->l==0){
1881 fprintf(stderr,",row=%e\n",term);
1882 }else{
1883 fprintf(stderr,",row=%e\n,",term*exp(-dell));
1884 }
1885 #endif
1886 oldl = pt->l;
1887 ++pt;
1888 if(i+1==n || oldl != pt->l){
1889 if(oldl == 0){
1890 #ifdef RESID_DEBUG
1891 fprintf(stderr,"linear ");
1892 #endif
1893 res += sum;
1894 }else{
1895 #ifdef RESID_DEBUG
1896 fprintf(stderr,"exp dell=%f, exp(-dell)=%f sum=%f: ",dell,exp(-dell),sum);
1897 #endif
1898 res += sum * exp(-dell);
1899 }
1900 #ifdef RESID_DEBUG
1901 fprintf(stderr,"i = %d, res = %f\n",i,res);
1902 #endif
1903 sum = 0;
1904 if(i+1<n){
1905 dell = ipow(delta,pt->l);
1906 //ldell = pt->l*dell;
1907 }
1908 }
1909 }
1910
1911 /* gaussian terms */
1912 n = HD->ng;
1913 //fprintf(stderr,"THERE ARE %d GAUSSIAN TERMS\n",n);
1914 gt = &(HD->gt[0]);
1915 for(i=0; i<n; ++i){
1916 #ifdef RESID_DEBUG
1917 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);
1918 #endif
1919 double d1 = delta - gt->epsilon;
1920 double t1 = tau - gt->gamma;
1921 double f1 = gt->t*(gt->t - 1) + 4. * gt->beta * tau * (tau * (gt->beta*SQ(t1) - 0.5) - t1*gt->t);
1922 double e1 = -gt->alpha*SQ(d1) - gt->beta*SQ(t1);
1923 sum = gt->n * f1 * pow(tau,gt->t - 2) * pow(delta,gt->d) * exp(e1);
1924 //fprintf(stderr,"sum = %f\n",sum);
1925 res += sum;
1926 ++gt;
1927 }
1928
1929 /* critical terms */
1930 n = HD->nc;
1931 ct = &(HD->ct[0]);
1932 for(i=0; i<n; ++i){
1933 #ifdef RESID_DEBUG
1934 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);
1935 #endif
1936 DEFINE_DELTA;
1937 DEFINE_DELB;
1938 DEFINE_DDELBDTAU;
1939 DEFINE_DPSIDTAU;
1940
1941 double d2DELbdtau2 = 2. * ct->b * (DELB/DELTA) + 4. * SQ(theta) * ct->b * (ct->b - 1) * (DELB/SQ(DELTA));
1942
1943 double d2PSIdtau2 = 2. * ct->D * PSI * (2. * ct->D * SQ(t1) -1.);
1944
1945 sum = ct->n * delta * (d2DELbdtau2 * PSI + 2 * dDELbdtau*dPSIdtau + DELB * d2PSIdtau2);
1946 res += sum;
1947 ++ct;
1948 }
1949
1950 #ifdef RESID_DEBUG
1951 fprintf(stderr,"phir_tautau = %f\n",res);
1952 #endif
1953 return res;
1954 }
1955
1956 /* === THIRD DERIVATIVES (this is getting boring now) === */
1957
1958 #ifdef INCLUDE_THIRD_DERIV_CODE
1959 /**
1960 Third derivative of helmholtz residual function, with respect to
1961 delta (thrice).
1962 */
1963 double helm_resid_deldeldel(double tau,double delta,const HelmholtzRunData *HD){
1964 double sum = 0, res = 0;
1965 double D;
1966 unsigned n, i;
1967 const HelmholtzPowTerm *pt;
1968 const HelmholtzGausTerm *gt;
1969 const HelmholtzCritTerm *ct;
1970
1971 #ifdef RESID_DEBUG
1972 fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
1973 #endif
1974
1975 #if 1
1976 /* major shortcut, but not efficient */
1977 double ddel = 0.0000000001;
1978 return (helm_resid_deldel(tau,delta+ddel,HD) - helm_resid_deldel(tau,delta,HD))/ddel;
1979 #endif
1980
1981 /* seem to be errors in the following, still haven't tracked them all down. */
1982
1983 #if 1
1984 /* wxmaxima code:
1985 a*delta^d*%e^(-delta^l)*tau^t
1986 diff(%,delta,3);
1987 */
1988 /* power terms */
1989 n = HD->np;
1990 pt = &(HD->pt[0]);
1991 D = ipow(delta,pt->l);
1992 unsigned oldl;
1993 for(i=0; i<n; ++i){
1994 double d = pt->d;
1995 double l = pt->l;
1996 double lpart = pt->l
1997 ? D*((D-1)*(D-2)-1) * l*SQ(l)
1998 + 3*D*(1-d)*(D-1) * SQ(l)
1999 + D*(3*SQ(d-1)-1) * l
2000 : 0;
2001 sum += pt->a * pow(tau,pt->t) * ipow(delta, d-3) * (d*(d-1)*(d-2) + lpart);
2002 oldl = pt->l;
2003 ++pt;
2004 if(i+1==n || oldl != pt->l){
2005 if(oldl == 0){
2006 res += sum; // note special meaning of l==0 case: no exponential
2007 }else{
2008 res += sum * exp(-D);
2009 }
2010 sum = 0;
2011 D = ipow(delta,pt->l);
2012 }
2013 }
2014
2015 //fprintf(stderr,"DELDELDEL fiff = %f, sum = %f ",fdiff, res);
2016 #endif
2017
2018 #if 1
2019 /* gaussian terms */
2020 n = HD->ng;
2021 //fprintf(stderr,"THERE ARE %d GAUSSIAN TERMS\n",n);
2022 gt = &(HD->gt[0]);
2023 for(i=0; i<n; ++i){
2024 double D = delta - gt->epsilon;
2025 double D2 = SQ(D);
2026 double T2 = SQ(tau - gt->gamma);
2027 double A = gt->alpha * delta;
2028 double A2 = SQ(A);
2029 double d = gt->d;
2030 double d2 = SQ(d);
2031
2032 // this expression calculated from wxMaxima using subsitutions for
2033 // D=delta-epsilon and A=alpha*delta.
2034 double f1 =
2035 - (8*A*A2) * D*D2
2036 + (12*d * A2) * D2
2037 + (12 * delta * A2 + (6*d - 6*d2)*A) * D
2038 - (6 * d * delta * A + d*d2 - 3*d2 + 2*d);
2039
2040 res += gt->n * pow(tau,gt->t) * pow(delta, d - 3.)
2041 * f1
2042 * exp(-(gt->alpha * D2 + gt->beta * T2));
2043 ++gt;
2044 }
2045 #endif
2046
2047 #if 1
2048 /* critical terms */
2049 n = HD->nc;
2050 ct = &(HD->ct[0]);
2051 for(i=0; i<n; ++i){
2052 #ifdef RESID_DEBUG
2053 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);
2054 #endif
2055
2056 DEFINE_DELTA;
2057 DEFINE_DELB;
2058 DEFINE_DPSIDDELTA;
2059 DEFINE_DDELDDELTA;
2060 DEFINE_DDELBDDELTA;
2061
2062 DEFINE_D2PSIDDELTA2;
2063 DEFINE_D2DELDDELTA2;
2064 DEFINE_D2DELBDDELTA2;
2065
2066 DEFINE_D3PSIDDELTA3;
2067 DEFINE_D3DELDDELTA3;
2068 DEFINE_D3DELBDDELTA3;
2069
2070 sum = ct->n * (
2071 delta * (DELB*d3PSIddelta3 + 3 * dDELbddelta * d2PSIddelta2 + 3 * d2DELbddelta2 * dPSIddelta + PSI * d3DELbddelta3)
2072 + 3 * (DELB*d2PSIddelta2 + 2 * dDELbddelta * dPSIddelta + d2DELbddelta2 * PSI)
2073 );
2074
2075 res += sum;
2076 ++ct;
2077 }
2078 #endif
2079
2080 return res;
2081 }
2082
2083 /**
2084 Third derivative of helmholtz residual function, with respect to
2085 delta once and tau (twice).
2086 */
2087 double helm_resid_deltautau(double tau,double delta,const HelmholtzRunData *HD){
2088
2089 double ddel = 0.0000000001;
2090 return (helm_resid_tautau(tau,delta+ddel,HD) - helm_resid_tautau(tau,delta,HD))/ddel;
2091
2092 }
2093
2094
2095 /**
2096 Third derivative of helmholtz residual function, with respect to
2097 delta (thrice).
2098 */
2099 double helm_resid_deldeltau(double tau,double delta,const HelmholtzRunData *HD){
2100
2101 double ddel = 0.0000000001;
2102 return (helm_resid_deltau(tau,delta+ddel,HD) - helm_resid_deltau(tau,delta,HD))/ddel;
2103
2104 }
2105
2106
2107
2108 /**
2109 Third derivative of helmholtz residual function, with respect to
2110 tau (thrice).
2111 */
2112 double helm_resid_tautautau(double tau,double delta,const HelmholtzRunData *HD){
2113
2114 double dtau = 0.0000000001;
2115 return (helm_resid_tautau(tau+dtau,delta,HD) - helm_resid_tautau(tau,delta,HD))/dtau;
2116 }
2117
2118 #endif
2119

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