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

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

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

revision 1864 by jpye, Sun Sep 14 08:39:51 2008 UTC revision 1865 by jpye, Mon Sep 15 08:40:14 2008 UTC
# Line 50  double helm_cp0(double T, const IdealDat Line 50  double helm_cp0(double T, const IdealDat
50      fprintf(stderr,"np = %d\n",data->np);      fprintf(stderr,"np = %d\n",data->np);
51  #endif  #endif
52      for(i = 0; i<data->np; ++i, ++pt){      for(i = 0; i<data->np; ++i, ++pt){
53          term = pt->a0 * pow(T, pt->t0);          term = pt->c * pow(T, pt->t);
54  #if 0  #if 0
55          fprintf(stderr,"i = %d: ",i);          fprintf(stderr,"i = %d: ",i);
56          fprintf(stderr,"power term, a = %f, t = %f, val = %f\n",pt->a0, pt->t0, term);          fprintf(stderr,"power term, c = %f, t = %f, val = %f\n",pt->c, pt->t, term);
57  #endif  #endif
58          sum += term;          sum += term;
59      }      }
# Line 105  double helm_ideal(double tau, double del Line 105  double helm_ideal(double tau, double del
105      /* power terms */      /* power terms */
106      pt = &(data->pt[0]);      pt = &(data->pt[0]);
107      for(i = 0; i<data->np; ++i, ++pt){      for(i = 0; i<data->np; ++i, ++pt){
108          double a = pt->a0;          double c = pt->c;
109          double t = pt->t0;          double t = pt->t;
110          if(t == 0){          if(t == 0){
111              term = a*log(tau);              term = c*log(tau);
112          }else{          }else{
113  #ifdef TEST  #ifdef TEST
114              assert(t!=-1);              assert(t!=-1);
115  #endif  #endif
116              term = -a / (t*(t+1)) * pow(Tstar_on_tau,t);              term = -c / (t*(t+1)) * pow(Tstar_on_tau,t);
117                fprintf(stderr,"i = %d, c = %f, t = %f, term = %f\n",i,c,t,term);
118          }          }
119          sum += pt->a0 * pow(tau, pt->t0);          sum += term;
120      }      }
121    
122      /* 'exponential' terms */      /* 'exponential' terms */
123      et = &(data->et[0]);      et = &(data->et[0]);
124      for(i=0; i<data->ne; ++i, ++et){      for(i=0; i<data->ne; ++i, ++et){
125          sum += et->b * log(1 - exp(-et->beta / Tstar_on_tau));          term = et->b * log(1 - exp(-et->beta / Tstar_on_tau));
126            fprintf(stderr,"exp i=%d, b=%f, beta=%f, term = %f\n",i,et->b, et->beta, term);
127            sum += term;
128      }      }
129    
130    #ifdef TEST
131        fprintf(stderr,"phi0 = %f\n",sum);
132    #endif
133    
134      return sum;      return sum;
135  }  }
136    
# Line 142  double helm_ideal_tau(double tau, double Line 149  double helm_ideal_tau(double tau, double
149    
150      pt = &(data->pt[0]);      pt = &(data->pt[0]);
151      for(i = 0; i<data->np; ++i, ++pt){      for(i = 0; i<data->np; ++i, ++pt){
152          double a = pt->a0;          double c = pt->c;
153          double t = pt->t0;          double t = pt->t;
154          if(t==0){          if(t==0){
155              term = a / tau;              term = c / tau;
156          }else{          }else{
157              // term = -a / (t*(t+1)) * pow(Tstar_on_tau,t);              // term = -c / (t*(t+1)) * pow(Tstar_on_tau,t);
158              term = a/(t+1)*pow(Tstar_on_tau,t)/tau;              term = c/(t+1)*pow(Tstar_on_tau,t)/tau;
159          }          }
160  #ifdef TEST  #ifdef TEST
161          if(isinf(term)){          if(isinf(term)){
162              fprintf(stderr,"Error with infinite-valued term with i = %d, a = %f, t = %f\n", i,a ,t);              fprintf(stderr,"Error with infinite-valued term with i = %d, c = %f, t = %f\n", i,c ,t);
163              abort();              abort();
164          }          }
165  #endif  #endif

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

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