/[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 1903 by jpye, Thu Sep 25 07:34:56 2008 UTC revision 1904 by jpye, Sat Sep 27 03:32:21 2008 UTC
# Line 238  double helm_resid(double tau, double del Line 238  double helm_resid(double tau, double del
238      double dell,ldell, term, sum, res = 0;      double dell,ldell, term, sum, res = 0;
239      unsigned n, i;      unsigned n, i;
240      const HelmholtzPowTerm *pt;      const HelmholtzPowTerm *pt;
     const HelmholtzExpTerm *et;  
241      const HelmholtzGausTerm *gt;      const HelmholtzGausTerm *gt;
242    
243      n = data->np;      n = data->np;
# Line 287  double helm_resid(double tau, double del Line 286  double helm_resid(double tau, double del
286          }          }
287      }      }
288    
     /* now the exponential terms */  
     n = data->ne;  
     //fprintf(stderr,"THERE ARE %d EXPONENTIAL TERMS at %p\n",n, data->et);  
     et = &(data->et[0]);  
     for(i=0; i< n; ++i){  
 #ifdef RESID_DEBUG  
         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);  
 #endif        
         double e1 = -et->phi * delta*delta  
                      + 2 * et->phi * delta  
                      - et->beta * tau * tau  
                      + 2 * et->beta * et->gamma * tau  
                      - et->phi  
                      - et->beta * et->gamma * et->gamma;  
         sum = et->a * pow(tau,et->t) * ipow(delta,et->d) * exp(e1);  
         //fprintf(stderr,"sum = %f\n",sum);  
         res += sum;  
         ++et;  
     }  
   
289  #if 1  #if 1
290      /* gaussian terms */      /* gaussian terms */
291      n = data->ng;      n = data->ng;
# Line 341  double helm_resid_del(double tau,double Line 320  double helm_resid_del(double tau,double
320      double dell, ldell;      double dell, ldell;
321      unsigned n, i;      unsigned n, i;
322      const HelmholtzPowTerm *pt;      const HelmholtzPowTerm *pt;
     const HelmholtzExpTerm *et;  
323      const HelmholtzGausTerm *gt;      const HelmholtzGausTerm *gt;
324    
325      n = data->np;      n = data->np;
# Line 371  double helm_resid_del(double tau,double Line 349  double helm_resid_del(double tau,double
349          }          }
350      }      }
351    
     /* exponential terms */  
     n = data->ne;  
     et = &(data->et[0]);  
     for(i=0; i< n; ++i){  
         //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);  
           
         double del2 = delta*delta;  
         double tau2 = tau*tau;  
         double gam2 = et->gamma * et->gamma;  
         double e1 = -et->phi * del2  
                      + 2 * et->phi * delta  
                      - et->beta * tau2  
                      + 2 * et->beta * et->gamma * tau  
                      - et->phi  
                      - et->beta * gam2;  
         sum = -et->a * pow(tau,et->t) * ipow(delta,et->d-1)  
             * (2 * et->phi * del2 - 2 * et->phi * delta - et->d)  
             * exp(e1);  
         //fprintf(stderr,"sum = %f\n",sum);  
         res += sum;  
         ++et;  
     }  
   
352  #if 1  #if 1
353      /* gaussian terms */      /* gaussian terms */
354      n = data->ng;      n = data->ng;
# Line 432  double helm_resid_tau(double tau,double Line 387  double helm_resid_tau(double tau,double
387      unsigned l;      unsigned l;
388      unsigned n, i;      unsigned n, i;
389      const HelmholtzPowTerm *pt;      const HelmholtzPowTerm *pt;
     const HelmholtzExpTerm *et;  
390      const HelmholtzGausTerm *gt;      const HelmholtzGausTerm *gt;
391    
392      n = data->np;      n = data->np;
# Line 467  double helm_resid_tau(double tau,double Line 421  double helm_resid_tau(double tau,double
421          }          }
422      }      }
423    
 #if 1  
     /* now the exponential terms */  
     n = data->ne;  
     et = &(data->et[0]);  
     for(i=0; i< n; ++i){  
         //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);  
           
         double tau2 = tau*tau;  
         double del2 = delta*delta;  
         double gam2 = et->gamma * et->gamma;  
         double e1 = -et->phi * del2  
                      + 2 * et->phi * delta  
                      - et->beta * tau2  
                      + 2 * et->beta * et->gamma * tau  
                      - et->phi  
                      - et->beta * gam2;  
         sum = -et->a * pow(tau,et->t - 1) * ipow(delta,et->d)  
             * (2 * et->beta * tau2 - 2 * et->beta * et->gamma * tau - et->t)  
             * exp(e1);  
         //fprintf(stderr,"sum = %f\n",sum);  
         res += sum;  
         ++et;  
     }  
 #endif  
   
424  //#define RESID_DEBUG  //#define RESID_DEBUG
425      /* gaussian terms */      /* gaussian terms */
426      n = data->ng;      n = data->ng;
Legend:
Removed from v.1903  
changed lines
  Added in v.1904
john.pye@anu.edu.au
 ViewVC Help
Powered by ViewVC 1.1.22