content: # -*- coding: utf-8 -*-
  """
  Created on Mon Apr 24 14:13:47 2017
  
  @author: ishort
  """
  
  import math
  import Useful
  import ToolBox
  #import numpy
  #JB#
  #from matplotlib.pyplot import plot, title, show, scatter
  #storage for fits (not all may be used)
  uw = []
  uwa = []
  uwb = []
  uwStage = []
  uwbStage = []
  uwu = []
  uwl = []
  uua=[]
  uub=[]
  """
  #a function to create a cubic function fit extrapolation
  def cubicFit(x,y):
      coeffs = numpy.polyfit(x,y,3)
      #returns an array of coefficents for the cubic fit of the form
      #Ax^3 + Bx^2 + Cx + D as [A,B,C,D]
      return coeffs
  #this will work for any number of data points!
  def valueFromFit(fit,x):
      #return the value y for a given fit, at point x
      return (fit[0]*(x**3)+fit[1]*(x**2)+fit[2]*x+fit[3])
  #holds the five temperature at which we have partition function data
  """
  masterTemp = [130, 500, 3000, 8000, 10000]
  #JB#
  
  #def levelPops(lam0In, logNStage, chiL, log10UwStage, gwL, numDeps, temp): 
  def levelPops(lam0In, logNStage, chiL, logUw, gwL, numDeps, temp):    
      """ Returns depth distribution of occupation numbers in lower level of b-b transition,
  
  // Input parameters:
  // lam0 - line centre wavelength in nm
  // logNStage - log_e density of absorbers in relevent ion stage (cm^-3)
  // logFlu - log_10 oscillator strength (unitless)
  // chiL - energy of lower atomic E-level of b-b transition in eV
  // Also needs atsmopheric structure information:
  // numDeps
  // temp structure """
      
      c = Useful.c()
      logC = Useful.logC()
      k = Useful.k()
      logK = Useful.logK()
      logH = Useful.logH()
      logEe = Useful.logEe()
      logMe = Useful.logMe()
  
      ln10 = math.log(10.0)
      logE = math.log10(math.e); #// for debug output
      log2pi = math.log(2.0 * math.pi)
      log2 = math.log(2.0)
  
      #//double logNl = logNlIn * ln10;  // Convert to base e
  
  
      #// Parition functions passed in are 2-element vectore with remperature-dependent base 10 log Us
      #// Convert to natural logs:
      #double thisLogUw, Ttheta;
      thisLogUw = 0.0 # //default initialization
      #logUw = [ 0.0 for i in range(5) ]
      logE10 = math.log(10.0)
      #print("log10UwStage ", log10UwStage)
      #for kk in range(len(logUw)):
      #    logUw[kk] = logE10*log10UwStage[kk] #// lburns new loop
          
      logGwL = math.log(gwL)
      
      
      #//System.out.println("chiL before: " + chiL);
      #// If we need to subtract chiI from chiL, do so *before* converting to tiny numbers in ergs!
      #////For testing with Ca II lines using gS3 internal line list only:
      #//boolean ionized = true;
      #//if (ionized) {
      #//    //System.out.println("ionized, doing chiL - chiI: " + ionized);
      #//    //         chiL = chiL - chiI;
      #//             chiL = chiL - 6.113;
      #//          }
      #//   //
  
      #//Log of line-center wavelength in cm
      logLam0 = math.log(lam0In) #// * 1.0e-7);
  
      #// energy of b-b transition
      logTransE = logH + logC - logLam0 #//ergs
      
      if (chiL <= 0.0):
          chiL = 1.0e-49
      logChiL = math.log(chiL) + Useful.logEv() #// Convert lower E-level from eV to ergs
      
      logBoltzFacL = logChiL - Useful.logK() #// Pre-factor for exponent of excitation Boltzmann factor
      boltzFacL = math.exp(logBoltzFacL)
  
      boltzFacGround = 0.0 / k #//I know - its zero, but let's do it this way anyway'
  
  
      #// return a 1D numDeps array of logarithmic number densities
      #// level population of lower level of bb transition (could be in either stage I or II!) 
          
      logNums = [ 0.0 for i in range(numDeps)]
  
      #double num, logNum, expFac;
      #JB#
      #print("thisLogUw:",numpy.shape(logUw))
      logUwFit = ToolBox.cubicFit(masterTemp,logUw)#u(T) fit
      uw.append(logUwFit)
      #JB#
  
      for id in range(numDeps):
  
      #//Determine temperature dependenet partition functions Uw:
          
          #Ttheta = 5040.0 / temp[0][id]
  #//NEW Determine temperature dependent partition functions Uw: lburns
          thisTemp = temp[0][id]
          
          """
          if (Ttheta >= 1.0):
              thisLogUw = logUw[0]
          
          if (Ttheta <= 0.5):
              thisLogUw = logUw[1]
          
          if (Ttheta > 0.5 and Ttheta < 1.0):
              thisLogUw = ( logUw[1] * (Ttheta - 0.5)/(1.0 - 0.5) ) \
                        + ( logUw[0] * (1.0 - Ttheta)/(1.0 - 0.5) )
          """             
      #JB#
          thisLogUw = ToolBox.valueFromFit(logUwFit,thisTemp)#u(T) value extrapolated
      #JB#
  
          if (thisTemp >= 10000.0):
              thisLogUw = logUw[4]
          
          if (thisTemp <= 130.0):
              thisLogUw = logUw[0]
              
          """
          
          if (thisTemp > 130 and thisTemp <= 500):
              thisLogUw = logUw[1] * (thisTemp - 130)/(500 - 130) \
                        + logUw[0] * (500 - thisTemp)/(500 - 130)
          
          if (thisTemp > 500 and thisTemp <= 3000):
              thisLogUw = logUw[2] * (thisTemp - 500)/(3000 - 500) \
                        + logUw[1] * (3000 - thisTemp)/(3000 - 500)
          
          if (thisTemp > 3000 and thisTemp <= 8000):
              thisLogUw = logUw[3] * (thisTemp - 3000)/(8000 - 3000) \
                        + logUw[2] * (8000 - thisTemp)/(8000 - 3000)
          
          if (thisTemp > 8000 and thisTemp < 10000):
              thisLogUw = logUw[4] * (thisTemp - 8000)/(10000 - 8000) \
                        + logUw[3] * (10000 - thisTemp)/(10000 - 8000)
          
          """
          #print("logUw ", logUw, " thisLogUw ", thisLogUw)
                             
  
          #//System.out.println("LevPops: ionized branch taken, ionized =  " + ionized);
          #// Take stat weight of ground state as partition function:
          logNums[id] = logNStage[id] - boltzFacL / temp[0][id] + logGwL - thisLogUw #// lower level of b-b transition
          #print("LevelPopsServer.stagePops id ", id, " logNStage[id] ", logNStage[id], " boltzFacL ", boltzFacL, " temp[0][id] ", temp[0][id], " logGwL ", logGwL, " thisLogUw ", thisLogUw, " logNums[id] ", logNums[id]);
  
          #// System.out.println("LevelPops: id, logNums[0][id], logNums[1][id], logNums[2][id], logNums[3][id]: " + id + " "
          #//          + Math.exp(logNums[0][id]) + " "
          #//         + Math.exp(logNums[1][id]) + " "
          #//          + Math.exp(logNums[2][id]) + " "
          #//        + Math.exp(logNums[3][id]));
          #//System.out.println("LevelPops: id, logNums[0][id], logNums[1][id], logNums[2][id], logNums[3][id], logNums[4][id]: " + id + " "
          #//        + logE * (logNums[0][id]) + " "
          #//        + logE * (logNums[1][id]) + " "
          #//        + logE * (logNums[2][id]) + " "
          # //        + logE * (logNums[3][id]) + " "
          #//        + logE * (logNums[4][id]) );
          #//System.out.println("LevelPops: id, logIonFracI, logIonFracII: " + id + " " + logE*logIonFracI + " " + logE*logIonFracII
          #//        + "logNum, logNumI, logNums[0][id], logNums[1][id] "
          #//        + logE*logNum + " " + logE*logNumI + " " + logE*logNums[0][id] + " " + logE*logNums[1][id]);
          #//System.out.println("LevelPops: id, logIonFracI: " + id + " " + logE*logIonFracI
          #//        + "logNums[0][id], boltzFacL/temp[0][id], logNums[2][id]: " 
          #//        + logNums[0][id] + " " + boltzFacL/temp[0][id] + " " + logNums[2][id]);
      #//id loop
      #stop
      #print (uw)
      return logNums
  
  #//This version - ionization equilibrium *WITHOUT* molecules - logNum is TOTAL element population
  #def stagePops2(logNum, Ne, chiIArr, log10UwAArr,  \
  #               numMols, logNumB, dissEArr, log10UwBArr, logQwABArr, logMuABArr, \
  #               numDeps, temp):
  def stagePops(logNum, Ne, chiIArr, logUw,  \
                 numDeps, temp):
      #line 1: //species A data - ionization equilibrium of A
      #line 2: //data for set of species "B" - molecular equlibrium for set {AB}
      """Ionization equilibrium routine WITHOUT molecule formation:
      // Returns depth distribution of ionization stage populations 
  
      // Input parameters:
      // logNum - array with depth-dependent total element number densities (cm^-3) 
      // chiI1 - ground state ionization energy of neutral stage 
      // chiI2 - ground state ionization energy of singly ionized stage 
      // Also needs atsmopheric structure information:
      // numDeps
      // temp structure 
      // rho structure
      // Atomic element A is the one whose ionization fractions are being computed
      //  
      """
  
      ln10 = math.log(10.0)
      logE = math.log10(math.e) #// for debug output
      log2pi = math.log(2.0 * math.pi)
      log2 = math.log(2.0)
  
      numStages = len(chiIArr)  #// + 1; //need one more stage above the highest stage to be populated
  
  #//    var numMols = dissEArr.length;
  
  
  #// Parition functions passed in are 2-element vectore with remperature-dependent base 10 log Us
  #// Convert to natural logs:
          #double Ttheta, thisTemp;
  #//Default initializations:
  #//We need one more stage in size of saha factor than number of stages we're actualy populating
      thisLogUw = [ 0.0 for i in range(numStages+1) ]
      
      for i in range(numStages+1):
          thisLogUw[i] = 0.0
          
  
      logE10 = math.log(10.0)
  
  #//atomic ionization stage Boltzmann factors:
          #double logChiI, logBoltzFacI;
      boltzFacI = [ 0.0 for i in range(numStages) ]
      #print("numStages ", numStages, " Useful.logEv ", Useful.logEv())
      for i in range(numStages):
          #print("i ", i, " chiIArr ", chiIArr[i])
          logChiI = math.log(chiIArr[i]) + Useful.logEv() 
          logBoltzFacI = logChiI  - Useful.logK()
          boltzFacI[i] = math.exp(logBoltzFacI)
          
  
      logSahaFac = log2 + (3.0 / 2.0) * (log2pi + Useful.logMe() + Useful.logK() - 2.0 * Useful.logH())
  
      #// return a 2D 5 x numDeps array of logarithmic number densities
      #// Row 0: neutral stage ground state population
      #// Row 1: singly ionized stage ground state population
      #// Row 2: doubly ionized stage ground state population        
      #// Row 3: triply ionized stage ground state population        
      #// Row 4: quadruply ionized stage ground state population        
          #double[][] logNums = new double[numStages][numDeps];
      logNums = [ [ 0.0 for i in range(numDeps)] for j in range(numStages) ]
  
      #//We need one more stage in size of saha factor than number of stages we're actualy populating
      #//   for index accounting pirposes
      #//   For atomic ionization stages:
  
      logSaha = [ [ 0.0 for i in range(numStages+1)] for j in range(numStages+1) ]
      saha = [ [ 0.0 for i in range(numStages+1)] for j in range(numStages+1) ]
  #//
          
      logIonFrac = [ 0.0 for i in range(numStages) ]
          #double expFac, logNe;
  
  #// Now - molecular variables:
  
  
      thisLogUwA = 0.0 #// element A 
      #thisLogQwAB = math.log(300.0) 
      
  #//For clarity: neutral stage of atom whose ionization equilibrium is being computed is element A
  #// for molecule formation:
      logUwA = [ 0.0 for i in range(5) ]
  
      
      #JB#
      uua=[]
      #uub=[]
      #qwab=[]
      for iStg in range(numStages):
          currentUwArr=list(logUw[iStg])#u(T) determined values
          UwFit = ToolBox.cubicFit(masterTemp,currentUwArr)#u(T) fit
          uua.append(UwFit)
          #print(logUw[iStg])
  
      
      for id in range(numDeps):
  
          #//// reduce or enhance number density by over-all Rosseland opcity scale parameter
          #//
          #//Row 1 of Ne is log_e Ne in cm^-3
          logNe = Ne[1][id]
  
          #//Determine temperature dependent partition functions Uw:
          thisTemp = temp[0][id]
          #Ttheta = 5040.0 / thisTemp
  
      #JB#
      #use temps and partition values to create a function
      #then use said function to extrapolate values for all points
          thisLogUw[numStages] = 0.0
          for iStg in range(numStages):
              thisLogUw[iStg] = ToolBox.valueFromFit(uua[iStg],thisTemp)#u(T) value extrapolated
  
      #JB#
           
  #// NEW Determine temperature dependent partition functions Uw: lburns
          if (thisTemp <= 130.0):
              for iStg in range(numStages):
                  thisLogUw[iStg] = logUw[iStg][0]            
              #for iMol in range(numMols):
              #    thisLogUwB[iMol] = logUwB[iMol][0]
  
          if (thisTemp >= 10000.0):
              for iStg in range(numStages):
                  thisLogUw[iStg] = logUw[iStg][4]            
              #for iMol in range(numMols):
              #    thisLogUwB[iMol] = logUwB[iMol][4]
              
  
  #//For clarity: neutral stage of atom whose ionization equilibrium is being computed is element A
  #// for molecule formation:
          thisLogUwA = thisLogUw[0];
  
          #//Ionization stage Saha factors: 
          for iStg in range(numStages):
              #print("iStg ", iStg)
              logSaha[iStg+1][iStg] = logSahaFac - logNe - (boltzFacI[iStg] /temp[0][id]) + (3.0 * temp[1][id] / 2.0) + thisLogUw[iStg+1] - thisLogUw[iStg]
              saha[iStg+1][iStg] = math.exp(logSaha[iStg+1][iStg])
              
  
  #//Compute log of denominator is ionization fraction, f_stage 
          denominator = 1.0 #//default initialization - leading term is always unity 
  #//ion stage contributions:
          for jStg in range(1, numStages+1):
              addend = 1.0 #//default initialization for product series
              for iStg in range(jStg):
                  #//console.log("jStg " + jStg + " saha[][] indices " + (iStg+1) + " " + iStg); 
                  addend = addend * saha[iStg+1][iStg] 
                 
              denominator = denominator + addend 
              
  
  #// 
          logDenominator = math.log(denominator) 
  
          logIonFrac[0] = -1.0 * logDenominator     #// log ionization fraction in stage I
  
          for jStg in range(1, numStages):
              addend = 0.0 #//default initialization for product series
              for iStg in range(jStg):
                  #//console.log("jStg " + jStg + " saha[][] indices " + (iStg+1) + " " + iStg); 
                  addend = addend + logSaha[iStg+1][iStg]
                 
              logIonFrac[jStg] = addend - logDenominator
  
          for iStg in range(numStages):
              logNums[iStg][id] = logNum[id] + logIonFrac[iStg]
                
      #//id loop
  
      return logNums;
      #//end method stagePops
      
  
  #end method levelPops
  
  #def stagePops2(logNum, Ne, chiIArr, log10UwAArr,  \
  #               numMols, logNumB, dissEArr, log10UwBArr, logQwABArr, logMuABArr, \
  #               numDeps, temp):
  def stagePops2(logNum, Ne, chiIArr, logUw,  \
                 numMols, logNumB, dissEArr, logUwB, logQwABArr, logMuABArr, \
                 numDeps, temp):
      #line 1: //species A data - ionization equilibrium of A
      #line 2: //data for set of species "B" - molecular equlibrium for set {AB}
      """Ionization equilibrium routine that accounts for molecule formation:
      // Returns depth distribution of ionization stage populations 
  
      // Input parameters:
      // logNum - array with depth-dependent total element number densities (cm^-3) 
      // chiI1 - ground state ionization energy of neutral stage 
      // chiI2 - ground state ionization energy of singly ionized stage 
      // Also needs atsmopheric structure information:
      // numDeps
      // temp structure 
      // rho structure
      // Atomic element A is the one whose ionization fractions are being computed
      //  Element B refers to array of other species with which A forms molecules AB """
  
      ln10 = math.log(10.0)
      logE = math.log10(math.e) #// for debug output
      log2pi = math.log(2.0 * math.pi)
      log2 = math.log(2.0)
  
      numStages = len(chiIArr)  #// + 1; //need one more stage above the highest stage to be populated
  
  #//    var numMols = dissEArr.length;
  
  
  #// Parition functions passed in are 2-element vectore with remperature-dependent base 10 log Us
  #// Convert to natural logs:
          #double Ttheta, thisTemp;
  #//Default initializations:
  #//We need one more stage in size of saha factor than number of stages we're actualy populating
      thisLogUw = [ 0.0 for i in range(numStages+1) ]
      
      for i in range(numStages+1):
          thisLogUw[i] = 0.0
          
  
      logE10 = math.log(10.0)
  
  #//atomic ionization stage Boltzmann factors:
          #double logChiI, logBoltzFacI;
      boltzFacI = [ 0.0 for i in range(numStages) ]
      #print("numStages ", numStages, " Useful.logEv ", Useful.logEv())
      for i in range(numStages):
          #print("i ", i, " chiIArr ", chiIArr[i])
          logChiI = math.log(chiIArr[i]) + Useful.logEv() 
          logBoltzFacI = logChiI  - Useful.logK()
          boltzFacI[i] = math.exp(logBoltzFacI)
          
  
      logSahaFac = log2 + (3.0 / 2.0) * (log2pi + Useful.logMe() + Useful.logK() - 2.0 * Useful.logH())
  
      #// return a 2D 5 x numDeps array of logarithmic number densities
      #// Row 0: neutral stage ground state population
      #// Row 1: singly ionized stage ground state population
      #// Row 2: doubly ionized stage ground state population        
      #// Row 3: triply ionized stage ground state population        
      #// Row 4: quadruply ionized stage ground state population        
          #double[][] logNums = new double[numStages][numDeps];
      logNums = [ [ 0.0 for i in range(numDeps)] for j in range(numStages) ]
  
      #//We need one more stage in size of saha factor than number of stages we're actualy populating
      #//   for index accounting pirposes
      #//   For atomic ionization stages:
  
      logSaha = [ [ 0.0 for i in range(numStages+1)] for j in range(numStages+1) ]
      saha = [ [ 0.0 for i in range(numStages+1)] for j in range(numStages+1) ]
  #//
          
      logIonFrac = [ 0.0 for i in range(numStages) ]
          #double expFac, logNe;
  
  #// Now - molecular variables:
  
  #//Treat at least one molecule - if there are really no molecules for an atomic species, 
  #//there will be one phantom molecule in the denominator of the ionization fraction
  #//with an impossibly high dissociation energy
      ifMols = True
      if (numMols == 0):
          ifMols = False
          numMols = 1
  #//This should be inherited, but let's make sure: 
          dissEArr[0] = 19.0 #//eV
     
  
  #//Molecular partition functions - default initialization:
         #double[] thisLogUwB = new double[numMols];
      
      thisLogUwB = [ 0.0 for i in range(numMols) ]
      
      for iMol in range(numMols):
          thisLogUwB[iMol] = 0.0 #// variable for temp-dependent computed partn fn of array element B 
         
      thisLogUwA = 0.0 #// element A 
      thisLogQwAB = math.log(300.0) 
  
  #//For clarity: neutral stage of atom whose ionization equilibrium is being computed is element A
  #// for molecule formation:
      logUwA = [ 0.0 for i in range(5) ]
      if (numMols > 0):
          for kk in range(len(logUwA)):
              logUwA[kk] = logUw[0][kk]
           #// lburns
  
          
        #//}
  #//// Molecular partition functions:
  
  #//Molecular dissociation Boltzmann factors:
      boltzFacIAB = [ 0.0 for i in range(numMols) ]
      logMolSahaFac = [ 0.0 for i in range(numMols) ]
      #//if (numMols > 0){
      #double logDissE, logBoltzFacIAB;
      for iMol in range(numMols):
          logDissE = math.log(dissEArr[iMol]) + Useful.logEv() 
          logBoltzFacIAB = logDissE  - Useful.logK()
          boltzFacIAB[iMol] = math.exp(logBoltzFacIAB)
          logMolSahaFac[iMol] = (3.0 / 2.0) * (log2pi + logMuABArr[iMol] + Useful.logK() - 2.0 * Useful.logH())
          #//console.log("iMol " + iMol + " dissEArr[iMol] " + dissEArr[iMol] + " logDissE " + logE*logDissE + " logBoltzFacIAB " + logE*logBoltzFacIAB + " boltzFacIAB[iMol] " + boltzFacIAB[iMol] + " logMuABArr " + logE*logMuABArr[iMol] + " logMolSahaFac " + logE*logMolSahaFac[iMol]);
          
         #//}
  #//   For molecular species:
      logSahaMol = [ 0.0 for i in range(numMols) ]
      invSahaMol = [ 0.0 for i in range(numMols) ]
      
      
      #JB#
      uua=[]
      uub=[]
      qwab=[]
      for iStg in range(numStages):
          currentUwArr=list(logUw[iStg])#u(T) determined values
          UwFit = ToolBox.cubicFit(masterTemp,currentUwArr)#u(T) fit
          uua.append(UwFit)
          #print(logUw[iStg])
  
      for iMol in range(numMols):
          currentUwBArr=list(logUwB[iMol])#u(T) determined values
          UwBFit = ToolBox.cubicFit(masterTemp,currentUwBArr)#u(T) fit
          uub.append(UwBFit)
          
      
      for id in range(numDeps):
  
          #//// reduce or enhance number density by over-all Rosseland opcity scale parameter
          #//
          #//Row 1 of Ne is log_e Ne in cm^-3
          logNe = Ne[1][id]
  
          #//Determine temperature dependent partition functions Uw:
          thisTemp = temp[0][id]
          #Ttheta = 5040.0 / thisTemp
  
      #JB#
      #use temps and partition values to create a function
      #then use said function to extrapolate values for all points
          thisLogUw[numStages] = 0.0
          for iStg in range(numStages):
              thisLogUw[iStg] = ToolBox.valueFromFit(uua[iStg],thisTemp)#u(T) value extrapolated
          for iMol in range(numMols):
              thisLogUwB[iMol] = ToolBox.valueFromFit(uub[iMol],thisTemp)#u(T) value extrapolated
      #JB#
           
  #// NEW Determine temperature dependent partition functions Uw: lburns
          if (thisTemp <= 130.0):
              for iStg in range(numStages):
                  thisLogUw[iStg] = logUw[iStg][0]            
              for iMol in range(numMols):
                  thisLogUwB[iMol] = logUwB[iMol][0]
  
          if (thisTemp >= 10000.0):
              for iStg in range(numStages):
                  thisLogUw[iStg] = logUw[iStg][4]            
              for iMol in range(numMols):
                  thisLogUwB[iMol] = logUwB[iMol][4]
              
  
          for iMol in range(numMols):
              if (thisTemp < 3000.0):
                  thisLogQwAB = ( logQwABArr[iMol][1] * (3000.0 - thisTemp)/(3000.0 - 500.0) ) \
                              + ( logQwABArr[iMol][2] * (thisTemp - 500.0)/(3000.0 - 500.0) )
           
              if ( (thisTemp >= 3000.0) and (thisTemp <= 8000.0) ):
                  thisLogQwAB = ( logQwABArr[iMol][2] * (8000.0 - thisTemp)/(8000.0 - 3000.0) ) \
                              + ( logQwABArr[iMol][3] * (thisTemp - 3000.0)/(8000.0 - 3000.0) )
           
              if ( thisTemp > 8000.0 ):
                  thisLogQwAB = ( logQwABArr[iMol][3] * (10000.0 - thisTemp)/(10000.0 - 8000.0) ) \
                              + ( logQwABArr[iMol][4] * (thisTemp - 8000.0)/(10000.0 - 8000.0) )
             
          #// iMol loop 
              
  #//For clarity: neutral stage of atom whose ionization equilibrium is being computed is element A
  #// for molecule formation:
          thisLogUwA = thisLogUw[0];
  
          #//Ionization stage Saha factors: 
          for iStg in range(numStages):
              #print("iStg ", iStg)
              logSaha[iStg+1][iStg] = logSahaFac - logNe - (boltzFacI[iStg] /temp[0][id]) + (3.0 * temp[1][id] / 2.0) + thisLogUw[iStg+1] - thisLogUw[iStg]
              saha[iStg+1][iStg] = math.exp(logSaha[iStg+1][iStg])
              
              
  #//Molecular Saha factors:
          for iMol in range(numMols):
              logSahaMol[iMol] = logMolSahaFac[iMol] - logNumB[iMol][id] - (boltzFacIAB[iMol] / temp[0][id]) + (3.0 * temp[1][id] / 2.0) + thisLogUwB[iMol] + thisLogUwA - thisLogQwAB
  #//For denominator of ionization fraction, we need *inverse* molecular Saha factors (N_AB/NI):
              logSahaMol[iMol] = -1.0 * logSahaMol[iMol]
              invSahaMol[iMol] = math.exp(logSahaMol[iMol])
  
  #//Compute log of denominator is ionization fraction, f_stage 
          denominator = 1.0 #//default initialization - leading term is always unity 
  #//ion stage contributions:
          for jStg in range(1, numStages+1):
              addend = 1.0 #//default initialization for product series
              for iStg in range(jStg):
                  #//console.log("jStg " + jStg + " saha[][] indices " + (iStg+1) + " " + iStg); 
                  addend = addend * saha[iStg+1][iStg] 
                 
              denominator = denominator + addend 
              
  #//molecular contribution
          if (ifMols == True):
              for iMol in range(numMols):
                  denominator = denominator + invSahaMol[iMol]
                         
  #// 
          logDenominator = math.log(denominator) 
  
          logIonFrac[0] = -1.0 * logDenominator     #// log ionization fraction in stage I
  
          for jStg in range(1, numStages):
              addend = 0.0 #//default initialization for product series
              for iStg in range(jStg):
                  #//console.log("jStg " + jStg + " saha[][] indices " + (iStg+1) + " " + iStg); 
                  addend = addend + logSaha[iStg+1][iStg]
                 
              logIonFrac[jStg] = addend - logDenominator
  
          for iStg in range(numStages):
              logNums[iStg][id] = logNum[id] + logIonFrac[iStg]
                
      #//id loop
  
      return logNums;
      #//end method stagePops
      
  def stagePops3(logNum, Ne, chiIArr, logUw, numDeps, temp):
      #Version for ChromaStarPyGas: logNum is now *neutral stage* population from Phil
      # Bennett's GAS package
      
      #line 1: //species A data - ionization equilibrium of A
      #line 2: //data for set of species "B" - molecular equlibrium for set {AB}
      """Ionization equilibrium routine that accounts for molecule formation:
      // Returns depth distribution of ionization stage populations 
  
      // Input parameters:
      // logNum - array with depth-dependent neutral stage  number densities (cm^-3) 
      // chiI1 - ground state ionization energy of neutral stage 
      // chiI2 - ground state ionization energy of singly ionized stage 
      // Also needs atsmopheric structure information:
      // numDeps
      // temp structure 
      // rho structure
      // Atomic element A is the one whose ionization fractions are being computed
      //  Element B refers to array of other species with which A forms molecules AB """
  
      ln10 = math.log(10.0)
      logE = math.log10(math.e) #// for debug output
      log2pi = math.log(2.0 * math.pi)
      log2 = math.log(2.0)
  
      numStages = len(chiIArr)  #// + 1; //need one more stage above the highest stage to be populated
  
  #//    var numMols = dissEArr.length;
  
  
  #// Parition functions passed in are 2-element vectore with remperature-dependent base 10 log Us
  #// Convert to natural logs:
          #double Ttheta, thisTemp;
  #//Default initializations:
  #//We need one more stage in size of saha factor than number of stages we're actualy populating
      thisLogUw = [ 0.0 for i in range(numStages+1) ]
      
      for i in range(numStages+1):
          thisLogUw[i] = 0.0
          
  
      logE10 = math.log(10.0)
  
  #//atomic ionization stage Boltzmann factors:
          #double logChiI, logBoltzFacI;
      boltzFacI = [ 0.0 for i in range(numStages) ]
      #print("numStages ", numStages, " Useful.logEv ", Useful.logEv())
      for i in range(numStages):
          #print("i ", i, " chiIArr ", chiIArr[i])
          logChiI = math.log(chiIArr[i]) + Useful.logEv() 
          logBoltzFacI = logChiI  - Useful.logK()
          boltzFacI[i] = math.exp(logBoltzFacI)
          
  
      logSahaFac = log2 + (3.0 / 2.0) * (log2pi + Useful.logMe() + Useful.logK() - 2.0 * Useful.logH())
  
      #// return a 2D 5 x numDeps array of logarithmic number densities
      #// Row 0: neutral stage ground state population
      #// Row 1: singly ionized stage ground state population
      #// Row 2: doubly ionized stage ground state population        
      #// Row 3: triply ionized stage ground state population        
      #// Row 4: quadruply ionized stage ground state population        
          #double[][] logNums = new double[numStages][numDeps];
      logNums = [ [ 0.0 for i in range(numDeps)] for j in range(numStages) ]
  
      #//We need one more stage in size of saha factor than number of stages we're actualy populating
      #//   for index accounting pirposes
      #//   For atomic ionization stages:
  
      #logSaha = [ [ 0.0 for i in range(numStages+1)] for j in range(numStages+1) ]
      #saha = [ [ 0.0 for i in range(numStages+1)] for j in range(numStages+1) ]
  #//
          
      #logIonFrac = [ 0.0 for i in range(numStages) ]
          #double expFac, logNe;
  
      
      #JB#
      uua=[]
      uub=[]
      qwab=[]
      for iStg in range(numStages):
          currentUwArr=list(logUw[iStg])#u(T) determined values
          UwFit = ToolBox.cubicFit(masterTemp,currentUwArr)#u(T) fit
          uua.append(UwFit)
          #print(logUw[iStg])
              
      for id in range(numDeps):
  
          #//// reduce or enhance number density by over-all Rosseland opcity scale parameter
          #//
          #//Row 1 of Ne is log_e Ne in cm^-3
          logNe = Ne[1][id]
  
          #//Determine temperature dependent partition functions Uw:
          thisTemp = temp[0][id]
          #Ttheta = 5040.0 / thisTemp
  
      #JB#
      #use temps and partition values to create a function
      #then use said function to extrapolate values for all points
          thisLogUw[numStages] = 0.0
          for iStg in range(numStages):
              thisLogUw[iStg] = ToolBox.valueFromFit(uua[iStg],thisTemp)#u(T) value extrapolated
      #JB#
           
  #// NEW Determine temperature dependent partition functions Uw: lburns
          if (thisTemp <= 130.0):
              for iStg in range(numStages):
                  thisLogUw[iStg] = logUw[iStg][0]            
  
          if (thisTemp >= 10000.0):
              for iStg in range(numStages):
                  thisLogUw[iStg] = logUw[iStg][4]
             
              
  #//For clarity: neutral stage of atom whose ionization equilibrium is being computed is element A
  #// for molecule formation:
          #thisLogUwA = thisLogUw[0];
  
          #//Ionization stage Saha factors:
          logNums[0][id] = logNum[id] 
          for iStg in range(1, numStages):
              #print("iStg ", iStg)
              thisLogSaha = logSahaFac - logNe - (boltzFacI[iStg-1] /temp[0][id]) + (3.0 * temp[1][id] / 2.0) + thisLogUw[iStg] - thisLogUw[iStg-1]
              #saha[iStg+1][iStg] = math.exp(logSaha[iStg+1][iStg])
              logNums[iStg][id] = logNums[iStg-1][id] + thisLogSaha
                
      #//id loop
  
      return logNums;
      #//end method stagePops
      
  #def sahaRHS(chiI, log10UwUArr, log10UwLArr, temp):
  def sahaRHS(chiI, logUwU, logUwL, temp):    
  
      """RHS of partial pressure formulation of Saha equation in standard form (N_U*P_e/N_L on LHS)
   // Returns depth distribution of LHS: Phi(T) === N_U*P_e/N_L (David Gray notation)
  
  // Input parameters:
  // chiI - ground state ionization energy of lower stage 
  // log10UwUArr, log10UwLArr - array of temperature-dependent partition function for upper and lower ionization stage
  // Also needs atsmopheric structure information:
  // numDeps
  // temp structure 
  //
  // Atomic element "A" is the one whose ionization fractions are being computed
  //  Element "B" refers to array of other species with which A forms molecules "AB" """
      
      ln10 = math.log(10.0)
      logE = math.log10(math.e) #// for debug output
      log2pi = math.log(2.0 * math.pi)
      log2 = math.log(2.0)
  
  #//    var numMols = dissEArr.length;
  
  #// Parition functions passed in are 2-element vectore with remperature-dependent base 10 log Us
  #// Convert to natural logs:
      #double Ttheta, thisTemp;
  #//Default initializations:
  #//We need one more stage in size of saha factor than number of stages we're actualy populating
      thisLogUwU = 0.0
      thisLogUwL = 0.0
  
      logE10 = math.log(10.0)
  #//We need one more stage in size of saha factor than number of stages we're actualy populating
      #logUwU = [0.0 for i in range(5)]
      #logUwL = [0.0 for i in range(5)]   
      for kk in range(len(logUwL)):
          logUwU[kk] = logUwL[kk]
      #   logUwL[kk] = logE10*log10UwLArr[kk]
           
  
      
  
      #//System.out.println("chiL before: " + chiL);
      #// If we need to subtract chiI from chiL, do so *before* converting to tiny numbers in ergs!
  
  #//atomic ionization stage Boltzmann factors:
      #double logChiI, logBoltzFacI;
      #double boltzFacI;
      logChiI = math.log(chiI) + Useful.logEv() 
      logBoltzFacI = logChiI  - Useful.logK()
      boltzFacI = math.exp(logBoltzFacI)
  
  #//Extra factor of k to get k^5/2 in the P_e formulation of Saha Eq.
      logSahaFac = log2 + (3.0 / 2.0) * (log2pi + Useful.logMe() + Useful.logK() - 2.0 * Useful.logH()) + Useful.logK()
  
      #//double[] logLHS = new double[numDeps];
      #double logLHS;
  
  #//   For atomic ionization stages:
      #double logSaha, saha, expFac;
  
  #//  for (int id = 0; id < numDeps; id++) {
  
  #//
  #//Determine temperature dependent partition functions Uw:
      thisTemp = temp[0]
      #Ttheta = 5040.0 / thisTemp
  
      """         
      if (Ttheta >= 1.0):
          thisLogUwU = logUwU[0]
          thisLogUwL = logUwL[0]
         
      if (Ttheta <= 0.5):
          thisLogUwU = logUwU[1]
          thisLogUwL = logUwL[1]
         
      if (Ttheta > 0.5 and Ttheta < 1.0):
          thisLogUwU = ( logUwU[1] * (Ttheta - 0.5)/(1.0 - 0.5) )
          + ( logUwU[0] * (1.0 - Ttheta)/(1.0 - 0.5) )
          thisLogUwL = ( logUwL[1] * (Ttheta - 0.5)/(1.0 - 0.5) )
          + ( logUwL[0] * (1.0 - Ttheta)/(1.0 - 0.5) )
      """
      #JB#
      currentUwUArr=list(logUwU)#u(T) determined values
      UwUFit = ToolBox.cubicFit(masterTemp,currentUwUArr)#u(T) fit
      thisLogUwU = ToolBox.valueFromFit(UwUFit,thisTemp)#u(T) value extrapolated
  
      currentUwLArr=list(logUwL)#u(T) determined values
      UwLFit = ToolBox.cubicFit(masterTemp,currentUwLArr)#u(T) fit
      thisLogUwL = ToolBox.valueFromFit(UwLFit,thisTemp)#u(T) value extrapolated
      #JB#
  
      #will need to do this one in Main as it goes through its own loop of temp
      #if thisTemp == superTemp[0][len(superTemp[0])]:
      #    uwu.append(UwUFit)
      #    uwl.append(UwLFit)
      #    
      #JB#
      
      if (thisTemp <= 130.0):
          thisLogUwU = logUwU[0]
          thisLogUwL = logUwL[0]
          
      if (thisTemp >= 10000.0):
          thisLogUwU = logUwU[4]
          thisLogUwL = logUwL[4]
          
      """   
      if (thisTemp > 130 and thisTemp <= 500):
          thisLogUwU = logUwU[1] * (thisTemp - 130)/(500 - 130) \
              + logUwU[0] * (500 - thisTemp)/(500 - 130)
          thisLogUwL = logUwL[1] * (thisTemp - 130)/(500 - 130) \
              + logUwL[0] * (500 - thisTemp)/(500 - 130)
          
      if (thisTemp > 500 and thisTemp <= 3000):
          thisLogUwU = logUwU[2] * (thisTemp - 500)/(3000 - 500) \
              + logUwU[1] * (3000 - thisTemp)/(3000 - 500)
          thisLogUwL = logUwL[2] * (thisTemp - 500)/(3000 - 500) \
              + logUwL[1] * (3000 - thisTemp)/(3000 - 500)
          
      if (thisTemp > 3000 and thisTemp <= 8000):
          thisLogUwU = logUwU[3] * (thisTemp - 3000)/(8000 - 3000) \
              + logUwU[2] * (8000 - thisTemp)/(8000 - 3000)
          thisLogUwL = logUwL[3] * (thisTemp - 3000)/(8000 - 3000) \
              + logUwL[2] * (8000 - thisTemp)/(8000 - 3000)
          
      if (thisTemp > 8000 and thisTemp < 10000): 
          thisLogUwU = logUwU[4] * (thisTemp - 8000)/(10000 - 8000) \
              + logUwU[3] * (10000 - thisTemp)/(10000 - 8000)
          thisLogUwL = logUwL[4] * (thisTemp - 8000)/(10000 - 8000) \
              + logUwL[3] * (10000 - thisTemp)/(10000 - 8000)
          
      if (thisTemp >= 10000):
          thisLogUwU = logUwU[4]
          thisLogUwL = logUwL[4]
      """
   
  
                    
  #//Ionization stage Saha factors: 
               
  #//Need T_kin^5/2 in the P_e formulation of Saha Eq.
      logSaha = logSahaFac - (boltzFacI /temp[0]) + (5.0 * temp[1] / 2.0) + thisLogUwU - thisLogUwL
      #// saha = Math.exp(logSaha);
  
      #//logLHS[id] = logSaha;
      logLHS = logSaha;
  #//    } //id loop
  
      return logLHS;
      #JB
      #return [logLHS,[[UwUFit,thisLogUwU],[UwLFit,thisLogUwL]]]
  #//
  #    } //end method sahaRHS    
      
  #def molPops(nmrtrLogNumB, nmrtrDissE, log10UwA, nmrtrLog10UwB, nmrtrLogQwAB, nmrtrLogMuAB, \
  #            numMolsB, logNumB, dissEArr, log10UwBArr, logQwABArr, logMuABArr,   \
  #            logGroundRatio, numDeps, temp):
  def molPops(nmrtrLogNumB, nmrtrDissE, logUwA, nmrtrLogUwB, nmrtrLogQwAB, nmrtrLogMuAB, \
              numMolsB, logNumB, dissEArr, logUwB, logQwABArr, logMuABArr,   \
              logGroundRatio, numDeps, temp):
      # line 1: //species A data - ionization equilibrium of A
      # //data for set of species "B" - molecular equlibrium for set {AB}
      
      """Diatomic molecular equilibrium routine that accounts for molecule formation:
   // Returns depth distribution of molecular population 
  
  // Input parameters:
  // logNum - array with depth-dependent total element number densities (cm^-3) 
  // chiI1 - ground state ionization energy of neutral stage 
  // chiI2 - ground state ionization energy of singly ionized stage 
  // Also needs atsmopheric structure information:
  // numDeps
  // temp structure 
  // rho structure
  //
  // Atomic element "A" is the one kept on the LHS of the master fraction, whose ionization fractions are included 
  //   in the denominator of the master fraction
  //  Element "B" refers to array of other sintpecies with which A forms molecules "AB" """
  
  
      logE = math.log10(math.e) #// for debug output
      #//System.out.println("molPops: nmrtrDissE " + nmrtrDissE + " log10UwA " + log10UwA[0] + " " + log10UwA[1] + " nmrtrLog10UwB " +
      #//     nmrtrLog10UwB[0] + " " + nmrtrLog10UwB[1] + " nmrtrLog10QwAB " + logE*nmrtrLogQwAB[2] + " nmrtrLogMuAB " + logE*nmrtrLogMuAB
      #//     + " numMolsB " + numMolsB + " dissEArr " + dissEArr[0] + " log10UwBArr " + log10UwBArr[0][0] + " " + log10UwBArr[0][1] + " log10QwABArr " +
      #//     logE*logQwABArr[0][2] + " logMuABArr " + logE*logMuABArr[0]);
      #//System.out.println("Line: nmrtrLog10UwB[0] " + logE*nmrtrLog10UwB[0] + " nmrtrLog10UwB[1] " + logE*nmrtrLog10UwB[1]);
  
      ln10 = math.log(10.0)
      log2pi = math.log(2.0 * math.pi)
      log2 = math.log(2.0)
  
      logE10 = math.log(10.0)
  #// Convert to natural logs:
      #double Ttheta, thisTemp;
  
  #//Treat at least one molecule - if there are really no molecules for an atomic species, 
  #//there will be one phantom molecule in the denominator of the ionization fraction
  #//with an impossibly high dissociation energy
      if (numMolsB == 0):
          numMolsB = 1
  #//This should be inherited, but let's make sure: 
          dissEArr[0] = 29.0 #//eV
  
      #//var molPops = function(logNum, numeratorLogNumB, numeratorDissE, numeratorLog10UwA, numeratorLog10QwAB, numeratorLogMuAB,  //species A data - ionization equilibrium of A
      #//Molecular partition functions - default initialization:
      thisLogUwB = [0.0 for i in range(numMolsB)]
      
      for iMol in range(numMolsB):
          thisLogUwB[iMol] = 0.0 #// variable for temp-dependent computed partn fn of array element B 
         
      thisLogUwA = 0.0 #// element A 
      nmrtrThisLogUwB = 0.0 #// element A 
      thisLogQwAB = math.log(300.0)
      nmrtrThisLogQwAB = math.log(300.0)
  
  #//For clarity: neutral stage of atom whose ionization equilibrium is being computed is element A
  #// for molecule formation:
      #logUwA = [0.0 for i in range(5)]
     
      #nmrtrLogUwB = [0.0 for i in range(5)]
      
      #for kk in range(len(logUwA)):
          #logUwA[kk] = logE10*log10UwA[kk]
          #nmrtrLogUwB[kk] = logE10*nmrtrLog10UwB[kk]
  
          #// lburns 
      
  #// Array of elements B for all molecular species AB:
      #double[][] logUwB = new double[numMolsB][2];
      #logUwB = [ [ 0.0 for i in range(5) ] for j in range(numMolsB) ]
      #//if (numMolsB > 0){
      #for iMol in range(numMolsB):
      #    for kk in range(5):
      #        logUwB[iMol][kk] = logE10*log10UwBArr[iMol][kk]
             # // lburns new loop
          
          
      #//}
  #// Molecular partition functions:
  #//       double nmrtrLogQwAB = logE10*nmrtrLog10QwAB;
  #//       double[] logQwAB = new double[numMolsB];
  #//      //if (numMolsB > 0){
  #//       for (int iMol = 0; iMol < numMolsB; iMol++){
  #//          logQwAB[iMol] = logE10*log10QwABArr[iMol];
  #//       }
  #      //}
  #//Molecular dissociation Boltzmann factors:
      nmrtrBoltzFacIAB = 0.0
      nmrtrLogMolSahaFac = 0.0
      logDissE = math.log(nmrtrDissE)  + Useful.logEv()
      #//System.out.println("logDissE " + logE*logDissE)
      logBoltzFacIAB = logDissE  - Useful.logK()
      #//System.out.println("logBoltzFacIAB " + logE*logBoltzFacIAB);
      nmrtrBoltzFacIAB = math.exp(logBoltzFacIAB)
      nmrtrLogMolSahaFac = (3.0 / 2.0) * (log2pi + nmrtrLogMuAB  + Useful.logK() - 2.0 * Useful.logH())
      #//System.out.println("nmrtrLogMolSahaFac " + logE*nmrtrLogMolSahaFac);
      #//System.out.println("nmrtrDissE " + nmrtrDissE + " logDissE " + logE*logDissE + " logBoltzFacIAB " + logE*logBoltzFacIAB + " nmrtrBoltzFacIAB " + nmrtrBoltzFacIAB + " nmrtrLogMuAB " + logE*nmrtrLogMuAB + " nmrtrLogMolSahaFac " + logE*nmrtrLogMolSahaFac);
      
      boltzFacIAB = [0.0 for i in range(numMolsB)]
      logMolSahaFac = [0.0 for i in range(numMolsB)]
      #//if (numMolsB > 0){
      for iMol in range(numMolsB):
          logDissE = math.log(dissEArr[iMol]) + Useful.logEv() 
          logBoltzFacIAB = logDissE  - Useful.logK()
          boltzFacIAB[iMol] = math.exp(logBoltzFacIAB)
          logMolSahaFac[iMol] = (3.0 / 2.0) * (log2pi + logMuABArr[iMol] + Useful.logK() - 2.0 * Useful.logH())
    #//System.out.println("logMolSahaFac[iMol] " + logE*logMolSahaFac[iMol]);
    #//System.out.println("iMol " + iMol + " dissEArr[iMol] " + dissEArr[iMol] + " logDissE " + logE*logDissE + " logBoltzFacIAB " + logE*logBoltzFacIAB + " boltzFacIAB[iMol] " + boltzFacIAB[iMol] + " logMuABArr " + logE*logMuABArr[iMol] + " logMolSahaFac " + logE*logMolSahaFac[iMol]);
                 
    #//double[] logNums = new double[numDeps]
   
    #//}
  #//   For molecular species:
      #double nmrtrSaha, nmrtrLogSahaMol, nmrtrLogInvSahaMol; //, nmrtrInvSahaMol;
      logMolFrac = [0.0 for i in range(numDeps)]
      logSahaMol = [0.0 for i in range(numMolsB)]
      invSahaMol = [0.0 for i in range(numMolsB)]
      
      
      #JB#
      currentUwAArr=list(logUwA)#u(T) determined values
      UwAFit = ToolBox.cubicFit(masterTemp, currentUwAArr)#u(T) fit    
  
      nmrtrLogUwBArr=list(nmrtrLogUwB)#u(T) determined values
      nmrtrLogUwBFit = ToolBox.cubicFit(masterTemp, nmrtrLogUwBArr)#u(T) fit
  
      #uwa.append(UwAFit)
      #uwb.append(nmrtrLogUwBFit)
      
      uwbFits=[]
      qwabFit = []
      for iMol in range(numMolsB):
            currentUwBArr=list(logUwB[iMol])
            UwBFit = ToolBox.cubicFit(masterTemp, currentUwBArr)
            uwbFits.append(UwBFit)
            
            currentLogQwABArr=list(logQwABArr[iMol])#u(T) determined values
            QwABFit = ToolBox.cubicFit(masterTemp, currentLogQwABArr)#u(T) fit
            qwabFit.append(QwABFit)
            
      #nmrtrQwABArr=list(nmrtrLogQwAB)#u(T) determined values
      #nmrtrQwABFit = ToolBox.cubicFit(masterTemp, nmrtrQwABArr)#u(T) fit
  
      #for Mols in range(numMolsB):
      #    currentLogUwBArr=list(logUwB[Mols])#u(T) determined values
      #    UwBFit=cubicFit(masterTemp,currentLogUwBArr)#u(T) fit
          
      #JB#
  #//
      
      temps=[]
      #valb=[]
      #vala=[]
      #valnb=[]
      #valqab=[]
      #valnmrtrqwb=[]
      #//  System.out.println("molPops: id      nmrtrLogNumB      logNumBArr[0]      logGroundRatio");
      for id in range(numDeps):
          
  
          #//System.out.format("%03d, %21.15f, %21.15f, %21.15f, %n", id, logE*nmrtrLogNumB[id], logE*logNumB[0][id], logE*logGroundRatio[id]);
  
          #//// reduce or enhance number density by over-all Rosseland opcity scale parameter
  
  #//Determine temparature dependent partition functions Uw:
          thisTemp = temp[0][id]
          temps.append(thisTemp)
          #Ttheta = 5040.0 / thisTemp
          """
          if (Ttheta >= 1.0):
              thisLogUwA = logUwA[0]
              nmrtrThisLogUwB = nmrtrLogUwB[0]
              for iMol in range(numMolsB):
                  thisLogUwB[iMol] = logUwB[iMol][0]
             
         
          if (Ttheta <= 0.5):
              thisLogUwA = logUwA[1]
              nmrtrThisLogUwB = nmrtrLogUwB[1]
              for iMol in range(numMolsB):
                  thisLogUwB[iMol] = logUwB[iMol][1]
             
          if (Ttheta > 0.5 and Ttheta < 1.0):
              thisLogUwA = ( logUwA[1] * ((Ttheta - 0.5)/(1.0 - 0.5)) ) \
              + ( logUwA[0] * ((1.0 - Ttheta)/(1.0 - 0.5)) )
              nmrtrThisLogUwB = ( nmrtrLogUwB[1] * ((Ttheta - 0.5)/(1.0 - 0.5)) ) \
              + ( nmrtrLogUwB[0] * ((1.0 - Ttheta)/(1.0 - 0.5)) )
              for iMol in range(numMolsB):
                  thisLogUwB[iMol] = ( logUwB[iMol][1] * ((Ttheta - 0.5)/(1.0 - 0.5)) ) \
                  + ( logUwB[iMol][0] * ((1.0 - Ttheta)/(1.0 - 0.5)) )
          """
      #JB#
          thisLogUwA = float(ToolBox.valueFromFit(UwAFit,thisTemp))#u(T) value extrapolated
          #vala.append(thisLogUwA)
          
          nmrtrThisLogUwB = float(ToolBox.valueFromFit(nmrtrLogUwBFit,thisTemp))#u(T) value extrapolated
          #valnb.append(nmrtrThisLogUwB)
          
      #for iMol in range(numMolsB):
      #    thisLogUwB[iMol]=logUwB[iMol]
      
          for iMol in range(numMolsB):
              thisLogUwB[iMol] = ToolBox.valueFromFit(uwbFits[iMol],thisTemp)#u(T) value extrapolated
              #valb.append(thisLogUwB[iMol])
              
          
  #// NEW Determine temperature dependent partition functions Uw: lburns
          thisTemp = temp[0][id]
          if (thisTemp <= 130.0):
              thisLogUwA = logUwA[0]
              nmrtrThisLogUwB = nmrtrLogUwB[0]
              for iMol in range(numMolsB):
                  thisLogUwB[iMol] = logUwB[iMol][0]
          if (thisTemp >= 10000.0):
              thisLogUwA = logUwA[4]
              nmrtrThisLogUwB = nmrtrLogUwB[4]
              for iMol in range(numMolsB):
                  thisLogUwB[iMol] = logUwB[iMol][4]           
          """
          if (thisTemp > 130 and thisTemp <= 500):
              thisLogUwA = logUwA[1] * (thisTemp - 130)/(500 - 130) \
                         + logUwA[0] * (500 - thisTemp)/(500 - 130)
              nmrtrThisLogUwB = nmrtrLogUwB[1] * (thisTemp - 130)/(500 - 130) \
                              + nmrtrLogUwB[0] * (500 - thisTemp)/(500 - 130)
              for iMol in range(numMolsB):
                  thisLogUwB[iMol] = logUwB[iMol][1] * (thisTemp - 130)/(500 - 130) \
                                   + logUwB[iMol][0] * (500 - thisTemp)/(500 - 130)
              
          
          if (thisTemp > 500 and thisTemp <= 3000):
              thisLogUwA = logUwA[2] * (thisTemp - 500)/(3000 - 500) \
                         + logUwA[1] * (3000 - thisTemp)/(3000 - 500)
              nmrtrThisLogUwB = nmrtrLogUwB[2] * (thisTemp - 500)/(3000 - 500) \
                              + nmrtrLogUwB[1] * (3000 - thisTemp)/(3000 - 500)
              for iMol in range(numMolsB):
                  thisLogUwB[iMol] = logUwB[iMol][2] * (thisTemp - 500)/(3000 - 500) \
                                   + logUwB[iMol][1] * (3000 - thisTemp)/(3000 - 500)
              
          
          if (thisTemp > 3000 and thisTemp <= 8000):
              thisLogUwA = logUwA[3] * (thisTemp - 3000)/(8000 - 3000) \
                         + logUwA[2] * (8000 - thisTemp)/(8000 - 3000)
              nmrtrThisLogUwB = nmrtrLogUwB[3] * (thisTemp - 3000)/(8000 - 3000) \
                              + nmrtrLogUwB[2] * (8000 - thisTemp)/(8000 - 3000)
              for iMol in range(numMolsB):
                  thisLogUwB[iMol] = logUwB[iMol][3] * (thisTemp - 3000)/(8000 - 3000) \
                                   + logUwB[iMol][2] * (8000 - thisTemp)/(8000 - 3000)
              
         
          if (thisTemp > 8000 and thisTemp < 10000):
              thisLogUwA = logUwA[4] * (thisTemp - 8000)/(10000 - 8000) \
                         + logUwA[3] * (10000 - thisTemp)/(10000 - 8000)
              nmrtrThisLogUwB = nmrtrLogUwB[4] * (thisTemp - 8000)/(10000 - 8000) \
                              + nmrtrLogUwB[3] * (10000 - thisTemp)/(10000 - 8000)
              for iMol in range(numMolsB):
                  thisLogUwB[iMol] = logUwB[iMol][4] * (thisTemp - 8000)/(10000 - 8000) \
                                  + logUwB[iMol][3] * (10000 - thisTemp)/(10000 - 8000)
              
          
          if (thisTemp >= 10000):
              thisLogUwA = logUwA[4]
              nmrtrThisLogUwB = nmrtrLogUwB[4]
              for iMol in range(numMolsB):
                  thisLogUwB[iMol] = logUwB[iMol][4]
  
          """
  
  
          #iMol loops for Q's
          
          for iMol in range(numMolsB):
              if (thisTemp < 3000.0):
                  thisLogQwAB = ( logQwABArr[iMol][1] * (3000.0 - thisTemp)/(3000.0 - 500.0) ) \
                  + ( logQwABArr[iMol][2] * (thisTemp - 500.0)/(3000.0 - 500.0) )         
              if ( (thisTemp >= 3000.0) and (thisTemp <= 8000.0) ):
                  thisLogQwAB = ( logQwABArr[iMol][2] * (8000.0 - thisTemp)/(8000.0 - 3000.0) ) \
                  + ( logQwABArr[iMol][3] * (thisTemp - 3000.0)/(8000.0 - 3000.0) )         
              if ( thisTemp > 8000.0 ):
                  thisLogQwAB = ( logQwABArr[iMol][3] * (10000.0 - thisTemp)/(10000.0 - 8000.0) ) \
                  + ( logQwABArr[iMol][4] * (thisTemp - 8000.0)/(10000.0 - 8000.0) )      
                  
              if (thisTemp < 3000.0):
                  nmrtrThisLogQwAB = ( nmrtrLogQwAB[1] * (3000.0 - thisTemp)/(3000.0 - 500.0) ) \
                  + ( nmrtrLogQwAB[2] * (thisTemp - 500.0)/(3000.0 - 500.0) )         
              if ( (thisTemp >= 3000.0) and (thisTemp <= 8000.0) ):
                  nmrtrThisLogQwAB = ( nmrtrLogQwAB[2] * (8000.0 - thisTemp)/(8000.0 - 3000.0) ) \
                  + ( nmrtrLogQwAB[3] * (thisTemp - 3000.0)/(8000.0 - 3000.0) )         
              if ( thisTemp > 8000.0 ):
                  nmrtrThisLogQwAB = ( nmrtrLogQwAB[3] * (10000.0 - thisTemp)/(10000.0 - 8000.0) ) \
                  + ( nmrtrLogQwAB[4] * (thisTemp - 8000.0)/(10000.0 - 8000.0) ) 
          
          
  #//For clarity: neutral stage of atom whose ionization equilibrium is being computed is element A
  #// for molecule formation:
  
  #   //Ionization stage Saha factors: 
  #//System.out.println("id " + id + " nmrtrLogNumB[id] " + logE*nmrtrLogNumB[id]);   
  #  // if (id == 16){ 
  #  //   System.out.println("id " + id + " nmrtrLogNumB[id] " + logE*nmrtrLogNumB[id] + " pp nmrtB " + (logE*(nmrtrLogNumB[id]+temp[1][id]+Useful.logK())) + " nmrtrThisLogUwB " + logE*nmrtrThisLogUwB + " thisLogUwA " + logE*thisLogUwA + " nmrtrLogQwAB " + logE*nmrtrThisLogQwAB);         
  #   //System.out.println("nmrtrThisLogUwB " + logE*nmrtrThisLogUwB + " thisLogUwA " + logE*thisLogUwA + " nmrtrThisLogQwAB " + logE*nmrtrThisLogQwAB);
  #   // }
          nmrtrLogSahaMol = nmrtrLogMolSahaFac - nmrtrLogNumB[id] - (nmrtrBoltzFacIAB / temp[0][id]) + (3.0 * temp[1][id] / 2.0) + nmrtrThisLogUwB + thisLogUwA - nmrtrThisLogQwAB
          nmrtrLogInvSahaMol = -1.0 * nmrtrLogSahaMol
          #//System.out.println("nmrtrLogInvSahaMol " + logE*nmrtrLogInvSahaMol);
          #//nmrtrInvSahaMol = Math.exp(nmrtrLogSahaMol);
          #//   if (id == 16){
          #//       System.out.println("nmrtrLogInvSahaMol " + logE*nmrtrLogInvSahaMol);
          #//   }
          #//   if (id == 16){
          #//        System.out.println("nmrtrBoltzFacIAB " + nmrtrBoltzFacIAB + " nmrtrThisLogUwB " + logE*nmrtrThisLogUwB + " thisLogUwA " + logE*thisLogUwA + " nmrtrThisLogQwAB " + nmrtrThisLogQwAB);   
          #//        System.out.println("nmrtrLogSahaMol " + logE*nmrtrLogSahaMol); // + " nmrtrInvSahaMol " + nmrtrInvSahaMol);
          #//   }
  
  #//Molecular Saha factors:
          for iMol in range(numMolsB):
  #//System.out.println("iMol " + iMol + " id " + id + " logNumB[iMol][id] " + logE*nmrtrLogNumB[id]);             
  #//System.out.println("iMol " + iMol + " thisLogUwB[iMol] " + logE*thisLogUwB[iMol] + " thisLogUwA " + logE*thisLogUwA + " thisLogQwAB " + logE*thisLogQwAB);
              logSahaMol[iMol] = logMolSahaFac[iMol] - logNumB[iMol][id] - (boltzFacIAB[iMol] / temp[0][id]) + (3.0 * temp[1][id] / 2.0) + float(thisLogUwB[iMol]) + thisLogUwA - thisLogQwAB
  #//For denominator of ionization fraction, we need *inverse* molecular Saha factors (N_AB/NI):
              logSahaMol[iMol] = -1.0 * logSahaMol[iMol]
              invSahaMol[iMol] = math.exp(logSahaMol[iMol])
              #//TEST invSahaMol[iMol] = 1.0e-99; //test
       #//     if (id == 16){
       #//         System.out.println("iMol " + iMol + " boltzFacIAB[iMol] " + boltzFacIAB[iMol] + " thisLogUwB[iMol] " + logE*thisLogUwB[iMol] + " logQwAB[iMol] " + logE*thisLogQwAB + " logNumB[iMol][id] " + logE*logNumB[iMol][id] + " logMolSahaFac[iMol] " + logE*logMolSahaFac[iMol]);   
       #//         System.out.println("iMol " + iMol + " logSahaMol " + logE*logSahaMol[iMol] + " invSahaMol[iMol] " + invSahaMol[iMol]);
       #//     }
           
  
  #//Compute log of denominator is ionization fraction, f_stage 
  #        //default initialization 
  #        //  - ratio of total atomic particles in all ionization stages to number in ground state: 
          denominator = math.exp(logGroundRatio[id]) #//default initialization - ratio of total atomic particles in all ionization stages to number in ground state 
  #//molecular contribution
          for iMol in range(numMolsB):
              #//  if (id == 16){
              #//   System.out.println("invSahaMol[iMol] " + invSahaMol[iMol] + " denominator " + denominator);
              #//  }
              denominator = denominator + invSahaMol[iMol]
             
  #// 
          logDenominator = math.log(denominator) 
          #//System.out.println("logGroundRatio[id] " + logE*logGroundRatio[id] + " logDenominator " + logE*logDenominator);
          #//  if (id == 16){
          #// System.out.println("id " + id + " logGroundRatio " + logGroundRatio[id] + " logDenominator " + logDenominator);
          #//  } 
          #//if (id == 36){
          #//     System.out.println("logDenominator " + logE*logDenominator);
          #// }
          #//var logDenominator = Math.log( 1.0 + saha21 + (saha32 * saha21) + (saha43 * saha32 * saha21) + (saha54 * saha43 * saha32 * saha21) );
  
          logMolFrac[id] = nmrtrLogInvSahaMol - logDenominator
          #//  if (id == 16){
          #//       System.out.println("id " + id + " logMolFrac[id] " + logE*logMolFrac[id]);
          #//  }
  
          #//logNums[id] = logNum[id] + logMolFrac;
      #} //id loop
      
      
      #JB - check (never used)#
      #print(uwa)
      #print(uwb)
      #title("logUwA")
      """
      plot(temps,vala)
      tempT=[]
      for t in masterTemp:
          tempT.append(valueFromFit(UwAFit,t))
      scatter(masterTemp,(tempT))
      show()
      
      #title("nmrtrlogUwB")
      plot(temps,valnb)
      tempT=[]
      for t in masterTemp:
          tempT.append(valueFromFit(nmrtrLogUwBFit,t))
      scatter(masterTemp,(tempT))
      show()
      
      #title("logUwB")
      plot(temps,valb)
      tempT=[]
      for t in masterTemp:
          tempT.append(valueFromFit(UwBFit,t))
      scatter(masterTemp,(tempT))
      show()
  
      #title("logQwAB")
      plot(temps,valqab)
      tempT=[]
      for t in masterTemp:
          tempT.append(valueFromFit(QwABFit,t))
      scatter(masterTemp,(tempT))
      show()
  
      #title("nmrtrlogQwAB")
      plot(temps,valnmrtrqwb)
      tempT=[]
      for t in masterTemp:
          tempT.append(valueFromFit(nmrtrQwABFit,t))
      scatter(masterTemp,(tempT))
      show()
      """
      #JB#
      return logMolFrac
      #//end method stagePops    