From 3662c314eef94f0e4c6f3a74684c8e3fb46650a9 Mon Sep 17 00:00:00 2001
From: Dave Goodwin <dgg@caltech.edu>
Date: Tue, 11 Jul 2006 15:25:57 +0000
Subject: [PATCH] fixed minor bug that caused problems in some cases when an
 element was excluded.

---
 Cantera/src/MultiPhaseEquil.cpp | 21 ++++++++++++++++++++-
 1 file changed, 20 insertions(+), 1 deletion(-)

diff --git a/Cantera/src/MultiPhaseEquil.cpp b/Cantera/src/MultiPhaseEquil.cpp
index 7927fed56..ad157d5b2 100644
--- a/Cantera/src/MultiPhaseEquil.cpp
+++ b/Cantera/src/MultiPhaseEquil.cpp
@@ -536,12 +536,17 @@ namespace Cantera {
 
     void MultiPhaseEquil::step(doublereal omega, vector_fp& deltaN) {
         index_t k, ik;
+        beginLogGroup("MultiPhaseEquil::step");
         if (omega < 0.0) 
             throw CanteraError("step","negative omega");
 
         for (ik = 0; ik < m_nel; ik++) {
             k = m_order[ik];
             m_lastmoles[k] = m_moles[k];
+            addLogEntry("component "+m_mix->speciesName(m_species[k])+" moles",
+                m_moles[k]);
+            addLogEntry("component "+m_mix->speciesName(m_species[k])+" step",
+                omega*deltaN[k]);
             m_moles[k] += omega * deltaN[k];
         }
 
@@ -557,6 +562,7 @@ namespace Cantera {
             }
         }
         updateMixMoles();
+        endLogGroup("MultiPhaseEquil::step");
     }
 
 
@@ -724,11 +730,15 @@ namespace Cantera {
                         for (k = 0; k < m_nsp; k++) {
                             kc = m_species[k];
                             if (m_mix->speciesPhaseIndex(kc) == ip) {
-                                stoich = nu[kc];
+                                stoich = nu[k]; // nu[kc];
                                 psum += stoich * stoich;
                             }
                         }
                         sum -= psum / (fabs(m_mix->phaseMoles(ip)) + TINY);
+                        if (ISNAN(sum)) {
+                            cout << " sum is nan. " << endl;
+                            cout << psum << " " << m_mix->phaseMoles(ip) << endl;
+                        }
                     }
                 }
                 rfctr = term1 + csum + sum;
@@ -736,8 +746,17 @@ namespace Cantera {
                     fctr = 1.0;
                 else
                     fctr = 1.0/(term1 + csum + sum);
+                if (ISNAN(fctr)) {
+                    cout << "fctr is nan." << endl;
+                    cout << term1 << " " << csum << " " << sum << " " << TINY << endl;
+                }
             }
             dxi[j] = -fctr*dg_rt;
+            if (ISNAN(dxi[j])) {
+                cout << "nan detected. " << endl;
+                cout << fctr << " " << dg_rt << endl;
+            }
+                
             index_t m;
              for (m = 0; m < m_nel; m++) {
                 if (m_moles[m_order[m]] <= 0.0 && (m_N(m, j)*dxi[j] < 0.0))