/[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 1985 by jpye, Tue Feb 3 01:15:17 2009 UTC revision 1986 by jpye, Tue Feb 3 05:19:15 2009 UTC
# Line 414  double helm_resid(double tau, double del Line 414  double helm_resid(double tau, double del
414          double d1 = delta - 1.;          double d1 = delta - 1.;
415          double t1 = tau - 1.;          double t1 = tau - 1.;
416          double theta = (1. - tau) + ct->A * pow(d1*d1, 0.5/ct->beta);          double theta = (1. - tau) + ct->A * pow(d1*d1, 0.5/ct->beta);
417          double psi = exp(-ct->A*d1*d1 - ct->D*t1*t1);          double psi = exp(-ct->C*d1*d1 - ct->D*t1*t1);
418          double DELTA = theta*theta + ct->B* pow(d1*d1, ct->a);          double DELTA = theta*theta + ct->B* pow(d1*d1, ct->a);
419          sum = ct->n * pow(DELTA, ct->b) * delta * psi;          sum = ct->n * pow(DELTA, ct->b) * delta * psi;
420          res += sum;          res += sum;
# Line 439  double helm_resid_del(double tau,double Line 439  double helm_resid_del(double tau,double
439      unsigned n, i;      unsigned n, i;
440      const HelmholtzPowTerm *pt;      const HelmholtzPowTerm *pt;
441      const HelmholtzGausTerm *gt;      const HelmholtzGausTerm *gt;
442        const HelmholtzCritTerm *ct;
443    
444  #ifdef RESID_DEBUG  #ifdef RESID_DEBUG
445          fprintf(stderr,"tau=%f, del=%f\n",tau,delta);          fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
# Line 469  double helm_resid_del(double tau,double Line 469  double helm_resid_del(double tau,double
469    
470      /* gaussian terms */      /* gaussian terms */
471      n = data->ng;      n = data->ng;
     //fprintf(stderr,"THERE ARE %d GAUSSIAN TERMS\n",n);  
472      gt = &(data->gt[0]);      gt = &(data->gt[0]);
473      for(i=0; i<n; ++i){      for(i=0; i<n; ++i){
474  #ifdef RESID_DEBUG  #ifdef RESID_DEBUG
475          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);
476  #endif  #endif
477          double val2;          sum = - gt->n * pow(tau,gt->t) * pow(delta, -1. + gt->d)
         val2 = - gt->n * pow(tau,gt->t) * pow(delta, -1. + gt->d)  
478              * (2. * gt->alpha * delta * (delta - gt->epsilon) - gt->d)              * (2. * gt->alpha * delta * (delta - gt->epsilon) - gt->d)
479              * 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)));
480          res += val2;          res += sum;
481  #ifdef RESID_DEBUG  #ifdef RESID_DEBUG
482          fprintf(stderr,"val2 = %f --> res = %f\n",val2,res);          fprintf(stderr,"sum = %f --> res = %f\n",sum,res);
483  #endif  #endif
484          ++gt;          ++gt;
485      }      }
486    
487      /* FIXME add critical terms calculation */      /* critical terms */
488        n = data->nc;
489        ct = &(data->ct[0]);
490        for(i=0; i<n; ++i){
491    #ifdef RESID_DEBUG
492            fprintf(stderr,"i = %d, CRITICAL, n = %e, a = %f, b = %f, beta = %f, A = %f, B = %f, C = %f, D = %f\n",i+1, ct->n, ct->a, ct->b, ct->beta, ct->A, ct->B, ct->C, ct->D);
493    #endif
494            double d1 = delta - 1.;
495            double t1 = tau - 1.;
496            double theta = (1. - tau) + ct->A * pow(d1*d1, 0.5/ct->beta);
497            double psi = exp(-ct->C*d1*d1 - ct->D*t1*t1);
498            double DELTA = theta*theta + ct->B* pow(d1*d1, ct->a);
499    
500            double dpsiddelta = -2. * ct->C * d1 * psi;
501    
502            double dDELddelta = d1 * (ct->A * theta * 2./ct->beta * pow(d1*d1, 0.5/ct->beta - 1) + 2* ct->B * ct->a * pow(d1*d1, ct->a - 1));
503    
504            double dDELbddelta = ct->b * pow(DELTA,ct->b - 1.) * dDELddelta;
505    
506            sum = ct->n * (pow(DELTA, ct->b) * (psi + delta * dpsiddelta) + dDELbddelta * delta * psi);
507            res += sum;
508            ++ct;
509        }
510    
511    
512      return res;      return res;
513  }  }
Legend:
Removed from v.1985  
changed lines
  Added in v.1986
john.pye@anu.edu.au
 ViewVC Help
Powered by ViewVC 1.1.22