diff --git a/Dipolar Gas Simulator/+Simulator/@Calculator/calculateChemicalPotential.m b/Dipolar Gas Simulator/+Simulator/@Calculator/calculateChemicalPotential.m
index 697ab38..a7c1c05 100644
--- a/Dipolar Gas Simulator/+Simulator/@Calculator/calculateChemicalPotential.m	
+++ b/Dipolar Gas Simulator/+Simulator/@Calculator/calculateChemicalPotential.m	
@@ -1,4 +1,4 @@
-function muchem = ChemicalPotential(psi,Params,Transf,VDk,V)
+function muchem = calculateChemicalPotential(psi,Params,Transf,VDk,V)
 
 %Parameters
 normfac = Params.Lx*Params.Ly*Params.Lz/numel(psi);
diff --git a/Dipolar Gas Simulator/+Simulator/@Calculator/calculateEnergyComponents.m b/Dipolar Gas Simulator/+Simulator/@Calculator/calculateEnergyComponents.m
index 49159ed..327bfde 100644
--- a/Dipolar Gas Simulator/+Simulator/@Calculator/calculateEnergyComponents.m	
+++ b/Dipolar Gas Simulator/+Simulator/@Calculator/calculateEnergyComponents.m	
@@ -1,4 +1,4 @@
-function E = EnergyComponents(psi,Params,Transf,VDk,V)
+function E = calculateEnergyComponents(psi,Params,Transf,VDk,V)
 
 %Parameters
 
diff --git a/Dipolar Gas Simulator/+Simulator/@Calculator/calculateNormalizedResiduals.m b/Dipolar Gas Simulator/+Simulator/@Calculator/calculateNormalizedResiduals.m
index 66e5316..a08e895 100644
--- a/Dipolar Gas Simulator/+Simulator/@Calculator/calculateNormalizedResiduals.m	
+++ b/Dipolar Gas Simulator/+Simulator/@Calculator/calculateNormalizedResiduals.m	
@@ -1,4 +1,4 @@
-function res = NormalizedResiduals(psi,Params,Transf,VDk,V,muchem)
+function res = calculateNormalizedResiduals(psi,Params,Transf,VDk,V,muchem)
 
 KEop= 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
 
diff --git a/Dipolar Gas Simulator/+Simulator/@Calculator/calculateNumericalHankelTransform.m b/Dipolar Gas Simulator/+Simulator/@Calculator/calculateNumericalHankelTransform.m
index 5e5be07..141de91 100644
--- a/Dipolar Gas Simulator/+Simulator/@Calculator/calculateNumericalHankelTransform.m	
+++ b/Dipolar Gas Simulator/+Simulator/@Calculator/calculateNumericalHankelTransform.m	
@@ -1,4 +1,4 @@
-function VDkSemi = NumericalHankelTransform(kr, kz, Rmax, Zmax, Nr)
+function VDkSemi = calculateNumericalHankelTransform(kr, kz, Rmax, Zmax, Nr)
 
 % accuracy inputs for numerical integration
 if(nargin==4)
diff --git a/Dipolar Gas Simulator/+Simulator/@Calculator/calculateOrderParameter.m b/Dipolar Gas Simulator/+Simulator/@Calculator/calculateOrderParameter.m
index 08b0fa4..a8ae9e4 100644
--- a/Dipolar Gas Simulator/+Simulator/@Calculator/calculateOrderParameter.m	
+++ b/Dipolar Gas Simulator/+Simulator/@Calculator/calculateOrderParameter.m	
@@ -1,4 +1,4 @@
-function [m_Order] = OrderParameter(psi,Transf,Params,VDk,V,T,muchem)
+function [m_Order] = calculateOrderParameter(psi,Transf,Params,VDk,V,T,muchem)
 
     NumRealiz = 100;
 
diff --git a/Dipolar Gas Simulator/+Simulator/@Calculator/calculatePhaseCoherence.m b/Dipolar Gas Simulator/+Simulator/@Calculator/calculatePhaseCoherence.m
index 9b3442e..0cf034f 100644
--- a/Dipolar Gas Simulator/+Simulator/@Calculator/calculatePhaseCoherence.m	
+++ b/Dipolar Gas Simulator/+Simulator/@Calculator/calculatePhaseCoherence.m	
@@ -1,4 +1,4 @@
-function [PhaseC] = PhaseCoherence(psi,Transf,Params)
+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);
diff --git a/Dipolar Gas Simulator/+Simulator/@Calculator/calculateTotalEnergy.m b/Dipolar Gas Simulator/+Simulator/@Calculator/calculateTotalEnergy.m
index 9d2c7ff..213f56d 100644
--- a/Dipolar Gas Simulator/+Simulator/@Calculator/calculateTotalEnergy.m	
+++ b/Dipolar Gas Simulator/+Simulator/@Calculator/calculateTotalEnergy.m	
@@ -1,4 +1,4 @@
-function E = TotalEnergy(psi,Params,Transf,VDk,V)
+function E = calculateTotalEnergy(psi,Params,Transf,VDk,V)
 
 %Parameters
 
diff --git a/Dipolar Gas Simulator/+Simulator/@Calculator/calculateVDCutoff.m b/Dipolar Gas Simulator/+Simulator/@Calculator/calculateVDCutoff.m
index dba1435..d5310c1 100644
--- a/Dipolar Gas Simulator/+Simulator/@Calculator/calculateVDCutoff.m	
+++ b/Dipolar Gas Simulator/+Simulator/@Calculator/calculateVDCutoff.m	
@@ -1,4 +1,4 @@
-function VDk = VDCutoff(Params, Transf)
+function VDk = calculateVDCutoff(Params, Transf)
 % 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
diff --git a/Dipolar Gas Simulator/+Simulator/@DipolarGas/initialize.m b/Dipolar Gas Simulator/+Simulator/@DipolarGas/initialize.m
index 63b7abb..2e4c9c2 100644
--- a/Dipolar Gas Simulator/+Simulator/@DipolarGas/initialize.m	
+++ b/Dipolar Gas Simulator/+Simulator/@DipolarGas/initialize.m	
@@ -1,4 +1,4 @@
-function [psi,V,VDk] = Initialize(Params,Transf)
+function [psi,V,VDk] = initialize(Params,Transf)
 
 format long
 X = Transf.X; Y = Transf.Y; Z = Transf.Z;
diff --git a/Dipolar Gas Simulator/+Simulator/@DipolarGas/runSimulation.m b/Dipolar Gas Simulator/+Simulator/@DipolarGas/runSimulation.m
index a99ddf8..57d966c 100644
--- a/Dipolar Gas Simulator/+Simulator/@DipolarGas/runSimulation.m	
+++ b/Dipolar Gas Simulator/+Simulator/@DipolarGas/runSimulation.m	
@@ -1,12 +1,11 @@
 %-% Run Simulation %-%
 clearvars
-cd(fileparts(matlab.desktop.editor.getActiveFilename));
-pwd
+
 % --- Obtain simulation parameters --- 
-[Params]    = SetupParameters();
+[Params]    = setupParameters();
 
 % --- Set up spatial grids and transforms --- 
-[Transf]    = SetupSpace(Params);
+[Transf]    = setupSpace(Params);
 
 % Plotter.visualizeSpace(Transf)
 
@@ -16,7 +15,7 @@ pwd
 
 mkdir(sprintf('./Data'))
 
-[psi,V,VDk] = Initialize(Params,Transf);
+[psi,V,VDk] = initialize(Params,Transf);
 
 Observ.EVec = []; Observ.NormVec = []; Observ.PCVec = []; Observ.tVecPlot = []; Observ.mucVec = [];
 t_idx = 1; %Start at t = 0;
diff --git a/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupParameters.m b/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupParameters.m
index d2e6932..feac996 100644
--- a/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupParameters.m	
+++ b/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupParameters.m	
@@ -1,4 +1,4 @@
-function [Params] = SetupParameters()
+function [Params] = setupParameters()
 
 %%--%% Parameters %%--%%
 
diff --git a/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupSpace.m b/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupSpace.m
index d970966..ea98ba4 100644
--- a/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupSpace.m	
+++ b/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupSpace.m	
@@ -1,4 +1,4 @@
-function [Transf] = SetupSpace(Params)
+function [Transf] = setupSpace(Params)
 Transf.Xmax = 0.5*Params.Lx;
 Transf.Ymax = 0.5*Params.Ly;
 Transf.Zmax = 0.5*Params.Lz;
diff --git a/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupSpaceRadial.m b/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupSpaceRadial.m
index e0a5d3c..ce2c96d 100644
--- a/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupSpaceRadial.m	
+++ b/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupSpaceRadial.m	
@@ -1,4 +1,4 @@
-function [Transf] = SetupSpaceRadial(Params,morder)
+function [Transf] = setupSpaceRadial(Params,morder)
 Zmax = 0.5*Params.Lz;
 Rmax = Params.Lr;
 Nz = Params.Nz;
diff --git a/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupWavefunction.m b/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupWavefunction.m
index 79ade2a..d1d90d3 100644
--- a/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupWavefunction.m	
+++ b/Dipolar Gas Simulator/+Simulator/@DipolarGas/setupWavefunction.m	
@@ -1,4 +1,4 @@
-function [psi] = SetupWavefunction()
+function [psi] = setupWavefunction()
 
 ellx = sqrt(Params.hbar/(Params.m*Params.wx))/Params.l0; 
 elly = sqrt(Params.hbar/(Params.m*Params.wy))/Params.l0;
diff --git a/Dipolar Gas Simulator/+Simulator/@ImaginaryTimeEvolution/splitStepFourier.m b/Dipolar Gas Simulator/+Simulator/@ImaginaryTimeEvolution/splitStepFourier.m
index c69ceb5..1397739 100644
--- a/Dipolar Gas Simulator/+Simulator/@ImaginaryTimeEvolution/splitStepFourier.m	
+++ b/Dipolar Gas Simulator/+Simulator/@ImaginaryTimeEvolution/splitStepFourier.m	
@@ -1,4 +1,4 @@
-function [psi] = SplitStepFourier(psi,Params,Transf,VDk,V,njob,t_idx,Observ)
+function [psi] = splitStepFourier(psi,Params,Transf,VDk,V,njob,t_idx,Observ)
 
 set(0,'defaulttextInterpreter','latex')
 set(groot, 'defaultAxesTickLabelInterpreter','latex'); set(groot, 'defaultLegendInterpreter','latex');
diff --git a/GaussianBeamABCD/BeamPropagation.py b/GaussianBeamABCD Calculator/BeamPropagation.py
similarity index 100%
rename from GaussianBeamABCD/BeamPropagation.py
rename to GaussianBeamABCD Calculator/BeamPropagation.py
diff --git a/GaussianBeamABCD/Example.ipynb b/GaussianBeamABCD Calculator/Example.ipynb
similarity index 100%
rename from GaussianBeamABCD/Example.ipynb
rename to GaussianBeamABCD Calculator/Example.ipynb
diff --git a/GaussianBeamABCD/DMD_Setup.py b/GaussianBeamABCD Calculator/calculateForDMDSetup.py
similarity index 100%
rename from GaussianBeamABCD/DMD_Setup.py
rename to GaussianBeamABCD Calculator/calculateForDMDSetup.py
diff --git a/GaussianBeamABCD/README.md b/GaussianBeamABCD/README.md
deleted file mode 100644
index a76085a..0000000
--- a/GaussianBeamABCD/README.md
+++ /dev/null
@@ -1,359 +0,0 @@
-# Gaussian Beam Propagation
-
-## Import files
-
-
-```python
-import BeamProp_Script as bs      # This is the script that handles the propagation
-import sympy as sym               # For Symbolic examples
-import numpy as np                # Handling of lists and for plotting
-import matplotlib.pyplot as plt   # Plotting
-```
-
-### Let's show what BeamProp_Script has
-
-
-```python
-help(bs)
-```
-
-    Help on module BeamProp_Script:
-    
-    NAME
-        BeamProp_Script - Created on Wed Feb 19 15:51:54 2020
-    
-    DESCRIPTION
-        @author: wrighta
-    
-    FUNCTIONS
-        W0(zr, lam)
-            Parameters
-            ----------
-            zr : float, integer, symbol
-                Rayleigh range in meters
-            lam : float, integer, symbol
-                Wavelength of light in meters
-            
-            Returns
-            -------
-            w0 : float, integer, symbol
-                Beam waist radius in meters
-        
-        Zr(wo, lam)
-            Parameters
-            ----------
-            wo : float, integer, or symbol
-                Beam waist radius in meters.
-            lam : float, integer, or symbol
-                Wavelength of light in meters.
-            
-            Returns
-            -------
-            zr : float, int, symbols
-                Rayleigh range for given beam waist and wavelength.
-        
-        lens(f)
-            Parameters
-            ----------
-            f : float or integer or sympy symbol in meters
-                Thin lens focal length in meters
-            
-            Returns
-            -------
-            mat : 2x2 matrix
-            [
-            [   1, 0],
-            [-1/f, 1]
-            ]
-        
-        mult(mat1, *argv)
-            Parameters
-            ----------
-            mat1 : 2x2 ABCD matrix
-                Last matrix light interacts with.
-            *argv : 2x2 ABCD matrices 
-                From left to right, the matrices should be entered such that the leftmost matrix interacts
-                with light temporally after the rightmost matrix.
-            
-            Returns
-            -------
-            Mat : 2x2 matrix
-                The ABCd matrix describing the whole optical system.
-        
-        plot(func, var, rang=array([0.  , 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.1 ,
-               0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.19, 0.2 , 0.21,
-               0.22, 0.23, 0.24, 0.25, 0.26, 0.27, 0.28, 0.29, 0.3 , 0.31, 0.32,
-               0.33, 0.34, 0.35, 0.36, 0.37, 0.38, 0.39, 0.4 , 0.41, 0.42, 0.43,
-               0.44, 0.45, 0.46, 0.47, 0.48, 0.49, 0.5 , 0.51, 0.52, 0.53, 0.54,
-               0.55, 0.56, 0.57, 0.58, 0.59, 0.6 , 0.61, 0.62, 0.63, 0.64, 0.65,
-               0.66, 0.67, 0.68, 0.69, 0.7 , 0.71, 0.72, 0.73, 0.74, 0.75, 0.76,
-               0.77, 0.78, 0.79, 0.8 , 0.81, 0.82, 0.83, 0.84, 0.85, 0.86, 0.87,
-               0.88, 0.89, 0.9 , 0.91, 0.92, 0.93, 0.94, 0.95, 0.96, 0.97, 0.98,
-               0.99, 1.  , 1.01, 1.02, 1.03, 1.04, 1.05, 1.06, 1.07, 1.08, 1.09,
-               1.1 , 1.11, 1.12, 1.13, 1.14, 1.15, 1.16, 1.17, 1.18, 1.19, 1.2 ,
-               1.21, 1.22, 1.23, 1.24, 1.25, 1.26, 1.27, 1.28, 1.29, 1.3 , 1.31,
-               1.32, 1.33, 1.34, 1.35, 1.36, 1.37, 1.38, 1.39, 1.4 , 1.41, 1.42,
-               1.43, 1.44, 1.45, 1.46, 1.47, 1.48, 1.49, 1.5 , 1.51, 1.52, 1.53,
-               1.54, 1.55, 1.56, 1.57, 1.58, 1.59, 1.6 , 1.61, 1.62, 1.63, 1.64,
-               1.65, 1.66, 1.67, 1.68, 1.69, 1.7 , 1.71, 1.72, 1.73, 1.74, 1.75,
-               1.76, 1.77, 1.78, 1.79, 1.8 , 1.81, 1.82, 1.83, 1.84, 1.85, 1.86,
-               1.87, 1.88, 1.89, 1.9 , 1.91, 1.92, 1.93, 1.94, 1.95, 1.96, 1.97,
-               1.98, 1.99, 2.  , 2.01, 2.02, 2.03, 2.04, 2.05, 2.06, 2.07, 2.08,
-               2.09, 2.1 , 2.11, 2.12, 2.13, 2.14, 2.15, 2.16, 2.17, 2.18, 2.19,
-               2.2 , 2.21, 2.22, 2.23, 2.24, 2.25, 2.26, 2.27, 2.28, 2.29, 2.3 ,
-               2.31, 2.32, 2.33, 2.34, 2.35, 2.36, 2.37, 2.38, 2.39, 2.4 , 2.41,
-               2.42, 2.43, 2.44, 2.45, 2.46, 2.47, 2.48, 2.49, 2.5 , 2.51, 2.52,
-               2.53, 2.54, 2.55, 2.56, 2.57, 2.58, 2.59, 2.6 , 2.61, 2.62, 2.63,
-               2.64, 2.65, 2.66, 2.67, 2.68, 2.69, 2.7 , 2.71, 2.72, 2.73, 2.74,
-               2.75, 2.76, 2.77, 2.78, 2.79, 2.8 , 2.81, 2.82, 2.83, 2.84, 2.85,
-               2.86, 2.87, 2.88, 2.89, 2.9 , 2.91, 2.92, 2.93, 2.94, 2.95, 2.96,
-               2.97, 2.98, 2.99]))
-            Parameters
-            ----------
-            func : Sympy function of one variable
-                Sympy function defining the beam width after the last optical element.
-            var : sympy variable
-                Variable in func that will be plotted.
-            rang : numpy array
-                Array of the values along the optical axis to be plotted
-            
-            Returns
-            -------
-            plot : matplotlib graph
-                Graph of the beam width of var
-        
-        prop(d)
-            Parameters
-            ----------
-            d : float or integer or sympy symbol
-                Distance light is propagating along the z-axis.
-            
-            Returns
-            -------
-            mat: 2x2 matrix
-            [
-            [1, d],
-            [0, 1]
-            ]
-        
-        q1_func(z, w0, lam, mat)
-            Parameters
-            ----------
-            z : float, int, symbol
-                Position of the beam waist in meters.
-            w0 : float, int, symbol
-                Radial waist size in meters (of the embedded Gaussian, i.e. W0/M).
-            lam : float, int, symbol
-                Wavelength of light in meters.
-            mat : float, int, symbol
-                The ABCD 2x2 matrix describing the optical system.
-            
-            Returns
-            -------
-            z: float, int, symbol
-                Position of the beam waist after the optical system
-            zr: float, int, symbol
-                Rayleigh range of the beam after the optical system
-        
-        q1_inv_func(z, w0, lam, mat)
-            Parameters
-            ----------
-            z : float, int, symbol
-                Position of the beam waist in meters.
-            w0 : float, int, symbol
-                Radial waist size in meters (of the embedded Gaussian, i.e. W0/M).
-            lam : float, int, symbol
-                Wavelength of light in meters.
-            mat : float, int, symbol
-                The ABCD 2x2 matrix describing the optical system.
-            
-            Returns
-            -------
-            R : float, int, symbol
-                Radius of curvature of the wavefront in meters.
-            w : float, int, symbol
-                Radius of the beam in meters.
-        
-        ray(y, theta)
-            Parameters
-            ----------
-            y : float or integer or sympy symbol in meters
-                The vertical height of a ray.
-            theta : float or integer in radians
-                The angle of divergence of the ray.
-            
-            Returns
-            -------
-            mat : 2x1 matrix
-                [
-                [y],
-                [teta]
-                ]
-    
-    DATA
-        oo = oo
-    
-    FILE
-        c:\users\wrighta\documents\beamprop\beamprop_script.py
-    
-    
-    
-
-## Let's first see how we define a beam and how we can visualize it propagating.
-
-### A Gaussian beam can be defined by it's (radial) waist, $w_0$, it's Rayleigh range, $z_R = \frac{\pi * w_0^2}{\lambda}$, and the location of its waist, $z_0$.
-
-
-```python
-w0 = 1E-3 # 1mm beam waist
-lam = 355E-9 # wavelength of 355 nm (UV)
-zR = bs.Zr(w0, lam) # Rayleigh range in m
-z0 = 0 # location of waist in m
-```
-
-### We now want to define our "optical system" using matrices. For this first example, we will just use a free space propagation matrix, and let the beam propagate a distance $d$ which we will define using a symbol.
-
-
-```python
-d = sym.symbols('d')
-M = bs.prop(d)
-```
-
-### We now use the bs script to do all the ABCD and q-parameter math, and return the waist and radius of curvature functions
-
-
-```python
-R, w = bs.q1_inv_func(0, w0, lam, M)
-```
-
-
-```python
-print('w = {}'.format(w))
-```
-
-    w = 0.001*(0.0127690021685256*d**2 + 1)**0.5
-    
-
-### And as simple as that, we have a function for our waist. Let's plot it and see what it looks like
-
-
-```python
-bs.plot(w, d, rang = np.arange(0,10))
-```
-
-
-![png](output_14_0.png)
-
-
-### Let's show what happens when a beam travels through a lens. We use the "mult" function to multiply multiple ABCD matrices together.
-
-
-```python
-w0 = 1E-3 # 1mm beam waist
-lam = 355E-9 # wavelength of 355 nm (UV)
-zR = bs.Zr(w0, lam) # Rayleigh range in m
-z0 = 0 # location of waist in m
-
-d = sym.symbols('d')
-M = bs.mult(bs.prop(d), bs.lens(.5), bs.prop(1))
-            
-R, w = bs.q1_inv_func(0, w0, lam, M)
-
-bs.plot(w, d, rang = np.arange(0,1,.01))
-```
-
-
-![png](output_16_0.png)
-
-
-### Lets look at how to expand and collimate a beam with a two lens system
-
-
-```python
-w0 = 1E-3 # 1mm beam waist
-lam = 355E-9 # wavelength of 355 nm (UV)
-zR = bs.Zr(w0, lam) # Rayleigh range in m
-z0 = 0 # location of waist in m
-
-d1, d2, d3, f1, f2 = sym.symbols('d1 d2 d3 f1 f2')
-
-M = bs.mult(bs.prop(d3),bs.lens(f2),bs.prop(d2), bs.lens(f1), bs.prop(d1))
-
-R, w = bs.q1_inv_func(0, w0, lam, M)
-```
-
-### For example, lets say the beam travels 1 m before hitting the first lens, and we want the beam to be 5x w0 after coming out of the second lens. We substitute d1 for 1 meter, since the beam propagates 1 meter, and we substitute d3 for 0, since we only care about the beam size right at the second lens. This gives us a relation between f1 and d2 (the separation between the lenses).
-
-
-```python
-w = w.subs(d1,1).subs(d3,0)
-f1_eq = sym.solve(w - 5*w0, f1)[0]
-print('f = {}'.format(f1_eq))
-```
-
-    f = 1.0084642216545e+15*d2*(1.12051580183833e+27*d2 - 4.41556446152598e+29*sqrt(1 - 0.000504320418227052*d2**2) + 8.88733242867719e+28)/(1.13000009595246e+42*d2**2 + 2.26000019190491e+42*d2 - 2.12276362486616e+45)
-    
-
-#### Suppose we wanted the distance between the lenses to be 1 meter, we could find what f1 we need.
-
-
-```python
-print('f1 = {:.2f} m, for a lens separation of 1 meter'.format(f1_eq.subs(d2, 1)))
-```
-
-    f1 = 0.17 m, for a lens separation of 1 meter
-    
-
-### Now we need to collimate the beam. Lets still assume the beam propagates 1 m, and f1 = .17 m.
-
-There are a couple different ways to think about collimation. One is that the beam size doesn't change over a long distance. The other is that the radius of curvature is infinite (i.e. a plane wave). Lets us the latter interpretation. Thus, we want to find the focal length f2 that makes R infinite, or that makes 1/R =0.
-
-
-```python
-R_coll = R.subs(d1,1).subs(d2,1).subs(f1,.17).subs(d3,0)
-f2_coll = sym.solve(1/R_coll,f2)[0]
-print('f2 = {:.2f}, for a collimated beam, 5x the original waist, after propagating 1m to the first lens of f1 = .17m, and propagating another 1m to the second lens'.format(f2_coll))
-```
-
-    f2 = 0.83, for a collimated beam, 5x the original waist, after propagating 1m to the first lens of f1 = .17m, and propagating another 1m to the second lens
-    
-
-### Lets plot the beam profile after the second lens, and see if it is collimated.
-
-
-```python
-M = bs.mult(bs.prop(d3),bs.lens(.83),bs.prop(1), bs.lens(.17), bs.prop(1))
-
-R, w = bs.q1_inv_func(0, w0, lam, M)
-
-bs.plot(w,d3)
-```
-
-
-![png](output_27_0.png)
-
-
-### Looks very collimated. Lets check the beam size (to make sure its 5* w0) and check the collimation
-
-
-```python
-expansion_factor = w.subs(d3,0)/ w0
-print('beam is w = {:.2f} x w0'.format(expansion_factor))
-```
-
-    beam is w = 4.90 x w0
-    
-
-
-```python
-beam_size_change = (w.subs(d3,10) - w.subs(d3,0)) / w.subs(d3,0) * 100
-print('Over 10 m after second lens, beam changes by {:.0f}%'.format(beam_size_change))
-```
-
-    Over 10 m after second lens, beam changes by 1%
-    
-
-
-```python
-
-```
diff --git a/IRF/doNoiseCorrelation.m b/IRF/doNoiseCorrelation.m
deleted file mode 100644
index e26e3dc..0000000
--- a/IRF/doNoiseCorrelation.m
+++ /dev/null
@@ -1,86 +0,0 @@
-function [result,avgpic]=doNoiseCorrelation(imgs,mask)
-%imgs:cell arrays of the Nim images to treat together. each element of
-%imgs should have the same size:(sy,sy)
-
-Nim=numel(imgs);
-
-%% initialize sizes
-enlarge=1;%use so that it works well whatever the images sizes are (even/odd)
-ROIsizey = size(imgs{1},1)-1;
-if enlarge; maxsizey = ROIsizey*2; else maxsizey = ROIsizey+1; end%*2;
-sy = ROIsizey+1;
-ROIsizex = size(imgs{1},2)-1;
-if enlarge; maxsizex = ROIsizex*2; else maxsizex = ROIsizex+1; end%*2
-sx = ROIsizex+1;
-
-xcors=(-ROIsizex/2:ROIsizex/2);%center=0
-ycors=(-ROIsizey/2:ROIsizey/2);%center=0
-avgpic = zeros(maxsizey,maxsizex);
-avgpicFFT = zeros(maxsizey,maxsizex);
-
-counter = 0;% not really need if all images of the table are used as then always equal to indim but can come useful is selection is done
-m0=0.01;%for plot
-
-
-
-for indim=1:Nim
-    
-    disp(['>>>>>>>>>>>>>>>>>>>>>> Image ' num2str(indim) '/' num2str(numel(Nim)) ' <<<<<<<<<<<<<<<<<<<<<'])
-    counter=counter+1;
-    
-    %% resized data for the NoiseCorrelations.
-    if nargin==1;ROI1=imgs{indim};else ROI1=imgs{indim}.*mask;end
-    text1='raw images';
-    
-    
-    pixToCalc = zeros(maxsizey,maxsizex);
-    pixToCalc(1:sy,1:sx) = ROI1;
-    pixToCalc = pixToCalc/sum(sum(pixToCalc));
-    
-    %Calculate correlation function
-    picAuxFFT = ifftshift(ifft2((abs(fft2(pixToCalc)).^2)));
-    
-    %for means
-    avgpic = avgpic+pixToCalc;
-    avgpicFFT = avgpicFFT+picAuxFFT;
-    
-    %temporary means:
-    avgpictemp = avgpic/counter;
-    figure(1); clf;
-    subplot(2,2,1);imagesc(avgpictemp); colorbar;hold all
-    
-    avgpicFFTtemp = avgpicFFT/counter;
-    subplot(2,2,2);imagesc(abs(avgpicFFTtemp)); colorbar;hold all
-    
-    avgpictemp = ifftshift(ifft2((abs(fft2(avgpictemp)).^2)));
-    subplot(2,2,3);imagesc(abs(avgpictemp)); colorbar;hold all
-    
-    
-    %temporary results:
-    result = (avgpicFFTtemp./avgpictemp-1);
-    
-    
-    
-    if enlarge result=result(ROIsizey/2+1:ROIsizey*3/2+1,ROIsizex/2+1:ROIsizex*3/2+1); end%/result(ROIsizey+1,ROIsizex+1);
-    subplot(2,2,4);imagesc(real(result),[-1 1]); colorbar;hold all
-    plot(ROIsizey/2+1,ROIsizex/2+1,'w+')
-    
-    normr=result(ROIsizey/2+1,ROIsizex/2+1);
-    disp(['Normalization:' num2str(normr)])
-    
-    %temporary plot
-    disp(['plot ' num2str(indim) '...'])
-    %plot correlation function
-    figure(100);clf
-    subplot(1,2,1); imagesc(ROI1);title(['Im' num2str(indim)])
-    subplot(1,2,2);  imagesc(xcors,ycors,real(result),[-m0/2,m0]);hold all; plot(0,0,'w+');colorbar;
-    nametitle = [text1 '- Nb averages: ',num2str(counter) ', norm: ' num2str(normr)];
-    title(nametitle);
-    
-    drawnow;
-    
-end
-avgpic=avgpic(1:sy,1:sx);
-end
-
-
diff --git a/IRF/fringeremoval.m b/IRF/fringeremoval.m
deleted file mode 100644
index 06b64b5..0000000
--- a/IRF/fringeremoval.m
+++ /dev/null
@@ -1,69 +0,0 @@
-function [optrefimages] = fringeremoval(absimages, refimages, bgmask)
-% FRINGEREMOVAL - Fringe removal and noise reduction from absorption images.
-% Creates an optimal reference image for each absorption image in a set as 
-% a linear combination of reference images, with coefficients chosen to 
-% minimize the least-squares residuals between each absorption image and 
-% the optimal reference image. The coefficients are obtained by solving a 
-% linear set of equations using matrix inverse by LU decomposition.
-%
-% Application of the algorithm is described in C. F. Ockeloen et al, Improved 
-% detection of small atom numbers through image processing, arXiv:1007.2136 (2010).
-%
-% Syntax:
-%    [optrefimages] = fringeremoval(absimages,refimages,bgmask);
-%
-% Required inputs:
-%    absimages    - Absorption image data, 
-%                   typically 16 bit grayscale images
-%    refimages    - Raw reference image data
-%       absimages and refimages are both cell arrays containing
-%       2D array data. The number of refimages can differ from the 
-%       number of absimages.
-%
-% Optional inputs:
-%    bgmask       - Array specifying background region used, 
-%                   1=background, 0=data. Defaults to all ones.
-% Outputs:
-%    optrefimages - Cell array of optimal reference images, 
-%                   equal in size to absimages.
-%
-
-% Dependencies: none
-%
-% Authors: Shannon Whitlock, Caspar Ockeloen
-% Reference: C. F. Ockeloen, A. F. Tauschinsky, R. J. C. Spreeuw, and
-%            S. Whitlock, Improved detection of small atom numbers through 
-%            image processing, arXiv:1007.2136
-% Email: 
-% May 2009; Last revision: 11 August 2010
-
-% Process inputs
-
-% Set variables, and flatten absorption and reference images
-nimgs  = size(absimages,3);
-nimgsR = size(refimages,3);
-xdim = size(absimages(:,:,1),2);
-ydim = size(absimages(:,:,1),1);
-
-R = single(reshape(refimages,xdim*ydim,nimgsR));
-A = single(reshape(absimages,xdim*ydim,nimgs));
-optrefimages=zeros(size(absimages)); % preallocate
-
-if not(exist('bgmask','var')); bgmask=ones(ydim,xdim); end
-k = find(bgmask(:)==1);  % Index k specifying background region
-
-% Ensure there are no duplicate reference images
-%R=unique(R','rows')'; % comment this line if you run out of memory
-
-% Decompose B = R*R' using singular value or LU decomposition
-[L,U,p] = lu(R(k,:)'*R(k,:),'vector');       % LU decomposition
-
-for j=1:nimgs
-    b=R(k,:)'*A(k,j);
-    % Obtain coefficients c which minimise least-square residuals  
-    lower.LT = true; upper.UT = true;
-    c = linsolve(U,linsolve(L,b(p,:),lower),upper);
-
-    % Compute optimised reference image
-    optrefimages(:,:,j)=reshape(R*c,[ydim xdim]);
-end
\ No newline at end of file
diff --git a/IRF/remove_blob_manually.m b/IRF/remove_blob_manually.m
deleted file mode 100644
index e449df6..0000000
--- a/IRF/remove_blob_manually.m
+++ /dev/null
@@ -1,179 +0,0 @@
-% Demo to threshold an image to find regions (blobs).
-% Then let user point to a blob that you want to eliminate.
-
-clc;    % Clear the command window.
-close all;  % Close all figures (except those of imtool.)
-imtool close all;  % Close all imtool figures if you have the Image Processing Toolbox.
-clearvars;  % Erase all existing variables.
-workspace;  % Make sure the workspace panel is showing.
-format long g;
-format compact;
-fontSize = 20;
-
-% Check that user has the Image Processing Toolbox installed.
-hasIPT = license('test', 'image_toolbox');
-if ~hasIPT
-	% User does not have the toolbox installed.
-	message = sprintf('Sorry, but you do not seem to have the Image Processing Toolbox.\nDo you want to try to continue anyway?');
-	reply = questdlg(message, 'Toolbox missing', 'Yes', 'No', 'Yes');
-	if strcmpi(reply, 'No')
-		% User said No, so exit.
-		return;
-	end
-end
-baseFileName = 'coins.png'; % Default
-
-% % Read in a standard MATLAB gray scale demo image.
-% folder = fullfile(matlabroot, '\toolbox\images\imdemos');
-% button = menu('Use which demo image?', 'CameraMan', 'Moon', 'Eight', 'Coins', 'Pout');
-% if button == 1
-% 	baseFileName = 'cameraman.tif';
-% elseif button == 2
-% 	baseFileName = 'moon.tif';
-% elseif button == 3
-% 	baseFileName = 'coins.png';
-% else
-% 	baseFileName = 'pout.tif';
-% end
-
-% Read in a standard MATLAB gray scale demo image.
-folder = fullfile(matlabroot, '\toolbox\images\imdemos');
-% Get the full filename, with path prepended.
-fullFileName = fullfile(folder, baseFileName);
-% Check if file exists.
-if ~exist(fullFileName, 'file')
-	% File doesn't exist -- didn't find it there.  Check the search path for it.
-	fullFileName = baseFileName; % No path this time.
-	if ~exist(fullFileName, 'file')
-		% Still didn't find it.  Alert user.
-		errorMessage = sprintf('Error: %s does not exist in the search path folders.', fullFileName);
-		uiwait(warndlg(errorMessage));
-		return;
-	end
-end
-grayImage = imread(fullFileName);
-% Get the dimensions of the image.
-% numberOfColorBands should be = 1.
-[rows, columns, numberOfColorBands] = size(grayImage);
-if numberOfColorBands > 1
-	% It's not really gray scale like we expected - it's color.
-	% Convert it to gray scale by taking only the green channel.
-	grayImage = grayImage(:, :, 2); % Take green channel.
-end
-% Display the original gray scale image.
-subplot(2, 3, 1);
-imshow(grayImage, []);
-title('Original Grayscale Image', 'FontSize', fontSize);
-% Enlarge figure to full screen.
-set(gcf, 'Units', 'Normalized', 'OuterPosition', [0 0 1 1]);
-% Give a name to the title bar.
-set(gcf, 'Name', 'Demo by ImageAnalyst', 'NumberTitle', 'Off')
-
-% Let's compute and display the histogram.
-[pixelCount, grayLevels] = imhist(grayImage);
-subplot(2, 3, 2);
-bar(grayLevels, pixelCount);
-grid on;
-title('Histogram of original image', 'FontSize', fontSize);
-xlim([0 grayLevels(end)]); % Scale x axis manually.
-
-% Threshold the image.
-binaryImage = grayImage > 100;
-% Clean up a bit
-binaryImage = bwareaopen(binaryImage, 500);
-binaryImage = imfill(binaryImage, 'holes');
-% Display the binary image.
-subplot(2, 3, 3);
-imshow(binaryImage, []);
-title('Binary Image', 'FontSize', fontSize);
-
-doAnother = true;
-while doAnother
-	% Find pixels
-	[labeledImage, numberOfBlobs] = bwlabel(binaryImage);
-	measurements = regionprops(labeledImage, 'PixelIdxList', 'Centroid')
-	allCentroids = [measurements.Centroid];
-	centroidX = allCentroids(1:2:end);
-	centroidY = allCentroids(2:2:end);
-	% Plot the centroids over the blobs
-	hold on;
-	plot(centroidX, centroidY, 'bo', 'MarkerSize', 10);
-	axis on;
-	% Put text labels on them.
-	for k = 1 : numberOfBlobs
-		text(centroidX(k), centroidY(k)+10, num2str(k), 'Color', 'b', 'FontWeight', 'bold');
-	end
-	promptMessage = sprintf('On the binary image in the upper right,\nClick the region to remove,\nor Cancel to abort processing?');
-	titleBarCaption = 'Continue?';
-	subplot(2, 3, 3);
-	button = questdlg(promptMessage, titleBarCaption, 'Continue', 'Cancel', 'Continue');
-	if strcmpi(button, 'Cancel')
-		break;
-	end
-	[x,y] = ginput(1)
-	% Plot where they clicked.
-	plot(x, y, 'r+', 'MarkerSize', 20, 'LineWidth', 3);
-	
-	% Find which centroid this (x,y) is closest to
-	% First find out the distance from where user clicked to every other centroid.
-	xDistances = (centroidX - x);
-	yDistances = (centroidY - y);
-	distances = sqrt(xDistances .^ 2 + yDistances .^ 2);
-	% Find the closest one.
-	[minDistance, indexOfClosest] = min(distances)
-	% Plot an X over the closest blob.
-	plot(centroidX(indexOfClosest), centroidY(indexOfClosest), 'rx', 'MarkerSize', 40, 'LineWidth', 3);
-	% Draw a line between them.
-	line([x, centroidX(indexOfClosest)], [y, centroidY(indexOfClosest)], 'Color', 'r', 'LineWidth', 2);
-	
-	% Now remove this index.
-	keeperIndexes = 1 : numberOfBlobs; % All of them
-	keeperIndexes(indexOfClosest) = []; % Remove this particular blob from the list of blobs.
-	% Remove it from the labeled image.
-	newLabeledImage = ismember(labeledImage, keeperIndexes);
-	% Get new indexes in consequtive order since one if now missing.
-	newBinaryImage = newLabeledImage > 0; % All except selected blob.
-	% Display the binary image.
-	subplot(2, 3, 4);
-	imshow(newBinaryImage, []);
-	title('New Binary Image', 'FontSize', fontSize);
-	% Now make measurements all over again with the indicated blob removed (optional).
-	[labeledImage, numberOfBlobs] = bwlabel(binaryImage);
-	measurements = regionprops(labeledImage, 'Area');
-	
-	% Mask the image to make selected blob 0
-	% Get the selected blob alone
-	selectedBlob = binaryImage & ~newBinaryImage;
-	maskedImage1 = grayImage; % Initialize.
-	maskedImage1(selectedBlob) = 0;
-	% Display the masked image.
-	subplot(2, 3, 5);
-	imshow(maskedImage1, []);
-	title('Masked Image', 'FontSize', fontSize);
-	
-	% Fill the image with surrounding background.
-	% First enlarge blob
-	selectedBlob = imdilate(selectedBlob, ones(7));
-	% Now do the fill from the boundary.
-	maskedImage2 = roifill(grayImage, selectedBlob);
-	% Display the masked image.
-	subplot(2, 3, 6);
-	imshow(maskedImage2, []);
-	title('Filled Image', 'FontSize', fontSize);
-	
-	% If we've deleted the last blob, exit.
-	if numberOfBlobs <= 1
-		% Bail out if there are no more blobs.
-		break;
-	end
-	
-	cumulativeRemoval = true;
-	if cumulativeRemoval
-		% If you want the removal to be cumulative, set grayImage to be maskedImage2 or maskedImage1.
-		% Otherwise comment out the line below to start from the original gray image every time.
-		grayImage = maskedImage2;
-		binaryImage = newBinaryImage;
-	end
-end
-
-
diff --git a/IRF/Analysis.m b/Imaging Response Function Extractor/extractIRF.m
similarity index 98%
rename from IRF/Analysis.m
rename to Imaging Response Function Extractor/extractIRF.m
index 7aa3072..b98c5b0 100644
--- a/IRF/Analysis.m
+++ b/Imaging Response Function Extractor/extractIRF.m	
@@ -42,7 +42,7 @@ for k = 1 : length(files)
 end
 %% Fringe removal
 
-optrefimages               = fringeremoval(absimages, refimages);
+optrefimages               = removefringesInImage(absimages, refimages);
 absimages_fringe_removed   = absimages(:, :, :) - optrefimages(:, :, :);
 
 nimgs    = size(absimages_fringe_removed,3);
@@ -397,8 +397,8 @@ function [RadialResponseFunc] = RadialImagingResponseFunction(C, k, kmax)
     RadialResponseFunc = C(6)*1/2*A.*RadialResponseFunc;
 end
 
-function [optrefimages] = fringeremoval(absimages, refimages, bgmask)
-    % FRINGEREMOVAL - Fringe removal and noise reduction from absorption images.
+function [optrefimages] = removefringesInImage(absimages, refimages, bgmask)
+    % removefringesInImage - Fringe removal and noise reduction from absorption images.
     % Creates an optimal reference image for each absorption image in a set as 
     % a linear combination of reference images, with coefficients chosen to 
     % minimize the least-squares residuals between each absorption image and 
@@ -409,7 +409,7 @@ function [optrefimages] = fringeremoval(absimages, refimages, bgmask)
     % detection of small atom numbers through image processing, arXiv:1007.2136 (2010).
     %
     % Syntax:
-    %    [optrefimages] = fringeremoval(absimages,refimages,bgmask);
+    %    [optrefimages] = removefringesInImage(absimages,refimages,bgmask);
     %
     % Required inputs:
     %    absimages    - Absorption image data,