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

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

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

revision 2736 by jpye, Wed Dec 11 11:09:29 2013 UTC revision 2738 by jpye, Thu Dec 12 08:21:07 2013 UTC
# Line 31  Line 31 
31    
32  #define SQ(X) ((X)*(X))  #define SQ(X) ((X)*(X))
33    
34  //#define CP0_DEBUG  #define CP0_DEBUG
35  #ifdef CP0_DEBUG  #ifdef CP0_DEBUG
36  # include "color.h"  # include "color.h"
37  # define MSG(FMT, ...) \  # define MSG(FMT, ...) \
# Line 62  Line 62 
62      TODO check if ^^^ is still true      TODO check if ^^^ is still true
63  */  */
64  Phi0RunData *cp0_prepare(const IdealData *I, double R, double Tstar){  Phi0RunData *cp0_prepare(const IdealData *I, double R, double Tstar){
65    #ifdef CP0_DEBUG
66        if(I->type==IDEAL_CP0){
67            MSG("R=%f,Tstar=%f (cp0 data: cp0star=%f, Tstar=%f)",R,Tstar,I->data.cp0.cp0star, I->data.cp0.Tstar);
68        }else{
69            MSG("R=%f,Tstar=%f (cp0 data: Tstar=%f)",R,Tstar,I->data.phi0.Tstar);
70        }
71    #endif
72    
73      Phi0RunData *N = FPROPS_NEW(Phi0RunData);      Phi0RunData *N = FPROPS_NEW(Phi0RunData);
74      int i, add_const_term=1;      int i, add_const_term=1;
75      double Tred, cp0red, p;      double Tred, cp0red, p;
# Line 73  Phi0RunData *cp0_prepare(const IdealData Line 81  Phi0RunData *cp0_prepare(const IdealData
81          N->ne = I->data.cp0.ne;          N->ne = I->data.cp0.ne;
82          Tred = I->data.cp0.Tstar;          Tred = I->data.cp0.Tstar;
83          cp0red = I->data.cp0.cp0star;          cp0red = I->data.cp0.cp0star;
84          MSG("Preparing PHI0 data for ideal fluid (np = %d, ne = %d)",N->np, N->ne);          //MSG("Preparing PHI0 data for ideal fluid (np = %d, ne = %d)",N->np, N->ne);
85          MSG("Tred = %f, Tstar = %f", Tred, Tstar);          //MSG("Tred = %f, Tstar = %f", Tred, Tstar);
86          MSG("cp0red = %f, R = %f", cp0red, R);          //MSG("cp0red = %f, R = %f", cp0red, R);
87    
88          for(i=0; i < N->np; ++i)if(I->data.cp0.pt[i].t == 0)add_const_term = 0;          for(i=0; i < N->np; ++i)if(I->data.cp0.pt[i].t == 0)add_const_term = 0;
89          MSG("add_const_term = %d",add_const_term);          //MSG("add_const_term = %d",add_const_term);
90    
91          N->pt = FPROPS_NEW_ARRAY(Phi0RunPowTerm,N->np + add_const_term);          N->pt = FPROPS_NEW_ARRAY(Phi0RunPowTerm,N->np + add_const_term);
92          N->et = FPROPS_NEW_ARRAY(Phi0RunExpTerm, N->ne);          N->et = FPROPS_NEW_ARRAY(Phi0RunExpTerm, N->ne);
# Line 97  Phi0RunData *cp0_prepare(const IdealData Line 105  Phi0RunData *cp0_prepare(const IdealData
105          }          }
106    
107          if(add_const_term){          if(add_const_term){
108              MSG("WARNING: adding constant term %d in cp0, is that what you really want?",i);              //MSG("WARNING: adding constant term %d in cp0, is that what you really want?",i);
109              N->pt[i].a = -1;              N->pt[i].a = -1;
110              N->pt[i].p = 0;              N->pt[i].p = 0;
111              N->np++;              N->np++;
# Line 109  Phi0RunData *cp0_prepare(const IdealData Line 117  Phi0RunData *cp0_prepare(const IdealData
117          }          }
118    
119          if(cp0red != R){          if(cp0red != R){
120              MSG("WARNING: adjusting for R (=%f) != cp0red (=%f)...\n",R,cp0red);              //MSG("WARNING: adjusting for R (=%f) != cp0red (=%f)...\n",R,cp0red);
121              double X = cp0red / R;              double X = cp0red / R;
122              // scale for any differences in R and cpstar */              // scale for any differences in R and cpstar */
123              for(i=0; i < N->np; ++i)N->pt[i].a *= X;              for(i=0; i < N->np; ++i)N->pt[i].a *= X;
# Line 121  Phi0RunData *cp0_prepare(const IdealData Line 129  Phi0RunData *cp0_prepare(const IdealData
129          // TODO add checks for disallowed terms, eg p = 0 or p = 1?          // TODO add checks for disallowed terms, eg p = 0 or p = 1?
130          N->np = I->data.phi0.np;          N->np = I->data.phi0.np;
131          N->ne = I->data.phi0.ne;          N->ne = I->data.phi0.ne;
132          MSG("Preparing PHI0 data for ideal fluid (np = %d, ne = %d)",N->np, N->ne);          //MSG("Preparing PHI0 data for ideal fluid (np = %d, ne = %d)",N->np, N->ne);
133    
134          // power terms          // power terms
135          N->pt = FPROPS_NEW_ARRAY(Phi0RunPowTerm,N->np);          N->pt = FPROPS_NEW_ARRAY(Phi0RunPowTerm,N->np);
# Line 298  double ideal_phi_tautau(double tau, cons Line 306  double ideal_phi_tautau(double tau, cons
306      double sum = 0;      double sum = 0;
307      double term;      double term;
308    
309  #ifdef CP0_DEBUG  #ifdef IDEAL_DEBUG
310      fprintf(stderr,"\ttau = %f\n",tau);      fprintf(stderr,"\ttau = %f\n",tau);
311  #endif  #endif
312    
# Line 310  double ideal_phi_tautau(double tau, cons Line 318  double ideal_phi_tautau(double tau, cons
318          }else{          }else{
319              term = -pt->a * pt->p * (pt->p - 1) * pow(tau, pt->p);              term = -pt->a * pt->p * (pt->p - 1) * pow(tau, pt->p);
320          }          }
321  #ifdef CP0_DEBUG  #ifdef IDEAL_DEBUG
322          fprintf(stderr,"\tpt[%d] = ap(p-1)*tau^p (a = %e, p = %e) = %f\n",i,pt->a, pt->p, term);          fprintf(stderr,"\tpt[%d] = ap(p-1)*tau^p (a = %e, p = %e) = %f\n",i,pt->a, pt->p, term);
323  #endif  #endif
324          sum += term;          sum += term;
# Line 323  double ideal_phi_tautau(double tau, cons Line 331  double ideal_phi_tautau(double tau, cons
331          double e = exp(-x);          double e = exp(-x);
332          double d = (1-e)*(1-e);          double d = (1-e)*(1-e);
333          term = et->n * x*x * e / d;          term = et->n * x*x * e / d;
334  #ifdef CP0_DEBUG  #ifdef IDEAL_DEBUG
335          fprintf(stderr,"\tet[%d] = n x^2 exp(-x)/(1 - exp(-x))^2  (n = %e, x=gamma*tau, gamma = %e) = %f\n",i,et->n, et->gamma, term);          fprintf(stderr,"\tet[%d] = n x^2 exp(-x)/(1 - exp(-x))^2  (n = %e, x=gamma*tau, gamma = %e) = %f\n",i,et->n, et->gamma, term);
336  #endif  #endif
337          sum += term;          sum += term;
338      }      }
339      /* note, at this point, sum == cp0/R - 1 */      /* note, at this point, sum == cp0/R - 1 */
340    #ifdef IDEAL_DEBUG
341      MSG("sum = %f, phi_tautau = %f",sum, -sum/SQ(tau));      MSG("sum = %f, phi_tautau = %f",sum, -sum/SQ(tau));
342    #endif
343      return -sum/SQ(tau);      return -sum/SQ(tau);
344  }  }

Legend:
Removed from v.2736  
changed lines
  Added in v.2738

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