function square_wave
%built by keith regner 9/2016 with R2014b (8.4 OSX)
clear all
close all

%% description
% This demo shows how a periodic square wave can be built by summing sine 
% and cosine waves with specific frequencies and amplitudes (Fourier 
% series). Each term in the Fourier series is shown in the upper plot. 
% The sum of the terms (shown in the top plot) and the actual square wave 
% are shown in the bottom plot. Change the number of terms in the Fourier 
% series using the slider bar. As the number of terms increases, the 
% summation approaches the actual square wave. 

%% create wave
b.fund = 1;                         %fundamental frequency
b.terms = 1;                        %number of terms
b.time = linspace(0,2,5000)';
slide_min = 1;                      %slider min
slide_max = 101;                    %slider max

b.A = 1;                            %actual squaure wave
b.t0 = 0;
b.T = 1/2;
b.sq_act = b.A*2*1i*(atanh(exp(-1i*pi*(b.time-b.t0)/b.T)) -...
    atanh(exp(1i*pi*(b.time-b.t0)/b.T)))/pi;

%% find square wave components
[b.F,b.Fsum] = components(b,b.terms);

%% figure and axis formatting
%2x1 subplot with slider bar
b.fig = figure('Units','inches','Position',[5 .5 8 8.3],...
    'Name','square_wave','NumberTitle','off');
b.ax1 = axes('Units','inches','Position',[1 5 6 3],'FontSize',12,...
    'LineWidth',1.5,'XLim', [0 2],'YLim',[-1.75 1.75],'Box','on');
hold on
b.ax1.XLabel.String = 'time [sec]';
b.ax1.YLabel.String = 'amplitude [a.u.]';

b.ax2 = axes('Units','inches','Position',[1 1.3 6 3],'FontSize',12,...
    'LineWidth',1.5,'XLim', [0 2],'YLim',[-1.75 1.75],'Box','on');
hold on
b.ax2.XLabel.String = 'time [sec]';
b.ax2.YLabel.String = 'amplitude [a.u.]';

%% plot solution
b.h_comp = plot(b.ax1,b.time,b.F,'LineWidth',1.5);
b.h_act = plot(b.ax2,b.time,real(b.sq_act),'Color',[1 0 0],...
    'LineWidth',1.5);
b.h_sum = plot(b.ax2,b.time,b.Fsum,'Color',[0 0 0],...
    'LineWidth',1.5);
b.leg2 = legend(b.ax2,{'actual','sum of sines'},...
    'Location','NorthEast','Box','off','FontSize',12);

%% add user control
%slider
b.txtbox = annotation('textbox','Units','inches','Position',...
    [3.3 .1 4 .3],'String',...
    ['number of terms = ' num2str(b.terms)],'FontSize',12,...
    'LineStyle','none');
b.sl = uicontrol('Parent',b.fig,'Style','slide','Units','inches',...
    'Position',[1.5 .5 5 .2],'Min',slide_min,'Max',slide_max,...
    'Value',b.terms,'SliderStep',[1/100 1/20],...
    'Callback',{@slider_call,b},'BackgroundColor',b.fig.Color);
b.label1 = uicontrol('Parent',b.fig,'Units','inches','Style',...
    'text','Position',[.8 .4 .6 .3],'String',num2str(slide_min),...
    'FontSize',12,'BackgroundColor',b.fig.Color);
b.label2 = uicontrol('Parent',b.fig,'Units','inches','Style',...
    'text','Position',[6.5 .4 .8 .3],'String',num2str(slide_max),...
    'FontSize',12,'BackgroundColor',b.fig.Color);

set(findall(b.fig,'-property','FontUnits'),'FontUnits','normalized')
set(findall(b.fig,'-property','Units'),'Units','normalized')

end

%% callback functions
%slider
function [] = slider_call(varargin)
%callback for the slider
[H,b] = varargin{[1,3]};
slider_value = H.Value;

%calculate new solution
[F,Fsum] = components(b,slider_value);

%update plot and indicators
cla(b.ax1)
plot(b.ax1,b.time,F,'LineWidth',1.5);
set(b.h_sum,'YData',Fsum);
set(b.txtbox,'String',...
    ['number of terms = ' num2str(slider_value)])
end

%% complex solution function
function [F,Fsum] = components(b,num)
terms_mat = repmat(1:num,length(b.time),1);
time_mat = repmat(b.time,1,num);

F = (4/pi)*sin(2*pi*((2*terms_mat) - 1)*b.fund.*time_mat)./...
    ((2*terms_mat) - 1);
Fsum = sum(F,2);

end

