classdef HW3_Solution_NIntegrate < handle
    %NIntegrate Numerical integration
    %   Perfoms integration using mid-point method
    %
    % Properties:
    %         f       %function handle
    %         h_req   %requested step size
    %         h       %actual step size
    %         a       %start of integral
    %         b       %end of integral
    %         EndPoints %End points of sub-steps
    %         Int_Result %Result
    
    properties
        f       %function handle
        h_req   %requested step size
        h       %actual step size
        a       %start of integral
        b       %end of integral
        EndPoints %End points of sub-steps
        Int_Result %Result
        Int_Method %determine the Method
    end
    
    methods
        function find_h(obj)
            %This code is simply taken from our script
            %Note that a->obj.a, etc for working with class properties
            
            Nsteps=ceil((obj.b-obj.a)/obj.h_req)+1;
            obj.EndPoints=linspace(obj.a,obj.b,Nsteps);
            %calculate spacing
            obj.h=obj.EndPoints(2)-obj.EndPoints(1);
        end
        
        function IntFunction(obj)
            %This code is simply taken from our script
            %Note that Int_Result->obj.Int_Result, etc for working with class properties
            obj.Int_Result=0;
            switch obj.Int_Method
                %The first case is used MidPoint method
                case 'MidPoint'
                    for ii=1:length(obj.EndPoints)-1
                        x_mid=obj.EndPoints(ii)+obj.h/2;
                        obj.Int_Result=obj.Int_Result+obj.f(x_mid)*obj.h;
                    end
                    fprintf('This is the solution of MidPoint Method %d ',obj.Int_Result)
                    
                    %The second case is used Trapezoid method
                    % Area of Trapezoid with two bases of a,b and height of
                    % h is A=(a+b)*h/2
                case 'Trapezoid'
                    for ii=1:length(obj.EndPoints)-1    %note 'length' function.
                        % find two bases for each step and multiple by
                        % height, here height is obj.h/2
                        Int=(obj.h/2)*(obj.f(obj.EndPoints(ii))+obj.f(obj.EndPoints(ii+1)));
                        obj.Int_Result=obj.Int_Result+Int;   %calculate rectangle area and add
                    end
                    fprintf('This is the solution of Trapezoid Method %d ',obj.Int_Result)
                    
                    %The third case is used Simpson method
                    % Check the wiki to find Simpson rule
                    % Integration for each step is obtained by
                    % (h/3)*[f(a)+f(a+h)+4*f(h)]
                case 'Simpson'
                    for ii=1:length(obj.EndPoints)-1
                        x_mid=obj.EndPoints(ii)+obj.h/2;    %find mid-point of this section
                        Int_step=(obj.h/6)*(obj.f(obj.EndPoints(ii))+4*obj.f(x_mid) ...
                            +obj.f(obj.EndPoints(ii+1)));
                        obj.Int_Result=obj.Int_Result+Int_step;
                    end
                    fprintf('This is the solution of Simpson Method %d ',obj.Int_Result)
                    
                otherwise
                    fprintf('Please choose integration methods: MidPoint? Trapezoid? Simpson?\r')
            end
            Int_Analytical=-exp(-obj.b)+exp(-obj.a)
            Int_Diff=Int_Analytical-obj.Int_Result;
            fprintf('This is the difference of the result from the exact known value %d \n',Int_Diff)
            
        end
        
        
        % Plot the function in exact value and Trapezoid approximation
        function PlotTrapezoid(obj)
            figure
            % First I plot f(x)= exp(-x) in 100 times smaller than h.
            X= obj.a:obj.h/100:obj.b;
            Y = obj.f(X);
            plot(X,Y,'b','linewidth',4);
            hold on
            % To create each trapezoid I need 5 points and put all of them
            % in X and Y
            for kk = 1:length(obj.EndPoints)-1;
                x_left = obj.EndPoints(kk);
                x_right = obj.EndPoints(kk+1);
                X = [x_left, x_right,x_right,x_left,x_left];
                Y = [0,0,obj.f(x_right),obj.f(x_left),0];
                plot(X,Y,'r','linewidth',1);
            end
            
            % You can also use stem command, but I recommended to use the
            % above method that is more straightforward.
            % XT= obj.a:obj.h:obj.b;
            % YT = obj.f(XT);
            % plot(XT,YT)
            % stem(XT,YT)
            xlabel('x')
            ylabel('exp(-x)')
            title('Trapezoid Method','fontsize',12)
            legend('Exp(-x)', 'Trapezoid Method')
            hold off
        end
    end
    
end

