/[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 1864 by jpye, Mon Sep 15 03:36:16 2008 UTC revision 1865 by jpye, Mon Sep 15 08:40:14 2008 UTC
# Line 184  double helmholtz_a(double T, double rho, Line 184  double helmholtz_a(double T, double rho,
184    
185  #ifdef TEST  #ifdef TEST
186      fprintf(stderr,"helmholtz_a: T = %f, rho = %f\n",T,rho);      fprintf(stderr,"helmholtz_a: T = %f, rho = %f\n",T,rho);
187        fprintf(stderr,"multiplying by RT = %f\n",data->R*T);
188  #endif  #endif
189    
190      return data->R * T * (helm_ideal(tau,delta,data->ideal) - helm_resid(tau,delta,data));      return data->R * T * (helm_ideal(tau,delta,data->ideal) - helm_resid(tau,delta,data));
# Line 232  static double ipow(double x, int n){ Line 233  static double ipow(double x, int n){
233      Residual part of helmholtz function.      Residual part of helmholtz function.
234  */  */
235  double helm_resid(double tau, double delta, const HelmholtzData *data){  double helm_resid(double tau, double delta, const HelmholtzData *data){
236  #if 1      double dell,ldell, sum, res = 0;
     double sum, res = 0;  
237      unsigned n, i;      unsigned n, i;
238      const HelmholtzPowTerm *pt;      const HelmholtzPowTerm *pt;
239      const HelmholtzExpTerm *et;      const HelmholtzExpTerm *et;
# Line 243  double helm_resid(double tau, double del Line 243  double helm_resid(double tau, double del
243    
244      /* power terms */      /* power terms */
245      sum = 0;      sum = 0;
246        dell = ipow(delta,pt->l);
247        ldell = pt->l * dell;
248      unsigned oldl;      unsigned oldl;
249      for(i=0; i<n; ++i){      for(i=0; i<n; ++i){
250          sum += pt->a * pow(tau, pt->t) * ipow(delta, pt->d);          sum += pt->a * pow(tau, pt->t) * ipow(delta, pt->d);
251            fprintf(stderr,"i = %d,               sum = %f\n",i,sum);
252          oldl = pt->l;          oldl = pt->l;
253          ++pt;          ++pt;
254          if(i+1==n || oldl != pt->l){          if(i+1==n || oldl != pt->l){
255              if(oldl == 0){              if(oldl == 0){
256                    fprintf(stderr,"linear ");
257                  res += sum;                  res += sum;
258              }else{              }else{
259                  res += sum * exp(-ipow(delta,pt->l));                  fprintf(stderr,"exp dell=%f, exp(-dell)=%f sum=%f: ",dell,exp(-dell),sum);
260                    res += sum * exp(-dell);
261              }              }
262                fprintf(stderr,"i = %d, res = %f\n",i,res);
263              sum = 0;              sum = 0;
264                dell = ipow(delta,pt->l);
265                ldell = pt->l*dell;
266          }          }
267      }      }
268    
269    #if 0
270      /* now the exponential terms */      /* now the exponential terms */
271      n = data->ne;      n = data->ne;
272      et = &(data->et[0]);      et = &(data->et[0]);
# Line 275  double helm_resid(double tau, double del Line 284  double helm_resid(double tau, double del
284          res += sum;          res += sum;
285          ++et;          ++et;
286      }      }
287    #endif
288    
289  #ifdef TEST  #ifdef TEST
290      fprintf(stderr,"phir = %f\n",res);      fprintf(stderr,"phir = %f\n",res);
# Line 282  double helm_resid(double tau, double del Line 292  double helm_resid(double tau, double del
292      return res;      return res;
293  }  }
294    
 #else  
     double sum;  
     double res = 0;  
     double delX;  
     unsigned l;  
     unsigned n, i;  
     const HelmholtzPowTerm *pt;  
     const HelmholtzExpTerm *et;  
   
     n = data->np;  
     pt = &(data->pt[0]);  
   
     delX = 1;  
   
     l = 0;  
     sum = 0;  
     for(i=0; i<n; ++i){  
         //fprintf(stderr,"i = %d, a = %e, t = %f, d = %d, l = %d\n",i+1, pt->a, pt->t, pt->d, pt->l);  
         sum += pt->a * pow(tau, pt->t) * ipow(delta, pt->d);  
         ++pt;  
         //fprintf(stderr,"l = %d\n",l);  
         if(i+1==n || l != pt->l){  
             if(l==0){  
                 //fprintf(stderr,"Adding non-exp term\n");  
                 res += sum;  
             }else{  
                 //fprintf(stderr,"Adding exp term with l = %d, delX = %e\n",l,delX);  
                 res += sum * exp(-delX);  
             }  
             /* set l to new value */  
             if(i+1!=n){  
                 l = pt->l;  
                 //fprintf(stderr,"New l = %d\n",l);  
                 delX = ipow(delta,l);  
                 sum = 0;  
             }  
         }  
     }  
   
     /* now the exponential terms */  
     n = data->ne;  
     et = &(data->et[0]);  
     for(i=0; i< n; ++i){  
 #if 0  
         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);  
 #endif  
           
         double e1 = -et->phi * delta*delta  
                      + 2 * et->phi * delta  
                      - et->beta * tau * tau  
                      + 2 * et->beta * et->gamma * tau  
                      - et->phi  
                      - et->beta * et->gamma * et->gamma;  
         sum = et->a * pow(tau,et->t) * ipow(delta,et->d) * exp(e1);  
         //fprintf(stderr,"sum = %f\n",sum);  
         res += sum;  
         ++et;  
     }  
   
     return res;  
 }  
 #endif  
   
295  /**  /**
296      Derivative of the helmholtz residual function with respect to      Derivative of the helmholtz residual function with respect to
297      delta.      delta.

Legend:
Removed from v.1864  
changed lines
  Added in v.1865

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