% -------------------------------------------------------------------------------- 
% This script plots a few solutions to HW2 
% 	in MAE105: Introduction to Mathematical Physics (spring 2015).
%
% April 26th, 2015
% Greg Wagner
%
% the first part plots the solution to problem 1, which is a solution of the 
%	heat equation for an initial condition which is a straight line
%	between x=0 and x=L.
%
% the second part plots the solution to Laplace's equation in a square, 
%	where [x,y] = [ [0,L] , [0,L] ], and with the boundary conditions
%	
%				u(x=0)=-1,	u(x=L)=1,	u(y=0)=0,	u(y=L)=0.
%
% -------------------------------------------------------------------------------- 
%% some preliminaries

% clear all variables
clear all

% plot parameters
lw = 2;
fs = 16;

% marker strings
markerstr{1} = 'k-';
markerstr{2} = 'b-';
markerstr{3} = 'r-';
markerstr{4} = 'g-';
markerstr{5} = 'c-';
markerstr{6} = 'm-';
markerstr{7} = 'y-';

% -------------------------------------------------------------------------------- 
%% solution to the heat equation

% physical parameters
a = 1;
k = 1;
L = 1;

% grid
x0 = 0;
xf = L;
nx = 1e3;
xv = linspace(x0,xf,nx)';

% time grid
t0 = 0;
tf = 1*L^2/k; 						% some characteristic time
nt = 5;
tv = linspace(t0,tf,nt);

% modes
nm = 1e2;

% initialize figure
figure(2), clf, set(gca,'FontSize',fs), grid, hold on, ylim([-1 0])

% plot solution at various times
for tt = 1:nt

	% initialize u(x,t(tt))
	u = zeros(nx,1);

	% add up the spatial modes
	for ii = 1:nm
		% exact form of coefficients for the linear initial condition
		Am = -8*a*L/( (2*ii-1)^2*pi^2 );
		u = u + Am*cos((ii-1/2)*pi*xv/L)*exp(-k*(ii-1/2)^2*pi^2/L^2*tv(tt));
	end

	% plot the solution
	figure(2)
	plot(xv,u,markerstr{tt},'LineWidth',lw)
	xlabel('x'), ylabel('u(x,t)')
	legendstr{tt} = ['t = ',num2str(tv(tt),'%0.2f')];

	%
	% note: can also make a movie using this loop...
	%

end

% add a legend
figure(2), legend(legendstr,'location','southeast')

% -------------------------------------------------------------------------------- 
%% laplace in a square
L = 1;

% grid
nx = 1e2; ny=nx;

% number of modes
np = 1e2;

% initialize u(x,y)
u = zeros(nx,ny);

% initialize grid
x = linspace(0,L,nx);
y = linspace(0,L,ny);
[X,Y]=meshgrid(x,y);

% construct u(x,y) from exact solution
for p = 1:np
	u = u + 4*sin((2*p-1)*pi*Y)/((2*p-1)*pi*sinh((2*p-1)*pi)) ...
			.* ( sinh((2*p-1)*pi*X)+sinh((2*p-1)*pi*(X-1)));
end

% plot u(x,y)
figure(3), clf, set(gca,'FontSize',16)
pcolor(X,Y,u), shading interp, colormap(parula), colorbar('eastoutside')
xlabel('x'), ylabel('y'), title('u(x,y)')
