diff --git a/Lecture32.tgz b/Lecture32.tgz
new file mode 100644
index 0000000..becc676
Binary files /dev/null and b/Lecture32.tgz differ
diff --git a/Lecture32/Makefile b/Lecture32/Makefile
new file mode 100644
index 0000000..3c85d2d
--- /dev/null
+++ b/Lecture32/Makefile
@@ -0,0 +1,19 @@
+# I am a comment, and I want to say that the variable CC will be
+# the compiler to use.
+CC=g++
+# Hey!, I am comment number 2. I want to say that CFLAGS will be the
+# options I'll pass to the compiler.
+LIBS=-I$(INCLUDEPATH) -L$(LD_LIBRARY_PATH) -lcpt
+
+
+all: ising
+
+ising: ising.cpp
+	$(CC) $^  $(LIBS) -o ising
+
+
+
+clean:
+	rm -rf *o ising
+
+
diff --git a/Lecture32/ising.cpp b/Lecture32/ising.cpp
new file mode 100644
index 0000000..e404e3d
--- /dev/null
+++ b/Lecture32/ising.cpp
@@ -0,0 +1,149 @@
+#include "cptstd.hpp"
+#include "linalg.hpp"
+#include "random.hpp"
+using namespace cpt;
+
+
+class Ising {
+public : 
+
+  Ising(double iJ=1.0, int iL=10, int iN=100, double iT=2.0, double iH=0.0) :
+    J(iJ), L(iL), N(iN), T(iT), H(iH), s(L,L), Lx(L), Ly(L)
+  {
+    s = Matrix<int,2>(Lx, Ly);
+    for (int i = 0; i < Lx; i++)
+      for (int j = 0; j < Ly; j++)
+	s[i][j] =  rng.rand() < 0.5 ? +1 : -1;   // hot start
+    compute_boltzmann_factors();
+    steps = 0;
+  }
+
+
+  
+  void compute_boltzmann_factors()
+  {
+    for (int i = -8; i <= 8; i += 4) {
+      w[i + 8][0] = exp( - (i * J - 2 * H) / T);
+      w[i + 8][2] = exp( - (i * J + 2 * H) / T);
+    }
+  }
+
+  bool metropolis_step()
+  {
+    // choose a random spin
+    int i = int(Lx * rng.rand());
+    int j = int(Ly * rng.rand());
+
+    // find its neighbors using periodic boundary conditions
+    int iPrev = i == 0 ? Lx-1 : i-1;
+    int iNext = i == Lx-1 ? 0 : i+1;
+    int jPrev = j == 0 ? Ly-1 : j-1;
+    int jNext = j == Ly-1 ? 0 : j+1;
+
+    // find sum of neighbors
+    int sumNeighbors = s[iPrev][j] + s[iNext][j] + s[i][jPrev] + s[i][jNext];
+    int delta_ss = 2*s[i][j]*sumNeighbors;
+
+    // ratio of Boltzmann factors
+    double ratio = w[delta_ss+8][1+s[i][j]];
+    if (rng.rand() < ratio) {
+      s[i][j] = -s[i][j];
+      return true;
+    } else return false;
+  }
+
+  double acceptanceRatio;
+
+  void one_monte_carlo_step_per_spin ( ) {
+    int accepts = 0;
+    for (int i = 0; i < N; i++)
+      if (metropolis_step())
+	++accepts;
+    acceptanceRatio = accepts/double(N);
+    ++steps;
+  }
+
+  double magnetizationPerSpin ( ) {
+    int sSum = 0;
+    for (int i = 0; i < Lx; i++)
+      for (int j = 0; j < Ly; j++) {
+	sSum += s[i][j];
+      }
+    return sSum / double(N);
+  }
+
+  double energyPerSpin ( ) {
+    int sSum = 0, ssSum = 0;
+    for (int i = 0; i < Lx; i++)
+      for (int j = 0; j < Ly; j++) {
+	sSum += s[i][j];
+	int iNext = i == Lx-1 ? 0 : i+1;
+	int jNext = j == Ly-1 ? 0 : j+1;
+	ssSum += s[i][j]*(s[iNext][j] + s[i][jNext]);
+      }
+    return -(J*ssSum + H*sSum)/N;
+  }
+
+
+protected : 
+
+  Random rng;                     // random number generator
+
+  double J;                       // ferromagnetic coupling
+  int L, Lx, Ly;                  // number of spins in x and y
+  int N;                          // number of spins
+  Matrix<int,2> s;                // the spins
+  double T;                       // temperature
+  double H;                       // magnetic field
+
+  double w[17][3];                // Boltzmann factors
+
+  int steps;                      // steps so far
+
+
+};
+
+int main (int argc, char *argv[]) {
+
+  int Lx, Ly, N;
+  double T, H; 
+     cout << " Two-dimensional Ising Model - Metropolis simulation\n"
+          << " ---------------------------------------------------\n"
+          << " Enter number of spins L in each direction: ";
+     cin >> Lx;
+     Ly = Lx;
+     N = Lx*Ly;
+     cout << " Enter temperature T: ";
+     cin >> T;
+     cout << " Enter magnetic field H: ";
+     cin >> H;
+     cout << " Enter number of Monte Carlo steps: ";
+     int MCSteps;
+     cin >> MCSteps;
+     Ising ising(1.0, Lx, N, T, H);
+
+     int thermSteps = int(0.2 * MCSteps);
+     cout << " Performing " << thermSteps
+          << " steps to thermalize the system ..." << flush;
+     for (int s = 0; s < thermSteps; s++)
+          ising.one_monte_carlo_step_per_spin();
+
+     cout << " Done\n Performing production steps ..." << flush;
+     double mAv = 0, m2Av = 0, eAv = 0, e2Av = 0;
+     ofstream dataFile("ising.data");
+     for (int s = 0; s < MCSteps; s++) {
+          ising.one_monte_carlo_step_per_spin();
+          double m = ising.magnetizationPerSpin();
+          double e = ising.energyPerSpin();
+          mAv += m; m2Av += m * m;
+          eAv += e; e2Av += e * e;
+          dataFile << m << '\t' << e << '\n';
+     }
+     dataFile.close();
+     mAv /= MCSteps; m2Av /= MCSteps;
+     eAv /= MCSteps; e2Av /= MCSteps;
+     cout << " \n Magnetization and energy per spin written in file "
+          << " \"ising.data\"" << endl;
+     cout << " <m> = " << mAv << " +/- " << sqrt(m2Av - mAv*mAv) << endl;
+     cout << " <e> = " << eAv << " +/- " << sqrt(e2Av - eAv*eAv) << endl;
+}
diff --git a/Lecture32/ising.py b/Lecture32/ising.py
new file mode 100644
index 0000000..b4a540f
--- /dev/null
+++ b/Lecture32/ising.py
@@ -0,0 +1,120 @@
+import math
+import random
+
+class Ising :
+    def __init__(self, J=1.0, L=10, N=100, T=2., H=0.) : 
+
+        self.J = J                       # spin-spin coupling += ferro, -= antiferro
+        self.L_x = L; self.L_y = L       # number of spins in x and y
+        self.N = N                       # total number of spins
+        self.s = []                      # L_x x L_y array of spin values
+        self.T = T                       # Temperature
+        self.H = H                       # magnetic field
+
+        self.w = []                      # Boltzmann factors at fixed T and H
+        self.steps = 0                   # Monte Carlo steps so far
+        self.acceptance_ratio = 0        # accepted steps / total number of steps
+
+        # create spin lattice and set spin randomly up or down (hot start)
+        for i in range(self.L_x):
+            self.s.append( [ ] )
+            for j in range(self.L_y):
+                self.s[i].append(random.choice( (-1, 1) ))
+        self.compute_Boltzmann_factors()
+        self.steps = 0
+
+
+
+    def compute_Boltzmann_factors(self):
+        self.w = []
+        for m in range(5):
+            self.w.append( [] )
+            sum_of_neighbors = -4 + 2 * m
+            for n in range(2):
+                s_i = -1 + 2 * n
+                factor = math.exp( -2.0 * (self.J * sum_of_neighbors + self.H) * s_i / self.T )
+                self.w[m].append(factor)
+
+
+
+    def Metropolis_step_accepted(self):
+
+        # choose a random spin
+        i = random.randrange(self.L_x)
+        j = random.randrange(self.L_y)
+
+        # find the sum of neighbors assuming periodic boundary conditions
+        sum_of_neighbors = ( self.s[(i-1)%self.L_x][j] + self.s[(i+1)%self.L_x][j] +
+                             self.s[i][(j-1)%self.L_y] + self.s[i][(j+1)%self.L_y] )
+
+        # access ratio of precomputed Boltzmann factors
+        ratio = self.w[2 + int(sum_of_neighbors/2)][int((1 + self.s[i][j])/2)]
+
+        # apply the Metropolis test
+        if ratio > 1.0 or ratio > random.random():
+            self.s[i][j] = -self.s[i][j]
+            return True
+        else:
+            return False
+
+    
+
+    def one_Monte_Carlo_step_per_spin(self):
+        accepts = 0
+        for n in range(self.N):
+            if self.Metropolis_step_accepted():
+                accepts += 1
+        self.acceptance_ratio = accepts / float(self.N)
+        self.steps += 1
+
+    def magnetization_per_spin(self):
+
+        s_sum = 0.0
+        for i in range(self.L_x):
+            for j in range(self.L_y):
+                s_sum += self.s[i][j]
+        return s_sum / float(self.N)
+
+    def energy_per_spin(self):
+
+        s_sum = 0.0
+        ss_sum = 0.0
+        for i in range(self.L_x):
+            for j in range(self.L_y):
+                s_sum += self.s[i][j]
+                ss_sum += self.s[i][j] * (self.s[(i+1)%self.L_x][j] + self.s[i][(j+1)%self.L_y])
+        return -(self.J * ss_sum + self.H * s_sum) / float(self.N)
+
+print " Two-dimensional Ising Model - Metropolis simulation"
+print " ---------------------------------------------------"
+L = int(input(" Enter number of spins L in each direction: "))
+T = float(input(" Enter temperature T: "))
+H = float(input(" Enter magnetic field H: "))
+MC_steps = int(input(" Enter number of Monte Carlo steps: "))
+ising = Ising(L=L, N=L*L, T=T, H=H)
+
+therm_steps = int(0.2 * MC_steps)
+print " Performing", therm_steps, "thermalization steps ..."
+for i in range(therm_steps):
+    ising.one_Monte_Carlo_step_per_spin()
+
+print " Done ... Performing production steps ..."
+m_av = 0.0; m2_av = 0.0; e_av = 0.0; e2_av = 0.0
+data_file = open("ising.data", "w")
+for i in range(MC_steps):
+    ising.one_Monte_Carlo_step_per_spin()
+    m = ising.magnetization_per_spin()
+    e = ising.energy_per_spin()
+    m_av += m
+    m2_av += m**2
+    e_av += e
+    e2_av += e**2
+    data_file.write(repr(m) + "\t" + repr(e) + "\n")
+data_file.close()
+print " M/spin and E/spin values written in ising.data"
+m_av /= float(MC_steps)
+m2_av /= float(MC_steps)
+e_av /= float(MC_steps)
+e2_av /= float(MC_steps)
+print " <m> =", m_av, "+/-", math.sqrt(m2_av - m_av**2)
+print " <e> =", e_av, "+/-", math.sqrt(e2_av - e_av**2)