function ezsurfspher(f,alpha,beta,h1,h2)
% ezsurfspher(f,alpha,beta,h1,h2)
%
% for sperical coordinates rho,phi,theta
% plot surface given by rho = f(phi,theta) over region given by
%   alpha <= theta <= beta
%   h1 <= phi <= h2
% where f is symbolic expression of phi,theta
% h1,h2 are symbolic expressions of theta
% 
% Example: part of sphere with radius 1, center (1,0,0) with z <= x
%   syms phi theta
%   f = 2*sin(phi)*cos(theta)
%   h1 = atan(1/cos(theta))
%   h2 = pi
%   ezsurfspher(f,-pi/2,pi/2,h1,h2)
%   nice3d

f = sym(f); h1 = sym(h1); h2 = sym(h2);
alpha = double(alpha); beta = double(beta);
syms phi theta

syms PHT THT
PHR = (1-PHT)*h1+PHT*h2;
FT = subs(f,phi,PHR);
X = subs(FT*sin(PHR)*cos(theta),theta,THT);
Y = subs(FT*sin(PHR)*sin(theta),theta,THT);
Z = subs(FT*cos(PHR),theta,THT);
ezsurf(X,Y,Z+0*PHT+0*THT,[0,1,alpha,beta])
p = findobj(gca,'type','surface'); p = p(1);
set(p,'MeshStyle','row','FaceColor','interp');
title('')


