First, I took the ising2d.c file as a starting point. I implemented periodic boundary conditions by setting spins of the boundary to be the same. I then added a third integrated dimension to the Monte Carlo integration steps. By integrating over the spins in the lattice, the spontaneous magnetization was measured and plot as a function of temperature for $T < T_c$. I used both a hot and a cold start to do this to ensure that the system was working properly. Here is a plot showing both of the plots together, exhibiting this behavior:
From this graph, we can also see that the critical temperature is around 4.5. Here is a more detailed view of the graph:
We can verify here that the critical temperature is around 4.5. For both runs, I used 1000 Monte Carlo integration sweeps.
2. Critical Temperature via Spin Correlator
We can also find the critical temperature via the spin correlator $G(x,y,z) = \left<\sigma(x_0+x,y_0+y,z_0+z)\sigma(x_0,y_0,z_0)\right>$ where $\sigma(x,y,z)$ is the spin at site $(x,y,z)$. This will gives us the correlation of a point $(x_0,y_0,z_0)$ and another point $(x,y,z)$ away from that point. We can define a more convenient function of one variable that is the correlation function integrated acrodd two dimensions, i.e. $\tilde{G}(z)=\displaystyle\sum_{x,y}G(x,y,z)$.
I modified the code from the previous section to calculate this function for a given temperature value. I then ran the program several times for $T > T_c$. It can be shown that $\tilde{G}(z) \sim \mathrm{const} \times e^{\frac{-z}{\xi}}$.
Solving for $M=\xi^{-1}$, we can solve for $M = -\frac{\mathrm{ln}(\tilde{G}(z))}{z}$. I therefore plotted the M value I found from running my program from 5.7 to 4.5 in steps of 0.1 at the longest offset (z=10). Here is the plot:
As you can see this agrees with the section 1 answer that $T_c \simeq 4.5$.
