/[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 3060 - (show annotations) (download) (as text)
Wed Aug 12 16:12:37 2015 UTC (3 years, 9 months ago) by sid
File MIME type: text/x-csrc
File size: 56984 byte(s)
update

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

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