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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1844 by jpye, Sat Aug 30 07:33:11 2008 UTC revision 1845 by jpye, Sat Aug 30 15:02:45 2008 UTC
# Line 152  double helmholtz_s(double T, double rho, Line 152  double helmholtz_s(double T, double rho,
152      assert(!isnan(data->R));      assert(!isnan(data->R));
153  #endif  #endif
154    
155        fprintf(stderr,"helm_ideal_tau = %f\n",helm_ideal_tau(tau,delta,data->ideal));
156        fprintf(stderr,"helm_resid_tau = %f\n",helm_resid_tau(tau,delta,data));
157        fprintf(stderr,"helm_ideal = %f\n",helm_ideal(tau,delta,data->ideal));
158        fprintf(stderr,"helm_resid = %f\n",helm_resid(tau,delta,data));
159      return data->R * (      return data->R * (
160          tau * (helm_ideal_tau(tau,delta,data->ideal) + helm_resid_tau(tau,delta,data))          tau * (helm_ideal_tau(tau,delta,data->ideal) + helm_resid_tau(tau,delta,data))
161          - helm_ideal(tau,delta,data->ideal) - helm_resid(tau,delta,data)          - helm_ideal(tau,delta,data->ideal) - helm_resid(tau,delta,data)
# Line 210  double helm_ideal(double tau, double del Line 214  double helm_ideal(double tau, double del
214          sum += pt->a0 * pow(tau, pt->t0);          sum += pt->a0 * pow(tau, pt->t0);
215      }      }
216    
217  #if 0  #if 1
218      et = &(data->et[0]);      et = &(data->et[0]);
219      for(i=0; i<data->ne; ++i, ++et){      for(i=0; i<data->ne; ++i, ++et){
220          sum += et->b * log( 1 - exp(- et->B * tau));          sum += et->b * log( 1 - exp(- et->B * tau));
# Line 236  double helm_ideal_tau(double tau, double Line 240  double helm_ideal_tau(double tau, double
240          sum += pt->a0 * pt->t0 * pow(tau, pt->t0 - 1);          sum += pt->a0 * pt->t0 * pow(tau, pt->t0 - 1);
241      }      }
242    
243  #if 0  #if 1
244        /* FIXME we're missing the '2.5' in the alpha0 expression for Nitrogen... */
245      et = &(data->et[0]);      et = &(data->et[0]);
246      for(i=0; i<data->ne; ++i, ++et){      for(i=0; i<data->ne; ++i, ++et){
247          sum += et->b * log( 1 - exp(- et->B * tau));          sum += et->b * log( 1 - exp(- et->B * tau));
# Line 246  double helm_ideal_tau(double tau, double Line 251  double helm_ideal_tau(double tau, double
251  }    }  
252    
253  /**  /**
254      Residual part of helmholtz function. Note: we have NOT prematurely      Residual part of helmholtz function.
     optimised here ;-)  
255  */  */
256  double helm_resid(double tau, double delta, const HelmholtzData *data){  double helm_resid(double tau, double delta, const HelmholtzData *data){
257    #if 0
258        double sum, res = 0;
259        unsigned n, i;
260        const HelmholtzPowTerm *pt;
261        const HelmholtzExpTerm *et;
262    
263        n = data->np;
264        pt = &(data->pt[0]);
265    
266        /* power terms */
267        sum = 0;
268        unsigned oldl;
269        for(i=0; i<n; ++i){
270            sum += pt->a * pow(tau, pt->t) * ipow(delta, pt->d);
271            oldl = pt->l;
272            ++pt;
273            if(i+1==n || oldl != pt->l){
274                if(oldl == 0){
275                    res += sum;
276                }else{
277                    res += sum * exp(-ipow(delta,pt->l));
278                }
279                sum = 0;
280            }
281        }
282    
283        /* now the exponential terms */
284        n = data->ne;
285        et = &(data->et[0]);
286        for(i=0; i< n; ++i){
287            fprintf(stderr,"i = %d, a = %e, t = %f, d = %d, phi = %d, beta = %d, gamma = %f\n",i+1, et->a, et->t, et->d, et->phi, et->beta, et->gamma);
288            
289            double e1 = -et->phi * delta*delta
290                         + 2 * et->phi * delta
291                         - et->beta * tau * tau
292                         + 2 * et->beta * et->gamma * tau
293                         - et->phi
294                         - et->beta * et->gamma * et->gamma;
295            sum = et->a * pow(tau,et->t) * ipow(delta,et->d) * exp(e1);
296            //fprintf(stderr,"sum = %f\n",sum);
297            res += sum;
298            ++et;
299        }
300    
301        return res;
302    }
303    
304    #else
305      double sum;      double sum;
306      double res = 0;      double res = 0;
307      double delX;      double delX;
308      unsigned l;      unsigned l;
309      unsigned np, i;      unsigned n, i;
310      const HelmholtzPowTerm *pt;      const HelmholtzPowTerm *pt;
311        const HelmholtzExpTerm *et;
312    
313      np = data->np;      n = data->np;
314      pt = &(data->pt[0]);      pt = &(data->pt[0]);
315    
316      delX = 1;      delX = 1;
317    
318      l = 0;      l = 0;
319      sum = 0;      sum = 0;
320      for(i=0; i<np; ++i){      for(i=0; i<n; ++i){
321          //fprintf(stderr,"i = %d, a = %e, t = %f, d = %d, l = %d\n",i+1, pt->a, pt->t, pt->d, pt->l);          //fprintf(stderr,"i = %d, a = %e, t = %f, d = %d, l = %d\n",i+1, pt->a, pt->t, pt->d, pt->l);
322          sum += pt->a * pow(tau, pt->t) * ipow(delta, pt->d);          sum += pt->a * pow(tau, pt->t) * ipow(delta, pt->d);
323          ++pt;          ++pt;
324          //fprintf(stderr,"l = %d\n",l);          //fprintf(stderr,"l = %d\n",l);
325          if(i+1==np || l != pt->l){          if(i+1==n || l != pt->l){
326              if(l==0){              if(l==0){
327                  //fprintf(stderr,"Adding non-exp term\n");                  //fprintf(stderr,"Adding non-exp term\n");
328                  res += sum;                  res += sum;
# Line 278  double helm_resid(double tau, double del Line 331  double helm_resid(double tau, double del
331                  res += sum * exp(-delX);                  res += sum * exp(-delX);
332              }              }
333              /* set l to new value */              /* set l to new value */
334              if(i+1!=np){              if(i+1!=n){
335                  l = pt->l;                  l = pt->l;
336                  //fprintf(stderr,"New l = %d\n",l);                  //fprintf(stderr,"New l = %d\n",l);
337                  delX = ipow(delta,l);                  delX = ipow(delta,l);
# Line 287  double helm_resid(double tau, double del Line 340  double helm_resid(double tau, double del
340          }          }
341      }      }
342    
343        /* now the exponential terms */
344        n = data->ne;
345        et = &(data->et[0]);
346        for(i=0; i< n; ++i){
347            fprintf(stderr,"i = %d, a = %e, t = %f, d = %d, phi = %d, beta = %d, gamma = %f\n",i+1, et->a, et->t, et->d, et->phi, et->beta, et->gamma);
348            
349            double e1 = -et->phi * delta*delta
350                         + 2 * et->phi * delta
351                         - et->beta * tau * tau
352                         + 2 * et->beta * et->gamma * tau
353                         - et->phi
354                         - et->beta * et->gamma * et->gamma;
355            sum = et->a * pow(tau,et->t) * ipow(delta,et->d) * exp(e1);
356            //fprintf(stderr,"sum = %f\n",sum);
357            res += sum;
358            ++et;
359        }
360    
361      return res;      return res;
362  }  }
363    #endif
364    
365  /**  /**
366      Derivative of the helmholtz residual function with respect to      Derivative of the helmholtz residual function with respect to
# Line 335  double helm_resid_del(double tau,double Line 407  double helm_resid_del(double tau,double
407          double del2 = delta*delta;          double del2 = delta*delta;
408          double tau2 = tau*tau;          double tau2 = tau*tau;
409          double gam2 = et->gamma * et->gamma;          double gam2 = et->gamma * et->gamma;
410          sum = -et->a * pow(tau,et->t) * ipow(delta,et->d-1)          double e1 = -et->phi * del2
             * (2 * et->phi * del2 - 2 * et->phi * delta - et->d)  
             * exp(-et->phi * del2  
411                       + 2 * et->phi * delta                       + 2 * et->phi * delta
412                       - et->beta * tau2                       - et->beta * tau2
413                       + 2 * et->beta * et->gamma * tau                       + 2 * et->beta * et->gamma * tau
414                       - et->phi                       - et->phi
415                       - et->beta * gam2                       - et->beta * gam2;
416                 );          sum = -et->a * pow(tau,et->t) * ipow(delta,et->d-1)
417                * (2 * et->phi * del2 - 2 * et->phi * delta - et->d)
418                * exp(e1);
419          //fprintf(stderr,"sum = %f\n",sum);          //fprintf(stderr,"sum = %f\n",sum);
420          res += sum;          res += sum;
421          ++et;          ++et;
# Line 363  double helm_resid_tau(double tau,double Line 435  double helm_resid_tau(double tau,double
435      double res = 0;      double res = 0;
436      double delX;      double delX;
437      unsigned l;      unsigned l;
438      unsigned np, i;      unsigned n, i;
439      const HelmholtzPowTerm *pt;      const HelmholtzPowTerm *pt;
440        const HelmholtzExpTerm *et;
441    
442      np = data->np;      n = data->np;
443      pt = &(data->pt[0]);      pt = &(data->pt[0]);
444    
445      delX = 1;      delX = 1;
446    
447      l = 0;      l = 0;
448      sum = 0;      sum = 0;
449      for(i=0; i<np; ++i){      for(i=0; i<n; ++i){
450          if(pt->t){          if(pt->t){
451              //fprintf(stderr,"i = %d, a = %e, t = %f, d = %d, l = %d\n",i+1, pt->a, pt->t, pt->d, pt->l);              //fprintf(stderr,"i = %d, a = %e, t = %f, d = %d, l = %d\n",i+1, pt->a, pt->t, pt->d, pt->l);
452              sum += pt->a * pow(tau, pt->t - 1) * ipow(delta, pt->d) * pt->t;              sum += pt->a * pow(tau, pt->t - 1) * ipow(delta, pt->d) * pt->t;
453          }          }
454          ++pt;          ++pt;
455          //fprintf(stderr,"l = %d\n",l);          //fprintf(stderr,"l = %d\n",l);
456          if(i+1==np || l != pt->l){          if(i+1==n || l != pt->l){
457              if(l==0){              if(l==0){
458                  //fprintf(stderr,"Adding non-exp term\n");                  //fprintf(stderr,"Adding non-exp term\n");
459                  res += sum;                  res += sum;
# Line 389  double helm_resid_tau(double tau,double Line 462  double helm_resid_tau(double tau,double
462                  res += sum * exp(-delX);                  res += sum * exp(-delX);
463              }              }
464              /* set l to new value */              /* set l to new value */
465              if(i+1!=np){              if(i+1!=n){
466                  l = pt->l;                  l = pt->l;
467                  //fprintf(stderr,"New l = %d\n",l);                  //fprintf(stderr,"New l = %d\n",l);
468                  delX = ipow(delta,l);                  delX = ipow(delta,l);
# Line 398  double helm_resid_tau(double tau,double Line 471  double helm_resid_tau(double tau,double
471          }          }
472      }      }
473    
474    #if 1
475        /* now the exponential terms */
476        n = data->ne;
477        et = &(data->et[0]);
478        for(i=0; i< n; ++i){
479            //fprintf(stderr,"i = %d, a = %e, t = %f, d = %d, phi = %d, beta = %d, gamma = %f\n",i+1, et->a, et->t, et->d, et->phi, et->beta, et->gamma);
480            
481            double tau2 = tau*tau;
482            double del2 = delta*delta;
483            double gam2 = et->gamma * et->gamma;
484            double e1 = -et->phi * del2
485                         + 2 * et->phi * delta
486                         - et->beta * tau2
487                         + 2 * et->beta * et->gamma * tau
488                         - et->phi
489                         - et->beta * gam2;
490            sum = -et->a * pow(tau,et->t - 1) * ipow(delta,et->d)
491                * (2 * et->beta * tau2 - 2 * et->beta * et->gamma * tau - et->t)
492                * exp(e1);
493            //fprintf(stderr,"sum = %f\n",sum);
494            res += sum;
495            ++et;
496        }
497    #endif
498    
499        return res;
500    
501      return res;      return res;
502  }    }  
503    
Legend:
Removed from v.1844  
changed lines
  Added in v.1845
john.pye@anu.edu.au
 ViewVC Help
Powered by ViewVC 1.1.22