/[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 1835 by jpye, Mon Aug 25 08:18:23 2008 UTC revision 1836 by jpye, Mon Aug 25 08:32:38 2008 UTC
# Line 209  double helm_ideal_tau(double tau, double Line 209  double helm_ideal_tau(double tau, double
209      optimised here ;-)      optimised here ;-)
210  */  */
211  double helm_resid(double tau, double delta, const HelmholtzData *data){  double helm_resid(double tau, double delta, const HelmholtzData *data){
       
212      double sum;      double sum;
213      double phir = 0;      double res = 0;
214      unsigned i;      double delX;
215        unsigned l;
216        unsigned nr, i;
217        const HelmholtzATDL *atdl;
218    
219      const HelmholtzATDL *atdl = &(data->atdl[0]);      nr = data->nr;
220            atdl = &(data->atdl[0]);
     for(i=0; i<5; ++i){  
         phir += atdl->a * pow(tau, atdl->t) * ipow(delta, atdl->d);  
         ++atdl;  
     }  
221    
222      sum = 0;      delX = 1;
     for(i=5; i<10; ++i){  
         sum += atdl->a * pow(tau, atdl->t) * ipow(delta, atdl->d);  
         ++atdl;  
     }  
     phir += exp(-delta) * sum;  
   
     sum = 0;  
     for(i=10; i<17; ++i){  
         sum += atdl->a * pow(tau, atdl->t) * ipow(delta, atdl->d);  
         ++atdl;  
     }  
     phir += exp(-delta*delta) * sum;  
223    
224        l = 0;
225      sum = 0;      sum = 0;
226      for(i=17; i<21; ++i){      for(i=0; i<nr; ++i){
227            //fprintf(stderr,"i = %d, a = %e, t = %f, d = %d, l = %d\n",i+1, atdl->a, atdl->t, atdl->d, atdl->l);
228          sum += atdl->a * pow(tau, atdl->t) * ipow(delta, atdl->d);          sum += atdl->a * pow(tau, atdl->t) * ipow(delta, atdl->d);
229          ++atdl;          ++atdl;
230            //fprintf(stderr,"l = %d\n",l);
231            if(i+1==nr || l != atdl->l){
232                if(l==0){
233                    //fprintf(stderr,"Adding non-exp term\n");
234                    res += sum;
235                }else{
236                    //fprintf(stderr,"Adding exp term with l = %d, delX = %e\n",l,delX);
237                    res += sum * exp(-delX);
238                }
239                /* set l to new value */
240                if(i+1!=nr){
241                    l = atdl->l;
242                    //fprintf(stderr,"New l = %d\n",l);
243                    delX = ipow(delta,l);
244                    sum = 0;
245                }
246            }
247      }      }
     phir += exp(-delta*delta*delta) * sum;  
248    
249      return phir;      return res;
250  }  }
251    
252  /**  /**
253      Derivative of the helmholtz residual function with respect to      Derivative of the helmholtz residual function with respect to
254      delta.      delta.
255    
256      THERE IS AN ERROR IN THIS FUNCTION.      THERE APPEARS TO BE AN ERROR IN THIS FUNCTION.
257  */    */  
258  double helm_resid_del(double tau,double delta, const HelmholtzData *data){  double helm_resid_del(double tau,double delta, const HelmholtzData *data){
       
259      double sum;      double sum;
260      double phir = 0;      double res = 0;
261      unsigned i;      double delX, XdelX;
262      double X;      unsigned l;
263      double delX;      unsigned nr, i;
264      double XdelX;      const HelmholtzATDL *atdl;
265    
266      const HelmholtzATDL *atdl = &(data->atdl[0]);      nr = data->nr;
267            atdl = &(data->atdl[0]);
     for(i=0; i<5; ++i){  
         //fprintf(stderr,"i = %d, a = %e, t = %f, d = %d\n",i+1, atdl->a, atdl->t, atdl->d);  
         phir += atdl->a * pow(tau, atdl->t) * pow(delta, atdl->d - 1) * atdl->d;  
         ++atdl;  
     }  
   
     sum = 0;  
     X = 1;  
     delX = ipow(delta,X);  
     XdelX = X*delX;  
     for(i=5; i<10; ++i){  
         //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) * pow(delta, atdl->d - 1) * (XdelX - atdl->d);  
         ++atdl;  
     }  
     //fprintf(stderr,"sum = %f\n",sum);  
     phir += exp(-delX) * sum;  
268    
269      sum = 0;      delX = 1;
     X = 2;  
     delX = ipow(delta,X);  
     XdelX = X*delX;  
     for(i=10; i<17; ++i){  
         //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) * pow(delta, atdl->d - 1) * (XdelX - atdl->d);  
         ++atdl;  
     }  
     //fprintf(stderr,"sum = %f\n",sum);  
     phir += exp(-delX) * sum;  
270    
271        l = 0;
272      sum = 0;      sum = 0;
273      X = 3;      XdelX = 0;
274      delX = ipow(delta,X);      for(i=0; i<nr; ++i){
275      XdelX = X*delX;          //fprintf(stderr,"i = %d, a = %e, t = %f, d = %d, l = %d\n",i+1, atdl->a, atdl->t, atdl->d, atdl->l);
276      for(i=17; i<21; ++i){          sum += atdl->a * pow(tau, atdl->t) * ipow(delta, atdl->d - 1) * (atdl->d - XdelX);
277          //fprintf(stderr,"i = %d, a = %e, t = %f, d = %d\n",i+1, atdl->a, atdl->t, atdl->d);          ++atdl;
278          sum -= atdl->a * pow(tau, atdl->t) * pow(delta, atdl->d - 1) * (XdelX - atdl->d);          //fprintf(stderr,"l = %d\n",l);
279          ++atdl;          if(i+1==nr || l != atdl->l){
280                if(l==0){
281                    //fprintf(stderr,"Adding non-exp term\n");
282                    //fprintf(stderr,"sum = %f\n",sum);
283                    res += sum;
284                }else{
285                    //fprintf(stderr,"Adding exp term with l = %d, delX = %e\n",l,delX);
286                    //fprintf(stderr,"sum = %f\n",sum);
287                    res += sum * exp(-delX);
288                }
289                /* set l to new value */
290                if(i+1!=nr){
291                    l = atdl->l;
292                    delX = ipow(delta,l);
293                    XdelX = l * delX;
294                    //fprintf(stderr,"New l = %d, XdelX = %f\n",l,XdelX);
295                    sum = 0;
296                }
297            }
298      }      }
     //fprintf(stderr,"sum = %f\n",sum);  
     phir += exp(-delX) * sum;  
299    
300      return phir;      return res;
301  }  }
302    
303  /**  /**

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

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