diff --git a/Images/MCS_diif_L.png b/Images/MCS_diif_L.png
new file mode 100644
index 0000000..6872cb0
Binary files /dev/null and b/Images/MCS_diif_L.png differ
diff --git a/InitEnergy.py b/InitEnergy.py
new file mode 100644
index 0000000..52678c1
--- /dev/null
+++ b/InitEnergy.py
@@ -0,0 +1,18 @@
+def Ener(ising,L,J):
+    E=0
+    M=0
+    for i in range (L):
+        for j in range (L):
+            a,b,c,d=i+1,i-1,j+1,j-1
+            if i==L-1:
+                a=0
+            if i==0:
+                b=L-1
+            if j==L-1:
+                c=0
+            if j==0:
+                d=L-1
+        
+            M+=ising[i][j]
+            E=E-J*float(ising[i][j]*(ising[a][j]+ising[b][j]+ising[i][c]+ising[i][d]))
+    return (M,E)
\ No newline at end of file
diff --git a/LatticeLen.py b/LatticeLen.py
new file mode 100644
index 0000000..a8da6b6
--- /dev/null
+++ b/LatticeLen.py
@@ -0,0 +1,35 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from RandomSpin import randomIsing
+from InitEnergy import Ener
+from Montecarlo import MCS
+T=1.4
+J=1
+Len=[]
+for L in range(10,40,10):
+    ising=randomIsing(L)         
+    M,E=Ener(ising,L,J)
+    print(M,E)
+    E*=0.5
+    N=100000
+    # ener=[]
+    mag=[]
+    for n in range(N):
+        ising,E,M=MCS(ising,L,T,J,E,M)
+        # enspin=E/(L*L)
+        Mspin=M/(L*L)
+    
+        # ener.append(enspin)
+        mag.append(Mspin)
+
+    Len.append(mag)
+
+num=np.arange(0,n+1,1)
+plt.plot(num,Len[0],color='r',label='L=10')
+plt.plot(num,Len[1],color='g',label='L=20')
+plt.plot(num,Len[2],color='b',label='L=30')
+plt.title('M/N vs MCS')
+plt.xlabel('MCS')
+plt.ylabel('M/N')
+plt.legend()
+plt.show()
\ No newline at end of file
diff --git a/Montecarlo.py b/Montecarlo.py
new file mode 100644
index 0000000..4881f1b
--- /dev/null
+++ b/Montecarlo.py
@@ -0,0 +1,36 @@
+import numpy as np
+def MCS(ising,L,T,J,E,M):
+    for k in range(L):
+        for l in range (L):
+            r=np.random.uniform(0,1)
+            i=int(r*(L-1))
+            t=np.random.uniform(0,1)
+            j=int(t*(L-1))
+            a,b,c,d=i+1,i-1,j+1,j-1
+            if i==L-1:
+                a=0
+            if i==0:
+                b=L-1
+            if j==L-1:
+                c=0
+            if j==0:
+                d=L-1
+            Ei=-J*float(ising[i][j]*(ising[a][j]+ising[b][j]+ising[i][c]+ising[i][d]))
+            ising[i][j]=-ising[i][j]
+            Ef=-J*float(ising[i][j]*(ising[a][j]+ising[b][j]+ising[i][c]+ising[i][d]))
+            dE=Ef-Ei
+            
+            if(dE<=0.0):
+                E+=dE
+                M+=(2.0*float(ising[i][j]))
+            else:
+                u=np.exp(-dE/T)
+                h=np.random.uniform(0,1)
+                if(h<u):
+                    E+=dE
+                    M+=(2.0*float(ising[i][j]))
+
+                else:
+                    ising[i][j]=-ising[i][j]
+        
+    return(ising,E,M)
\ No newline at end of file
diff --git a/README.md b/README.md
index 0981234..79480a0 100644
--- a/README.md
+++ b/README.md
@@ -1,4 +1,29 @@
 # Ising Model
 Ising model is a mathematical model to used to understant Ferromagnetism an its phase transition to paramagnetic behavior above critical temperature. This model consists 
-of a N  fixed lattice sites arranged to form a n(n=1,2,3) dimensional array. Each of the lattice sites has a spin $S_i$ associated with it, which can take values +1 or -1.
+of a N fixed lattice sites arranged to form a n(n=1,2,3) dimensional array. Each of the lattice sites has a spin $S_i$ associated with it, which can take values +1 or -1.
 The entire configuration of spin { $S_i$ } can be used to determine the energy of the system which is given by the hamiltonian:
+
+![Alt text](Images/Screenshot%20from%202024-02-06%2017-21-36-modified.png)
+
+Here $\epsilon_{ij}$ is the interaction energy between nearest neighbour spin and H is the external magnetic field. We can simplify the problem by considering an isotropic interaction with a constant $\epsilon$. The hamiltonian becomes:
+
+![Alt text](Images/Screenshot%20from%202024-02-06%2017-21-57-modified.png)
+
+Given the energy we can find all the thermodynamic variables by finding the partition function given by:
+
+![Alt text](Images/Screenshot%20from%202024-02-06%2017-21-57-modified%20(1).png)
+
+The following thermodynamic quantities are to be found to understand the phase transitions from ferromagnetic to paramagnetic transition:
+
+![Alt text](Images/Screenshot%20from%202024-02-06%2017-21-57-modified%20(3).png)
+
+# Phase Transitions
+
+ The ferromagnetic to paramagnetic phase transition is a second order phase transition. When the temperature T
+ approaches the Curie temperature T,the magnetization M(T),which is the order parameter of the transition, continuously goes to zero. Above Tc, the spins point in different directions and their magnetic fields are canceled. This disordered configuration is
+caused by the random thermal movement of the spins. The higher the temperature, the more difficult it is for any orderly arrangement of spins to be maintained. 
+But there is still interactions with the neighbours unlike a paramagnetic material. But when the temperature is lowered below critical temperature ( $T_c$ )the interactions become stronger then thermal 
+fluctuations and the spins align
+spontaneously. Instead of canceling each other out, the
+individual magnetic fields are added, producing a macroscopic
+magnetic field.
diff --git a/RandomSpin.py b/RandomSpin.py
new file mode 100644
index 0000000..2387702
--- /dev/null
+++ b/RandomSpin.py
@@ -0,0 +1,13 @@
+import numpy as np
+
+def randomIsing(L):
+    ising=np.zeros((L,L),dtype=float)
+    for i in range(L):
+        for j in range(L):
+            p=np.random.uniform(0,1)
+            if (p<0.5):
+                ising[i][j]=-1
+            else:
+                ising[i][j]=1
+    return (ising)
+
diff --git a/__pycache__/InitEnergy.cpython-38.pyc b/__pycache__/InitEnergy.cpython-38.pyc
new file mode 100644
index 0000000..7f2b0dd
Binary files /dev/null and b/__pycache__/InitEnergy.cpython-38.pyc differ
diff --git a/__pycache__/Montecarlo.cpython-38.pyc b/__pycache__/Montecarlo.cpython-38.pyc
new file mode 100644
index 0000000..c939f6f
Binary files /dev/null and b/__pycache__/Montecarlo.cpython-38.pyc differ
diff --git a/__pycache__/RandomSpin.cpython-38.pyc b/__pycache__/RandomSpin.cpython-38.pyc
new file mode 100644
index 0000000..2e76c11
Binary files /dev/null and b/__pycache__/RandomSpin.cpython-38.pyc differ
diff --git a/montecarlo.py b/montecarlo.py
index 4c3c08b..d7acc74 100644
--- a/montecarlo.py
+++ b/montecarlo.py
@@ -1,78 +1,20 @@
 import numpy as np
 import matplotlib.pyplot as plt
-L=5
-T=1
-ising=np.zeros((L,L),dtype=float)
-
-for i in range(L):
-    for j in range(L):
-        p=np.random.uniform(0,1)
-        if (p<0.5):
-            ising[i][j]=-1
-        else:
-            ising[i][j]=1
-            
+from RandomSpin import randomIsing
+from InitEnergy import Ener
+from Montecarlo import MCS
+L=20
+T=2.1
+ising=randomIsing(L)         
 J=1
-E=0
-M=0
-for i in range (L):
-    for j in range (L):
-        a,b,c,d=i+1,i-1,j+1,j-1
-        if i==L-1:
-            a=0
-        if i==0:
-            b=L-1
-        if j==L-1:
-            c=0
-        if j==0:
-            d=L-1
-    
-        M+=ising[i][j]
-        E=E-J*float(ising[i][j]*(ising[a][j]+ising[b][j]+ising[i][c]+ising[i][d]))
+M,E=Ener(ising,L,J)
 print(M,E)
-
 E*=0.5
-
 N=100000
 ener=[]
 mag=[]
 for n in range(N):
-   
-
-    for k in range(L):
-        for l in range (L):
-            
-            r=np.random.uniform(0,1)
-            i=int(r*(L-1))
-            t=np.random.uniform(0,1)
-            j=int(t*(L-1))
-            a,b,c,d=i+1,i-1,j+1,j-1
-            if i==L-1:
-                a=0
-            if i==0:
-                b=L-1
-            if j==L-1:
-                c=0
-            if j==0:
-                d=L-1
-            Ei=-J*float(ising[i][j]*(ising[a][j]+ising[b][j]+ising[i][c]+ising[i][d]))
-            ising[i][j]=-ising[i][j]
-            Ef=-J*float(ising[i][j]*(ising[a][j]+ising[b][j]+ising[i][c]+ising[i][d]))
-            dE=Ef-Ei
-            
-            if(dE<=0.0):
-                E+=dE
-                M+=(2.0*float(ising[i][j]))
-            else:
-                u=np.exp(-dE/T)
-                h=np.random.uniform(0,1)
-                if(h<u):
-                    E+=dE
-                    M+=(2.0*float(ising[i][j]))
-
-                else:
-                    ising[i][j]=-ising[i][j]
-        
+    ising,E,M=MCS(ising,L,T,J,E,M)
     enspin=E/(L*L)
     Mspin=M/(L*L)
    
@@ -80,7 +22,6 @@
     mag.append(Mspin)
 num=np.arange(0,n+1,1)
 plt.plot(num,mag,linewidth=1)
-# plt.plot(num,ener,linewidth=1)
 plt.show()
 
 
diff --git a/pi.py b/pi.py
index e67e03f..48374be 100644
--- a/pi.py
+++ b/pi.py
@@ -9,7 +9,6 @@
 count=0
 pi_value=[]
 x_value=np.arange(0,iterations,1)
-
 for i in range(iterations):  
     x1=x[i]
     y1=y[i]
diff --git a/test.py b/test.py
deleted file mode 100644
index 55a8194..0000000
--- a/test.py
+++ /dev/null
@@ -1,60 +0,0 @@
-
-import numpy as np
-import matplotlib.pyplot as plt
-t=10000
-h=1.05457*pow(10,-34)
-Kb=1.380649*pow(10,-23)
-w=pow(10,9)
-# h=1
-# Kb=1
-# w=1
-P=[]
-for j in range(1,500):
-    T=j
-
-    def energy(n):
-        y= ((n+(1/2))*h*w)
-        return(y)
-
-    def prob(n):
-        y=np.exp(-energy(n)/(Kb*T))
-        return(y)
-
-
-    n1=0
-    total=[]
-    
-    for i in range (t):
-        n=np.random.randint(0,100)
-        if energy(n)<=energy(n1):
-            total.append(energy(n))
-          
-            n1=n
-        else:
-            k=np.random.uniform(0,1)
-            R=prob(n)/prob(n1)
-            if R>=k:
-                total.append(energy(n))
-                
-                n1=n
-            else:
-                #total.append(energy(n1))
-                
-                n1=n1
-
-
-    # E=total[3000:-1]
-    avgenergy=np.mean(total)
-    P.append(avgenergy) 
-
-temp=np.arange(1,499,1)
-# plt.plot(temp,P)
-# plt.show()
-f=np.copy(P[1:len(P)])
-r=np.copy(P[0:len(P)-1])
-
-d=r-f
-plt.plot(temp,d)
-plt.show()
-    
-
diff --git a/test2.py b/test2.py
deleted file mode 100644
index e2f1828..0000000
--- a/test2.py
+++ /dev/null
@@ -1,43 +0,0 @@
-import numpy as np
-import matplotlib.pyplot as plt
-t=100000
-hbar = 1.0545718e-34  # Planck's constant / 2*pi
-k_B = 1.380649e-23   # Boltzmann constant
-m = 1e-26            # mass of oscillator
-omega = 1e9          # frequency of oscillator
-def energy(n):
-    y= ((n+(1/2))*hbar*omega)
-    return(y)
-
-def prob(n):
-    y=np.exp(-energy(n)/(k_B*300))
-    return(y)
-
-
-n1=0
-total=[]
-ener=0
-for i in range (1,t):
-    n=np.random.randint(0,1000)
-    if energy(n)<=energy(n1):
-        ener+=energy(n)/i
-        total.append(energy(n))
-        
-        n1=n
-    else:
-        k=np.random.uniform(0,1)
-        R=prob(n)/prob(n1)
-        if R>=k:
-            ener+=energy(n)/i
-            total.append(energy(n))
-            
-            n1=n
-        else:
-            ener+=energy(n1)/i
-            total.append(energy(n1))
-            
-            n1=n1
-
-num=np.arange(1,t,1)
-plt.plot(num,total)
-plt.show()
diff --git a/test3.py b/test3.py
deleted file mode 100644
index b9bec43..0000000
--- a/test3.py
+++ /dev/null
@@ -1,49 +0,0 @@
-import numpy as np
-import matplotlib.pyplot as plt
-N=10000
-hbar = 1.0545718e-34  # Planck's constant / 2*pi
-k_B = 1.380649e-23   # Boltzmann constant
-m = 1e-26            # mass of oscillator
-omega = 1e9          # frequency of oscillator
-def energy(n):
-        y= ((n+(1/2))*hbar*omega)
-        return(y)
-tot=[]
-for j in range(1,300):
-    T=j
-    n1=0
-    t=0
-    E=np.zeros(N)
-    for i in range (N):
-        n=np.random.randint(0,100)
-        if energy(n)<=energy(n1):
-            E[i]=energy(n)
-            
-            
-            n1=n
-        else:
-            k=np.random.uniform(0,1)
-            R= np.exp( -(n - n1) / k_B*T )
-            if R>=k:
-                E[i]=energy(n)
-                n1=n
-                
-                
-            else:
-                
-                #E[i]=energy(n)
-                n1=n1
-
-    total=np.mean(E)
-    tot.append(total)
-
-temp=np.arange(1,300,1)
-plt.plot(temp,tot)
-plt.show()
-
-
-
-
-
-   
-   
\ No newline at end of file