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

Contents of /trunk/models/johnpye/fprops/helmholtz.c

Parent Directory Parent Directory | Revision Log Revision Log


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

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