/[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 1906 by jpye, Sat Sep 27 08:56:29 2008 UTC revision 1907 by jpye, Sat Sep 27 09:45:06 2008 UTC
# Line 484  double helm_resid_tau(double tau,double Line 484  double helm_resid_tau(double tau,double
484    
485  /**  /**
486      Mixed derivative of the helmholtz residual function with respect to      Mixed derivative of the helmholtz residual function with respect to
487      delta and tau      delta and tau.
   
     FIXME this function is WRONG.  
488  */  */
489  double helm_resid_deltau(double tau,double delta,const HelmholtzData *data){  double helm_resid_deltau(double tau,double delta,const HelmholtzData *data){
490      double dell,ldell, term, sum = 0, res = 0;      double dell,ldell, term, sum = 0, res = 0;
# Line 557  double helm_resid_deltau(double tau,doub Line 555  double helm_resid_deltau(double tau,doub
555      FIXME this function is WRONG.      FIXME this function is WRONG.
556  */  */
557  double helm_resid_deldel(double tau,double delta,const HelmholtzData *data){  double helm_resid_deldel(double tau,double delta,const HelmholtzData *data){
558            double sum = 0, res = 0;
559      double sum;      double dell, ldell;
560      double phir = 0;      unsigned n, i;
561      unsigned i;      const HelmholtzPowTerm *pt;
562      unsigned X;      const HelmholtzGausTerm *gt;
     double XdelX;  
   
     const HelmholtzPowTerm *pt = &(data->pt[0]);  
       
     for(i=0; i<5; ++i){  
         phir += pt->a * pow(tau, pt->t) * ipow(delta, pt->d - 2) * (SQ(pt->d) - X);  
         ++pt;  
     }  
563    
     sum = 0;  
     X = 1;  
     XdelX = delta;  
     for(i=5; i<10; ++i){  
         sum += pt->a * pow(tau, pt->t) * ipow(delta, pt->d - 2) * (SQ(XdelX) - X*XdelX - 2*pt->d*XdelX + XdelX + SQ(pt->d) - pt->d);  
         ++pt;  
     }  
     phir += exp(-delta) * sum;  
564    
565      sum = 0;  #ifdef RESID_DEBUG
566      X = 2;          fprintf(stderr,"tau=%f, del=%f\n",tau,delta);
567      XdelX = 2*delta*delta;  #endif
568      for(i=10; i<17; ++i){  
569          sum += pt->a * pow(tau, pt->t) * ipow(delta, pt->d - 2) * (SQ(XdelX) - X*XdelX - 2*pt->d*XdelX + XdelX + SQ(pt->d) - pt->d);      /* power terms */
570        n = data->np;
571        pt = &(data->pt[0]);
572        dell = ipow(delta,pt->l);
573        ldell = pt->l * dell;
574        unsigned oldl;
575        for(i=0; i<n; ++i){
576            double lpart = pt->l ? SQ(ldell) + ldell*(1. - 2*pt->d - pt->l) : 0;
577            sum += pt->a * pow(tau, pt->t) * ipow(delta, pt->d - 2) * (pt->d*(pt->d - 1) + lpart);
578            oldl = pt->l;
579          ++pt;          ++pt;
580            if(i+1==n || oldl != pt->l){
581                if(oldl == 0){
582                    res += sum;
583                }else{
584                    res += sum * exp(-dell);
585                }
586                sum = 0;
587                dell = ipow(delta,pt->l);
588                ldell = pt->l*dell;
589            }
590      }      }
     phir += exp(-delta*delta) * sum;  
591    
592      sum = 0;      /* gaussian terms */
593      X = 3;      n = data->ng;
594      XdelX = 3*delta*delta*delta;      //fprintf(stderr,"THERE ARE %d GAUSSIAN TERMS\n",n);
595      for(i=17; i<21; ++i){      gt = &(data->gt[0]);
596          sum += pt->a * pow(tau, pt->t) * ipow(delta, pt->d - 2) * (SQ(XdelX) - X*XdelX - 2*pt->d*XdelX + XdelX + SQ(pt->d) - pt->d);      for(i=0; i<n; ++i){
597          ++pt;          double s1 = SQ(delta - gt->epsilon);
598            double f1 = 2*delta*gt->alpha *(2*gt->d*gt->epsilon
599                - delta * (2*gt->d + 1 - 2 * gt->alpha * s1));
600            res += - gt->n * pow(tau,gt->t) * pow(delta, -1. + gt->d)
601                * f1
602                * exp(-(gt->alpha * s1 + gt->beta*SQ(tau-gt->gamma)));
603            ++gt;
604      }      }
     phir += exp(-delta*delta*delta) * sum;  
605    
606      return phir;      return res;
607  }  }
608    

Legend:
Removed from v.1906  
changed lines
  Added in v.1907

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