diff --git a/Dipolar Gas Simulator/+Helper/ProgressBar.m b/Dipolar Gas Simulator/+Helper/ProgressBar.m
new file mode 100644
index 0000000..e30a4e0
--- /dev/null
+++ b/Dipolar Gas Simulator/+Helper/ProgressBar.m	
@@ -0,0 +1,68 @@
+classdef ProgressBar < handle
+% class for command-line progress-bar notification.
+    properties
+        strPercentageLength;
+        strDotsMaximum;
+    end
+    methods
+        %--- constructor
+        function this = ProgressBar()
+             %% Initialization
+            % Vizualization parameters
+            this.strPercentageLength = 10;   %   Length of percentage string (must be >5)
+            this.strDotsMaximum      = 10;   %   The total number of dots in a progress bar
+        end
+        %--- print method
+        function run(this, msg)
+            % This function creates a text progress bar. It should be called with a 
+            % STRING argument to initialize and terminate. Otherwise the number corresponding 
+            % to progress in % should be supplied.
+            % INPUTS:   C   Either: Text string to initialize or terminate 
+            %                       Percentage number to show progress 
+            % OUTPUTS:  N/A
+            % Example:  Please refer to demo_textprogressbar.m
+            % Author: Paul Proteus (e-mail: proteus.paul (at) yahoo (dot) com)
+            % Version: 1.0
+            % Changes tracker:  29.06.2010  - First version
+            % Inspired by: http://blogs.mathworks.com/loren/2007/08/01/monitoring-progress-of-a-calculation/
+            %% Main 
+            persistent strCR;                %   Carriage return pesistent variable
+            if isempty(strCR) && ~ischar(msg)
+                % Progress bar must be initialized with a string
+                error('The text progress must be initialized with a string!');
+            elseif isempty(strCR) && ischar(msg)
+                % Progress bar - initialization
+                fprintf('%s',msg);
+                strCR = -1;
+            elseif ~isempty(strCR) && ischar(msg)
+                % Progress bar  - termination
+                strCR = [];  
+                fprintf([msg '\n']);
+            elseif isnumeric(msg)
+                % Progress bar - normal progress
+                msg = floor(msg);
+                percentageOut = [num2str(msg) '%%'];
+                percentageOut = [percentageOut repmat(' ',1,this.strPercentageLength-length(percentageOut)-1)];
+                nDots = floor(msg/100*this.strDotsMaximum);
+                dotOut = ['[' repmat('.',1,nDots) repmat(' ',1,this.strDotsMaximum-nDots) ']'];
+                strOut = [percentageOut dotOut];
+                
+                % Print it on the screen
+                if strCR == -1
+                    % Don't do carriage return during first run
+                    fprintf(strOut);
+                else
+                    % Do it during all the other runs
+                    fprintf([strCR strOut]);
+                end
+                
+                % Update carriage return
+                strCR = repmat('\b',1,length(strOut)-1);
+                
+            else
+                % Any other unexpected input
+                error('Unsupported argument type');
+            end
+        end
+    end
+end
\ No newline at end of file
diff --git a/Dipolar Gas Simulator/+Helper/parforNotifications.m b/Dipolar Gas Simulator/+Helper/parforNotifications.m
deleted file mode 100644
index 4ad3af4..0000000
--- a/Dipolar Gas Simulator/+Helper/parforNotifications.m	
+++ /dev/null
@@ -1,148 +0,0 @@
-% Copyright (c) 2019 Andrea Alberti
-%
-% All rights reserved.
-classdef parforNotifications < handle
-    properties
-        N;   % number of iterations
-        text = 'Please wait ...';   % text to show
-        width = 50;
-        showWarning = true;
-    end
-    properties (GetAccess = public, SetAccess = private)
-        n;
-    end
-    properties (Access = private)
-        inProgress = false;
-        percent;
-        DataQueue;
-        usePercent;
-        Nstr;
-        NstrL;
-        lastComment;
-    end
-    methods
-        function this = parforNotifications()
-            this.DataQueue = parallel.pool.DataQueue;
-            afterEach(this.DataQueue, @this.updateStatus);
-        end
-        % Start progress bar
-        function PB_start(this,N,varargin)
-            assert(isscalar(N) && isnumeric(N) && N == floor(N) && N>0, 'Error: ''N'' must be a scalar positive integer.');
-            
-            this.N = N;
-            
-            p = inputParser;
-            addParameter(p,'message','Please wait: ');
-            addParameter(p,'usePercentage',true);
-            
-            parse(p,varargin{:});
-            
-            this.text = p.Results.message;
-            assert(ischar(this.text), 'Error: ''Message'' must be a string.');
-            
-            this.usePercent = p.Results.usePercentage;
-            assert(isscalar(this.usePercent) && islogical(this.usePercent), 'Error: ''usePercentage'' must be a logical scalar.');
-            
-            this.percent = 0;
-            this.n = 0;
-            this.lastComment = '';
-            if this.usePercent
-                fprintf('%s [%s]: %3d%%\n',this.text, char(32*ones(1,this.width)),0);
-            else
-                this.Nstr = sprintf('%d',this.N);
-                this.NstrL = numel(this.Nstr);
-                fprintf('%s [%s]: %s/%s\n',this.text, char(32*ones(1,this.width)),[char(32*ones(1,this.NstrL-1)),'0'],this.Nstr);
-            end
-            
-            this.inProgress = true;
-        end
-        % Iterate progress bar
-        function PB_iterate(this,str)
-            if nargin == 1
-                send(this.DataQueue,'');
-            else
-                send(this.DataQueue,str);
-            end
-        end
-        function warning(this,warn_id,msg)
-            if this.showWarning
-                msg = struct('Action','Warning','Id',warn_id,'Message',msg);
-                send(this.DataQueue,msg);
-            end
-        end
-        function PB_reprint(this)
-            p = round(100*this.n/this.N);
-            
-            this.percent = p;
-            
-            cursor_pos=1+round((this.width-1)*p/100);
-            
-            if p < 100
-                sep_char = '|';
-            else
-                sep_char = '.';
-            end
-            
-            if this.usePercent
-                fprintf('%s [%s%s%s]: %3d%%\n', this.text, char(46*ones(1,cursor_pos-1)), sep_char, char(32*ones(1,this.width-cursor_pos)),p);
-            else
-                nstr=sprintf('%d',this.n);
-                fprintf('%s [%s%s%s]: %s/%s\n', this.text, char(46*ones(1,cursor_pos-1)), sep_char, char(32*ones(1,this.width-cursor_pos)),[char(32*ones(1,this.NstrL-numel(nstr))),nstr],this.Nstr);
-            end
-        end
-        function updateStatus(this,data)
-            
-            if ischar(data)
-                
-                this.n = this.n + 1;
-                
-                p = round(100*this.n/this.N);
-                
-                if p >= this.percent+1 || this.n == this.N
-                    this.percent = p;
-                    
-                    cursor_pos=1+round((this.width-1)*p/100);
-                    
-                    if p < 100
-                        sep_char = '|';
-                    else
-                        sep_char = '.';
-                    end
-                    
-                    if ~isempty(data)
-                        comment = [' (',data,')'];
-                    else
-                        comment = '';
-                    end
-                    
-                    if this.usePercent
-                        fprintf('%s%s%s%s]: %3d%%%s\n',char(8*ones(1,58+numel(this.lastComment))), char(46*ones(1,cursor_pos-1)), sep_char, char(32*ones(1,this.width-cursor_pos)),p,comment);
-                    else
-                        nstr=sprintf('%d',this.n);
-                        fprintf('%s%s%s%s]: %s/%s%s\n',char(8*ones(1,55+2*numel(this.Nstr)+numel(this.lastComment))), char(46*ones(1,cursor_pos-1)), sep_char, char(32*ones(1,this.width-cursor_pos)),[char(32*ones(1,this.NstrL-numel(nstr))),nstr],this.Nstr,comment)
-                    end
-                    
-                    this.lastComment = comment;
-                    
-                    
-                    if p == 100
-                        this.inProgress = false;
-                    end
-                end
-                
-            else
-                switch data.Action
-                    case 'Warning'
-                        warning(data.Id,[data.Message,newline]);
-                        if this.inProgress
-                            this.PB_reprint();
-                        end
-                end
-                
-            end
-            
-        end
-    end
-end
-
-
diff --git a/Dipolar Gas Simulator/+Scripts/test.m b/Dipolar Gas Simulator/+Scripts/test.m
index 835e7a9..a828fe5 100644
--- a/Dipolar Gas Simulator/+Scripts/test.m	
+++ b/Dipolar Gas Simulator/+Scripts/test.m	
@@ -21,7 +21,7 @@ OptionsStruct.Dimensions              = [40, 40, 20];
 OptionsStruct.CutoffType              = 'Cylindrical';
 OptionsStruct.SimulationMode          = 'ImaginaryTimeEvolution'; % 'ImaginaryTimeEvolution' | 'RealTimeEvolution'
 OptionsStruct.TimeStepSize            = 50E-6;                    % in s
-OptionsStruct.NumberOfTimeSteps       = 2E6;                      % in s
+OptionsStruct.NumberOfTimeSteps       = 10;                      % in s
 OptionsStruct.EnergyTolerance         = 5E-10;
 
 OptionsStruct.SaveData                = true;
diff --git a/Dipolar Gas Simulator/+Simulator/@Calculator/Calculator.m b/Dipolar Gas Simulator/+Simulator/@Calculator/Calculator.m
index aae055e..a80c05c 100644
--- a/Dipolar Gas Simulator/+Simulator/@Calculator/Calculator.m	
+++ b/Dipolar Gas Simulator/+Simulator/@Calculator/Calculator.m	
@@ -43,7 +43,7 @@ classdef Calculator < handle & matlab.mixin.Copyable
             this.TotalEnergy                   = this.CalculatorDefaults.TotalEnergy;
             this.CutoffType                    = p.Results.CutoffType;
         end
-        
+
         function restoreDefaults(this)
             this.ChemicalPotential          = this.CalculatorDefaults.ChemicalPotential;
             this.EnergyComponents           = this.CalculatorDefaults.EnergyComponents;
diff --git a/Dipolar Gas Simulator/+Simulator/@Calculator/calculateChemicalPotential.m b/Dipolar Gas Simulator/+Simulator/@Calculator/calculateChemicalPotential.m
deleted file mode 100644
index 96a6a84..0000000
--- a/Dipolar Gas Simulator/+Simulator/@Calculator/calculateChemicalPotential.m	
+++ /dev/null
@@ -1,29 +0,0 @@
-function muchem = calculateChemicalPotential(~,psi,Params,Transf,VDk,V)
-
-    %Parameters
-    normfac = Params.Lx*Params.Ly*Params.Lz/numel(psi);
-    KEop= 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
-    
-    % DDIs
-    frho=fftn(abs(psi).^2);
-    Phi=real(ifftn(frho.*VDk));
-    
-    Eddi = (Params.gdd*Phi.*abs(psi).^2);
-    
-    %Kinetic energy
-    Ekin = KEop.*abs(fftn(psi)*normfac).^2;
-    Ekin = trapz(Ekin(:))*Transf.dkx*Transf.dky*Transf.dkz/(2*pi)^3;
-    
-    %Potential energy
-    Epot = V.*abs(psi).^2;
-    
-    %Contact interactions
-    Eint = Params.gs*abs(psi).^4;
-    
-    %Quantum fluctuations
-    Eqf = Params.gammaQF*abs(psi).^5;
-    
-    %Total energy
-    muchem = Ekin + trapz(Epot(:) + Eint(:) + Eddi(:) + Eqf(:))*Transf.dx*Transf.dy*Transf.dz; %
-    muchem = muchem / Params.N;
-end
\ No newline at end of file
diff --git a/Dipolar Gas Simulator/+Simulator/@Calculator/calculateEnergyComponents.m b/Dipolar Gas Simulator/+Simulator/@Calculator/calculateEnergyComponents.m
deleted file mode 100644
index 92dd32c..0000000
--- a/Dipolar Gas Simulator/+Simulator/@Calculator/calculateEnergyComponents.m	
+++ /dev/null
@@ -1,35 +0,0 @@
-function E = calculateEnergyComponents(~,psi,Params,Transf,VDk,V)
-
-    %Parameters
-    
-    KEop= 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
-    normfac = Params.Lx*Params.Ly*Params.Lz/numel(psi);
-    
-    % DDIs
-    frho = fftn(abs(psi).^2);
-    Phi = real(ifftn(frho.*VDk));
-    
-    Eddi = 0.5*Params.gdd*Phi.*abs(psi).^2;
-    E.Eddi = trapz(Eddi(:))*Transf.dx*Transf.dy*Transf.dz;
-    
-    % EddiTot = trapz(Eddi(:))*Transf.dx*Transf.dy*Transf.dz;
-    
-    %Kinetic energy
-    % psik = ifftshift(fftn(fftshift(psi)))*normfac;
-    
-    Ekin = KEop.*abs(fftn(psi)*normfac).^2;
-    E.Ekin = trapz(Ekin(:))*Transf.dkx*Transf.dky*Transf.dkz/(2*pi)^3;
-    
-    % Potential energy
-    Epot = V.*abs(psi).^2;
-    E.Epot = trapz(Epot(:))*Transf.dx*Transf.dy*Transf.dz;
-    
-    %Contact interactions
-    Eint = 0.5*Params.gs*abs(psi).^4;
-    E.Eint = trapz(Eint(:))*Transf.dx*Transf.dy*Transf.dz;
-    
-    %Quantum fluctuations
-    Eqf = 0.4*Params.gammaQF*abs(psi).^5;
-    E.Eqf = trapz(Eqf(:))*Transf.dx*Transf.dy*Transf.dz;
-    
-end
diff --git a/Dipolar Gas Simulator/+Simulator/@Calculator/calculateNormalizedResiduals.m b/Dipolar Gas Simulator/+Simulator/@Calculator/calculateNormalizedResiduals.m
deleted file mode 100644
index 134e426..0000000
--- a/Dipolar Gas Simulator/+Simulator/@Calculator/calculateNormalizedResiduals.m	
+++ /dev/null
@@ -1,25 +0,0 @@
-function res = calculateNormalizedResiduals(~,psi,Params,Transf,VDk,V,muchem)
-
-    KEop= 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
-    
-    % DDIs
-    frho=fftn(abs(psi).^2);
-    Phi=real(ifftn(frho.*VDk));
-    
-    Eddi = Params.gdd*Phi.*psi;
-    
-    %Kinetic energy
-    Ekin = ifftn(KEop.*fftn(psi));
-    
-    %Potential energy
-    Epot = V.*psi;
-    
-    %Contact interactions
-    Eint = Params.gs*abs(psi).^2.*psi;
-    
-    %Quantum fluctuations
-    Eqf = Params.gammaQF*abs(psi).^3.*psi;
-    
-    %Total energy
-    res =  trapz(abs(Ekin(:) + Epot(:) + Eint(:) + Eddi(:) + Eqf(:) - muchem*psi(:))*Transf.dx*Transf.dy*Transf.dz)/trapz(abs(muchem*psi(:))*Transf.dx*Transf.dy*Transf.dz);
-end
\ No newline at end of file
diff --git a/Dipolar Gas Simulator/+Simulator/@Calculator/calculateNumericalHankelTransform.m b/Dipolar Gas Simulator/+Simulator/@Calculator/calculateNumericalHankelTransform.m
deleted file mode 100644
index b3d96c6..0000000
--- a/Dipolar Gas Simulator/+Simulator/@Calculator/calculateNumericalHankelTransform.m	
+++ /dev/null
@@ -1,39 +0,0 @@
-function VDkSemi = calculateNumericalHankelTransform(~,kr,kz,Rmax,Zmax,Nr)
-
-    % accuracy inputs for numerical integration
-    if(nargin==5)
-        Nr = 5e4;
-    end
-    
-    Nz = 64;
-    farRmultiple = 2000;
-    
-    % midpoint grids for the integration over 0<z<Zmax, Rmax<r<farRmultiple*Rmax (i.e. starts at Rmax)
-    dr=(farRmultiple-1)*Rmax/Nr;
-    r = ((1:Nr)'-0.5)*dr+Rmax;
-    dz=Zmax/Nz;
-    z = ((1:Nz)-0.5)*dz;
-    [R, Z] = ndgrid(r,z);
-    Rsq = R.^2 + Z.^2;
-    
-    % real space interaction to be transformed
-    igrandbase = (1 - 3*Z.^2./Rsq)./Rsq.^(3/2);
-    
-    % do the Hankel/Fourier-Bessel transform numerically
-    % prestore to ensure each besselj is only calculated once
-    % cell is faster than (:,:,krn) slicing
-    Nkr = numel(kr);
-    besselr = cell(Nkr,1);
-    for krn = 1:Nkr
-        besselr{krn} = repmat(r.*besselj(0,kr(krn)*r),1,Nz);
-    end
-    
-    VDkSemi = zeros([Nkr,numel(kz)]);
-    for kzn = 1:numel(kz)
-        igrandbasez = repmat(cos(kz(kzn)*z),Nr,1) .* igrandbase;
-        for krn = 1:Nkr
-            igrand = igrandbasez.*besselr{krn};
-            VDkSemi(krn,kzn) = VDkSemi(krn,kzn) - sum(igrand(:))*dz*dr;
-        end
-    end
-end
\ No newline at end of file
diff --git a/Dipolar Gas Simulator/+Simulator/@Calculator/calculateOrderParameter.m b/Dipolar Gas Simulator/+Simulator/@Calculator/calculateOrderParameter.m
deleted file mode 100644
index 15a28f4..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;
-
-    avgpsi = 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));
-    
-        % --- %
-        avgpsi = avgpsi + abs(sum(psi_j(:)))/NumRealiz;
-        avgpsi2 = avgpsi2 + sum(abs(psi_j(:)).^2)/NumRealiz;
-
-    end
-
-    m_Order = 1/sqrt(Mx*My*Mz)*avgpsi/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 449cf36..0000000
--- a/Dipolar Gas Simulator/+Simulator/@Calculator/calculatePhaseCoherence.m	
+++ /dev/null
@@ -1,19 +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/calculateTotalEnergy.m b/Dipolar Gas Simulator/+Simulator/@Calculator/calculateTotalEnergy.m
deleted file mode 100644
index f53f7a2..0000000
--- a/Dipolar Gas Simulator/+Simulator/@Calculator/calculateTotalEnergy.m	
+++ /dev/null
@@ -1,32 +0,0 @@
-function E = calculateTotalEnergy(~,psi,Params,Transf,VDk,V)
-
-    %Parameters
-    
-    KEop= 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
-    normfac = Params.Lx*Params.Ly*Params.Lz/numel(psi);
-    
-    % DDIs
-    frho = fftn(abs(psi).^2);
-    Phi = real(ifftn(frho.*VDk));
-    
-    Eddi = 0.5*Params.gdd*Phi.*abs(psi).^2;
-    
-    % EddiTot = trapz(Eddi(:))*Transf.dx*Transf.dy*Transf.dz;
-    
-    %Kinetic energy
-    % psik = ifftshift(fftn(fftshift(psi)))*normfac;
-    
-    Ekin = KEop.*abs(fftn(psi)*normfac).^2;
-    Ekin = trapz(Ekin(:))*Transf.dkx*Transf.dky*Transf.dkz/(2*pi)^3;
-    
-    % Potential energy
-    Epot = V.*abs(psi).^2;
-    
-    %Contact interactions
-    Eint = 0.5*Params.gs*abs(psi).^4;
-    
-    %Quantum fluctuations
-    Eqf = 0.4*Params.gammaQF*abs(psi).^5;
-    
-    E = Ekin + trapz(Epot(:) + Eint(:) + Eddi(:) + Eqf(:))*Transf.dx*Transf.dy*Transf.dz;
-end
\ No newline at end of file
diff --git a/Dipolar Gas Simulator/+Simulator/@Calculator/calculateVDCutoff.m b/Dipolar Gas Simulator/+Simulator/@Calculator/calculateVDCutoff.m
deleted file mode 100644
index 3c49b4e..0000000
--- a/Dipolar Gas Simulator/+Simulator/@Calculator/calculateVDCutoff.m	
+++ /dev/null
@@ -1,64 +0,0 @@
-function VDk = calculateVDCutoff(this,Params,Transf,TransfRad)
-% makes the dipolar interaction matrix, size numel(Params.kr) * numel(Params.kz)
-% Rmax and Zmax are the interaction cutoffs
-% VDk needs to be multiplied by Cdd
-% approach is that of Lu, PRA 82, 023622 (2010)
-
-% == Calulating the DDI potential in Fourier space with appropriate cutoff == %
-% Cylindrical (semianalytic)
-% Cylindrical infinite Z, polarized along x (analytic)
-% Spherical
-
-    switch this.CutoffType
-	    case 'Cylindrical' %Cylindrical (semianalytic)
-            Zcutoff = Params.Lz/2;            
-            alph    = acos((Transf.KX*sin(Params.theta)*cos(Params.phi)+Transf.KY*sin(Params.theta)*sin(Params.phi)+Transf.KZ*cos(Params.theta))./sqrt(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2));
-            alph(1) = pi/2;
-            
-            % Analytic part of cutoff for slice 0<z<Zmax, 0<r<Inf Ronen, PRL 98, 030406 (2007)
-            cossq = cos(alph).^2;
-            VDk   = cossq-1/3;
-            sinsq = 1 - cossq;
-            VDk   = VDk + exp(-Zcutoff*sqrt(Transf.KX.^2+Transf.KY.^2)).*( sinsq .* cos(Zcutoff * Transf.KZ) - sqrt(sinsq.*cossq).*sin(Zcutoff * Transf.KZ) );
-            
-            % Nonanalytic part
-            % For a cylindrical cutoff, we need to construct a kr grid based on the 3D parameters using Bessel quadrature
-            VDkNon      = this.calculateNumericalHankelTransform(TransfRad.kr, TransfRad.kz, TransfRad.Rmax, Zcutoff);
-            
-            % Interpolating the nonanalytic part onto 3D grid
-            fullkr = [-flip(TransfRad.kr)',TransfRad.kr'];
-            [KR,KZ] = ndgrid(fullkr,TransfRad.kz);
-            [KX3D,KY3D,KZ3D] = ndgrid(ifftshift(Transf.kx),ifftshift(Transf.ky),ifftshift(Transf.kz));
-            KR3D = sqrt(KX3D.^2 + KY3D.^2);
-            fullVDK = [flip(VDkNon',2),VDkNon']';
-            VDkNon = interpn(KR,KZ,fullVDK,KR3D,KZ3D,'spline',0); %Last argument is -1/3 for full VDk. 0 for nonanalytic piece?
-            VDkNon = fftshift(VDkNon);
-            
-            VDk = VDk + VDkNon;
-    
-        case 'CylindricalInfiniteZ' %Cylindrical infinite Z, polarized along x -- PRA 107, 033301 (2023)
-	        alph = acos((Transf.KX*sin(Params.theta)*cos(Params.phi)+Transf.KY*sin(Params.theta)*sin(Params.phi)+Transf.KZ*cos(Params.theta))./sqrt(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2));
-	        alph(1) = pi/2;
-	        rhoc = max([abs(Transf.x),abs(Transf.y)]);
-	        KR = sqrt(Transf.KX.^2+Transf.KY.^2);
-	        func = @(n,u,v) v.^2./(u.^2+v.^2).*(v.*besselj(n,u).*besselk(n+1,v) - u.*besselj(n+1,u).*besselk(n,v));    
-	        VDk = -0.5*func(0,KR*rhoc,abs(Transf.KZ)*rhoc) + (Transf.KX.^2./KR.^2 - 0.5).*func(2,KR*rhoc,abs(Transf.KZ)*rhoc);
-	        VDk = (1/3)*(3*(cos(alph).^2)-1) - VDk;
-    
-            VDk(KR==0) = -1/3 + 1/2*abs(Transf.KZ(KR==0))*rhoc.*besselk(1,abs(Transf.KZ(KR==0))*rhoc);
-	        VDk(Transf.KZ==0) = 1/6 + (Transf.KX(Transf.KZ==0).^2-Transf.KY(Transf.KZ==0).^2)./(KR(Transf.KZ==0).^2).*(1/2 - besselj(1,KR(Transf.KZ==0)*rhoc)./(KR(Transf.KZ==0)*rhoc));
-	        VDk(1,1,1) = 1/6;            
-        
-        case 'Spherical' %Spherical
-            Rcut = min(Params.Lx/2,Params.Ly/2,Params.Lz/2);
-            alph = acos((Transf.KX*sin(Params.theta)*cos(Params.phi)+Transf.KY*sin(Params.theta)*sin(Params.phi)+Transf.KZ*cos(Params.theta))./sqrt(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2));
-	        alph(1) = pi/2;
-    
-            K = sqrt(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
-            VDk = (cos(alph).^2-1/3).*(1 + 3*cos(Rcut*K)./(Rcut^2.*K.^2) - 3*sin(Rcut*K)./(Rcut^3.*K.^3));
-        
-        otherwise
-	        disp('Choose a valid DDI cutoff type!')
-            return
-    end
-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 1cc2d12..b2907de 100644
--- a/Dipolar Gas Simulator/+Simulator/@DipolarGas/propagateWavefunction.m	
+++ b/Dipolar Gas Simulator/+Simulator/@DipolarGas/propagateWavefunction.m	
@@ -1,4 +1,4 @@
-function [psi] = propagateWavefunction(this,psi,Params,Transf,VDk,V,njob,t_idx,Observ)
+function [psi] = propagateWavefunction(this,psi,Params,Transf,VDk,V,t_idx,Observ)
     set(0,'defaulttextInterpreter','latex')
     set(groot, 'defaultAxesTickLabelInterpreter','latex'); set(groot, 'defaultLegendInterpreter','latex');
     
@@ -9,10 +9,16 @@ function [psi] = propagateWavefunction(this,psi,Params,Transf,VDk,V,njob,t_idx,O
             KEop= 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
             Observ.residual = 1; Observ.res = 1; 
             
-            muchem = this.Calculator.ChemicalPotential(psi,Params,Transf,VDk,V);
+            muchem = this.Calculator.calculateChemicalPotential(psi,Params,Transf,VDk,V);
             AdaptIdx = 0;
             
+            pb = Helper.ProgressBar();
+            pb.run('Running evolution in imaginary time: ');
+
             while t_idx < Params.sim_time_cut_off
+
+                pb.run(t_idx);
+
                 %kin
                 psi = fftn(psi);
                 psi = psi.*exp(-0.5*1i*dt*KEop);
@@ -34,12 +40,12 @@ function [psi] = propagateWavefunction(this,psi,Params,Transf,VDk,V,njob,t_idx,O
                 Norm = trapz(abs(psi(:)).^2)*Transf.dx*Transf.dy*Transf.dz;
                 psi = sqrt(Params.N)*psi/sqrt(Norm);
             
-                muchem = this.Calculator.ChemicalPotential(psi,Params,Transf,VDk,V);
+                muchem = this.Calculator.calculateChemicalPotential(psi,Params,Transf,VDk,V);
             
                 if mod(t_idx,1000) == 0        
             
                     %Change in Energy
-                    E = this.Calculator.TotalEnergy(psi,Params,Transf,VDk,V);
+                    E = this.Calculator.calculateTotalEnergy(psi,Params,Transf,VDk,V);
                     E = E/Norm;
                     Observ.EVec = [Observ.EVec E];
             
@@ -47,12 +53,12 @@ function [psi] = propagateWavefunction(this,psi,Params,Transf,VDk,V,njob,t_idx,O
                     Observ.mucVec = [Observ.mucVec muchem];
             
                     %Normalized residuals
-                    res = this.Calculator.NormalizedResiduals(psi,Params,Transf,VDk,V,muchem);
+                    res = this.Calculator.calculateNormalizedResiduals(psi,Params,Transf,VDk,V,muchem);
                     Observ.residual = [Observ.residual res];
                     
                     Observ.res_idx = Observ.res_idx + 1;
             
-                    save(sprintf('./Data/Run_%03i/psi_gs.mat',njob),'psi','muchem','Observ','t_idx','Transf','Params','VDk','V');
+                    save(sprintf('./Data/Run_%03i/psi_gs.mat',Params.njob),'psi','muchem','Observ','t_idx','Transf','Params','VDk','V');
                     
                     %Adaptive time step -- Careful, this can quickly get out of control
                     relres = abs(Observ.residual(Observ.res_idx)-Observ.residual(Observ.res_idx-1))/Observ.residual(Observ.res_idx);
@@ -78,16 +84,17 @@ function [psi] = propagateWavefunction(this,psi,Params,Transf,VDk,V,njob,t_idx,O
             end
             
             %Change in Energy
-            E = this.Calculator.TotalEnergy(psi,Params,Transf,VDk,V);
+            E = this.Calculator.calculateTotalEnergy(psi,Params,Transf,VDk,V);
             E = E/Norm;
             Observ.EVec = [Observ.EVec E];
             
             % Phase coherence
-            [PhaseC] = this.Calculator.PhaseCoherence(psi,Transf,Params);
+            [PhaseC] = this.Calculator.calculatePhaseCoherence(psi,Transf,Params);
             Observ.PCVec = [Observ.PCVec PhaseC];
             
             Observ.res_idx = Observ.res_idx + 1;
-            save(sprintf('./Data/Run_%03i/psi_gs.mat',njob),'psi','muchem','Observ','t_idx','Transf','Params','VDk','V');
+            save(sprintf('./Data/Run_%03i/psi_gs.mat',Params.njob),'psi','muchem','Observ','t_idx','Transf','Params','VDk','V');
+            pb.run(' - Job Completed!\n');
         
         case 'RealTimeEvolution'
 
@@ -95,9 +102,15 @@ function [psi] = propagateWavefunction(this,psi,Params,Transf,VDk,V,njob,t_idx,O
         
             KEop= 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
             
-            muchem = chemicalpotential(psi,Params,Transf,VDk,V);
+            muchem = this.Calculator.calculateChemicalPotential(psi,Params,Transf,VDk,V);
             
+            pb = Helper.ProgressBar();
+            pb.run('Running evolution in real time: ');
+
             while t_idx < Params.sim_time_cut_off
+
+                pb.run(t_idx);
+
                 % Parameters at time t
                 
                 %kin
@@ -117,7 +130,7 @@ function [psi] = propagateWavefunction(this,psi,Params,Transf,VDk,V,njob,t_idx,O
                 psi = psi.*exp(-0.5*1i*dt*KEop);
                 psi = ifftn(psi);
                 
-                muchem = chemicalpotential(psi,Params,Transf,VDk,V);
+                muchem = this.Calculator.calculateChemicalPotential(psi,Params,Transf,VDk,V);
             
                 if mod(t_idx,1000)==0        
                     %Change in Normalization
@@ -125,18 +138,18 @@ function [psi] = propagateWavefunction(this,psi,Params,Transf,VDk,V,njob,t_idx,O
                     Observ.NormVec = [Observ.NormVec Norm];
             
                     %Change in Energy
-                    E = energytotal(psi,Params,Transf,VDk,V);
+                    E = this.Calculator.calculateTotalEnergy(psi,Params,Transf,VDk,V);
                     E = E/Norm;
                     Observ.EVec = [Observ.EVec E];
             
                     % Phase coherence
-                    [PhaseC] = PhaseCoherence(psi,Transf);
+                    [PhaseC] = this.Calculator.calculatePhaseCoherence(psi,Transf);
                     Observ.PCVec = [Observ.PCVec PhaseC];
             
                     Observ.tVecPlot = [Observ.tVecPlot tVal];
                     Observ.res_idx = Observ.res_idx + 1;
                     
-                    save(sprintf('./Data/Run_%03i/TimeEvolution/psi_%i.mat',njob,Observ.res_idx),'psi','muchem','Observ','t_idx');        
+                    save(sprintf('./Data/Run_%03i/TimeEvolution/psi_%i.mat',Params.njob,Observ.res_idx),'psi','muchem','Observ','t_idx');        
                 end
                 if any(isnan(psi(:)))
                     disp('NaNs encountered!')
@@ -150,18 +163,19 @@ function [psi] = propagateWavefunction(this,psi,Params,Transf,VDk,V,njob,t_idx,O
             Observ.NormVec = [Observ.NormVec Norm];
             
             %Change in Energy
-            E = energytotal(psi,Params,Transf,VDk,V);
+            E = this.Calculator.calculateTotalEnergy(psi,Params,Transf,VDk,V);
             E = E/Norm;
             Observ.EVec = [Observ.EVec E];
             
             % Phase coherence
-            [PhaseC] = PhaseCoherence(psi,Transf);
+            [PhaseC] = this.Calculator.calculatePhaseCoherence(psi,Transf);
             Observ.PCVec = [Observ.PCVec PhaseC];
             
             Observ.tVecPlot = [Observ.tVecPlot tVal];
             
             Observ.res_idx = Observ.res_idx + 1;
-            save(sprintf('./Data/Run_%03i/TimeEvolution/psi_%i.mat',njob,Observ.res_idx),'psi','muchem','Observ','t_idx');
+            save(sprintf('./Data/Run_%03i/TimeEvolution/psi_%i.mat',Params.njob,Observ.res_idx),'psi','muchem','Observ','t_idx');
+            pb.run(' - Job Completed!\n');
         
         otherwise
 	        disp('Choose a valid DDI cutoff type!')
diff --git a/Dipolar Gas Simulator/+Simulator/@DipolarGas/run.m b/Dipolar Gas Simulator/+Simulator/@DipolarGas/run.m
index 4db3dda..1d4adfb 100644
--- a/Dipolar Gas Simulator/+Simulator/@DipolarGas/run.m	
+++ b/Dipolar Gas Simulator/+Simulator/@DipolarGas/run.m	
@@ -13,10 +13,10 @@ function [Params, Transf, psi,V,VDk] = run(this)
     [psi,V,VDk] = this.initialize(Params,Transf,TransfRad);
     
     Observ.EVec = []; Observ.NormVec = []; Observ.PCVec = []; Observ.tVecPlot = []; Observ.mucVec = [];
-    t_idx = 1; %Start at t = 0;
+    t_idx = 1; % Start at t = 0;
     Observ.res_idx = 1;
 
     % --- Run Simulation --- 
-    % mkdir(sprintf('./Data/Run_%03i',Params.njob))
-    % [psi] = this.propagateWavefunction(psi,Params,Transf,VDk,V,njob,t_idx,Observ);
+    mkdir(sprintf('./Data/Run_%03i',Params.njob))
+    [psi] = this.propagateWavefunction(psi,Params,Transf,VDk,V,t_idx,Observ);
 end