/[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 1897 by jpye, Thu Sep 25 03:46:29 2008 UTC revision 1898 by jpye, Thu Sep 25 07:01:05 2008 UTC
# Line 332  double helm_resid(double tau, double del Line 332  double helm_resid(double tau, double del
332      return res;      return res;
333  }  }
334    
335    #define RESID_DEBUG
336  /**  /**
337      Derivative of the helmholtz residual function with respect to      Derivative of the helmholtz residual function with respect to
338      delta.      delta.
# Line 347  double helm_resid_del(double tau,double Line 348  double helm_resid_del(double tau,double
348      n = data->np;      n = data->np;
349      pt = &(data->pt[0]);      pt = &(data->pt[0]);
350    
351    #ifdef RESID_DEBUG
352            fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
353    #endif
354    
355      sum = 0;      sum = 0;
356      dell = ipow(delta,pt->l);      dell = ipow(delta,pt->l);
357      ldell = pt->l * dell;      ldell = pt->l * dell;
# Line 399  double helm_resid_del(double tau,double Line 404  double helm_resid_del(double tau,double
404  #ifdef RESID_DEBUG  #ifdef RESID_DEBUG
405          fprintf(stderr,"i = %d, GAUSSIAN, n = %e, t = %f, d = %f, alpha = %f, beta = %f, gamma = %f, epsilon = %f\n",i+1, gt->n, gt->t, gt->d, gt->alpha, gt->beta, gt->gamma, gt->epsilon);          fprintf(stderr,"i = %d, GAUSSIAN, n = %e, t = %f, d = %f, alpha = %f, beta = %f, gamma = %f, epsilon = %f\n",i+1, gt->n, gt->t, gt->d, gt->alpha, gt->beta, gt->gamma, gt->epsilon);
406  #endif  #endif
407          double d1 = delta - gt->epsilon;  #define SQ(X) ((X)*(X))
408          double t1 = tau - gt->gamma;          double val2;
409          double e1 = -gt->alpha*d1*d1 - gt->beta*t1*t1;          val2 = - gt->n * pow(tau,gt->t) * pow(delta, -1. + gt->d)
410          double m1 = gt->n * pow(tau,gt->t) * pow(delta,gt->d - 1);              * (2. * gt->alpha * delta * (delta - gt->epsilon) - gt->d)
411          double f1 = -(2.*gt->alpha*delta*delta - 2.*gt->alpha*gt->epsilon*delta - gt->d);              * exp(-(gt->alpha * SQ(delta-gt->epsilon) + gt->beta*SQ(tau-gt->gamma)));
412  #ifdef RESID_DEBUG          res += val2;
         fprintf(stderr,"t1 = %f\n",t1);  
         fprintf(stderr,"d1 = %f\n",d1);  
         fprintf(stderr,"e1 = %f\n",e1);  
         fprintf(stderr,"n = %f, m1 = %f\n", gt->n, m1);  
         fprintf(stderr,"f1 = %f\n",f1);  
 #endif  
         sum = m1 * f1 * exp(e1);  
         res += sum;  
413  #ifdef RESID_DEBUG  #ifdef RESID_DEBUG
414          fprintf(stderr,"sum = %f, res = %f\n",sum,res);          fprintf(stderr,"val2 = %f --> res = %f\n",val2,res);
415  #endif  #endif
416    #undef SQ
417          ++gt;          ++gt;
418      }      }
419  #endif  #endif
# Line 506  double helm_resid_tau(double tau,double Line 504  double helm_resid_tau(double tau,double
504    
505  #define SQ(X) ((X)*(X))  #define SQ(X) ((X)*(X))
506          double val2;          double val2;
507          val2 = -gt->n * pow(tau,gt->t - 1.) * pow(delta, gt->d)          val2 = gt->n * pow(tau,gt->t - 1.) * pow(delta, gt->d)
508              * (2. * gt->beta * SQ(tau) - 2. * gt->beta * gt->gamma * tau - gt->t)              * (2. * gt->beta * SQ(tau) - 2. * gt->beta * gt->gamma * tau + gt->t)
509              * exp(-gt->alpha * SQ(delta-gt->epsilon) + -gt->beta*SQ(tau-gt->gamma));              * exp(-(gt->alpha * SQ(delta-gt->epsilon) + gt->beta*SQ(tau-gt->gamma)));
510            res += val2;
511  #ifdef RESID_DEBUG  #ifdef RESID_DEBUG
512          fprintf(stderr,"res = %f\n",res);          fprintf(stderr,"res = %f\n",res);
513  #endif  #endif

Legend:
Removed from v.1897  
changed lines
  Added in v.1898

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