%                       Program fseries
%
%       This program sums the real Fourier series on [-2,2] for several
%   different functions.  User enters the choice of function f (and 
%   sometimes the choice of sine or cosine series), and the number of
%   terms 'N' to be summed up.  The program lists the coefficients in 
%   expansion, and plots the Nth partial sum against the function f.
%   It also computes the mean square error sigma = (1/2L) ||f - s_N||^2.
%   In our case, L = 2.
%       The choices of functions are 
%
%   1) f(x) given on [0,2] is f(x) = x for 0 < x < 1, f(x) = 2-x for 1 < x < 2.
%      User then has choice of making an even or odd extension of f to [-2,2].
%   2) f(x) given on [0,2] is f(x) = 1-x/2.  User has choice of making even
%      or odd extension of f to [-2,2].
%   3) f(x) given on [-2,2] is f(x) = 1 for |x| < 1, f = 0 for 1 < |x| <2.
%   4) f(x) given on [-2,2] is f(x) = x(x+2)(x-2).


x =  -2:.005:2;
m = input('Enter choice of data 1,2,3,or 4.  ')



%  first choice of data
if m == 1
   disp(' 1 for the odd extension to [-2,2] (sine series)  ')
   disp(' 2 for the even extension to [-2,2] (cosine series)'  )
   mm = input('Enter 1 or 2     ')
N = input('enter the number of terms to be summed   ')
n = 1:N;

%  odd extension of first choice of data
   if mm  == 1
       sumsig = 0
       b = zeros(size(n));
       for k = 1:N
          b(k) = (8.0/(k*pi)^2)*.5*(1-(-1)^k)*(-1)^((k-1)/2);
          sumsig = sumsig +b(k)^2;
       end


        fprintf( 'These are the coefficients b(k) \n\n')
        fprintf( '         k            b(k)       \n\n')
        fprintf( '       %2.0f       %.8f  \n', [n;b]  );

       sigma2 = 1/3 - .5*sumsig;
   fprintf( ' \n')
   fprintf(' For N = %2.0f, the mean square error sigma^2 = %1.8f\n', N, sigma2)
                
        partialsum = 0;
        for k = 1:N
            partialsum = partialsum + b(k)*sin(k*x*pi/2);
        end

        f = (x <-1).*(-2-x) +( (x >= -1) -(x >1)).*x +(x> 1).*(2-x);
        z = zeros(size(x));
        plot(x,f,x,partialsum,'g',x,z,'r')

%    even extension of first choice of data 
     elseif mm == 2
       sumsig = .5;
       a = zeros(size(n));
       for k = 1:N
          a(k) = (8.0/(k*pi)^2)*((-1)^(k/2)-1)*.5*(1+(-1)^k);
          sumsig = sumsig + a(k)^2;
       end
    
        aa = [1,a];

        n = [0,n];

        fprintf('These are the coefficients a(k) \n\n')
        fprintf( '         k            a(k)       \n\n')
        fprintf( '       %2.0f       %.8f  \n', [n;aa]  )
        fprintf( ' \n')
       

        sigma2 = 1/3 -.5*sumsig;

   fprintf(' For N = %2.0f, the mean square error sigma^2 = %1.8f\n', N,sigma2)

       partialsum = .5*ones(size(x));
       for k = 1:N
          partialsum = partialsum +a(k)*cos(k*x*pi/2);
       end

       f1 = (x< -1).*(x+2) +((x< 0) - (x < -1)).*(-x);
       f2 = ((x < 1) -(x < 0)).*x + (x >= 1).*(2-x);
       f = f1 + f2;
       plot(x,f,x,partialsum,'g')

   end

% begin second choice of data

elseif m == 2 

   disp(' 1 for the odd extension on [-2,2] (sine series) ' )
   disp(' 2 for even extension [-2,2] (cosine series) ')
   mm = input(' Enter 1 or 2  ')
   N = input('enter the number of terms to be summed  ')
   n = 1:N;

%  odd extension of second choice of data 
   if mm == 1 
      sumsig = 0;
      b = zeros(size(n));
      for k = 1:N
         b(k) = 2.0/(k*pi);
         sumsig = sumsig +b(k)^2;
      end

      fprintf('These are the coefficients b(k)\n \n')
      fprintf('       k           b(k)          \n\n')
      fprintf('     %2.0f        %.8f    \n', [n;b]  )
      fprintf('\n')

      sigma2 = 1/3 - .5*sumsig;

      fprintf('For N = %2.0f, the mean square error sigma = %1.8f\n', N,sigma2)


      partialsum = 0;
      for k = 1:N
         partialsum = partialsum + b(k)*sin(k*x*pi/2);
      end
     
      f = (x < 0).*(-1 -x/2) + (x >= 0).*(1-x/2);
      plot(x,f,x, partialsum, 'g')

%  even extension of second choice of data
   elseif mm == 2
       sumsig = .5;
       a = zeros(size(n));
       for k = 1:N
          a(k) = 2*(1-(-1)^k)/(k*pi)^2;
          sumsig = sumsig + a(k)^2;
       end 

       aa = [1.0, a];
       n = [0,n];

       fprintf('These are the coefficients a(k)  \n \n')
       fprintf('       k              a(k)        \n \n')
       fprintf('     %2.0f          %.8f  \n', [n;aa] )
       fprintf('\n')


       sigma2 = 1/3 - .5*sumsig;

       fprintf('For N = %2.0f the mean square error sigma^2 = %.8f\n', N, sigma2)

       partialsum = .5*ones(size(x));
       for k = 1:N
           partialsum = partialsum +a(k)*cos(k*x*pi/2);
       end
       
       f = (x< 0).* (x/2 +1) + (x >= 0).* (1- x/2);
       plot(x,f,x,partialsum,'g')
% end the logical choice of mm = 1 or mm = 2 for second choice of data
   end

% third choice of data (an even function )
elseif m == 3
N = input('enter the number to terms to be summed   ')
n = 1:N;
    sumsig = .5;
    a = zeros(size(n));
    for k = 1:N
        a(k) = (1/(k*pi))* (-1)^((k-1)/2)*(1-(-1)^k);
        sumsig = sumsig + a(k)^2;
    end


    n = [0,n];
    aa = [1,a];

    fprintf('These are the coefficients a(k)   \n\n')
    fprintf('          k             a(k)      \n\n')
    fprintf('        %2.0f          %.8f       \n', [n;aa])
    fprintf('\n')


    sigma2 = .5 - .5*sumsig;

    fprintf('For N = %2.0f, the mean square error sigma^2 = %.8f \n', N, sigma2)

    partialsum = .5*ones(size(x));
    for k = 1:N
        partialsum = partialsum + a(k)*cos(k*x*pi/2);
    end

    f = (x >= -1) - (x > 1);
    plot(x,f,x,partialsum,'g')

% begin fourth choice of data (an odd function)

else 
N = input('enter the number of terms to be summed  ')
n = 1:N;

    b = zeros(size(n));
    sumsig = 0;
    for k = 1:N
        b(k) = -12*(2/(k*pi))^3 *(-1)^k;
        sumsig = sumsig + b(k)^2;
    end

    fprintf('These are the coefficients b(k)  \n\n')
    fprintf('         k             b(k)          \n\n')
    fprintf('       %2.0f          %.8f          \n', [n;b])
    fprintf('\n')

    sigma2 = 512/105 - .5*sumsig;

    fprintf('For N = %2.0f, the mean square error sigma^2 = %.8f\n', N, sigma2)
    partialsum = 0;
    for k = 1:N
       partialsum = partialsum + b(k)*sin(k*x*pi/2);
    end

    f = x.*(x+2).*(2-x);

    plot(x,f,x,partialsum,'g')

end