> restart:with(plots):
h:=10:
r:=4:
alpha:=arctan(r/h):
beta:=0.4:
n:=20:
hs1:=7.4: # centre première sphère
rs1:=(h-hs1)*sin(alpha):# son rayon
hc1:=hs1+rs1*sin(alpha): # cote cercle de contact et directrice 1
xd1:=solve(hs1-rs1/cos(beta)+tan(beta)*x=hc1,x):
xaa:=solve(hs1-rs1/cos(beta)+tan(beta)*x=cot(alpha)*x+h,x):zaa:=cot(alpha)*xaa+h:
xa:=solve(hs1-rs1/cos(beta)+tan(beta)*x=-cot(alpha)*x+h,x):za:=-cot(alpha)*xa+h:
aa:=sqrt((xaa-xa)**2+(zaa-za)**2):
hs2:=hs1-aa/cos(alpha):
rs2:=(h-hs2)*sin(alpha):
hc2:=hs2+rs2*sin(alpha):
xd2:=solve(hs1-rs1/cos(beta)+tan(beta)*x=hc2,x):
plan1:=plot3d(hc1,x=-r..r,y=-r..r,grid=[12,12],style=wireframe,color=blue):
plan2:=plot3d(hs1-rs1/cos(beta)+tan(beta)*x,x=-r..r,y=-r..r,style=wireframe,grid=[20,13],color=grey):
plan3:=plot3d(hc2,x=-r..r,y=-r..r,grid=[12,12],style=wireframe,color=blue):
demico:=plot3d(h-sqrt(x**2+y**2)/tan(alpha),x=-r..r,y=-r..0,style=patchnogrid,color=yellow):
coco:=seq(spacecurve([[0,0,h],[r*cos(t),r*sin(t),0]],numpoints=2),t=seq(k*2*Pi/n,k=1..n)),spacecurve([r*cos(t),r*sin(t),0],t=0..2*Pi):
sphere1:=plot3d([rs1*cos(phi)*cos(theta),rs1*cos(phi)*sin(theta),hs1+rs1*sin(phi)],theta=0..2*Pi,phi=-Pi/2..Pi/2,style=patchnogrid):
sphere2:=plot3d([rs2*cos(phi)*cos(theta),rs2*cos(phi)*sin(theta),hs2+rs2*sin(phi)],theta=0..2*Pi,phi=-Pi/2..Pi/2,style=patchnogrid):
directrice1:=spacecurve([[xd1,-4,hc1],[xd1,4,hc1]],color=black,thickness=3):
foyer1:=pointplot3d([rs1*sin(beta),0,hs1-rs1*cos(beta)],style=point,symbol=circle,thickness=3,color=black):
directrice1:=spacecurve([[xd1,-4,hc1],[xd1,4,hc1]],color=black,thickness=3):
foyer2:=pointplot3d([-rs2*sin(beta),0,hs2+rs2*cos(beta)],style=point,symbol=circle,thickness=3,color=black):
directrice2:=spacecurve([[xd2,-4,hc2],[xd2,4,hc2]],color=black,thickness=3):
xp:=unapply(r*cos(t)*(1-zz/h),t):
yp:=unapply(r*sin(t)*(1-zz/h),t):
zp:=unapply(solve(zz=hs1-rs1/cos(beta)+tan(beta)*xp(t),zz),t):
a:=2:
ptcourant:=spacecurve([[rs1*sin(beta),0,hs1-rs1*cos(beta)],[subs(zz=zp(a),xp(a)),subs(zz=zp(a),yp(a)),zp(a)],[xd1,subs(zz=zp(a),yp(a)),hc1]],color=black),spacecurve([[-rs2*sin(beta),0,hs2+rs2*cos(beta)],[subs(zz=zp(a),xp(a)),subs(zz=zp(a),yp(a)),zp(a)],[xd2,subs(zz=zp(a),yp(a)),hc2]],color=black):
ellipse:=spacecurve([subs(zz=zp(t),xp(t)),subs(zz=zp(t),yp(t)),zp(t)],t=0..2*Pi,thickness=3,color=red):
rc1:=rs1*cos(alpha):rc2:=rs2*cos(alpha):
CercleContact1:=spacecurve([rc1*cos(theta),rc1*sin(theta),hs1+rs1*sin(alpha)],theta=0..2*Pi,thickness=2,color=black):
CercleContact2:=spacecurve([rc2*cos(theta),rc2*sin(theta),hs2+rs2*sin(alpha)],theta=0..2*Pi,thickness=2,color=black):

> display([coco,demico,plan1,plan2,plan3,sphere1,sphere2,ellipse,foyer1,foyer2,ptcourant,directrice1,directrice2,CercleContact1,CercleContact2],scaling=constrained,view=[-4..4,-4..4,1..10]);

>