classdef HW8_wave_decay < handle
    %this class takes some info on a wave-packet and simulates its dacay
    %over time
    
    properties
        N %number of samples
        DeltaX %the length of the interval in spacial domain in Angstrom
        X %spacial smples
        Sigma %the half width of the Gaussian wave-packet
        M %mass of the particle
        T %length of the time interval
        A %Gaussian wave-packet
        SigmaX %the variance of the spreading function
        Tvec
    end
    
    methods
        function obj = gauss_wave(obj)
            %this function generate a Gaussian wave-packet
            obj.X = linspace(-obj.DeltaX/2,obj.DeltaX/2,obj.N); %Angstrom
            DX = obj.X(2)-obj.X(1);
            obj.A = normpdf(obj.X,0,obj.Sigma);
            obj.A = obj.A/sqrt(obj.A*obj.A'*DX); %Normalize 
            %figure
            %plot(obj.X,obj.A)
            %xlabel('X');ylabel('\Psi(x)')
            
        end
        
        function obj = dec_wave(obj)
           %this function simulate a decaying wave-packet
           Y=fftshift(fft(obj.A));
           DX = obj.X(2)-obj.X(1);
           K=(0:obj.N-1);
           K(K>obj.N/2)=K(K>obj.N/2)-obj.N;
           K=fftshift(K);

           obj.Tvec=linspace(0,obj.T,100); %time where 1.5*10^-16 s = 1
           obj.Sigma = zeros(1,length(obj.Tvec));
           figure
           for tt=1:length(obj.Tvec)
               Wave=0;
               t=obj.Tvec(tt);
                for kk=1:obj.N
                   KS=K(kk);           %Spatial frequency
                   KWV=2*pi*KS;        %Wave number
                   W=1/(2*obj.M)*KWV^2;    %Omega
                   PlaneWave=Y(kk)/obj.N*exp(1i* (KWV/obj.N*(1:obj.N)-W*t) );
                   Wave=Wave+PlaneWave;
                end

                plot(obj.X,abs(Wave))
                axis([-100 100 0 max(obj.A)])
                xlabel('X(A)');ylabel('|\Psi(x)|')

                P=abs(Wave); 
                P=P/max(P); %Spatial probability density
                obj.SigmaX(tt)=length(P(P>max(P)/2))/2*DX; %Find sigma of Gaussian
                s1 = sprintf('time: %d s',obj.Tvec(tt)*1.54e-11);
                s2=sprintf('Sigma: %d A',obj.SigmaX(tt));
                text(20,max(obj.A)/1.8,s1)
                text(20,max(obj.A)/2,s2);
                pause(.1)
            end


        end
    end
    
end

