Changeset 476 for trunk/matml

Show
Ignore:
Timestamp:
10/05/2009 07:44:57 PM (3 years ago)
Author:
powell
Message:

New macro eliminates singular returns of pow() for exponents below 2.

Files:
1 modified

Legend:

Unmodified
Added
Removed

  • trunk/matml/src/ternary/freenergy.c

    r475 r476  
    1717  (((flags) == NO_INFINITY || (flags) == WITH_DERIVATIVES_NO_INFINITY) ? \ 
    1818   (-1e100) : TERNARY_NEGATIVE_INFINITY) /*+ "Negative infinity" definition for energy derivatives +*/ 
     19 
     20#define POWER2(A,B) \ 
     21  (((B)<0 && (A) <= 0.) ? LOCAL_INFINITY : pow ((A),(B))) \ 
     22  /*+ Contingency in pow for polynomial second derivatives +*/ 
    1923 
    2024/*++++++++++++++++++++++++++++++++++++++ 
    … …  
    281285  double ret = 
    282286    eparams->R*T*(((eparams->S1==0.) ? 0. : 
    283                    ((C2+C3>=1) ? LOCAL_INFINITY : eparams->S1/C1)) 
     287                   ((C2+C3>=1.) ? LOCAL_INFINITY : eparams->S1/C1)) 
    284288                  + ((C2<=0.) ? LOCAL_INFINITY : eparams->S2/C2) 
    285289                  + ((C2>=1.) ? LOCAL_INFINITY : eparams->S4/(1.-C2))) 
    … …  
    291295  for (i=0; i<eparams->n_poly; i++) 
    292296    ret += ((eparams->poly) [i].coeff + (eparams->poly) [i].coeff * T) * 
    293       pow (C1, (eparams->poly) [i].n1-2) * 
    294       pow (C2, (eparams->poly) [i].n2-2) * pow (C3, (eparams->poly) [i].n3) * 
     297      POWER2(C1, (eparams->poly) [i].n1-2) * 
     298      POWER2(C2, (eparams->poly) [i].n2-2) * pow (C3, (eparams->poly) [i].n3) * 
    295299      ((eparams->poly) [i].n2 * ((eparams->poly) [i].n2 - 1) * C1 * C1 
    296300       - 2. * (eparams->poly) [i].n1 * (eparams->poly) [i].n2 * C1 * C2 
    … …  
    360364  for (i=0; i<eparams->n_poly; i++) 
    361365    ret += ((eparams->poly) [i].coeff + (eparams->poly) [i].coeff * T) * 
    362       pow (C1, (eparams->poly) [i].n1-2) * 
    363       pow (C2, (eparams->poly) [i].n2) * pow (C3, (eparams->poly) [i].n3-2) * 
     366      POWER2(C1, (eparams->poly) [i].n1-2) * 
     367      pow (C2, (eparams->poly) [i].n2) * POWER2(C3, (eparams->poly) [i].n3-2) * 
    364368      ((eparams->poly) [i].n3 * ((eparams->poly) [i].n3 - 1) * C1 * C1 
    365369       - 2. * (eparams->poly) [i].n1 * (eparams->poly) [i].n3 * C1 * C3 
    … …  
    429433  for (i=0; i<eparams->n_poly; i++) 
    430434    ret += ((eparams->poly) [i].coeff + (eparams->poly) [i].coeff * T) * 
    431       pow (C1, (eparams->poly) [i].n1-2) * 
     435      POWER2(C1, (eparams->poly) [i].n1-2) * 
    432436      pow (C2, (eparams->poly) [i].n2-1) * pow (C3, (eparams->poly) [i].n3-1) * 
    433437      ((eparams->poly) [i].n2 * (eparams->poly) [i].n3 * C1 * C1