
% Generating solutions of the Laplace equation
% in a finite rectangular plate by superimposing two
% solutions.

format compact;
W = 7;				%width of plate
H = 5;				%Height of plate
Tl = 150;			%temperature of plate's left boundary
Tb = 100;			%temperature of plate's bottom's boundary

Npts=101;
Nplots=80;
Nrows=5;			% number of columns
Ncols=Nplots/Nrows;
Nskip=2;			% use every Nskip coefs
x=linspace(0,W,Npts);
y=linspace(0,H,Npts);
[X,Y]=meshgrid(x,y);

fs=8;
plotsols=0;
 if(plotsols==1)
 figure(2)
 clf
 figure(3)
 clf
end

bb=zeros(1,Nplots);		% bottom
bl=zeros(1,Nplots);		% left

% Fourier coefs for Tb temperature
for jj=1:Nskip:Nplots
 if(jj/2==floor(jj/2))
  bn = 0;
 else    
  bn = (2*Tb)*2/(jj*pi*sinh(jj*pi*H/W));
 end
 bb(jj) = bn;
end

% Fourier coefs for Tl temperature
for jj=1:Nskip:Nplots
 if(jj/2==floor(jj/2))
  bn = 0;
 else    
  bn = (2*Tl)*2/(jj*pi*sinh(jj*pi*W/H));
 end
 bl(jj) = bn;
end

%% Fourier coefs for Tb=tent temperature (see exercise 7.9.23)
%hb=20;			%height tent
%for jj=1:Nskip:Nplots
% if(jj/2==floor(jj/2))
%  bn = 0;
% else    
%  bn = (8*hb/pi^2)*(-1)^((jj-1)/2)/(jj^2*sinh(jj*pi*H/W));
% end
% bb(jj) = bn;
%end
%
%hl=25;			%height tent
%for jj=1:Nskip:Nplots
% if(jj/2==floor(jj/2))
%  bn = 0;
% else    
%  bn = (8*hl/pi^2)*(-1)^((jj-1)/2)/(jj^2*sinh(jj*pi*W/H));
% end
% bl(jj) = bn;
%end

bb(1:Nskip:Nplots)
bl(1:Nskip:Nplots)

Ub=zeros(Npts,Npts);
Ul=zeros(Npts,Npts);

for jj=1:Nskip:Nplots		% build Fourier series
  Unb=bb(jj)*sinh(jj*pi*(H-Y)/W).*sin(jj*pi*X/W);
  Unl=bl(jj)*sinh(jj*pi*(W-X)/H).*sin(jj*pi*Y/H);
  Ub=Ub+Unb;
  Ul=Ul+Unl;
  if(plotsols==1)
   figure(2)
   subplot(2*Nrows,Ncols,jj+Ncols*floor((jj-1)/Ncols))
   set(gca,'FontSize',[fs]);
   surf(X,Y,Unb);
   shading interp
   %colormap(gray)
   zl = ['z(b_{',num2str(jj),'})'];
   xlabel('x','Fontsize',fs); ylabel('y','Fontsize',fs); zlabel(zl,'Fontsize',fs);axis tight
   view([12 46])
   if(jj==1)
     ax=axis;
     ax(5)=-ax(6)/3; 
   end
   axis(ax);
   set(gca,'FontSize',[fs]);
  
   subplot(2*Nrows,Ncols,jj+Ncols*floor((jj-1)/Ncols)+Ncols)
   set(gca,'FontSize',[fs]);
   surf(X,Y,Unb);
   shading interp
   colormap(jet)
   view([0 90])
   xlabel('x','Fontsize',fs); ylabel('y','Fontsize',fs); zlabel('z','Fontsize',fs);axis tight
   set(gca,'FontSize',[fs]);
   %fprintf('Press any key to continue...\n');

   figure(3)
   subplot(2*Nrows,Ncols,jj+Ncols*floor((jj-1)/Ncols))
   set(gca,'FontSize',[fs]);
   surf(X,Y,Unl);
   shading interp
   %colormap(gray)
   zl = ['z(b_{',num2str(jj),'})'];
   xlabel('x','Fontsize',fs); ylabel('y','Fontsize',fs); zlabel(zl,'Fontsize',fs);axis tight
   view([12 46])
   if(jj==1)
     ax=axis;
     ax(5)=-ax(6)/3; 
   end
   axis(ax);
   set(gca,'FontSize',[fs]);
  
   subplot(2*Nrows,Ncols,jj+Ncols*floor((jj-1)/Ncols)+Ncols)
   set(gca,'FontSize',[fs]);
   surf(X,Y,Unl);
   shading interp
   colormap(jet)
   view([0 90])
   xlabel('x','Fontsize',fs); ylabel('y','Fontsize',fs); zlabel('z','Fontsize',fs);axis tight
   set(gca,'FontSize',[fs]);
   %fprintf('Press any key to continue...\n');
   %pause
   pause(0.0001)
  end
end

% Now plot the Fourier combo:
fs=16;
figure(4)
clf

  subplot(2,2,1)
  set(gca,'FontSize',[fs]);
  surf(X,Y,Ub);
  shading interp
  %colormap(gray)
  xlabel('x','Fontsize',fs); ylabel('y','Fontsize',fs); zlabel('z','Fontsize',fs);axis tight
  view([60 32])
  set(gca,'FontSize',[fs]);

  subplot(2,2,3)
  set(gca,'FontSize',[fs]);
  surf(X,Y,Ub);
  shading interp
  colormap(jet)
  view([0 90])
  xlabel('x','Fontsize',fs); ylabel('y','Fontsize',fs); zlabel('z','Fontsize',fs);axis tight
  colorbar
  set(gca,'FontSize',[fs]);
  axis equal

  subplot(2,2,2)
  set(gca,'FontSize',[fs]);
  surf(X,Y,Ul);
  shading interp
  %colormap(gray)
  xlabel('x','Fontsize',fs); ylabel('y','Fontsize',fs); zlabel('z','Fontsize',fs);axis tight
  view([60 32])
  set(gca,'FontSize',[fs]);

  subplot(2,2,4)
  set(gca,'FontSize',[fs]);
  surf(X,Y,Ul);
  shading interp
  colormap(jet)
  view([0 90])
  xlabel('x','Fontsize',fs); ylabel('y','Fontsize',fs); zlabel('z','Fontsize',fs);axis tight
  colorbar
  set(gca,'FontSize',[fs]);
  axis equal

fs=16;
figure(5)
clf

  subplot(2,1,1)
  set(gca,'FontSize',[fs]);
  surf(X,Y,Ub+Ul);
  shading interp
  %colormap(gray)
  xlabel('x','Fontsize',fs); ylabel('y','Fontsize',fs); zlabel('z','Fontsize',fs);axis tight
  view([60 32])
  set(gca,'FontSize',[fs]);
  axis tight

  subplot(2,1,2)
  set(gca,'FontSize',[fs]);
  surf(X,Y,Ub+Ul);
  shading interp
  colormap(jet)
  view([0 90])
  xlabel('x','Fontsize',fs); ylabel('y','Fontsize',fs); zlabel('z','Fontsize',fs);axis tight
  colorbar
  set(gca,'FontSize',[fs]);
  axis equal
  axis tight
