tip =   10.0;       (* angle in degrees the axis makes with the normal *)

density = 0.81;

tippy = tip*Pi/180;

b = Tan[tippy];            (* y = b z *)	

r = 1;
  
  
p = r/Sqrt[b^2+1];

Print[" top is at p ",p];

hz = -3/8 r/Sqrt[1+1/b^2] ;     (* positive h means h below *)

hy = -1/b hz;

Print["check" Sqrt[hz^2+hy^2]," should be ",3.0/8.0 r]

Print["p is now ",N[p]," cap starts at ",N[-r/Sqrt[b^2+1]]];

ans1 = Integrate[1,{x,-Sqrt[r^2 - z^2 - y^2],+Sqrt[r^2 - z^2 -y^2]}];

(*
ans2 = Integrate[ans1,{y,b z,Sqrt[r^2- z^2]}];
*)

ans2 = Pi(r^2-z^2) -(r^2-z^2)ArcCos[- b z/Sqrt[r^2-z^2]] + b z Sqrt[r^2-z^2(1+b^2)]


cap = Integrate[Pi Sqrt[r^2 - z^2]^2,{z,-r,-p}];

Print["cap," ,cap];

volume[w_] := NIntegrate[ans2,{z,-p,w}]+cap - 2/3 density Pi r^3;

delta = 0.000001;

w = 0;

Do[ w = w - volume[w] delta/(volume[w+delta]-volume[w]),{i,1,30}];


w = Re[w];

Mxz = 2/3 Integrate[ Sqrt[r^2 - (b^2+1) z^2]^3,{z,-p,w}]; 

ybar =  Mxz/(2/3*density*Pi r^3);

ybar = Re[ybar];

N1 = Integrate[z, {x,-Sqrt[r^2 - z^2 - y^2],+Sqrt[r^2 - z^2 -y^2]}];

(*
N2 = Integrate[N1, {y,b z,Sqrt[r^2- z^2]}];

Mxy =  -(p + 3/8(r-p))*cap + NIntegrate[ N2, {z,-p,w}];
*)

Mxy = -(p+3/8(r-p))*cap + NIntegrate[N1,{z,-p,w},{y,b z,Sqrt[r^2 - z^2]}]

zbar = Mxy/(2/3 density Pi r^3);

zbar = Re[zbar];

Print["(",ybar,",",zbar,")"];
Print[" full centroid is at  (",hy,",",hz,")"];


a1 = ParametricPlot[{r Cos[t],r Sin[t]},{t,(Pi/2-tippy)-Pi,(Pi/2-tippy)}];

a2 = ParametricPlot[{b z,z},{z,-p,p}];

a3 = Plot[w,{y, b w,r},PlotStyle->{Thickness[0.008]}];

a4 = Graphics[Point[{hy,hz}]];

a5 = Graphics[Point[Re[{ybar,zbar}]]];

Show[a1,a2,a3,a4,a5,AspectRatio->Automatic];

Print[ybar," ",hy," bottom, full centroid "];