classdef MOTSimulator < handle & matlab.mixin.Copyable

    properties (Access = public)
        SimulationMode; % MOT type
        TimeStep;
        SimulationTime;
        NumberOfAtoms;
        
        InitialPositions;
        InitialVelocities;
        
        NozzleLength;
        NozzleRadius;
        Beta;
        ApertureCut;
        OvenDistance;
        OvenTemperature;
        MagneticGradient;
        NozzleExitDivergence;
        MOTExitDivergence;
        MOTDistance;
        
        BluePower;
        BlueDetuning;
        BlueBeamRadius;
        BlueBeamWaist;
        BlueWaveVector;
        BlueSaturationIntensity;
        
        OrangePower;
        OrangeDetuning;
        OrangeBeamRadius;
        OrangeBeamWaist;
        OrangeWaveVector;
        OrangeSaturationIntensity;
        
        CoolingBeamPower;
        CoolingBeamWaveVector;
        CoolingBeamLinewidth;
        CoolingBeamDetuning;
        CoolingBeamRadius;
        CoolingBeamWaist;
        CoolingBeamSaturationIntensity;
        
        SidebandPower;
        SidebandDetuning;
        SidebandBeamRadius;
        SidebandBeamWaist;
        SidebandBeamSaturationIntensity;
                
        PushBeamPower;
        PushBeamWaveVector;
        PushBeamLinewidth;
        PushBeamDetuning;
        PushBeamRadius;
        PushBeamWaist;
        PushBeamDistance;
        DistanceBetweenPushBeamAnd3DMOTCenter;
        PushBeamSaturationIntensity;
        
        ZeemanSlowerBeamPower;
        ZeemanSlowerBeamDetuning;
        ZeemanSlowerBeamRadius;
        ZeemanSlowerBeamWaist;
        ZeemanSlowerBeamSaturationIntensity;
        
        TotalPower;
        LandegFactor;
        MagneticSubLevel;
        
        CaptureVelocity;
        VelocityCutoff;
        ClausingFactor;
        ThetaArray;
        AngularDistribution;
        NormalizationConstantForAngularDistribution;
        
        %Flags
        SpontaneousEmission;
        Sideband;
        ZeemanSlowerBeam;
        Gravity;
        AtomicBeamCollision;
        
        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, 'ZeemanSlowerBeam',     false,...
                @islogical);
            addParameter(p, 'Gravity',              false,...
                @islogical);
            addParameter(p, 'AtomicBeamCollision',  true,...
                @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.ZeemanSlowerBeam     = p.Results.ZeemanSlowerBeam;
            s.Gravity              = p.Results.Gravity;
            s.AtomicBeamCollision  = p.Results.AtomicBeamCollision;
            
            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 <= 10000, 'Not time efficient to compute for atom numbers larger than 10,000!');
            this.NumberOfAtoms = val;
        end
        function ret = get.NumberOfAtoms(this)
            ret = this.NumberOfAtoms;
        end
        
        function set.InitialPositions(this,val)
            this.InitialPositions = val;
        end
        function ret = get.InitialPositions(this)
            ret = this.InitialPositions;
        end
        function set.InitialVelocities(this,val)
            this.InitialVelocities = val;
        end
        function ret = get.InitialVelocities(this)
            ret = this.InitialVelocities;
        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.BlueWaveVector(this, val)
            this.BlueWaveVector = val;
        end
        function ret = get.BlueWaveVector(this)
            ret = this.BlueWaveVector;
        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.OrangeWaveVector(this, val)
            this.OrangeWaveVector = val;
        end
        function ret = get.OrangeWaveVector(this)
            ret = this.OrangeWaveVector;
        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.CoolingBeamWaveVector(this, val)
            this.CoolingBeamWaveVector = val;
        end
        function ret = get.CoolingBeamWaveVector(this)
            ret = this.CoolingBeamWaveVector;
        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.PushBeamWaveVector(this, val)
            this.PushBeamWaveVector = val;
        end
        function ret = get.PushBeamWaveVector(this)
            ret = this.PushBeamWaveVector;
        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.AngularDistribution(this,val)
            this.AngularDistribution = val;
        end
        function ret = get.AngularDistribution(this)
            ret = this.AngularDistribution;
        end
        function set.ThetaArray(this,val)
            this.ThetaArray = val;
        end
        function ret = get.ThetaArray(this)
            ret = this.ThetaArray;
        end
        function set.NormalizationConstantForAngularDistribution(this,val)
            this.NormalizationConstantForAngularDistribution = val;
        end
        function ret = get.NormalizationConstantForAngularDistribution(this)
            ret = this.NormalizationConstantForAngularDistribution;
        end
        
        function set.AtomicBeamCollision(this,val)
            this.AtomicBeamCollision=val;
        end
        
        function ret=get.AtomicBeamCollision(this)
            ret=this.AtomicBeamCollision;
        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 = 4* this.CoolingBeamPower/(pi*this.CoolingBeamWaist^2)/this.CoolingBeamSaturationIntensity/10; % two beams are reflected
        end
        
        function ret       = get.SidebandSaturationParameter(this)
            ret = 4*this.SidebandPower/(pi*this.SidebandBeamWaist^2)/this.SidebandBeamSaturationIntensity/10;
        end
        
        function ret       = get.PushBeamSaturationParameter(this)
            ret = this.PushBeamPower/(pi*this.PushBeamWaist^2)/this.PushBeamSaturationIntensity/10;
        end
        
        function ret       = get.ZeemanSlowerBeamSaturationParameter(this)
            ret = this.ZeemanSlowerBeamPower/(pi*this.ZeemanSlowerBeamWaist^2)/this.ZeemanSlowerBeamSaturationIntensity/10;
        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)/9) * ((Helper.PhysicsConstants.BoltzmannConstant * this.OvenTemperatureinKelvin)/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
    
end