/[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 1834 by jpye, Fri Aug 22 07:23:58 2008 UTC revision 1835 by jpye, Mon Aug 25 08:18:23 2008 UTC
# Line 36  Line 36 
36    
37  /* forward decls */  /* forward decls */
38    
39  static double helm_ideal(double tau, double delta, const HelmholtzData *data);  static double helm_ideal(double tau, double delta, const IdealData *data);
40  static double helm_ideal_tau(double tau, double delta, const HelmholtzData *data);  static double helm_ideal_tau(double tau, double delta, const IdealData *data);
41  static double helm_resid(double tau, double delta, const HelmholtzData *data);  static double helm_resid(double tau, double delta, const HelmholtzData *data);
42  static double helm_resid_del(double tau, double delta, const HelmholtzData *data);  static double helm_resid_del(double tau, double delta, const HelmholtzData *data);
43  static double helm_resid_tau(double tau, double delta, const HelmholtzData *data);  static double helm_resid_tau(double tau, double delta, const HelmholtzData *data);
# Line 74  double helmholtz_u(double T, double rho, Line 74  double helmholtz_u(double T, double rho,
74      double delta = rho / data->rho_star;      double delta = rho / data->rho_star;
75    
76  #ifdef TEST  #ifdef TEST
77      fprintf(stderr,"ideal_tau = %f\n",helm_ideal_tau(tau,delta,data));      fprintf(stderr,"ideal_tau = %f\n",helm_ideal_tau(tau,delta,data->ideal));
78      fprintf(stderr,"resid_tau = %f\n",helm_resid_tau(tau,delta,data));      fprintf(stderr,"resid_tau = %f\n",helm_resid_tau(tau,delta,data));
79      fprintf(stderr,"R T = %f\n",data->R * data->T_star);      fprintf(stderr,"R T = %f\n",data->R * data->T_star);
80  #endif  #endif
81    
82      return data->R * data->T_star * (helm_ideal_tau(tau,delta,data) + helm_resid_tau(tau,delta,data));      return data->R * data->T_star * (helm_ideal_tau(tau,delta,data->ideal) + helm_resid_tau(tau,delta,data));
83  }  }
84    
85  /**  /**
# Line 95  double helmholtz_h(double T, double rho, Line 95  double helmholtz_h(double T, double rho,
95      double tau = data->T_star / T;      double tau = data->T_star / T;
96      double delta = rho / data->rho_star;      double delta = rho / data->rho_star;
97    
98      return data->R * T * (1 + tau * (helm_ideal_tau(tau,delta,data) + helm_resid_tau(tau,delta,data)) + delta*helm_resid_del(tau,delta,data));      return 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));
99  }  }
100    
101  /**  /**
# Line 112  double helmholtz_s(double T, double rho, Line 112  double helmholtz_s(double T, double rho,
112      double delta = rho / data->rho_star;      double delta = rho / data->rho_star;
113    
114      return data->R * (      return data->R * (
115          tau * (helm_ideal_tau(tau,delta,data) + helm_resid_tau(tau,delta,data))          tau * (helm_ideal_tau(tau,delta,data->ideal) + helm_resid_tau(tau,delta,data))
116          - helm_ideal(tau,delta,data) - helm_resid(tau,delta,data)          - helm_ideal(tau,delta,data->ideal) - helm_resid(tau,delta,data)
117      );      );
118  }  }
119    
# Line 151  static double ipow(double x, int n){ Line 151  static double ipow(double x, int n){
151  /**  /**
152      Ideal component of helmholtz function      Ideal component of helmholtz function
153  */    */  
154  double helm_ideal(double tau, double delta, const HelmholtzData *data){  double helm_ideal(double tau, double delta, const IdealData *data){
     const double *a0;  
155    
156      double tau13 = pow(tau,1./3.);      const IdealPowTerm *pt;
157      double taum32 = pow(tau,-3./2.);      const IdealExpTerm *et;
     double taum74 = pow(tau,-7./4.);  
158    
159      a0 = &(data->a0[0]);      unsigned i;
160        double sum = log(delta) + data->c + data->m * tau - log(tau);
161        double term;
162    
163        //fprintf(stderr,"constant = %f, linear = %f", data->c, data->m);
164        //fprintf(stderr,"initial terms = %f\n",sum);
165        pt = &(data->pt[0]);
166        for(i = 0; i<data->np; ++i, ++pt){
167            term = pt->a0 * pow(tau, pt->t0);
168            //fprintf(stderr,"i = %d: a0 = %f, t0 = %f, term = %f\n",i,pt->a0, pt->t0, term);
169            sum += pt->a0 * pow(tau, pt->t0);
170        }
171    
172    #if 0
173        et = &(data->et[0]);
174        for(i=0; i<data->ne; ++i, ++et){
175            sum += et->b * log( 1 - exp(- et->B * tau));
176        }
177    #endif
178    
179      return log(delta) + a0[0] + a0[1] * tau - log(tau) + a0[2] * tau13 + a0[3] * taum32 + a0[4] * taum74;      return sum;
180  }  }
181    
182  /**  /**
183      Partial dervivative of ideal component of helmholtz residual function with      Partial dervivative of ideal component of helmholtz residual function with
184      respect to tau.      respect to tau.
185  */    */  
186  double helm_ideal_tau(double tau, double delta, const HelmholtzData *data){  double helm_ideal_tau(double tau, double delta, const IdealData *data){
187      const double *a0;      const IdealPowTerm *pt;
188        const IdealExpTerm *et;
     double taum114 = pow(tau,-11./4.);  
     double taum52 = pow(tau,-5./2.);  
     double taum23 = pow(tau,-2./3.);  
   
     //fprintf(stderr,"tau = %f, taum23 = %f\n",tau,taum23);  
     //fprintf(stderr,"tau = %f, taum52 = %f\n",tau,taum52);  
     //fprintf(stderr,"tau = %f, taum114 = %f\n",tau,taum114);  
189    
190      a0 = &(data->a0[0]);      unsigned i;
191        double sum = -1./tau + data->m;
192    
193      //unsigned i;   for(i=0;i<5;++i)fprintf(stderr,"a0[%u] = %f\n",i,a0[i]);      pt = &(data->pt[0]);
194        for(i = 0; i<data->np; ++i, ++pt){
195            sum += pt->a0 * pt->t0 * pow(tau, pt->t0 - 1);
196        }
197    
198      double res;  #if 0
199      res = a0[2]/3.*taum23 - 1./tau - 3./2.*a0[3]*taum52 - 7./4.*a0[4]*taum114 + a0[1];      et = &(data->et[0]);
200      //fprintf(stderr,"res = %f\n",res);      for(i=0; i<data->ne; ++i, ++et){
201      return res;          sum += et->b * log( 1 - exp(- et->B * tau));
202        }
203    #endif
204        return sum;
205  }    }  
206    
207  /**  /**
# Line 309  double helm_resid_tau(double tau,double Line 325  double helm_resid_tau(double tau,double
325    
326      delX = 1;      delX = 1;
327    
 #if 1  
328      l = 0;      l = 0;
329      sum = 0;      sum = 0;
330      for(i=0; i<nr; ++i){      for(i=0; i<nr; ++i){
# Line 337  double helm_resid_tau(double tau,double Line 352  double helm_resid_tau(double tau,double
352          }          }
353      }      }
354    
 #else  
     /* old code, not so flexible */  
   
     delX = 1;  
   
     for(i=0; i<5; ++i){  
         if(atdl->t){  
             //fprintf(stderr,"i = %d, a = %e, t = %f, d = %d\n",i+1, atdl->a, atdl->t, atdl->d);  
             res += atdl->a * pow(tau, atdl->t - 1) * ipow(delta, atdl->d) * atdl->t;  
         }  
         ++atdl;  
     }  
   
     sum = 0;  
     for(i=5; i<10; ++i){  
         if(atdl->t){  
             //fprintf(stderr,"i = %d, a = %e, t = %f, d = %d\n",i+1, atdl->a, atdl->t, atdl->d);  
             sum += atdl->a * pow(tau, atdl->t - 1) * ipow(delta, atdl->d) * atdl->t;  
         }  
         ++atdl;  
     }  
     res += exp(-delta) * sum;  
   
     sum = 0;  
     for(i=10; i<17; ++i){  
         if(atdl->t){  
             //fprintf(stderr,"i = %d, a = %e, t = %f, d = %d\n",i+1, atdl->a, atdl->t, atdl->d);  
             sum += atdl->a * pow(tau, atdl->t - 1) * ipow(delta, atdl->d) * atdl->t;  
         }  
         ++atdl;  
     }  
     res += exp(-delta*delta) * sum;  
   
     sum = 0;  
     for(i=17; i<21; ++i){  
         if(atdl->t){  
             //fprintf(stderr,"i = %d, a = %e, t = %f, d = %d\n",i+1, atdl->a, atdl->t, atdl->d);  
             sum += atdl->a * pow(tau, atdl->t - 1) * ipow(delta, atdl->d) * atdl->t;  
         }  
         ++atdl;  
     }  
     res += exp(-delta*delta*delta) * sum;  
 #endif  
   
355      return res;      return res;
356  }    }  
357    

Legend:
Removed from v.1834  
changed lines
  Added in v.1835

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