function fourier(N,L)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% find and plot first N terms of 
%%% Fourier series for a function fun() assumed
%%% to be periodic on (-L,L)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% M. Peszynska, 11/1/2008
%% � Copyright Department of Mathematics, Oregon State University
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


acof = zeros(N,1);
bcof = zeros(N,1);


myfc = @(x)(fun(x));
ccof0 = 1/L*quad(myfc,-L,L);
for m=1:N
    myfs = @(x)(fun(x).*sin(m*x*pi/L));
    bcof(m) = 1/L*quad(myfs,-L,L);

    myfc = @(x)(fun(x).*cos(m*x*pi/L));
    acof(m) = 1/L*quad(myfc,-L,L);
end

plot(0:N,[ccof0; acof],1:N,bcof),
title('Fourier coefficients of the function','FontSize',14);
axis( [0 N min([ccof0;acof;bcof])*1.2 max([ccof0;acof;bcof])*1.2]);
legend('a_n','b_n');
pause

set(0,'DefaultLineLineWidth',2);  
xplot=linspace(-L,L,1000);
myv = myfourier(xplot,L,N,ccof0,acof,bcof);

for m=1:N
    plot(xplot,fun(xplot),xplot,myfourier(xplot,L,m,ccof0,acof,bcof));
    legend('Function','Fourier series',2);
    axis([-L L min(myv)*1.2 max(myv)*1.2 ]);
    title(sprintf('Function and its Fourier series with m=%d terms',m),'FontSize',14);
    pause
end
end


function v = myfourier(x,L,nterms,ccof0,acof,bcof)
    v = ccof0/2 + 0*x;
    if nterms > length(acof), return; end;
    for j=1:nterms
        v = v+acof(j)*cos(j*pi/L*x) + bcof(j)*sin(j*pi/L*x);
    end
end

function v = fun(x)
%%%% function periodic
v=fun2(x);
end

function v = fun1(x)
%%%% function periodic
v=0*x;
for j=1:length(x)
    if x(j) <=0
        v(j)=-1;
    else
        v(j)=1;
    end
        
end
end

function v = fun2(x)
%%%% function periodic
v=0*x;
for j=1:length(x)
    if x(j) <=0
        v(j)=pi+x(j);
    else
        v(j)=pi-x(j);
    end
        
end
end
