| 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 |
|
Parent Directory
| 
Revision Log