#disjoint.maple # # Frank Sottile # College Station # 22 August 2006 #################################################################### # # The copyright to this file, and the images it draws are owned by # Frank Sottile. We allow the free use of this, but you must first # ask permission of either Frank Sottile or Thorsten Theobald, and # acknowledge our creation of these images. # #################################################################### # # This file computes and draws an example of four disjoint spheres # in space having 12 common tangents. Three have centres at the # corners of an equilateral triangle of side length sqrt(3), # # [0,1,-35/100], [-sqrt(3)/2,-1/2,-35/100], [sqrt(3)/2,-1/2,-35/100], # # and have radius 4/5. # # The fourth is centred at [0,0,0] with radius 1/4. # # Since sqrt(1+z^2) > 4/5 + 1/4, they are disjoint. # #################################################################### interface(quiet=true): with(plots): plotsetup(x11,plotoptions=`width=3in,height=2.5in`): Check := proc(p,v,c,ra) # # This computes the distance from c of the line through p # with direction v, and compares it to ra # local f,i: f := sum( (p[i]+t*v[i]-c[i])^2,i=1..3 ): lprint(evalf(sqrt(minimize(f))-ra,8)); end proc: ##################################################### Centres:=[[0,1,z],[-sqrt(3)/2,-1/2,z],[sqrt(3)/2,-1/2,z]]: Radii:= [4/5,4/5,4/5]: r:=1/4: rr:=4/5: z:=-35/100: # # Checks that the spheres are disjoint # #evalf(sqrt(1+z^2)),evalf(r),evalf(rr),evalf(sqrt(1+z^2)-r-rr); # # Here are the equations for a line through the point p with direction # v to be tangent to the spheres. # # We first tries to set up the derivation of p using linear algebra # M:=linalg[matrix](Centres): v:=linalg[vector]([v1,v2,v3]): c:=[seq(linalg[vector](Centres[i]),i=1..3)]: vsq:=linalg[dotprod](v,v,orthogonal): A:=linalg[multiply](linalg[inverse](M), linalg[vector]([seq((linalg[dotprod](c[i],v,orthogonal))^2,i=1..3)])): B:= linalg[multiply](linalg[inverse](M), linalg[vector]([seq(linalg[dotprod](c[i],c[i],orthogonal)+r^2-Radii[i]^2,i=1..3)])): # # But it fails (@#@$%^&* MAPLE !), so we do a work-around # p:=linalg[vector]([seq(simplify(-A[i]/2/vsq + B[i]/2), i=1..3)]): Quartic:=simplify(linalg[dotprod](p,p,orthogonal)-r^2): Cubic:= simplify(linalg[dotprod](p,v,orthogonal)): # # Simplify is not really working properly. I have introduced a major hack # to get around this. Warning: it may fail sometimes # De:=op(1,denom(Quartic)): Quartic:=numer(Quartic)/De: De:=op(1,op(1,denom(Cubic))): Cubic:=numer(Cubic)/De: ######################################################################## Digits:=20: #G:=Groebner[Basis](subs(v1=1,[Quartic,Cubic]),plex(v2,v3)): G:=Groebner[gbasis](subs(v1=1,[Quartic,Cubic]),plex(v2,v3)): Vs3:=fsolve(G[1]=0,v3): Vs2:=[seq(solve(subs(v3=Vs3[i],G[2])=0,v2),i=1..12)]: P:=[seq(evalf(subs(v1=1, v2=Vs2[i], v3=Vs3[i],evalm(p))), i=1..12)]: V1:=[seq(1/evalf(sqrt(1+Vs2[i]^2+Vs3[i]^2)),i=1..12)]: V2:=[seq(evalf(V1[i]*Vs2[i]),i=1..12)]: V3:=[seq(evalf(V1[i]*Vs3[i]),i=1..12)]: ################################################################### # # Checks that the 12 solution lines are tangent to the 4 spheres # #for j from 1 to 12 do # for k from 1 to 3 do # Check(P[j], [V1[j],V2[j],V3[j]], Centres[k], rr); # end do: # Check(P[j], [V1[j],V2[j],V3[j]], [0,0,0], r); #end do: # ################################################################### # The Four Spheres S4:=plot3d([r*sin(ph)*cos(th),r*sin(ph)*sin(th),r*cos(ph)],ph=0..Pi, th=-Pi..Pi,grid=[10,30],color=magenta): S1:=plot3d([Centres[1][1]+Radii[1]*sin(ph)*cos(th), Centres[1][2]+Radii[1]*sin(ph)*sin(th), Centres[1][3]+Radii[1]*cos(ph)],ph=0..Pi, th=-Pi..Pi,grid=[14,30],color=cyan): S2:=plot3d([Centres[2][1]+Radii[2]*sin(ph)*cos(th), Centres[2][2]+Radii[2]*sin(ph)*sin(th), Centres[2][3]+Radii[2]*cos(ph)],ph=0..Pi, th=-Pi..Pi,grid=[14,30],color=pink): S3:=plot3d([Centres[3][1]+Radii[3]*sin(ph)*cos(th), Centres[3][2]+Radii[3]*sin(ph)*sin(th), Centres[3][3]+Radii[3]*cos(ph)],ph=0..Pi, th=-Pi..Pi,grid=[14,30],color=yellow): ################################################################### # aquamarine black blue brown coral cyan # gold green gray grey khaki magenta # maroon orange navy pink plum red # sienna tan turquoise violet wheat white # yellow # 1 2 3 4 5 6 7 8 9 10 11 12 CO:=[black,green, tan,aquamarine,brown,coral,gray,violet, red ,turquoise,blue,plum]: L :=[ -3.5, -2.5, -3 , -2.5 , -2.5, -2.5,-4.5, -5 , -3 , -4.5 ,-3.5 , -5 ]: U :=[ 3.5, 2.5 , 3 , 2.5 , 2.5, 2.5 , 5 , 5.5 , 3 , 5 , 3.5 , 5.5 ]: # r ul r ul ur ur r ul # l br l br bl bl l br for i from 1 to 12 do T||i:=spacecurve([P[i][1]+t*V1[i],P[i][2]+t*V2[i],P[i][3]+t*V3[i]] ,t=L[i]..U[i], color=CO[i],thickness=2,numpoints=2): od: M:=[0,0,0,0, .6,.6,.5,.4,0,0,0]: for i from 5 to 8 do T||i||1:=spacecurve([P[i][1]+t*V1[i],P[i][2]+t*V2[i],P[i][3]+t*V3[i]] ,t=L[i]..M[i], color=CO[i],thickness=2,numpoints=2): T||i||2:=spacecurve([P[i][1]+t*V1[i],P[i][2]+t*V2[i],P[i][3]+t*V3[i]] ,t=M[i]..U[i], color=CO[i],thickness=2,numpoints=2): od: AA:=[ S2,####### pink # T1,T3,T9,T11, # black# tan# red# blue # S4,####### cyan # T51,T61, T71,T81,# brown # coral # gray# violet # S1,####### magenta # T2,T4,T10, T12,# green# aquamarine # turquoise # plum # S3,####### yellow # T52,T62,T72,T82, # brown # coral# gray# violet T82, T12, NULL]: # # Draws pictures # #plotsetup(ps,plotoutput=`1.eps`, # plotoptions=`color,portrait,width=3in,height=3in,border`); #plotsetup(gif,plotoutput="disjoint.gif",plotoptions=`height=2500,width=2500`): plotsetup(gif,plotoutput="disjoint.gif",plotoptions=`height=1200,width=1200`): display3d(AA[],view=[-6..6,-6..6,-3..3],orientation=[30,60],scaling=CONSTRAINED); #display3d(AA[],view=[-3..3,-3..3,-2..2],orientation=[30,60],scaling=CONSTRAINED); #display3d(AA[],view=[-2..2,-2..2,-2..2],orientation=[30,60],scaling=CONSTRAINED);