/[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 1848 by jpye, Mon Sep 1 05:42:21 2008 UTC revision 1849 by jpye, Mon Sep 1 14:48:55 2008 UTC
# Line 1  Line 1 
1    /*  ASCEND modelling environment
2        Copyright (C) 2008 Carnegie Mellon University
3    
4        This program is free software; you can redistribute it and/or modify
5        it under the terms of the GNU General Public License as published by
6        the Free Software Foundation; either version 2, or (at your option)
7        any later version.
8    
9        This program is distributed in the hope that it will be useful,
10        but WITHOUT ANY WARRANTY; without even the implied warranty of
11        MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12        GNU General Public License for more details.
13    
14        You should have received a copy of the GNU General Public License
15        along with this program; if not, write to the Free Software
16        Foundation, Inc., 59 Temple Place - Suite 330,
17        Boston, MA 02111-1307, USA.
18    */
19    
20  #include <math.h>  #include <math.h>
21    
22  #include "ideal.h"  #include "ideal.h"
# Line 12  Line 31 
31    IDEAL COMPONENT RELATIONS    IDEAL COMPONENT RELATIONS
32  */  */
33    
34    /*
35        Hypothesis:
36        in calculating the ideal component relations, we have some integration
37        constants that ultimately allow the arbitrary scales of h and s to
38        be specified. Hence we can ignore components of helm_ideal that
39        are either constant or linear in tau, and work them out later when we
40        stick in the values of data->c and data->m.
41    */
42    
43  /**  /**
44      Ideal component of helmholtz function      Ideal component of helmholtz function
45  */    */  
# Line 21  double helm_ideal(double tau, double del Line 49  double helm_ideal(double tau, double del
49      const IdealExpTerm *et;      const IdealExpTerm *et;
50    
51      unsigned i;      unsigned i;
52      double sum = log(delta) + data->c + data->m * tau - log(tau);      double sum = log(delta) - log(tau) + data->c + data->m * tau;
53      double term;      double term;
54    
55      //fprintf(stderr,"constant = %f, linear = %f", data->c, data->m);      //fprintf(stderr,"constant = %f, linear = %f", data->c, data->m);
56      //fprintf(stderr,"initial terms = %f\n",sum);      //fprintf(stderr,"initial terms = %f\n",sum);
57      pt = &(data->pt[0]);      pt = &(data->pt[0]);
58      for(i = 0; i<data->np; ++i, ++pt){      for(i = 0; i<data->np; ++i, ++pt){
59          term = pt->a0 * pow(tau, pt->t0);          double a = pt->a0;
60            double t = pt->t0;
61            if(t == 0){
62                term = a*log(tau) /* + a */; /* term ignored, see above */
63            }else if(t == 1){
64                term = /* a * tau */ - a * tau * log(tau); /* term ignored, see above */
65            }else{
66                term = a * pow(tau,t) / ((t - 1) * t);
67            }
68          //fprintf(stderr,"i = %d: a0 = %f, t0 = %f, term = %f\n",i,pt->a0, pt->t0, term);          //fprintf(stderr,"i = %d: a0 = %f, t0 = %f, term = %f\n",i,pt->a0, pt->t0, term);
69          sum += pt->a0 * pow(tau, pt->t0);          sum += pt->a0 * pow(tau, pt->t0);
70      }      }
71    
72  #if 1      /* 'exponential' terms */
73      et = &(data->et[0]);      et = &(data->et[0]);
74      for(i=0; i<data->ne; ++i, ++et){      for(i=0; i<data->ne; ++i, ++et){
75          sum += et->b * log( 1 - exp(- et->B * tau));          sum += et->b * log( 1 - exp(-et->B * tau));
76      }      }
 #endif  
77    
78      return sum;      return sum;
79  }  }
# Line 52  double helm_ideal_tau(double tau, double Line 87  double helm_ideal_tau(double tau, double
87      const IdealExpTerm *et;      const IdealExpTerm *et;
88    
89      unsigned i;      unsigned i;
90        double term;
91      double sum = -1./tau + data->m;      double sum = -1./tau + data->m;
92    
93      pt = &(data->pt[0]);      pt = &(data->pt[0]);
94      for(i = 0; i<data->np; ++i, ++pt){      for(i = 0; i<data->np; ++i, ++pt){
95          sum += pt->a0 * pt->t0 * pow(tau, pt->t0 - 1);          double a = pt->a0;
96            double t = pt->t0;
97            if(t==0){
98                term = a / tau - a*t*pow(tau,t-1)/(t-1);
99            }else if(t==1){
100                term = -a*log(tau);
101            }else{
102                term = a*t*pow(tau,t-1)/(t-1);
103            }
104            sum += term;
105      }      }
106    
107  #if 1      /* 'exponential' terms */
     /* FIXME we're missing the '2.5' in the alpha0 expression for Nitrogen... */  
108      et = &(data->et[0]);      et = &(data->et[0]);
109      for(i=0; i<data->ne; ++i, ++et){      for(i=0; i<data->ne; ++i, ++et){
110          sum += et->b * log( 1 - exp(- et->B * tau));          sum += et->b * et->B  / (1 - exp(-tau*et->B));
111      }      }
 #endif  
112      return sum;      return sum;
113  }    }
Legend:
Removed from v.1848  
changed lines
  Added in v.1849
john.pye@anu.edu.au
 ViewVC Help
Powered by ViewVC 1.1.22