
% Steady state temperature of sphere of radius (top half@100C and bottom half@0C)

format compact;
Temp = 100;			% temperature of top sphere

NptsR=51;
NptsA=71;
Nskip=2;			% use every Nskip coefs
Nplots=16*Nskip;
Nrows=4;			% number of columns
Ncols=(Nplots/Nrows)/Nskip;
r=linspace(0,1,NptsR);
a=linspace(0,pi,NptsA);
[R,A]=meshgrid(r,a);


b=zeros(1,Nplots);

% Legendre coef for U0 temperature
for jj=1:Nskip:Nplots
  L=legendre(jj-1,0);
  b(jj)=L(1)*((2*jj+1)/(2*(jj+1)));
end

b(1:Nskip:Nplots)

U=zeros(NptsA,NptsR);
U=U+Temp/2;

fs=8;
figure(3)
clf
plotsols=1;
for jj=1:Nskip:Nplots		% build Fourier series
  L=legendre(jj,cos(A));
  L0 = L(1,:,:);
  L1 = real(reshape(L0,NptsA,NptsR));
  Un=Temp*b(jj)*R.^(jj).*L1;
  U=U+Un;
  if(plotsols==1)
   subplot(Nrows,Ncols,(jj+1)/Nskip)
   set(gca,'FontSize',[fs]);
   surf(R,A,Un);
   shading interp
   zl = ['u(b_{',num2str(jj),'})'];
   xlabel('r','Fontsize',fs); ylabel('\theta','Fontsize',fs); 
   view([30 60])
   set(gca,'FontSize',[fs]);
   if(jj==1)
     axis tight
     ax=axis;
     ax(5)=-ax(6); 
   end
   axis(ax);
   axis tight;
   ax=axis;
   text((ax(2)-ax(1))/2,ax(4),ax(6),zl,'Fontsize',fs);
   pause(0.0001)
  end
end

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

  subplot(2,1,1)
  set(gca,'FontSize',[fs]);
  surf(R,A,U);
  shading interp
  xlabel('r','Fontsize',fs); ylabel('\theta','Fontsize',fs); zlabel('u','Fontsize',fs);axis tight
  view([60 32])
  set(gca,'FontSize',[fs]);

  subplot(2,1,2)
  set(gca,'FontSize',[fs]);
  surf(R,A,U);
  shading interp
  colormap(jet)
  view([0 90])
  xlabel('r','Fontsize',fs); ylabel('\theta','Fontsize',fs); zlabel('u','Fontsize',fs);axis tight
  colorbar
  set(gca,'FontSize',[fs]);
  pause(0.0001)


fs=16;
figure(5)
clf
hold on
ncolors=128;
cc=jet(ncolors);
umin=min(min(U));
umax=max(max(U));

set(gca,'FontSize',[fs]);
%X=R.*cos(pi/2-A);
%Y=R.*sin(pi/2-A);
for jj=1:NptsR
  for ii=1:NptsA
  x=Temp*R(ii,jj).*cos(pi/2-A(ii,jj));
  y=Temp*R(ii,jj).*sin(pi/2-A(ii,jj));
  u=U(ii,jj);
  index=ceil((u-umin)*(ncolors-1)/(umax-umin))+1;
  color=cc(index,:);
  plot3(x,y,u,'o','MarkerFaceColor',color,'Color',color,'MarkerSize',6)
  plot3(-x,y,u,'o','MarkerFaceColor',color,'Color',color,'MarkerSize',6)
 end 
end
xlabel('x','Fontsize',fs); ylabel('y','Fontsize',fs); zlabel('u','Fontsize',fs);axis tight
view([0 90])
set(gca,'FontSize',[fs]);
axis square;

