
% Generating solutions of the Laplace equation
% in a finite rectangular plate

format compact;
W = 10;				%width of plate
H = 3;				%Height of plate
T = 100;			%temperature of plate's boundary for y=0

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

fs=8;
figure(3)
clf

b=zeros(1,Nplots);

BC=3;
switch BC

 case 1 % Fourier coefs for T(x,y=H)=T*x*(W-x) temperature
  for jj=1:Nskip:Nplots
   if(jj/2==floor(jj/2))
    bn = 0;
   else    
    bn = (4*T*W^1/(pi^3))/(jj^3*sinh(jj*pi*H/W));
   end
   b(jj) = bn;
  end

 case 2 % Fourier coefs for T(x,y=0)=100 temperature
  for jj=1:Nskip:Nplots
   if(jj/2==floor(jj/2))
    bn = 0;
   else    
    bn = (2*T)*2/(jj*pi*sinh(jj*pi*H/W));
   end
   b(jj) = bn;
  end

 case 3 % Fourier coefs for T(x,y=0)=tent temperature (see exercise 7.9.23)
  h=20;			%height tent
  for jj=1:Nskip:Nplots
   if(jj/2==floor(jj/2))
    bn = 0;
   else    
    bn = (8*h/pi^2)*(-1)^((jj-1)/2)/(jj^2*sinh(jj*pi*H/W));
   end
   b(jj) = bn;
  end

 otherwise
  disp('Unknown BC');

end

b(1:Nskip:Nplots)

U=zeros(Npts,Npts);

plotsols=1;
for jj=1:Nskip:Nplots		% build Fourier series
  Un=b(jj)*sinh(jj*pi*(H-Y)/W).*sin(jj*pi*X/W);
  U=U+Un;
  if(plotsols==1)
   subplot(2*Nrows,Ncols,jj+Ncols*floor((jj-1)/Ncols))
   set(gca,'FontSize',[fs]);
   surf(X,Y,Un);
   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,Un);
   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,1,1)
  set(gca,'FontSize',[fs]);
  surf(X,Y,U);
  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,1,2)
  set(gca,'FontSize',[fs]);
  surf(X,Y,U);
  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]);

