diff --git a/epod.m b/epod.m
index 0b615a53267f6ce6ad4f67802d6a3f1e2054c3d0..5bfee460b3dde79e1df1a4282cab11087f1bf5ed 100644
--- a/epod.m
+++ b/epod.m
@@ -24,6 +24,8 @@ end
 % Convert to single to conserve memory
 g_1 = single(g_1);
 g_2 = single(g_2);
+
+% Make NaNs zeros
  
 % Compute the temporal correlation matrices
 C_g1 = g_1*g_1';
@@ -51,28 +53,24 @@ Sig_inv1 = diag(1./diag(Sig1));
 Sig_inv2 = diag(1./diag(Sig2));
 
 % Compute the random coefficients
-a1 = psi1'*g_1;
-a2 = psi2'*g_2;
- 
-% Obtain the spatial modes
-phi1 = Sig_inv1*a1;
-phi2 = Sig_inv2*a2;
+a1 = psi1*Sig1;
+a2 = psi2*Sig2;
 
-% Calculate the extended coefficients
-b_1_2 = psi1'*g_2;
-b_2_1 = psi2'*g_1;
+% Calculate the extended modes (unscaled)
+P_1_2 = a1*Sig_inv1*g_2;
+P_2_1 = a2*Sig_inv2*g_1;
  
 % Obtain extended singular values
-Sig_1_2 = diag(sqrt(sum(b_1_2.^2,2)));
-Sig_2_1 = diag(sqrt(sum(b_2_1.^2,2)));
+Sig_1_2 = sqrt(diag(P_1_2*P_1_2'));
+Sig_2_1 = sqrt(diag(P_2_1*P_2_1'));
  
 % Obtain psuedo-inverse
 Sig_inv_1_2 = diag(1./diag(Sig_1_2));
 Sig_inv_2_1 = diag(1./diag(Sig_2_1));
- 
-% Obtain extended spatial modes
-phi_1_2 = Sig_inv_1_2*b_1_2;
-phi_2_1 = Sig_inv_2_1*b_2_1;
+
+% Get extended modes (scaled)
+phi_1_2 = Sig_inv_1_2*P_1_2;
+phi_2_1 = Sig_inv_2_1*P_2_1;
 
 % Obtain temporal correlations
 nt = size(g_1,1);