diff --git a/Dipolar-Gas-Simulator/+Scripts/run_locally.m b/Dipolar-Gas-Simulator/+Scripts/run_locally.m
index 6dad06f..1353782 100644
--- a/Dipolar-Gas-Simulator/+Scripts/run_locally.m
+++ b/Dipolar-Gas-Simulator/+Scripts/run_locally.m
@@ -22,9 +22,9 @@ OptionsStruct.Dimensions              = [18, 18, 18];
 OptionsStruct.IncludeDDICutOff        = true;
 OptionsStruct.CutoffType              = 'Cylindrical';
 OptionsStruct.SimulationMode          = 'ImaginaryTimeEvolution'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution'
-OptionsStruct.TimeStepSize            = 0.001;           % in s
-OptionsStruct.MinimumTimeStepSize     = 1E-5;            % in s
-OptionsStruct.TimeCutOff              = 1E6;             % in s
+OptionsStruct.TimeStepSize            = 0.002;           % in s
+OptionsStruct.MinimumTimeStepSize     = 1E-6;            % in s
+OptionsStruct.TimeCutOff              = 2E6;             % in s
 OptionsStruct.EnergyTolerance         = 5E-10;
 OptionsStruct.ResidualTolerance       = 1E-05;
 OptionsStruct.NoiseScaleFactor        = 0.01;
diff --git a/Dipolar-Gas-Simulator/+Simulator/@Calculator/calculateOrderParameter.m b/Dipolar-Gas-Simulator/+Simulator/@Calculator/calculateOrderParameter.m
deleted file mode 100644
index e58fd48..0000000
--- a/Dipolar-Gas-Simulator/+Simulator/@Calculator/calculateOrderParameter.m
+++ /dev/null
@@ -1,58 +0,0 @@
-function [m_Order] = calculateOrderParameter(~,psi,Transf,Params,VDk,V,T,muchem)
-
-    NumRealiz = 100;
-
-    Mx          = numel(Transf.x);
-    My          = numel(Transf.y);
-    Mz          = numel(Transf.z);
-
-    r           = normrnd(0,1,size(psi));
-    theta       = rand(size(psi));
-    noise       = r.*exp(2*pi*1i*theta);
-    
-    KEop        = 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
-    Gamma       = 1-1i*Params.gamma_S;
-    dt          = Params.dt;
-
-    avgpsi1     = 0;
-    avgpsi2     = 0;
-
-    for jj = 1:NumRealiz
-        % Generate initial state
-        xi      = sqrt(2*Params.gamma_S*Params.kbol*T*10^(-9)*dt/(Params.hbar*Params.w0*Transf.dx*Transf.dy*Transf.dz));
-        swapx   = randi(length(Transf.x),1,length(Transf.x));
-        swapy   = randi(length(Transf.y),1,length(Transf.y));
-        swapz   = randi(length(Transf.z),1,length(Transf.z));
-        psi_j   = psi + xi * noise(swapx,swapy,swapz);
-        
-        % --- % propagate forward in time 1 time step:
-        % Kin
-        psi_j   = fftn(psi_j);
-        psi_j   = psi_j.*exp(-0.5*1i*Gamma*dt*KEop);
-        psi_j   = ifftn(psi_j);
-    
-        % DDI
-        frho    = fftn(abs(psi_j).^2);
-        Phi     = real(ifftn(frho.*VDk));
-    
-        % Real-space
-        psi_j   = psi_j.*exp(-1i*Gamma*dt*(V + Params.gs*abs(psi_j).^2 + Params.gammaQF*abs(psi_j).^3 + Params.gdd*Phi - muchem));
-
-        % Kin
-        psi_j   = fftn(psi_j);
-        psi_j   = psi_j.*exp(-0.5*1i*Gamma*dt*KEop);
-        psi_j   = ifftn(psi_j);
-    
-        % Projection
-        kcut    = sqrt(2*Params.e_cut);
-        K       = (Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2)<kcut.^2;
-        psi_j   = ifftn(K.*fftn(psi_j));
-    
-        % --- %
-        avgpsi1 = avgpsi1 + abs(sum(psi_j(:)))   /NumRealiz;
-        avgpsi2 = avgpsi2 + sum(abs(psi_j(:)).^2)/NumRealiz;
-
-    end
-
-    m_Order = 1/sqrt(Mx*My*Mz)*avgpsi1/sqrt(avgpsi2);
-end
\ No newline at end of file
diff --git a/Dipolar-Gas-Simulator/+Simulator/@Calculator/calculatePhaseCoherence.m b/Dipolar-Gas-Simulator/+Simulator/@Calculator/calculatePhaseCoherence.m
deleted file mode 100644
index 22a6411..0000000
--- a/Dipolar-Gas-Simulator/+Simulator/@Calculator/calculatePhaseCoherence.m
+++ /dev/null
@@ -1,21 +0,0 @@
-function [PhaseC] = calculatePhaseCoherence(~,psi,Transf,Params)
-
-    norm             = sum(sum(sum(abs(psi).^2,1),2),3)*Transf.dx*Transf.dy*Transf.dz;
-    psi              = psi/sqrt(norm);
-    
-    NumGlobalShifts  = 800;
-    betas            = []; 
-    phishift         = [];
-
-    for jj = 1:NumGlobalShifts
-        phishift(jj) = -pi + 2*pi*(jj-1)/NumGlobalShifts;
-        betas(jj)    = sum(sum(sum(abs(angle(psi*exp(-1i*phishift(jj)))).*abs(psi).^2)));
-    end
-    [minbeta,minidx] = min(betas);
-    
-    psi              = psi*exp(-1i*phishift(minidx));
-    phi              = angle(psi);
-    
-    avgphi           = sum(sum(sum(phi.*abs(psi).^2,1),2),3)*Transf.dx*Transf.dy*Transf.dz;
-    PhaseC           = sum(sum(sum(abs(angle(psi)-avgphi).*abs(psi).^2)))*Transf.dx*Transf.dy*Transf.dz;
-end
\ No newline at end of file
diff --git a/Dipolar-Gas-Simulator/+Simulator/@Calculator/calculateVDk.m b/Dipolar-Gas-Simulator/+Simulator/@Calculator/calculateVDk.m
index 7a56837..266a6ad 100644
--- a/Dipolar-Gas-Simulator/+Simulator/@Calculator/calculateVDk.m
+++ b/Dipolar-Gas-Simulator/+Simulator/@Calculator/calculateVDk.m
@@ -8,6 +8,7 @@ function VDk = calculateVDk(this,Params,Transf,TransfRad,IncludeDDICutOff)
 % Cylindrical (semianalytic)
 % Cylindrical infinite Z, polarized along x (analytic)
 % Spherical
+% Custom Cylindrical
 
     if IncludeDDICutOff
         switch this.CutOffType
@@ -81,7 +82,4 @@ function VDk = calculateVDk(this,Params,Transf,TransfRad,IncludeDDICutOff)
         alph(1) = pi/2;
         VDk     = cos(alph).^2-1/3;
     end
-
-    VDk     = 3 * VDk;
-
 end
\ No newline at end of file
diff --git a/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/propagateWavefunction.m b/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/propagateWavefunction.m
index d1ee687..87198bf 100644
--- a/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/propagateWavefunction.m
+++ b/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/propagateWavefunction.m
@@ -1,9 +1,8 @@
 function [psi] = propagateWavefunction(this,psi,Params,Transf,VDk,V,t_idx,Observ)
-    set(0,'defaulttextInterpreter','latex')
-    set(groot, 'defaultAxesTickLabelInterpreter','latex'); set(groot, 'defaultLegendInterpreter','latex');
     
     switch this.SimulationMode
 	    case 'ImaginaryTimeEvolution'
+            
             dt              = -1j*abs(this.TimeStepSize); % Imaginary Time
             
             KEop            = 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
@@ -11,39 +10,40 @@ function [psi] = propagateWavefunction(this,psi,Params,Transf,VDk,V,t_idx,Observ
             Observ.res      = 1; 
             
             muchem          = this.Calculator.calculateChemicalPotential(psi,Params,Transf,VDk,V);
-            AdaptIdx        = 0;
             
             if this.PlotLive
                 Plotter.plotLive(psi,Params,Transf,Observ)
                 drawnow
             end
-
+            
+            AdaptIdx        = 0;
+            
             while t_idx < Params.sim_time_cut_off
-
+            
                 % kin
-                psi         = fftn(psi);
-                psi         = psi.*exp(-0.5*1i*dt*KEop);
-                psi         = ifftn(psi);
+                psi    = fftn(psi);
+                psi    = psi.*exp(-0.5*1i*dt*KEop);
+                psi    = ifftn(psi);
                 
                 % DDI
-                frho        = fftn(abs(psi).^2);
-                Phi         = real(ifftn(frho.*VDk));
+                frho   = fftn(abs(psi).^2);
+                Phi    = real(ifftn(frho.*VDk));
                 
                 % Real-space
-                psi         = psi.*exp(-1i*dt*(V + Params.gs*abs(psi).^2 + Params.gdd*Phi + Params.gammaQF*abs(psi).^3 - muchem));
+                psi    = psi.*exp(-1i*dt*(V + Params.gs*abs(psi).^2 + Params.gammaQF*abs(psi).^3 + Params.gdd*Phi - muchem));
             
                 % kin
-                psi         = fftn(psi);
-                psi         = psi.*exp(-0.5*1i*dt*KEop);
-                psi         = ifftn(psi);
+                psi    = fftn(psi);
+                psi    = psi.*exp(-0.5*1i*dt*KEop);
+                psi    = ifftn(psi);
                
                 % Renorm
-                Norm        = sum(abs(psi(:)).^2)*Transf.dx*Transf.dy*Transf.dz;
-                psi         = sqrt(Params.N)*psi/sqrt(Norm);
+                Norm   = sum(abs(psi(:)).^2)*Transf.dx*Transf.dy*Transf.dz;
+                psi    = sqrt(Params.N)*psi/sqrt(Norm);
             
-                muchem      = this.Calculator.calculateChemicalPotential(psi,Params,Transf,VDk,V);
-            
-                if mod(t_idx,500) == 0        
+                muchem = this.Calculator.calculateChemicalPotential(psi,Params,Transf,VDk,V);
+                
+                if mod(t_idx,1000) == 0        
             
                     % Change in Energy
                     E               = this.Calculator.calculateTotalEnergy(psi,Params,Transf,VDk,V);
@@ -93,10 +93,6 @@ function [psi] = propagateWavefunction(this,psi,Params,Transf,VDk,V,t_idx,Observ
             E              = E/Norm;
             Observ.EVec    = [Observ.EVec E];
             
-            % Phase coherence
-            [PhaseC]       = this.Calculator.calculatePhaseCoherence(psi,Transf,Params);
-            Observ.PCVec   = [Observ.PCVec PhaseC];
-            
             Observ.res_idx = Observ.res_idx + 1;
 
             disp('Run completed!');
@@ -145,10 +141,6 @@ function [psi] = propagateWavefunction(this,psi,Params,Transf,VDk,V,t_idx,Observ
                     E = E/Norm;
                     Observ.EVec = [Observ.EVec E];
             
-                    % Phase coherence
-                    [PhaseC] = this.Calculator.calculatePhaseCoherence(psi,Transf);
-                    Observ.PCVec = [Observ.PCVec PhaseC];
-            
                     Observ.tVecPlot = [Observ.tVecPlot tVal];
                     Observ.res_idx = Observ.res_idx + 1;
                     
diff --git a/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/setupParameters.m b/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/setupParameters.m
index 48b004d..e5f207f 100644
--- a/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/setupParameters.m
+++ b/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/setupParameters.m
@@ -63,6 +63,7 @@ function [Params] = setupParameters(this)
     
     % DDI strength
     Params.gdd     = 4*pi*Params.add/l0;
+    Params.gdd     = 3 * Params.gdd; % <------------- Necessary for the 3-D case
     
     % Trap gamma
     Params.gx      = (Params.wx/w0)^2;
diff --git a/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/setupSpaceRadial.m b/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/setupSpaceRadial.m
index 2686b0e..be4453b 100644
--- a/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/setupSpaceRadial.m
+++ b/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/setupSpaceRadial.m
@@ -3,31 +3,31 @@ function [Transf] = setupSpaceRadial(~,Params,morder)
     Params.Lr   = 0.5*min(Params.Lx,Params.Ly);
     Params.Nr   = max(Params.Nx,Params.Ny);
     
-    Zmax = 0.5*Params.Lz;
-    Rmax = Params.Lr;
-    Nz = Params.Nz;
-    Nr = Params.Nr;
+    Zmax        = 0.5*Params.Lz;
+    Rmax        = Params.Lr;
+    Nz          = Params.Nz;
+    Nr          = Params.Nr;
     
     if(nargin==2)
         morder=0; %only do Bessel J0
     end
     
     % Fourier grids
-    z=linspace(-Zmax,Zmax,Nz+1);
-    z=z(1:end-1);
-    dz=z(2)-z(1);
-    Kmax=Nz*2*pi/(4*Zmax);
-    kz=linspace(-Kmax,Kmax,Nz+1);
-    kz=kz(1:end-1);
+    z    = linspace(-Zmax,Zmax,Nz+1);
+    z    = z(1:end-1);
+    dz   = z(2)-z(1);
+    Kmax = Nz*2*pi/(4*Zmax);
+    kz   = linspace(-Kmax,Kmax,Nz+1);
+    kz   = kz(1:end-1);
     
     % Hankel grids and transform 
-    H = hankelmatrix(morder,Rmax,Nr);
-    r=H.r(:);
-    kr=H.kr(:);
-    T = diag(H.J/H.kmax)*H.T*diag(Rmax./H.J)*dz*(2*pi);
+    H    = hankelmatrix(morder,Rmax,Nr);
+    r    = H.r(:);
+    kr   = H.kr(:);
+    T    = diag(H.J/H.kmax)*H.T*diag(Rmax./H.J)*dz*(2*pi);
     Tinv = diag(H.J./Rmax)*H.T'*diag(H.kmax./H.J)/dz/(2*pi);
-    wr=H.wr;
-    wk=H.wk;
+    wr   = H.wr;
+    wk   = H.wk;
     % H.T'*diag(H.J/H.vmax)*H.T*diag(Rmax./H.J)
     
     [Transf.R,Transf.Z]=ndgrid(r,z);
diff --git a/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/setupWavefunction.m b/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/setupWavefunction.m
index ce42bb5..973e017 100644
--- a/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/setupWavefunction.m
+++ b/Dipolar-Gas-Simulator/+Simulator/@DipolarGas/setupWavefunction.m
@@ -8,9 +8,9 @@ function [psi] = setupWavefunction(~,Params,Transf)
     elly = sqrt(Params.hbar/(Params.m*Params.wy))/Params.l0;
     ellz = sqrt(Params.hbar/(Params.m*Params.wz))/Params.l0; 
     
-    Rx   = 4*ellx; 
-    Ry   = 4*elly; 
-    Rz   = 4*ellz;
+    Rx   = 8*ellx; 
+    Ry   = 8*elly; 
+    Rz   = 8*ellz;
     X0   = 0.0*Transf.Xmax; 
     Y0   = 0.0*Transf.Ymax; 
     Z0   = 0*Transf.Zmax;
@@ -28,5 +28,5 @@ function [psi] = setupWavefunction(~,Params,Transf)
     % Renormalize wavefunction
     Norm  = sum(abs(psi(:)).^2)*Transf.dx*Transf.dy*Transf.dz;
     psi   = sqrt(Params.N)*psi/sqrt(Norm);
-
+    
 end
\ No newline at end of file
diff --git a/Dipolar-Gas-Simulator/+Simulator/@Potentials/setupParameters.m b/Dipolar-Gas-Simulator/+Simulator/@Potentials/setupParameters.m
index c574d9b..90bafc9 100644
--- a/Dipolar-Gas-Simulator/+Simulator/@Potentials/setupParameters.m
+++ b/Dipolar-Gas-Simulator/+Simulator/@Potentials/setupParameters.m
@@ -1,7 +1,7 @@
 function [Params] = setupParameters(this)
     CONSTANTS      = Helper.PhysicsConstants;
     hbar           = CONSTANTS.PlanckConstantReduced; % [J.s]
-    w0             = 2*pi*100;                        % Angular frequency unit [s^-1]
+    w0             = 2*pi*61.658214297935530;         % Angular frequency unit [s^-1]
     
     % Mass, length scale
     Params.m       = CONSTANTS.Dy164Mass;