classdef Potentials < handle & matlab.mixin.Copyable

    properties (Access = private)
        
        PotentialDefaults     = struct('TrapPotentialType',           'Harmonic', ...
                                       'TrapFrequencies',        100 * ones(1,2));                                       
    end

    properties (Access = public)
        TrapPotentialType;
        TrapFrequencies;
        NumberOfGridPoints;
        Dimensions;
        SimulationParameters;
    end
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %- Methods
    
    methods
        function this  = Potentials(varargin)
            
            p = inputParser;
            p.KeepUnmatched = true;
            
            addParameter(p, 'TrapPotentialType', this.PotentialDefaults.TrapPotentialType,... 
                @(x) assert(any(strcmpi(x,{'None','Harmonic','SquareBox','RoundBox'}))));
            addParameter(p, 'TrapFrequencies',   this.PotentialDefaults.TrapFrequencies,...
                @(x) assert(isnumeric(x) && isvector(x) && all(x > 0)));
            addParameter(p, 'NumberOfGridPoints',    128 * ones(1,3),...
                @(x) assert(isnumeric(x) && isvector(x) && all(x > 0)));
            addParameter(p, 'Dimensions',        10 * ones(1,3),...
                @(x) assert(isnumeric(x) && isvector(x) && all(x > 0)));
            
            p.parse(varargin{:});

            this.TrapPotentialType    = p.Results.TrapPotentialType;
            this.TrapFrequencies      = p.Results.TrapFrequencies;
            this.NumberOfGridPoints   = p.Results.NumberOfGridPoints;
            this.Dimensions           = p.Results.Dimensions;

            this.SimulationParameters = this.setupParameters();
        end
        
        function [ret] = trap(this)
            format long
            switch this.TrapPotentialType
                case 'None'
                    x        = linspace(-0.5*this.SimulationParameters.Lx,0.5*this.SimulationParameters.Lx-this.SimulationParameters.Lx/this.SimulationParameters.Nx,this.SimulationParameters.Nx);
                    y        = linspace(-0.5*this.SimulationParameters.Ly,0.5*this.SimulationParameters.Ly-this.SimulationParameters.Ly/this.SimulationParameters.Ny,this.SimulationParameters.Ny);
                    [X,Y]    = ndgrid(x,y);

                    ret      = 0.0 * (X+Y);
                case 'Harmonic'
                    x        = linspace(-0.5*this.SimulationParameters.Lx,0.5*this.SimulationParameters.Lx-this.SimulationParameters.Lx/this.SimulationParameters.Nx,this.SimulationParameters.Nx);
                    y        = linspace(-0.5*this.SimulationParameters.Ly,0.5*this.SimulationParameters.Ly-this.SimulationParameters.Ly/this.SimulationParameters.Ny,this.SimulationParameters.Ny);
                    [X,Y]    = ndgrid(x,y);

                    ret      = 0.5*(this.SimulationParameters.gx.*X.^2+this.SimulationParameters.gy.*Y.^2); 
            end
        end
    
        function restoreDefaults(this)
            this.TrapPotentialType          = this.PotentialDefaults.TrapPotentialType;
            this.TrapFrequencies            = this.PotentialDefaults.TrapFrequencies;
        end
    end %

    methods
        function set.TrapFrequencies(this, val)
            assert(isnumeric(val) && isvector(val) && all(val > 0), 'Incorrectly defined trap frequencies!');
            this.TrapFrequencies = val;
        end
        function ret = get.TrapFrequencies(this)
            ret = this.TrapFrequencies;
        end
        function set.TrapPotentialType(this, val)
            assert(any(strcmpi(val,{'None','Harmonic'})), 'Trap potential of the specified type cannot be generated!');
            this.TrapPotentialType = val;
        end
        function ret = get.TrapPotentialType(this)
            ret = this.TrapPotentialType;
        end
        function set.NumberOfGridPoints(this, val)
            assert(isnumeric(val) && isvector(val) && all(val > 0), 'Incorrectly defined grid!');
            this.NumberOfGridPoints = val;
        end
        function ret = get.NumberOfGridPoints(this)
            ret = this.NumberOfGridPoints;
        end
        function set.Dimensions(this, val)
            assert(isnumeric(val) && isvector(val) && all(val > 0), 'Incorrectly defined dimensions!');
            this.Dimensions = val;
        end
        function ret = get.Dimensions(this)
            ret = this.Dimensions;
        end
    end % - setters and getters
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %- Methods
    
    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
            
            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)

        function y = custom_heaviside(x)
            % This function computes the Heaviside step function with a custom value for Heaviside(0).
            % x: input array
            % H0: value for Heaviside(0)
            
            % Use MATLAB's built-in heaviside function
            y = heaviside(x);
            
            % Replace the default value for Heaviside(0) with the custom value H0
            y(x == 0) = 1;
        end

        % 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 =  Simulator.Potentials(varargin{:});
            end
            singleObj = localObj;
        end
    end
    
end