From dc94ba0decdf6a26a46e49184c96abfa5e841d0d Mon Sep 17 00:00:00 2001 From: Dave Goodwin Date: Mon, 17 Dec 2007 07:51:37 +0000 Subject: [PATCH] *** empty log message *** --- Cantera/src/spectra/rotor.cpp | 74 ++++++++++++++++++---- Cantera/src/spectra/rotor.h | 21 +++++- test_problems/spectroscopy/Makefile.in | 9 ++- test_problems/spectroscopy/spectratest.cpp | 21 ++++-- 4 files changed, 104 insertions(+), 21 deletions(-) diff --git a/Cantera/src/spectra/rotor.cpp b/Cantera/src/spectra/rotor.cpp index 394c620de..af0b62c51 100644 --- a/Cantera/src/spectra/rotor.cpp +++ b/Cantera/src/spectra/rotor.cpp @@ -1,4 +1,7 @@ - +/** + * @file rotor.cpp + * +*/ #ifdef WIN32 #pragma warning(disable:4786) #pragma warning(disable:4503) @@ -9,10 +12,13 @@ namespace Cantera { - /* - * @param mu reduced mass in kg - * @param re bond length in meters - * @dipoleMoment permanent dipole moment in ... + /** + * @param Bv Rotational constant, wavenumbers + * @param Dv Coefficient describing centrifugal + * effects on the bond length. For a rigid rotor, Bv = 0. + * @param Hv Coefficient describing higher-order vibration-rotation + * interactions. For a rigid rotor, Hv = 0. + * @dipoleMoment permanent dipole moment. */ Rotor::Rotor(doublereal Bv, doublereal dipoleMoment, doublereal Dv, doublereal Hv ) : m_Bv(Bv), @@ -20,30 +26,72 @@ namespace Cantera { m_Hv(Hv), m_dipole(dipoleMoment) {} - // energy in wavenumbers + /** + * The energy of the level with rotational quantum number J, + * in wavenumber units. + * \f[ + * E(J) = J(J+1)B - [J(J+1)]^2 D + [J(J+1)]^3H + * \f] + * For a rigid rotor, only B is non-zero. The parameters B, D, and H + * are set in the constructor. + */ doublereal Rotor::energy_w(int J) { int jjp1 = J*(J + 1); return jjp1*(m_Bv + jjp1*(m_Hv*jjp1 - m_Dv)); } - doublereal Rotor::degeneracy(int J) { + /** + * The number of quantum states with the same J. + */ + int Rotor::degeneracy(int J) { return 2*J + 1; } + /** + * The rotational partition function. + * + * If T/Trot > 100, then the classical value (T/Trot) is + * is returned. Otherwise, it is computed as a sum + * \f[ + * z = \sum_{J=0}^{J_{max} (2J + 1) \exp(-E(J)/kT) + * \] + */ doublereal Rotor::partitionFunction(doublereal T, int cutoff) { int j = 0; - if (cutoff < 0) cutoff = 100; - doublereal dsum = 0.0, sum = 0.0; - for (j = 0; j < cutoff; j++) { - dsum = degeneracy(j)*exp(-Planck*energy_w(j)/(Boltzmann * T)); - sum += dsum; + doublereal T_Trot = wnum_to_J(m_Bv)/(Boltzmann*T); + if (T_Trot > 100.0) + return T_Trot; + else { + if (cutoff < 0) cutoff = 3.0*sqrt(T/m_Bv); + doublereal dsum = 0.0, sum = 0.0; + for (j = 0; j < cutoff; j++) { + dsum = degeneracy(j)*exp(-wnum_to_J(energy_w(j))/(Boltzmann * T)); + sum += dsum; + } + return sum; } - return sum; } + /** + * Ratio of the population of all states with rotational quantum + * number J to the ground state population. + */ + doublereal Rotor::relPopulation(int J, doublereal T) { + return degeneracy(J)*exp(-wnum_to_J(energy_w(J))/(Boltzmann*T)); + } + + /** The difference in the energies of an upper and a lower state. + * + */ doublereal Rotor::frequency(int J_lower, int J_upper) { return (energy_w(J_upper) - energy_w(J_lower)); } + + doublereal Rotor::intensity(int J_lower, int J_upper, doublereal T) { + int dJ = J_upper - J_lower; + if (dJ > 1 || dJ < -1) return 0; + return relPopulation(J_lower, T); + } } diff --git a/Cantera/src/spectra/rotor.h b/Cantera/src/spectra/rotor.h index 0aa699258..b037b030d 100644 --- a/Cantera/src/spectra/rotor.h +++ b/Cantera/src/spectra/rotor.h @@ -1,3 +1,8 @@ +/** + * @file rotor.h + * Header file for class Rotor. + */ + #ifndef CT_ROTOR #define CT_ROTOR @@ -5,6 +10,10 @@ namespace Cantera { + /** + * Class Rotor represents a non-rigid quantum-mechanical rotor. + * @ingroup Spectroscopy + */ class Rotor { public: Rotor() {} @@ -16,15 +25,23 @@ namespace Cantera { doublereal Dv = 0.0, doublereal Hv = 0.0); - // energy in wavenumbers + /// energy in wavenumbers doublereal energy_w(int J); - doublereal degeneracy(int J); + /// degeneracy + int degeneracy(int J); doublereal partitionFunction(doublereal T, int cutoff=-1); doublereal frequency(int J_lower, int J_upper); + doublereal relPopulation(int J, doublereal T); + + doublereal population(int J, doublereal T) { + return relPopulation(J,T)/partitionFunction(T); + } + doublereal intensity(int J_lower, int J_upper, doublereal T); + protected: doublereal m_Bv; diff --git a/test_problems/spectroscopy/Makefile.in b/test_problems/spectroscopy/Makefile.in index 0bfdbe0a8..d369747de 100644 --- a/test_problems/spectroscopy/Makefile.in +++ b/test_problems/spectroscopy/Makefile.in @@ -53,6 +53,13 @@ CANTERA_LIBDIR=@buildlib@ # required Cantera libraries CANTERA_LIBS = @LOCAL_LIBS@ -lctcxx +# Cantera library Dependencies +CANTERA_LIB_DEPS = $(CANTERA_LIBDIR)/libtransport.a \ + $(CANTERA_LIBDIR)/libthermo.a \ + $(CANTERA_LIBDIR)/libctspectra.a \ + $(CANTERA_LIBDIR)/libctnumerics.a \ + $(CANTERA_LIBDIR)/libctbase.a + # the directory where Cantera include files may be found. ifeq ($(src_dir_tree), 1) CANTERA_INCDIR=../../Cantera/src @@ -77,7 +84,7 @@ PROGRAM = $(PROG_NAME)$(EXE_EXT) all: $(PROGRAM) -$(PROGRAM): $(OBJS) $(CANTERA_LIBDIR)/libctbase.a +$(PROGRAM): $(OBJS) $(CANTERA_LIB_DEPS) $(CXX) -o $(PROGRAM) $(OBJS) $(LCXX_FLAGS) $(LINK_OPTIONS) \ $(CANTERA_LIBS) @LIBS@ $(FORT_LIBS) \ $(LCXX_END_LIBS) diff --git a/test_problems/spectroscopy/spectratest.cpp b/test_problems/spectroscopy/spectratest.cpp index 776470429..dc997ea0b 100644 --- a/test_problems/spectroscopy/spectratest.cpp +++ b/test_problems/spectroscopy/spectratest.cpp @@ -6,9 +6,20 @@ using namespace std; using namespace Cantera; int main() { - - Rotor* r = new Rotor(1.0); - double w = 8065.0; - cout << "eV: " << wnum_to_eV(w) << endl; - delete r; + double B; + double T; + cout << "enter B, T: "; + cin >> B >> T; + Rotor* r = new Rotor(B); + double theta = wnum_to_J(B)/Boltzmann; + cout << "theta = " << theta << endl; + double pop, cpop = 0.0; + for (int j = 0; j < 50; j++) { + pop = r->population(j, T); + cpop += pop; + if (cpop > 0.999) break; + cout << j << ", " << r->energy_w(j) << ", " + << r->frequency(j, j+1) << ", " + << pop << ", " << cpop << ", " << r->partitionFunction(T) << ", " << T/theta << endl; + } }