/[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 1840 by jpye, Fri Aug 29 03:52:04 2008 UTC revision 1841 by jpye, Fri Aug 29 07:27:59 2008 UTC
# Line 69  double helmholtz_p(double T, double rho, Line 69  double helmholtz_p(double T, double rho,
69      //fprintf(stderr,"p calc: rho = %f\n",rho);      //fprintf(stderr,"p calc: rho = %f\n",rho);
70      //fprintf(stderr,"p calc: delta = %f\n",delta);      //fprintf(stderr,"p calc: delta = %f\n",delta);
71      //fprintf(stderr,"p calc: R*T*rho = %f\n",data->R * T * rho);      //fprintf(stderr,"p calc: R*T*rho = %f\n",data->R * T * rho);
72    
73        //fprintf(stderr,"T = %f\n", T);
74        //fprintf(stderr,"rhob = %f, rhob* = %f, delta = %f\n", rho/data->M, data->rho_star/data->M, delta);
75  #endif  #endif
76            
77      return data->R * T * rho * (1. + delta * helm_resid_del(tau,delta,data));      return data->R * T * rho * (1 + delta * helm_resid_del(tau,delta,data));
78  }  }
79    
80  /**  /**
# Line 294  double helm_resid(double tau, double del Line 297  double helm_resid(double tau, double del
297      NOTE: POSSIBLY STILL AN ERROR IN THIS FUNCTION.      NOTE: POSSIBLY STILL AN ERROR IN THIS FUNCTION.
298  */    */  
299  double helm_resid_del(double tau,double delta, const HelmholtzData *data){  double helm_resid_del(double tau,double delta, const HelmholtzData *data){
300      double sum;      double term, x, sum;
301      double res = 0;      double res = 0;
302      double delX, XdelX;      double delX, XdelX;
303      unsigned l;      unsigned l;
# Line 305  double helm_resid_del(double tau,double Line 308  double helm_resid_del(double tau,double
308      n = data->np;      n = data->np;
309      pt = &(data->pt[0]);      pt = &(data->pt[0]);
310    
     delX = 1;  
311    
     l = 0;  
312      sum = 0;      sum = 0;
313    
314    
315        for(i=0; i<n; ++i){
316            term = pt->a * pow(tau, pt->t) * pow(delta, pt->d - 1) * (pt->d - pt->l*pow(delta,pt->l));
317            if(pt->l==0){  
318                x = 1;
319            }else{
320                x = exp(-pow(delta,pt->l));
321            }
322            //fprintf(stderr,"i = %d, a = %e, t = %f, d = %d, l = %d --> pow = %f, exp = %f, term = %f\n",i+1, pt->a, pt->t, pt->d, pt->l, term, x, term*x);
323            res += term * x;
324            ++pt;
325        }
326    #if 0
327        delX = 1;
328        l = 0;
329      XdelX = 0;      XdelX = 0;
330      for(i=0; i<n; ++i){      for(i=0; i<n; ++i){
331          //fprintf(stderr,"i = %d, a = %e, t = %f, d = %d, l = %d\n",i+1, pt->a, pt->t, pt->d, pt->l);          term = pt->a * pow(tau, pt->t) * ipow(delta, pt->d - 1) * (pt->d - XdelX);
332          sum += pt->a * pow(tau, pt->t) * pow(delta, pt->d - 1) * (pt->d - XdelX);          fprintf(stderr,"i = %d, a = %e, t = %f, d = %d, l = %d --> term = %f\n",i+1, pt->a, pt->t, pt->d, pt->l, term);
333            sum += term;
334          ++pt;          ++pt;
335          //fprintf(stderr,"l = %d\n",l);          //fprintf(stderr,"l = %d\n",l);
336          if(i+1==n || l != pt->l){          if(i+1==n || l != pt->l){
337              if(l==0){              if(l==0){
338                  //fprintf(stderr,"Adding non-exp term\n");                  //fprintf(stderr,"Adding non-exp term\n");
                 //fprintf(stderr,"sum = %f\n",sum);  
339                  res += sum;                  res += sum;
340              }else{              }else{
341                  //fprintf(stderr,"Adding exp term with l = %d, delX = %e\n",l,delX);                  //fprintf(stderr,"Adding exp term with l = %d, delX = %e\n",l,delX);
                 //fprintf(stderr,"sum = %f\n",sum);  
342                  res += sum * exp(-delX);                  res += sum * exp(-delX);
343              }              }
344                fprintf(stderr,"sum = %f, mult = %f, res = %f\n",sum,exp(-delX),res);
345              /* set l to new value */              /* set l to new value */
346              if(i+1!=n){              if(i+1!=n){
347                  l = pt->l;                  l = pt->l;
# Line 332  double helm_resid_del(double tau,double Line 349  double helm_resid_del(double tau,double
349                  XdelX = l * delX;                  XdelX = l * delX;
350                  //fprintf(stderr,"New l = %d, XdelX = %f\n",l,XdelX);                  //fprintf(stderr,"New l = %d, XdelX = %f\n",l,XdelX);
351                  sum = 0;                  sum = 0;
352                    fprintf(stderr,"sum zero\n");
353              }              }
354          }          }
355      }      }
356    #endif
357    
358  #if 0  #if 1
359      /* now the exponential terms */      /* now the exponential terms */
360      n = data->ne;      n = data->ne;
361      et = &(data->et[0]);      et = &(data->et[0]);
362      for(i=0; i< n; ++i){      for(i=0; i< n; ++i){
363          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);          //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);
364                    
365          double del2 = delta*delta;          double del2 = delta*delta;
366          double tau2 = tau*tau;          double tau2 = tau*tau;
# Line 355  double helm_resid_del(double tau,double Line 374  double helm_resid_del(double tau,double
374                       - et->phi                       - et->phi
375                       - et->beta * gam2                       - et->beta * gam2
376                 );                 );
377          fprintf(stderr,"sum = %f\n",sum);          //fprintf(stderr,"sum = %f\n",sum);
378          res += sum;          res += sum;
379          ++et;          ++et;
380      }      }
Legend:
Removed from v.1840  
changed lines
  Added in v.1841
john.pye@anu.edu.au
 ViewVC Help
Powered by ViewVC 1.1.22