%
%  w01plot.m
%
%  save file as "w01plot.m" (make sure that it is on the matlab path)
%  to run, type "w01plot" at the matlab prompt ">>"
%

%  initialize workspace & set some control variables
clear;  close all
Np = 200;  %  number of plot/integration points (this really matters!)

%  set x (row) vector (type "x" at the >>)
dx = pi/Np;  x = -pi:dx:pi;

%  constant C>0 
C = 0.5;

%  function y(x) at x-values (elementwise arithmetic, type "help arith")
y = sin(x) ./ (cosh(C)-cos(x));

%  choose figure & clear (type "help hold")
figure(1);  clf;  hold on

%  plot y(x) & set axis (type "help axis")
plot(x,y,'k','linewidth',2)
axis([-pi pi 1.1*min(y) 1.1*max(y)])

%  compute some fourier series (type "series(4)")
series = [1 2 3 4 5];  errors = zeros(size(series));
%series = 2.^(1:1:6);  errors = zeros(size(series));

%  this loop chooses different Ns-values as defined in 'series'
for js = 1:length(series)
    yf = 0;
    Ns = series(js);
    
    %  this loop sums the fourier series to the Ns-th term
    for j = 0:1:Ns
        Bj = 2*exp(-j*C);
        yf = yf + Bj*sin(j*x);
    end
    plot(x,yf,'b')
    
    %  trapezoidal rule integrated squared-error
    %  (don't double count endpoints)
    errors(js) = (sum((y-yf).^2) - (y(1)-yf(1))^2)*dx;
end

%  !!! put some info on the plot !!!
title('\bf fourier approximations')
xlabel('\bf x-axis')
ylabel('\bf y(x) in black;  y_N(x) in blue')

%  error plot (semi-log)
figure(2);  clf;  hold on
plot(series,log(errors),'kx')

title('\bf semi-log plot of error')
xlabel('\bf number of terms')
ylabel('\bf log integral of (y(x)-y_N(x))^2')

%  add line with slope -2C
plot([series(1) series(end)],...
    [series(end)-series(1) 0]*2*C+log(errors(end)),'--')