classdef Potentials < handle & matlab.mixin.Copyable

    properties (Access = private)
        
        PotentialDefaults     = struct('TrapPotentialType',           'Harmonic', ...
                                       'TrapFrequencies',        100 * ones(1,3), ...
                                       'TrapDepth',                           10, ...
                                       'BoxSize',                              5);
    end

    properties (Access = public)
        TrapPotentialType;
        TrapFrequencies;
        TrapDepth;
        BoxSize;
        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, 'TrapDepth',     this.PotentialDefaults.TrapDepth,...
                @(x) assert(isnumeric(x) && isscalar(x) && (x > 0)));
            addParameter(p, 'BoxSize', this.PotentialDefaults.BoxSize,...
                @(x) assert(isnumeric(x) && isscalar(x) && (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.TrapDepth            = p.Results.TrapDepth;
            this.BoxSize              = p.Results.BoxSize;
            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);
                    z        = linspace(-0.5*this.SimulationParameters.Lz,0.5*this.SimulationParameters.Lz-this.SimulationParameters.Lz/this.SimulationParameters.Nz,this.SimulationParameters.Nz);
                    [X,Y,Z]  = ndgrid(x,y,z);
                    ret      = 0.0 * (X+Y+Z);
                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);
                    z        = linspace(-0.5*this.SimulationParameters.Lz,0.5*this.SimulationParameters.Lz-this.SimulationParameters.Lz/this.SimulationParameters.Nz,this.SimulationParameters.Nz);
                    [X,Y,Z]  = ndgrid(x,y,z);
                    ret      = 0.5*(this.SimulationParameters.gx.*X.^2+this.SimulationParameters.gy.*Y.^2+this.SimulationParameters.gz*Z.^2); 
                case 'SquareBox'
                    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);
                    z        = linspace(-0.5*this.SimulationParameters.Lz,0.5*this.SimulationParameters.Lz-this.SimulationParameters.Lz/this.SimulationParameters.Nz,this.SimulationParameters.Nz);
                    [X,Y,Z]  = ndgrid(x,y,z);
                    
                    heaviside_X = this.custom_heaviside(-abs(X) + this.SimulationParameters.boxsize/2);
                    heaviside_Y = this.custom_heaviside(-abs(Y) + this.SimulationParameters.boxsize/2);
                    heaviside_Z = this.custom_heaviside(-abs(Z) + this.SimulationParameters.boxsize/2);
                    
                    ret = this.SimulationParameters.boxdepth * (1 - heaviside_X .* heaviside_Y .* heaviside_Z);
                case 'RoundBox'
                    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);
                    z        = linspace(-0.5*this.SimulationParameters.Lz,0.5*this.SimulationParameters.Lz-this.SimulationParameters.Lz/this.SimulationParameters.Nz,this.SimulationParameters.Nz);
                    [X,Y,Z]  = ndgrid(x,y,z);
                    
                    ret = this.SimulationParameters.boxdepth * (1 - this.custom_heaviside((this.SimulationParameters.boxsize/2)^2 - X.^2 - Y.^2 - Z.^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','SquareBox','RoundBox'})), '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