diff --git a/code/fwd_CEM_eltved_christensen/CEM_forward_FEniCS.py b/code/fwd_CEM_eltved_christensen/CEM_forward_FEniCS.py
new file mode 100644
index 0000000..582da71
--- /dev/null
+++ b/code/fwd_CEM_eltved_christensen/CEM_forward_FEniCS.py
@@ -0,0 +1,183 @@
+# %%
+"""
+Note by Amal Alghamdi: This code is copied from the project report: Depth 
+Dependency in Electrical Impedance Tomography with the Complete 
+Electrode Model by Anders Eltved and Nikolaj Vestbjerg Christensen (Appendix D.5). Some 
+modifications are made
+
+Created on Wed May 13 0 8 : 1 8 : 4 7 2015
+
+@author : ec0di
+"""
+from dolfin import *
+from matplotlib import pyplot as plt
+from utils import *
+import scipy.io as io
+
+# %% set up data
+case_name = 'case1'  # 'case1' , case3', 'case4', 'case_ref'
+KTC23_dir = './KTC23_data/'
+
+Imatr = io.loadmat(KTC23_dir+"ref.mat")["Injref"]
+
+if case_name == 'case1':
+    phantom_file_name = KTC23_dir+"true1.mat"
+    phantom = io.loadmat(phantom_file_name)["truth"]
+    Uel_ref = io.loadmat(KTC23_dir+"data1.mat")["Uel"].flatten()
+
+elif case_name == 'case2':
+    phantom_file_name = KTC23_dir+"true2.mat"
+    phantom = io.loadmat(phantom_file_name)["truth"]
+    Uel_ref = io.loadmat(KTC23_dir+"data2.mat")["Uel"].flatten()
+
+elif case_name == 'case3':
+    phantom_file_name = KTC23_dir+"true3.mat"
+    phantom = io.loadmat(phantom_file_name)["truth"]
+    Uel_ref = io.loadmat(KTC23_dir+"data3.mat")["Uel"].flatten()
+
+elif case_name == 'case4':
+    phantom_file_name = KTC23_dir+"true4.mat"
+    phantom = io.loadmat(phantom_file_name)["truth"]
+    Uel_ref = io.loadmat(KTC23_dir+"data4.mat")["Uel"].flatten()
+
+elif case_name == 'case_ref':
+    phantom_file_name = KTC23_dir+"true1.mat"
+    phantom = io.loadmat(phantom_file_name)["truth"]
+    phantom[:] = 0
+    Uel_ref = io.loadmat(KTC23_dir+"ref.mat")["Uelref"].flatten()
+
+else:
+    raise Exception("unknown case")
+
+# %% build eit-fenics model
+L = 32
+F = 50
+myeit = EITFenics(L, F)
+num_inj_tested = 76
+z = 1e-6 
+Uel_sim, Q, q_list = myeit.solve_forward(phantom, Imatr, num_inj_tested, z)
+
+#%%
+H = FunctionSpace(myeit.mesh, 'CG', 1)
+inc_pert_1 = interpolate(myeit.inclusion, H)
+sigma_perturb = Function(H)
+sigma_perturb.vector().set_local(inc_pert_1.vector().get_local() - 0.8)
+
+
+
+#project(myeit.inclusion - 0.8, myeit.V[L])
+#%%
+print("Solving P")
+sigma_background = Function(H)
+sigma_background.vector().set_local(inc_pert_1.vector().get_local()*0 + 0.8)
+
+w_list = myeit.solve_P(q_list, sigma_perturb, sigma_background, z)
+
+#%%
+# Solve forward for background phantom A
+phantomA = np.copy(phantom)
+phantomA[:] = 0.8
+Uel_sim_A, Q_A, q_list_A = myeit.solve_forward(phantomA, Imatr, num_inj_tested, z)
+
+#%%
+# Solve forward for background phantom AC
+phantomAC = np.copy(phantom)
+Uel_sim_AC, Q_AC, q_list_AC = myeit.solve_forward(phantomAC, Imatr, num_inj_tested, z)
+
+
+
+# %% Plot potential solution for each injection pattern
+H1 = FunctionSpace(myeit.mesh, 'CG', 1)
+h_func = Function(H1)
+for i, w in enumerate(w_list):
+     plt.figure()
+     h_func.vector().set_local(w.vector().get_local()[L:-1])
+     # h_func.vector().get_local()
+     im = plot(h_func)
+     plt.colorbar(im)
+#     plt.savefig(case_name+"_q_L_"+str(L)+"_inj_"+str(i)+".png")
+
+#%%
+
+mismatch_norm_list = []
+rel_error_list = []
+plot_diffs = False
+for j, q in enumerate(q_list):
+    # Compute the difference between the rhs and lhs in the paper, right before eq 5.1
+    lhs_diff = q_list_A[j].vector().get_local()[L:-1] - q_list_AC[j].vector().get_local()[L:-1]
+    
+    rhs_ = w_list[j].vector().get_local()[L:-1]
+    
+
+
+    rhs_diff_func = Function(H)
+    rhs_diff_func.vector().set_local(rhs_)
+    if plot_diffs:
+       plt.figure()
+       im = plot(rhs_diff_func)
+       plt.colorbar(im)
+       
+
+    lhs_diff_func = Function(H)
+    lhs_diff_func.vector().set_local(-lhs_diff)
+    if plot_diffs:       
+       plt.figure()
+       im = plot(lhs_diff_func)
+       plt.colorbar(im)
+       
+
+    mismatch_func = Function(H)
+    mismatch_func.vector().set_local(rhs_-(-lhs_diff))
+    if plot_diffs:     
+       plt.figure()
+       im = plot(mismatch_func)
+       plt.colorbar(im)
+       
+    
+    print(norm(lhs_diff_func))
+    print(norm(rhs_diff_func))
+    print(norm(mismatch_func))
+    mismatch_norm_list.append(norm(mismatch_func))
+    rel_error_list.append(norm(mismatch_func)/norm(rhs_diff_func))
+
+
+# %% postprocess
+plt.figure()
+plt.plot(Q.flatten(order='F'), label="Q")
+plt.legend()
+img_title = case_name+"_Q_L_"+str(L) +"_F_"+str(F)
+plt.title(img_title)
+plt.savefig(img_title+".png")
+
+plt.figure()
+plt.plot(Uel_sim[:num_inj_tested*31], label="Uel_sim")
+plt.legend()
+img_title = case_name+"_Uel_sim_L_"+str(L) +"_F_"+str(F)
+plt.title(img_title)
+plt.savefig(img_title+".png")
+
+plt.figure()
+plt.plot(Uel_ref[:num_inj_tested*31], label="Uel_ref")
+plt.legend()
+img_title = case_name+"_Uel_ref_L_"+str(L) +"_F_"+str(F)
+plt.title(img_title)
+plt.savefig(img_title+".png")
+
+plt.figure()
+plt.plot(Uel_ref[:num_inj_tested*31] -
+         Uel_sim[:num_inj_tested*31], label="Uel_ref - Uel_sim")
+plt.legend()
+img_title = case_name+"_error_L_"+str(L) +"_F_"+str(F)
+plt.title(img_title)
+plt.savefig(img_title+".png")
+
+plt.figure()
+plt.plot(Uel_sim[:num_inj_tested*31], label='U_sim')
+plt.plot(Uel_ref[:num_inj_tested*31] -
+         Uel_sim[:num_inj_tested*31], label='Uel_ref - Uel_sim')
+plt.legend()
+img_title = case_name+"_sim_and_error_L_"+str(L) +"_F_"+str(F)
+plt.title(img_title)
+plt.savefig(img_title+".png")
+
+# %%
diff --git a/code/fwd_CEM_eltved_christensen/KTC23_data/Mesh_sparse.mat b/code/fwd_CEM_eltved_christensen/KTC23_data/Mesh_sparse.mat
new file mode 100644
index 0000000..0e9b85f
Binary files /dev/null and b/code/fwd_CEM_eltved_christensen/KTC23_data/Mesh_sparse.mat differ
diff --git a/code/fwd_CEM_eltved_christensen/KTC23_data/data1.mat b/code/fwd_CEM_eltved_christensen/KTC23_data/data1.mat
new file mode 100644
index 0000000..734e37f
Binary files /dev/null and b/code/fwd_CEM_eltved_christensen/KTC23_data/data1.mat differ
diff --git a/code/fwd_CEM_eltved_christensen/KTC23_data/data2.mat b/code/fwd_CEM_eltved_christensen/KTC23_data/data2.mat
new file mode 100644
index 0000000..97e19ca
Binary files /dev/null and b/code/fwd_CEM_eltved_christensen/KTC23_data/data2.mat differ
diff --git a/code/fwd_CEM_eltved_christensen/KTC23_data/data3.mat b/code/fwd_CEM_eltved_christensen/KTC23_data/data3.mat
new file mode 100644
index 0000000..ffc2ef5
Binary files /dev/null and b/code/fwd_CEM_eltved_christensen/KTC23_data/data3.mat differ
diff --git a/code/fwd_CEM_eltved_christensen/KTC23_data/data4.mat b/code/fwd_CEM_eltved_christensen/KTC23_data/data4.mat
new file mode 100644
index 0000000..f3d0385
Binary files /dev/null and b/code/fwd_CEM_eltved_christensen/KTC23_data/data4.mat differ
diff --git a/code/fwd_CEM_eltved_christensen/KTC23_data/ref.mat b/code/fwd_CEM_eltved_christensen/KTC23_data/ref.mat
new file mode 100644
index 0000000..7124008
Binary files /dev/null and b/code/fwd_CEM_eltved_christensen/KTC23_data/ref.mat differ
diff --git a/code/fwd_CEM_eltved_christensen/KTC23_data/true1.mat b/code/fwd_CEM_eltved_christensen/KTC23_data/true1.mat
new file mode 100644
index 0000000..572be02
Binary files /dev/null and b/code/fwd_CEM_eltved_christensen/KTC23_data/true1.mat differ
diff --git a/code/fwd_CEM_eltved_christensen/KTC23_data/true2.mat b/code/fwd_CEM_eltved_christensen/KTC23_data/true2.mat
new file mode 100644
index 0000000..c2ac5dd
Binary files /dev/null and b/code/fwd_CEM_eltved_christensen/KTC23_data/true2.mat differ
diff --git a/code/fwd_CEM_eltved_christensen/KTC23_data/true3.mat b/code/fwd_CEM_eltved_christensen/KTC23_data/true3.mat
new file mode 100644
index 0000000..0dc2137
Binary files /dev/null and b/code/fwd_CEM_eltved_christensen/KTC23_data/true3.mat differ
diff --git a/code/fwd_CEM_eltved_christensen/KTC23_data/true4.mat b/code/fwd_CEM_eltved_christensen/KTC23_data/true4.mat
new file mode 100644
index 0000000..a1d5a55
Binary files /dev/null and b/code/fwd_CEM_eltved_christensen/KTC23_data/true4.mat differ
diff --git a/code/fwd_CEM_eltved_christensen/squaredomainCEM.py b/code/fwd_CEM_eltved_christensen/squaredomainCEM.py
new file mode 100644
index 0000000..dd7787b
--- /dev/null
+++ b/code/fwd_CEM_eltved_christensen/squaredomainCEM.py
@@ -0,0 +1,86 @@
+
+#%% 
+"""
+Note by Amal Alghamdi: This code is copied from the project report: Depth 
+Dependency in Electrical Impedance Tomography with the Complete 
+Electrode Model by Anders Eltved and Nikolaj Vestbjerg Christensen (Appendix D.5). Some 
+modifications are made to make it run against FEniCS 2016.2.
+
+Created on Mon Mar 9 05:17:48 2015
+@author : ec0di
+"""
+import numpy as np 
+from dolfin import *
+from mshr import *
+import dolfin as dl
+import matplotlib.pyplot as plt
+
+Z=0.1; 
+I=[2,-2]
+x0 =0; y0 =0; x1 =1; y1=1 
+L=2
+mesh = RectangleMesh(Point(x0 , y0) , Point(x1 , y1), 10, 10) 
+
+class sigma_fun(Expression):
+    def eval(self ,values ,x): 
+        values [0]=1
+
+sigma = sigma_fun(degree=2)
+
+class Left (SubDomain) :
+    def inside(self , x, on_boundary): 
+        return near(x[0] , 0.0)
+class Right(SubDomain):
+    def inside(self , x, on_boundary): 
+        return near(x[0] , 1.0)
+
+left = Left() 
+right = Right()
+
+boundaries = FacetFunction("size_t", mesh) 
+boundaries.set_all (0)
+left.mark(boundaries , 1)
+right.mark(boundaries , 2)
+
+dS = Measure('ds', domain=mesh)[boundaries]
+
+R = FunctionSpace (mesh , "R" ,0)
+H1 = FunctionSpace(mesh,"CG",1)
+
+#R = FiniteElement("R", mesh.ufl_cell(), 0)
+#H1 = FiniteElement("CG", mesh.ufl_cell(), 1)
+
+mixedspaces=R.ufl_element()*R.ufl_element()*H1.ufl_element()*R.ufl_element()
+
+V = FunctionSpace(mesh, mixedspaces)
+u = TrialFunction(V)
+v = TestFunction(V)
+
+f = 0*dS(1)
+B = sigma*inner(nabla_grad(u[L]) , nabla_grad(v[L]))*dx
+for i in range(L):
+    B += 1/Z*(u[L]-u[i])*(v[L]-v[i])*dS(i+1)
+    B += (v[L+1]*u[i]/assemble(1*dS(i+1)))*dS(i+1) 
+    B += (u[L+1]*v[i]/assemble(1*dS(i+1)))*dS(i+1)
+    f += (I[i]*v[i]/assemble(1*dS(i+1)))*dS(i+1)
+          
+q = Function(V)
+solve(B==f, q)
+
+Q=q.split(deepcopy=True)
+U=[]
+
+for i in range(2):
+    U.append(Q[ i ])
+
+u = Q[1]#Q[2]
+
+for i in range(2):
+    print("U"+str(i+1)+": "+str(U[i].compute_vertex_values()[0]))
+
+#%%
+im = plot(q[2])
+plt.colorbar(im)
+interactive()
+plt.show()
+plt.savefig("squaredomainCEM.png")
diff --git a/code/fwd_CEM_eltved_christensen/utils.py b/code/fwd_CEM_eltved_christensen/utils.py
new file mode 100644
index 0000000..ff91584
--- /dev/null
+++ b/code/fwd_CEM_eltved_christensen/utils.py
@@ -0,0 +1,260 @@
+"""
+Note by Amal Alghamdi: This code is copied from the project report: Depth 
+Dependency in Electrical Impedance Tomography with the Complete 
+Electrode Model by Anders Eltved and Nikolaj Vestbjerg Christensen (Appendix D.5). Some 
+modifications are made.
+"""
+
+import numpy as np
+from scipy.interpolate import RegularGridInterpolator
+from matplotlib import pyplot as plt
+from dolfin import *
+from mshr import *
+
+class  EITFenics:
+    def __init__(self, L=32, F=50):
+        self.L = L
+        self.F = F
+        self._create_mesh()
+        self._build_subdomains()
+        self.V, self.dS = build_spaces(self.mesh, L, self.subdomains)
+
+    def _create_mesh(self):
+        R = 1  # radius of circle
+        n = 300 # number of polygons to approximate circle
+        self.mesh = generate_mesh(Circle(Point(0, 0), R, n), self.F)  # generate mesh
+        print("Mesh created with %d elements" % self.mesh.num_entities(2))
+    def _build_subdomains(self):
+        L = self.L
+        e_l = np.pi / L
+        d_e = 2*np.pi / L - e_l
+
+        # Define subdomain mesh
+        self.subdomains = MeshFunction("size_t", self.mesh, self.mesh.topology().dim()-1, 0)
+
+        # Define subdomains
+        def twopiarctan(x):
+            val = np.arctan2(x[1], x[0])
+            if val < 0:
+                val = val+2*np.pi
+            return val
+
+        class e(SubDomain):
+            def inside(self, x, on_boundary):
+                theta = twopiarctan(x)
+                # print "theta inside", theta
+                if theta1 > theta2:
+                    return on_boundary and ((theta >= 0
+                                            and theta <= theta2) or (theta >= theta1
+                                                                    and theta <= 2*np.pi))
+                return on_boundary and theta >= theta1 \
+                    and theta <= theta2
+                # return  theta>=theta1 and theta<=theta2
+
+        for i in range(1, L+1):
+            shift_theta = np.pi/2 - np.pi/(2*L)
+            # print "shift_theta", shift_theta
+            # print L
+            theta1 = np.mod((i - 1) * (e_l+d_e) + shift_theta, 2*np.pi)
+            theta2 = np.mod(theta1+e_l, 2*np.pi)
+            # print i
+            # print theta1
+            # print theta2
+            e1 = e()  # create instance
+            e1 .mark(self.subdomains, i)  # mark subdomains
+            xdmf = XDMFFile("subdomains.xdmf")
+            xdmf.write(self.subdomains)
+    def _create_inclusion(self, phantom):
+        high_conductivity = 1e1
+        low_conductivity = 1e-2
+        background_conductivity = 0.8
+        # Define conductivity
+        phantom_float = np.zeros(phantom.shape)
+        phantom_float[phantom == 0] = background_conductivity
+        phantom_float[np.isclose(phantom, 1, rtol=0.01)] = low_conductivity
+        phantom_float[phantom == 2] = high_conductivity
+
+        plt.figure()
+        im = plt.imshow(phantom_float)
+        plt.colorbar(im)  # norm= 'log'
+        plt.savefig("phantom.png")
+
+        self.inclusion = Inclusion(phantom_float, degree=0)
+    def solve_forward(self, phantom, injection_patterns, num_inj_tested, z=1e-6):
+        self._create_inclusion(phantom)
+        L = self.L
+
+        # Define vector of contact impedance
+        # z = 10e-6  # 0.1 # Impedence
+        Z = []
+        for i in range(self.L):
+            Z.append(z)
+
+        # Define H1 room
+        H1 = FunctionSpace(self.mesh, 'CG', 1)
+
+        # Loop over current patterns
+        num_inj = 76  # Number of injection pattern
+        # num_inj_tested = 76
+
+        B = build_b(self.inclusion, Z, self.V, self.dS, L)
+        Q = np.zeros((L, num_inj))
+        Diff = np.zeros((L-1, num_inj))
+        q_list = []
+
+        for i in range(num_inj)[:num_inj_tested]:
+            print("injection pattern"+str(i))
+            Q_i, q = solver(injection_patterns[:, i], B, self.V, self.dS, L)
+            q_list.append(q)
+
+            Q[:, i] = Q_i
+            Diff[:, i] = np.diff(Q_i)
+
+        Uel_sim = -Diff.flatten(order='F')
+        return Uel_sim, Q, q_list
+
+    def solve_P(self, y_list, sigma_perturb, background_sigma, z=1e-6):
+        L = self.L
+
+        # Define vector of contact impedance
+        # z = 10e-6  # 0.1 # Impedence
+        Z = []
+        for i in range(self.L):
+            Z.append(z)
+
+        # Define H1 room
+        H1 = FunctionSpace(self.mesh, 'CG', 1)
+
+        B = build_b(background_sigma, Z, self.V, self.dS, L)
+
+        v = TestFunction(self.V)
+
+        w_list = []
+        for y in y_list:
+
+            f = -sigma_perturb * inner(nabla_grad(y[L]), nabla_grad(v[L])) * dx 
+
+            rhs = assemble(f)
+
+            # Compute solution
+            w = Function(self.V)
+            solve(B, w.vector(), rhs)
+            w_list.append(w)
+
+        return w_list
+
+
+class Inclusion(UserExpression):
+    def __init__(self, phantom, **kwargs):
+        super().__init__(**kwargs)
+        x_grid = np.linspace(-1, 1, 256)
+        y_grid = np.linspace(-1, 1, 256)
+        self._interpolater = RegularGridInterpolator(
+            (x_grid, y_grid), phantom, method="nearest")
+
+    def eval(self, values, x):
+        values[0] = self._interpolater([x[0], x[1]])
+
+def build_subdomains(L, mesh):
+    def twopiarctan(x):
+        val = np.arctan2(x[1], x[0])
+        if val < 0:
+            val = val+2*np.pi
+        return val
+
+    e_l = np.pi / L
+    d_e = 2*np.pi / L - e_l
+
+    # Define subdomain mesh
+    subdomains = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
+
+   # Define subdomains
+    class e(SubDomain):
+        def inside(self, x, on_boundary):
+            theta = twopiarctan(x)
+            # print "theta inside", theta
+            if theta1 > theta2:
+                return on_boundary and ((theta >= 0
+                                         and theta <= theta2) or (theta >= theta1
+                                                                  and theta <= 2*np.pi))
+            return on_boundary and theta >= theta1 \
+                and theta <= theta2
+            # return  theta>=theta1 and theta<=theta2
+
+    for i in range(1, L+1):
+        shift_theta = np.pi/2 - np.pi/(2*L)
+        # print "shift_theta", shift_theta
+        # print L
+        theta1 = np.mod((i - 1) * (e_l+d_e) + shift_theta, 2*np.pi)
+        theta2 = np.mod(theta1+e_l, 2*np.pi)
+        # print i
+        # print theta1
+        # print theta2
+        e1 = e()  # create instance
+        e1 .mark(subdomains, i)  # mark subdomain
+
+    return subdomains
+
+def build_spaces(mesh, L, subdomains):
+    R = FunctionSpace(mesh, "R", 0)
+    H1 = FunctionSpace(mesh, "CG", 1)
+
+    spacelist = None
+
+    for i in range(1, L+1):
+
+        if i == 1:
+            spacelist = R.ufl_element()
+        else:
+            spacelist *= R.ufl_element()
+
+    spacelist *= H1.ufl_element()
+    spacelist *= R.ufl_element()
+
+    # Create function space
+    V = FunctionSpace(mesh, spacelist)
+
+    # Define new measures associated with the boundaries
+    dS = Measure('ds', domain=mesh)[subdomains]
+
+    return V, dS
+
+
+def build_b(sigma, Z, V, dS, L):
+
+    # Define trial and test functions
+    u = TrialFunction(V)
+    v = TestFunction(V)
+
+    B = sigma * inner(nabla_grad(u[L]), nabla_grad(v[L])) * dx
+
+    for i in range(L):
+        B += 1/Z[i] * (u[L]-u[i])*(v[L]-v[i]) * dS(i + 1)
+        #TODO: check if this is correct for P operator
+        B += (v[L+1]*u[i] / assemble(1*dS(i+1))) * dS(i+1)
+        B += (u[L+1]*v[i] / assemble(1*dS(i+1))) * dS(i+1)
+
+    return assemble(B)
+
+
+def solver(I, B, V, dS, L):  # sigma ,L, I , Z ,mesh, subdomains )
+   # def 2 pi function
+
+    # Define trial and test functions
+    u = TrialFunction(V)
+    v = TestFunction(V)
+
+    f = 0*dS(1)
+
+    for i in range(L):
+        f += (I[i] * v[i] / assemble(1*dS(i+1))) * dS(i+1)
+
+    rhs = assemble(f)
+
+    # Compute solution
+    q = Function(V)
+    solve(B, q.vector(), rhs)
+
+    Q = q.vector().get_local()[:L]
+
+    return Q, q