% shape22.m  surface of shape22 U
%
% |\____/|    |
% | ____ |    |   2D
% |/    \|    |
function shape22
  clear
  diary('shape22_m.out')
  format compact
  disp('shape22.m running')
  n = 17;
  m = 7;
  % [xg,yg] = meshgrid(1:1:n,1:1:m); % then fill with my data
  xg=zeros(n,m);
  yg=zeros(n,m);
  zg=zeros(n,m);
  r = 1.0;      % can scale
  yc1 = 1.25*r; % can scale
  yc2 = -yc1;
  xc = 0.0;
  ang1 = 5.0*pi/4.0; % to 7 pi/4
  ang2 = 3.0*pi/4.0; % to pi/4 - dang
  dang = pi/(2.0*(n-1));

  disp('generate points');
  ang = ang1;
  for i=1:n % top
    x = xc+cos(ang);
    y = yc1+sin(ang);
    xg(i,1) = x;
    yg(i,1) = y;
    zg(i,1) = U(x,y);
    ang = ang + dang;
  end %i
  ang = ang2;
  for i=1:n; % bottom
    x = xc+cos(ang);
    y = yc2+sin(ang);
    xg(i,m) = x;
    yg(i,m) = y;
    zg(i,m) = U(x,y);
    ang = ang - dang;
  end % i
  disp(' ');
  disp('generate grid');
  for i=1:n; % grid
    yd = (yg(i,1)-yg(i,m))/(m-1);
    disp('yd');
    disp(yd);
    for j=2:6;  %grid 1+1, m-1
      x = xg(i,1);
      y = yg(i,1)-(j-1)*yd;
      xg(i,j) = x;
      yg(i,j) = y;
      zg(i,j) = U(x,y);
    end % j
  end % i
  disp('xg');
  disp(xg);
  disp('yg');
  disp(yg);
  disp('zg');
  disp(zg);
  disp('shape22.m end printing')
  diary off

  % plot surface
  surf(xg,yg,zg);
  title('shape22 surf')
  xlabel('X');
  ylabel('Y');
  zlabel('Z');

  % plot
  plot3(xg,yg,zg);
  title('shape22 plot2')
  xlabel('X');
  ylabel('Y');
  zlabel('Z');
  return
  
  function z = U(x,y)
    z=x*x*x+2.0*y*y*y;
  end % U

end % shape22.m
