classdef MOTSimulator < handle & matlab.mixin.Copyable

    properties (Access = public)
        SimulationMode; % MOT type
        TimeStep;
        SimulationTime;
        NumberOfAtoms;
        
        ParticleDynamicalQuantities;
        
        NozzleLength;
        NozzleRadius;
        Beta;
        ApertureCut;
        OvenDistance;
        OvenTemperature;
        MagneticGradient;
        NozzleExitDivergence;
        MOTExitDivergence;
        MOTDistance;
        
        BluePower;
        BlueDetuning;
        BlueBeamRadius;
        BlueBeamWaist;
        BlueWaveNumber;
        BlueSaturationIntensity;
        
        OrangePower;
        OrangeDetuning;
        OrangeBeamRadius;
        OrangeBeamWaist;
        OrangeWaveNumber;
        OrangeSaturationIntensity;
        
        CoolingBeamPower;
        CoolingBeamWaveNumber;
        CoolingBeamLinewidth;
        CoolingBeamDetuning;
        CoolingBeamRadius;
        CoolingBeamWaist;
        CoolingBeamSaturationIntensity;
        
        SidebandPower;
        SidebandDetuning;
        SidebandBeamRadius;
        SidebandBeamWaist;
        SidebandBeamSaturationIntensity;
                
        PushBeamPower;
        PushBeamWaveNumber;
        PushBeamLinewidth;
        PushBeamDetuning;
        PushBeamRadius;
        PushBeamWaist;
        PushBeamDistance;
        DistanceBetweenPushBeamAnd3DMOTCenter;
        PushBeamSaturationIntensity;
        
        ZeemanSlowerBeamPower;
        ZeemanSlowerBeamDetuning;
        ZeemanSlowerBeamRadius;
        ZeemanSlowerBeamWaist;
        ZeemanSlowerBeamSaturationIntensity;
        
        TotalPower;
        LandegFactor;
        MagneticSubLevel;
        
        CaptureVelocity;
        VelocityCutoff;
        ClausingFactor;
        ReducedClausingFactor;
        ReducedFlux;
        TimeSpentInInteractionRegion;
        
        %Flags
        SpontaneousEmission;
        Sideband;
        PushBeam;
        ZeemanSlowerBeam;
        Gravity;
        BackgroundCollision;
        
        DebugMode;
        DoSave;
        SaveDirectory;
        
        Results;
    end
    
    properties (SetAccess = private, GetAccess = public)
        SimulationParameters
    end
    
    properties (Dependent, SetAccess = private)
        CoolingBeamSaturationParameter;
        SidebandSaturationParameter;
        PushBeamSaturationParameter;
        ZeemanSlowerBeamSaturationParameter;
        OvenTemperatureinKelvin;
        AverageVelocity;
        AtomicBeamDensity;
        MeanFreePath;
        CollisionTime;
    end
        
    methods
        function s  = MOTSimulator(varargin)
            
            p = inputParser;
            p.KeepUnmatched = true;
            addParameter(p, 'SimulationMode', '2D',... 
                @(x) any(strcmpi(x,{'2D','3D'})));
            addParameter(p, 'TimeStep',       10e-06,...
                @(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));  
            addParameter(p, 'SimulationTime',       3e-03,...
                @(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));  
            addParameter(p, 'SpontaneousEmission',  false,...
                @islogical);
            addParameter(p, 'Sideband',             false,...
                @islogical);
            addParameter(p, 'PushBeam',             false,...
                @islogical);
            addParameter(p, 'ZeemanSlowerBeam',     false,...
                @islogical);
            addParameter(p, 'Gravity',              false,...
                @islogical);
            addParameter(p, 'BackgroundCollision',  false,...
                @islogical);
            addParameter(p, 'DebugMode',            false,...
                @islogical);
            addParameter(p, 'SaveData',             false,...
                @islogical);
            addParameter(p, 'SaveDirectory',          pwd,...
                @ischar);
            
            p.parse(varargin{:});
            
            s.SimulationMode       = p.Results.SimulationMode;
            
            s.TimeStep             = p.Results.TimeStep;
            s.SimulationTime       = p.Results.SimulationTime;  
            
            s.SpontaneousEmission  = p.Results.SpontaneousEmission;
            s.Sideband             = p.Results.Sideband;
            s.PushBeam             = p.Results.PushBeam;
            s.ZeemanSlowerBeam     = p.Results.ZeemanSlowerBeam;
            s.Gravity              = p.Results.Gravity;
            s.BackgroundCollision  = p.Results.BackgroundCollision;
            
            s.DebugMode            = p.Results.DebugMode;
            s.DoSave               = p.Results.SaveData;
            s.SaveDirectory        = p.Results.SaveDirectory;
            
            
            s.reinitializeSimulator();
            
            poolobj = gcp('nocreate'); % Check if pool is open
            if isempty(poolobj)
                parpool;
            end
            
        end
    end % - lifecycle    

    methods
        function set.TimeStep(this, val)
            assert(val > 1e-06, 'Not time efficient to compute for time steps smaller than 1 microsecond!');
            this.TimeStep = val;
        end
        function ret = get.TimeStep(this)
            ret = this.TimeStep;
        end
        function set.SimulationTime(this, val)
%             assert(val <= 5e-03, 'Not time efficient to compute for time spans longer than 5 milliseconds!');
            this.SimulationTime = val;
        end
        function ret = get.SimulationTime(this)
            ret = this.SimulationTime;
        end
        function set.NumberOfAtoms(this, val)
            assert(val <= 20000, 'Not time efficient to compute for atom numbers larger than 20,000!');
            this.NumberOfAtoms = val;
        end
        function ret = get.NumberOfAtoms(this)
            ret = this.NumberOfAtoms;
        end
        
        function set.ParticleDynamicalQuantities(this,val)
            this.ParticleDynamicalQuantities = val;
        end
        function ret = get.ParticleDynamicalQuantities(this)
            ret = this.ParticleDynamicalQuantities;
        end
        
        function set.NozzleLength(this,val)
            this.NozzleLength = val;
        end
        function ret = get.NozzleLength(this)
            ret = this.NozzleLength;
        end
        function set.NozzleRadius(this,val)
            this.NozzleRadius = val;
        end
        function ret = get.NozzleRadius(this)
            ret = this.NozzleRadius;
        end
        function set.Beta(this,val)
            this.Beta = val;
        end
        function ret = get.Beta(this)
            ret = this.Beta;
        end
        function set.ApertureCut(this,val)
            this.ApertureCut = val;
        end
        function ret = get.ApertureCut(this)
            ret = this.ApertureCut;
        end
        function set.OvenDistance(this,val)
            this.OvenDistance = val;
        end
        function ret = get.OvenDistance(this)
            ret = this.OvenDistance;
        end
        function set.OvenTemperature(this,val)
            this.OvenTemperature = val;
        end
        function ret = get.OvenTemperature(this)
            ret = this.OvenTemperature;
        end
        function set.MagneticGradient(this,val)
            this.MagneticGradient = val;
        end
        function ret = get.MagneticGradient(this)
            ret = this.MagneticGradient;
        end
        function set.NozzleExitDivergence(this,val)
            this.NozzleExitDivergence = val;
        end
        function ret = get.NozzleExitDivergence(this)
            ret = this.NozzleExitDivergence;
        end
        function set.MOTExitDivergence(this,val)
            this.MOTExitDivergence = val;
        end
        function ret = get.MOTExitDivergence(this)
            ret = this.MOTExitDivergence;
        end
        function set.MOTDistance(this,val)
            this.MOTDistance = val;
        end
        function ret = get.MOTDistance(this)
            ret = this.MOTDistance;
        end
        
        function set.BluePower(this,val)
            this.BluePower = val;
        end
        function ret = get.BluePower(this)
            ret = this.BluePower;
        end
        function set.BlueDetuning(this, val)
            this.BlueDetuning = val;
        end
        function ret = get.BlueDetuning(this)
            ret = this.BlueDetuning;
        end
        function set.BlueBeamRadius(this, val)
            this.BlueBeamRadius = val;
        end
        function ret = get.BlueBeamRadius(this)
            ret = this.BlueBeamRadius;
        end
        function set.BlueBeamWaist(this, val)
            this.BlueBeamWaist = val;
        end
        function ret = get.BlueBeamWaist(this)
            ret = this.BlueBeamWaist;
        end
        function set.BlueWaveNumber(this, val)
            this.BlueWaveNumber = val;
        end
        function ret = get.BlueWaveNumber(this)
            ret = this.BlueWaveNumber;
        end
        function set.BlueSaturationIntensity(this, val)
            this.BlueSaturationIntensity = val;
        end
        function ret = get.BlueSaturationIntensity(this)
            ret = this.BlueSaturationIntensity;
        end
        
        function set.OrangePower(this,val)
            this.OrangePower = val;
        end
        function ret = get.OrangePower(this)
            ret = this.OrangePower;
        end
        function set.OrangeDetuning(this, val)
            this.OrangeDetuning = val;
        end
        function ret = get.OrangeDetuning(this)
            ret = this.OrangeDetuning;
        end
        function set.OrangeBeamRadius(this, val)
            this.OrangeBeamRadius = val;
        end
        function ret = get.OrangeBeamRadius(this)
            ret = this.OrangeBeamRadius;
        end
        function set.OrangeBeamWaist(this, val)
            this.OrangeBeamWaist = val;
        end
        function ret = get.OrangeBeamWaist(this)
            ret = this.OrangeBeamWaist;
        end
        function set.OrangeWaveNumber(this, val)
            this.OrangeWaveNumber = val;
        end
        function ret = get.OrangeWaveNumber(this)
            ret = this.OrangeWaveNumber;
        end
        function set.OrangeSaturationIntensity(this, val)
            this.OrangeSaturationIntensity = val;
        end
        function ret = get.OrangeSaturationIntensity(this)
            ret = this.OrangeSaturationIntensity;
        end
        
        function set.CoolingBeamPower(this,val)
            this.CoolingBeamPower = val;
        end
        function ret = get.CoolingBeamPower(this)
            ret = this.CoolingBeamPower;
        end
        function set.CoolingBeamDetuning(this, val)
            this.CoolingBeamDetuning = val;
        end
        function ret = get.CoolingBeamDetuning(this)
            ret = this.CoolingBeamDetuning;
        end
        function set.CoolingBeamRadius(this, val)
            this.CoolingBeamRadius = val;
        end
        function ret = get.CoolingBeamRadius(this)
            ret = this.CoolingBeamRadius;
        end
        function set.CoolingBeamWaist(this, val)
            this.CoolingBeamWaist = val;
        end
        function ret = get.CoolingBeamWaist(this)
            ret = this.CoolingBeamWaist;
        end
        function set.CoolingBeamWaveNumber(this, val)
            this.CoolingBeamWaveNumber = val;
        end
        function ret = get.CoolingBeamWaveNumber(this)
            ret = this.CoolingBeamWaveNumber;
        end
        function set.CoolingBeamLinewidth(this, val)
            this.CoolingBeamLinewidth = val;
        end
        function ret = get.CoolingBeamLinewidth(this)
            ret = this.CoolingBeamLinewidth;
        end
        function set.CoolingBeamSaturationIntensity(this, val)
            this.CoolingBeamSaturationIntensity = val;
        end
        function ret = get.CoolingBeamSaturationIntensity(this)
            ret = this.CoolingBeamSaturationIntensity;
        end
        
        function set.SidebandPower(this,val)
            this.SidebandPower = val;
        end
        function ret = get.SidebandPower(this)
            ret = this.SidebandPower;
        end
        function set.SidebandDetuning(this, val)
            this.SidebandDetuning = val;
        end
        function ret = get.SidebandDetuning(this)
            ret = this.SidebandDetuning;
        end
        function set.SidebandBeamRadius(this, val)
            this.SidebandBeamRadius = val;
        end
        function ret = get.SidebandBeamRadius(this)
            ret = this.SidebandBeamRadius;
        end
        function set.SidebandBeamWaist(this, val)
            this.SidebandBeamWaist = val;
        end
        function ret = get.SidebandBeamWaist(this)
            ret = this.SidebandBeamWaist;
        end
        function set.SidebandBeamSaturationIntensity(this, val)
            this.SidebandBeamSaturationIntensity = val;
        end
        function ret = get.SidebandBeamSaturationIntensity(this)
            ret = this.SidebandBeamSaturationIntensity;
        end
        
        function set.PushBeamPower(this,val)
            this.PushBeamPower = val;
        end
        function ret = get.PushBeamPower(this)
            ret = this.PushBeamPower;
        end
        function set.PushBeamDetuning(this, val)
            this.PushBeamDetuning = val;
        end
        function ret = get.PushBeamDetuning(this)
            ret = this.PushBeamDetuning;
        end
        function set.PushBeamRadius(this, val)
            this.PushBeamRadius = val;
        end
        function ret = get.PushBeamRadius(this)
            ret = this.PushBeamRadius;
        end
        function set.PushBeamWaist(this, val)
            this.PushBeamWaist = val;
        end
        function ret = get.PushBeamWaist(this)
            ret = this.PushBeamWaist;
        end
        function set.PushBeamWaveNumber(this, val)
            this.PushBeamWaveNumber= val;
        end
        function ret = get.PushBeamWaveNumber(this)
            ret = this.PushBeamWaveNumber;
        end
        function set.PushBeamLinewidth(this, val)
            this.PushBeamLinewidth = val;
        end
        function ret = get.PushBeamLinewidth(this)
            ret = this.PushBeamLinewidth;
        end
        function set.PushBeamDistance(this, val)
            this.PushBeamDistance = val;
        end
        function ret = get.PushBeamDistance(this)
            ret = this.PushBeamDistance;
        end
        function set.DistanceBetweenPushBeamAnd3DMOTCenter(this, val)
            this.DistanceBetweenPushBeamAnd3DMOTCenter = val;
        end
        function ret = get.DistanceBetweenPushBeamAnd3DMOTCenter(this)
            ret = this.DistanceBetweenPushBeamAnd3DMOTCenter;
        end
        function set.PushBeamSaturationIntensity(this, val)
            this.PushBeamSaturationIntensity = val;
        end
        function ret = get.PushBeamSaturationIntensity(this)
            ret = this.PushBeamSaturationIntensity;
        end
        
        function set.ZeemanSlowerBeamPower(this,val)
            this.ZeemanSlowerBeamPower = val;
        end
        function ret = get.ZeemanSlowerBeamPower(this)
            ret = this.ZeemanSlowerBeamPower;
        end
        function set.ZeemanSlowerBeamDetuning(this, val)
            this.ZeemanSlowerBeamDetuning = val;
        end
        function ret = get.ZeemanSlowerBeamDetuning(this)
            ret = this.ZeemanSlowerBeamDetuning;
        end
        function set.ZeemanSlowerBeamRadius(this, val)
            this.ZeemanSlowerBeamRadius = val;
        end
        function ret = get.ZeemanSlowerBeamRadius(this)
            ret = this.ZeemanSlowerBeamRadius;
        end
        function set.ZeemanSlowerBeamWaist(this, val)
            this.ZeemanSlowerBeamWaist = val;
        end
        function ret = get.ZeemanSlowerBeamWaist(this)
            ret = this.ZeemanSlowerBeamWaist;
        end
        function set.ZeemanSlowerBeamSaturationIntensity(this, val)
            this.ZeemanSlowerBeamSaturationIntensity = val;
        end
        function ret = get.ZeemanSlowerBeamSaturationIntensity(this)
            ret = this.ZeemanSlowerBeamSaturationIntensity;
        end
        
        function set.TotalPower(this,val)
            this.TotalPower = val;
        end
        function ret = get.TotalPower(this)
            ret = this.TotalPower;
        end
        function set.LandegFactor(this,val)
            this.LandegFactor = val;
        end
        function ret = get.LandegFactor(this)
            ret = this.LandegFactor;
        end
        function set.MagneticSubLevel(this,val)
            this.MagneticSubLevel = val;
        end
        function ret = get.MagneticSubLevel(this)
            ret = this.MagneticSubLevel;
        end
        
        function set.CaptureVelocity(this,val)
            this.CaptureVelocity = val;
        end
        function ret = get.CaptureVelocity(this)
            ret = this.CaptureVelocity;
        end
        function set.VelocityCutoff(this,val)
            this.VelocityCutoff = val;
        end
        function ret = get.VelocityCutoff(this)
            ret = this.VelocityCutoff;
        end
        function set.ClausingFactor(this,val)
            this.ClausingFactor = val;
        end
        function ret = get.ClausingFactor(this)
            ret = this.ClausingFactor;
        end
        function set.ReducedClausingFactor(this,val)
            this.ReducedClausingFactor = val;
        end
        function ret = get.ReducedClausingFactor(this)
            ret = this.ReducedClausingFactor;
        end
        function set.ReducedFlux(this,val)
            this.ReducedFlux = val;
        end
        function ret = get.ReducedFlux(this)
            ret = this.ReducedFlux;
        end
        function set.TimeSpentInInteractionRegion(this,val)
            this.TimeSpentInInteractionRegion = val;
        end
        function ret = get.TimeSpentInInteractionRegion(this)
            ret = this.TimeSpentInInteractionRegion;
        end
        
        function set.DebugMode(this, val)
            this.DebugMode = val;
        end
        function ret = get.DebugMode(this)
            ret = this.DebugMode;
        end
        function set.DoSave(this, val)
            this.DoSave = val;
        end
        function ret = get.DoSave(this)
            ret = this.DoSave;
        end
        function set.SaveDirectory(this, val)
            this.SaveDirectory = val;
        end
        function ret = get.SaveDirectory(this)
            ret = this.SaveDirectory;
        end
        
        function set.Results(this, val)
            this.Results = val;
        end
        function ret = get.Results(this)
            ret = this.Results;
        end
        
    end % - setters and getters
    
    methods
        function ret       = get.CoolingBeamSaturationParameter(this)
            ret = 0.1 * (4 * this.CoolingBeamPower) / (pi*this.CoolingBeamWaist^2 * this.CoolingBeamSaturationIntensity); % two beams are reflected
        end
        
        function ret       = get.SidebandSaturationParameter(this)
            ret = 0.1 * (4 * this.SidebandPower) / (pi*this.SidebandBeamWaist^2 * this.SidebandBeamSaturationIntensity);
        end
        
        function ret       = get.PushBeamSaturationParameter(this)
            ret = 0.1 * this.PushBeamPower/(pi * this.PushBeamWaist^2 * this.PushBeamSaturationIntensity);
        end
        
        function ret       = get.ZeemanSlowerBeamSaturationParameter(this)
            ret = 0.1 * this.ZeemanSlowerBeamPower / (pi * this.ZeemanSlowerBeamWaist^2 * this.ZeemanSlowerBeamSaturationIntensity);
        end
        
        function ret       = get.OvenTemperatureinKelvin(this)
            ret = this.OvenTemperature + Helper.PhysicsConstants.ZeroKelvin;
        end
        
        function ret       = get.AverageVelocity(this)
            %See Background collision probability section in Barbiero
            ret = sqrt((8 * pi *Helper.PhysicsConstants.BoltzmannConstant*this.OvenTemperatureinKelvin)/ (9 * Helper.PhysicsConstants.Dy164Mass));
        end
        
        function ret       = get.AtomicBeamDensity(this)
           %See Background collision probability section in Barbiero 
           ret  = this.calculateFreeMolecularRegimeFlux() / (this.AverageVelocity * pi * (this.Beta*this.NozzleLength/2)^2);  
        end
            
        function ret       = get.MeanFreePath(this)
            % Cross section  = pi ( 2 * Van-der-waals radius of Dy)^2;
            % Van-der-waals radius of Dy = 281e-12 
            %See Expected atomic flux section and Background collision probability section in Barbiero
            ret = 1/(sqrt(2) * ( pi * (2*281e-12)^2) * this.AtomicBeamDensity);  
        end
        
        function ret       = get.CollisionTime(this)
            ret = 3 * this.MeanFreePath/this.AverageVelocity; %See Background collision probability section in Barbiero 
        end
        
    end % - getters for dependent properties
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    methods(Access = protected)
        function cp = copyElement(this)
            % Shallow copy object
            cp = copyElement@matlab.mixin.Copyable(this);
            
            % Forces the setter to redefine the function handles to the new copied object
            % cp.potentialType = this.potentialType;
            
            pl = properties(this);
            for k = 1:length(pl)
                sc = superclasses(this.(pl{k}));
                if any(contains(sc,{'matlab.mixin.Copyable'}))
                    cp.(pl{k}) = this.(pl{k}).copy();
                end
            end
        end
    end
    
    methods (Static)
        
        % Creates an Instance of Class, ensures singleton behaviour (that there
        % can only be one Instance of this class
        function singleObj = getInstance(varargin)
            % Creates an Instance of Class, ensures singleton behaviour
            persistent localObj;
            if isempty(localObj) || ~isvalid(localObj)
                localObj =  MOTSimulator(varargin{:});
            end
            singleObj = localObj;
        end
    end
    
end