restart;<Text-field layout="Heading 1" style="Heading 1">1. Plotting the Quadric and the Septic and Computing a Solution </Text-field>Load coefficient list and build quadric and septicCoeffsQS:=readdata(cat(prefix,"qscoeffsliftedmaple.m2"),integer,42):
CoeffsQ:=Matrix([seq(CoeffsQS[1][i] / CoeffsQS[2][i], i = 1.. 6)]):
MonomsQ:=Matrix([[1],[y],[x],[y^2],[x*y],[x^2]]):
CoeffsS:=Matrix([seq(CoeffsQS[1][i] / CoeffsQS[2][i], i = 7..42)]):
MonomsS:=Transpose(Matrix([[1, y, x, y^2, x*y, x^2, y^3, x*y^2, x^2*y, x^3, y^4, x*y^3, x^2*y^2, x^3*y, x^4, y^5, x*y^4, x^2*y^3, x^3*y^2, x^4*y, x^5, y^6, x*y^5, x^2*y^4, x^3*y^3, x^4*y^2, x^5*y, x^6, y^7, x*y^6, x^2*y^5, x^3*y^4, x^4*y^3, x^5*y^2, x^6*y, x^7]])):
Q:=(CoeffsQ.MonomsQ)[1,1];
S:=(CoeffsS.MonomsS)[1,1];Plot the Quadric and the Septic ...Digits:=8:
xran := -30..20:
yran := -20..20:
with(plots):
plotQS:=implicitplot([Q,S], x=xran, y=yran, numpoints = 50000, color = [red, blue], axes = boxed, thickness = 2):
display(plotQS);
Compute a solution Digits:=10:
SOL := fsolve({Q,S},{x,y});
filenamex:=cat(prefix,"solx.m2"):
filenamey:=cat(prefix,"soly.m2"):
SOLS:=subs(SOL,[x,y]): solsx:=[SOLS[1]]: solsy:=[SOLS[2]]:
writedata(filenamex, solsx); writedata(filenamey, solsy);
Z:=(sscanf(readline(cat(prefix,"g2value.m2")), %a))[1];xx := solve(S,x):
h:=numer(simplify(expand(subs(x =xx, Q)))):
factor(h);Ysol:=[fsolve(h, y)];
Sols:=[seq([subs(y = Ysol[i], xx),Ysol[i],Z], i = 1..nops(Ysol))];SXTC:=readline(cat(prefix,"SEXTICmaple.m2")):
sextic:=(sscanf(SXTC, %a))[1]:n = 4 looks goodsextsub:=subs({X = Sols[1][1], Y = Sols[1][2], Z = Sols[1][3]}, sextic):
implicitplot(subs(U =1, sextsub), V = -10..10, W = -10..10, numpoints = 10000);
implicitplot(subs(V =1, sextsub), U = -10..10, W = -10..10, numpoints = 10000);
implicitplot(subs(W =1, sextsub), U = -10..10, V = -10..10, numpoints = 10000);implicitplot(subs(V =1, sextsub), U = 0..2, W = 1..5, numpoints = 100000, axes = boxed);<Text-field layout="Heading 1" style="Heading 1">2. Discussing and Sampling the Plane Modell of the Coppler Curve</Text-field>Load coefficient list and build sextic...Digits:=20;
CoeffsSxtc := Matrix(readdata(cat(prefix,"sxtcmaple.m2"),28)):
nrmlz:=max( seq(abs(CoeffsSxtc[1,i]), i = 1..28));
CoeffsSxtc := (1/nrmlz)*CoeffsSxtc:
CoeffsLine := Matrix(readdata(cat(prefix,"linemaple.m2"),3)):
MonomsSxtc:= Transpose(Matrix([[ w^6, v*w^5, u*w^5, v^2*w^4, u*v*w^4, u^2*w^4, v^3*w^3, u*v^2*w^3, u^2*v*w^3, u^3*w^3, v^4*w^2, u*v^3*w^2, u^2*v^2*w^2, u^3*v*w^2, u^4*w^2, v^5*w, u*v^4*w, u^2*v^3*w, u^3*v^2*w, u^4*v*w, u^5*w, v^6, u*v^5, u^2*v^4, u^3*v^3, u^4*v^2, u^5*v, u^6]])):
MonomsLine:= Matrix([[w],[v],[u]]):
Sxtc := (CoeffsSxtc.MonomsSxtc)[1,1];
Line := (CoeffsLine.MonomsLine)[1,1];plotsetup(gif, plotoutput=`quadricseptic.gif`,
plotoptions=`height=1024,width=1024`);Lets take a look at our sextic...Digits:=8:
graphv:=implicitplot([subs(v=1,Sxtc),subs(v=1,Line)], u = -8..6, w =-2..30, numpoints=100000, color=[red,black], axes = boxed, thickness=2 ):
display(graphv);Now lets sample the curve...euclnorm := proc(M)
return(sqrt(M[1,1]^2 + M[1,2]^2));
end proc:newton:=proc(F,x,y, P0, kmax, tol)
local X, Xold, Xnew, DFpsinv,k, DfpiX, DeltaXsim, DeltaXnat, Done;
Xold := P0:
DFpsinv := MatrixInverse(Matrix([diff(F,x),diff(F,y)]));
for k from 0 to kmax do
DfpiX := Transpose(evalf(subs({x = Xold[1,1], y = Xold[1,2]},DFpsinv))):
DeltaXsim := -evalf(DfpiX*evalf(subs({x = Xold[1,1], y = Xold[1,2]},F))):
Xnew := Xold + DeltaXsim:
X[k] := Xnew:
DeltaXnat := DfpiX*evalf(subs({x = Xnew[1,1], y = Xnew[1,2]},F)):
Xold := Xnew;
if evalb(euclnorm(DeltaXsim) < tol) then
Done:=true:
break:
end if:
if evalb(euclnorm(DeltaXnat) > euclnorm(DeltaXsim)) then
Done:=false:
break:
end if:
od:
return(Done,k,Xnew):
end proc:
continuation:=proc(F,x,y,P0,Pf,kmax,tol,imax,s0,sMin,tol2)
local newPt, tngt, Xtmp,i , s, Xtmp2, Xs, tolMin;
tolMin:= 100;
tngt := - Matrix([-diff(F,y), diff(F,x)])* (1/sqrt(diff(F,x)^2 + diff(F,y)^2));
newPt := newton(F,x,y,P0,kmax, tol):
Xtmp := newPt[3];
if evalb(newPt[1] = true) then
for i from 0 to imax do
s := s0;
Xtmp2 := newPt[3] + s*evalf(subs({x= Xtmp[1,1], y = Xtmp[1,2]}, tngt));
newPt := newton(F,x,y,Xtmp2, kmax, tol);
Xtmp := newPt[3];
Xs[i] := newPt[3];
if newPt[1] = false then
print(cat("newtons method failed at i = ", i));
break
end if:
if tolMin > evalf(euclnorm(newPt[3] - Pf)) then
tolMin:=evalf(euclnorm(newPt[3] - Pf));
fi:
if tolMin < tol2 then
print(cat("arrived at Pf at i = ", i));
break
end if:
od:
end if:
print(tolMin):
return([seq([Xs[k][1,1], Xs[k][1,2]], k = 0.. i-1)] );
end proc:
Find the singularities (the cubic is copied and pasted from m2)Digits:=10;
sing1:=fsolve({subs(v=1, Line), subs(v=1, Sxtc)}, {u,w}, {u = 1..2, w = 3..7});
sing2:=fsolve({subs(v=1, Line), subs(v=1, Sxtc)}, {u,w}, {u = 1..2, w = 1..3});
sing3:=fsolve({subs(v=1, Line), subs(v=1, Sxtc)}, {u,w}, {u = 1..2, w = -1..1});fsolve(subs({u=1.61,v=1}, Sxtc), w , w= 2..4);Sampling the loopsF:= subs(u = x, v = 1, w = y, Sxtc):Digits:=20;
P0 := Matrix([1.8064, max(fsolve(subs({u=1.806,v=1}, Sxtc), w))]);
Pf := subs(sing1, Matrix([u,w]));
RES1 := continuation(F,x,y,P0, Pf, 10, 0.001, 2000, 0.05, 10, 10^(-3)):nops(RES1);Digits:=20;
P0 := Matrix([1.81, max(fsolve(subs({u=1.81,v=1}, Sxtc), w, w = 4..6))]):
Pf := subs(sing2, Matrix([u,w]));
RES2 := continuation(F,x,y,P0, Pf, 10, 0.0001, 400, -0.01, 10, 10^(-2)):
nops(RES2);Digits:=20:
P0 := Matrix([1.61, max(fsolve(subs({u=1.61,v=1}, Sxtc), w, w = 2..4))]):
Pf := subs(sing2, Matrix([u,w]));
RES3 := continuation(F,x,y,P0, Pf, 10, 0.0001, 400, 0.01, 10, 10^(-2)):nops(RES3);Digits:=20;
P0 := Matrix([1.63, max(fsolve(subs({u=1.63,v=1}, Sxtc), w, w = 2..4))]):
Pf := subs(sing1, Matrix([u,w]));
RES4 := continuation(F, x, y, P0, Pf, 10, 0.0001, 400, - 0.01, 10, 10^(-2)):nops(RES4);Digits:=20;
P0 := Matrix([1.5, min(fsolve(subs({u=1.5, v=1}, Sxtc), w, w = -1..1))]);
Pf := subs(sing3, Matrix([u,w]));
RES5 := continuation(F, x, y, P0, Pf, 10, 0.00001, 300, 0.05, 10, 10^(-2)):nops(RES5);Digits:=20;
P0 := Matrix([1.49, min(fsolve(subs({u = 1.49, v = 1}, Sxtc), w, w = -1..1))]);
Pf := subs(sing3,Matrix([u,w]));
RES6 := continuation(F, x, y, P0, Pf, 10, 0.00001, 300, - 0.05, 10, 10^(-2)):nops(RES6);nops(RES1) + nops(RES2) + nops(RES3) + nops(RES4); p0 :=plot([[P0[1,1], P0[1,2]]],color = black, style = point):
pts1:=plot(RES1, color = blue , style = line, thickness = 3):
pts2:=plot(RES2, color = magenta , style = line, thickness = 3):
pts3:=plot(RES3, color = green , style = line, thickness = 3):
pts4:=plot(RES4, color = yellow , style = line, thickness = 3):
display(pts1,pts2, pts3, pts4);Write Data to file... coppler := fopen(cat(prefix,"coppler.m2"), WRITE):
fprintf(coppler,"{",0):
for i from 1 to nops(RES1) do
fprintf(coppler,"matrix{{%.50f , %.50f , %.50f}} , \n", RES1[i][1], 1, RES1[i][2]);
od:
for i from 1 to nops(RES2) do
fprintf(coppler,"matrix{{%.50f , %.50f , %.50f}} , \n", RES2[i][1], 1, RES2[i][2]);
od:
for i from 1 to nops(RES3) do
fprintf(coppler,"matrix{{%.50f , %.50f , %.50f}} , \n", RES3[i][1], 1, RES3[i][2]);
od:
for i from 1 to nops(RES4) do
fprintf(coppler,"matrix{{%.50f , %.50f , %.50f}} , \n", RES4[i][1], 1, RES4[i][2]);
od:
fprintf(coppler,"null }",0);
fclose(coppler);coppler2 := fopen(cat(prefix,"coppler2.m2"), WRITE):
fprintf(coppler2,"{",0):
for i from 1 to nops(RES5) do
fprintf(coppler2,"matrix{{%.50f , %.50f , %.50f}} , \n", RES5[i][1], 1, RES5[i][2]);
od:
for i from 1 to nops(RES6) do
fprintf(coppler2,"matrix{{%.50f , %.50f , %.50f}} , \n", RES6[i][1], 1, RES6[i][2]);
od:
fprintf(coppler2,"null }",0);
fclose(coppler2);<Text-field layout="Heading 1" style="Heading 1">3. Plotting the Mechanism</Text-field>with(plottools): with(VectorCalculus):Radius:=0.1:Ls:=readdata(cat(prefix,"copplercoord.m2"), 18): nops(Ls);
GLs := readdata(cat(prefix,"statlegs.m2"), 18):GLs[1] - Ls[1];Qs:=[seq( Ls[1][3*i+1..3*(i+1)], i = 0..5)];
Qcenter := 1/6*(Qs[1] + Qs[2] + Qs[3] + Qs[4] + Qs[5] + Qs[6]);
Qsnew := [ Qs[1] - Qcenter, Qs[2] - Qcenter, Qs[3] - Qcenter, Qs[4] - Qcenter, Qs[5] - Qcenter, Qs[6] - Qcenter];
plane := implicitplot3d( z + 4, x = -6..4, y = -2.1..3, z = -6..4, color = gray, scaling = constrained, style = patchnogrid):
L1:= line(Qsnew[1], [Qsnew[1][1], Qsnew[1][2], -4], color = blue, thickness = 3):
L2:= line(Qsnew[2], [Qsnew[2][1], Qsnew[2][2], -4], color = blue, thickness = 3):
L3:= line(Qsnew[3], [Qsnew[3][1], Qsnew[3][2], -4], color = blue, thickness = 3):
L4:= line(Qsnew[4], [Qsnew[4][1], Qsnew[4][2], -4], color = blue, thickness = 3):
L5:= line(Qsnew[5], [Qsnew[5][1], Qsnew[5][2], -4], color = blue, thickness = 3):
L6:= line(Qsnew[6], [Qsnew[6][1], Qsnew[6][2], -4], color = blue, thickness = 3):
Q1:= sphere(Qsnew[1],0.1):
Q2:= sphere(Qsnew[2],0.1):
Q3:= sphere(Qsnew[3],0.1):
Q4:= sphere(Qsnew[4],0.1):
Q5:= sphere(Qsnew[5],0.1):
Q6:= sphere(Qsnew[6],0.1):
baseplatf := plane,L1,L2,L3,L4,L5,L6,Q1,Q2,Q3,Q4,Q5,Q6:
plotPlatform:=proc(Ps,Qs,basecol, topcol, legcol)
local Pp1, Pp2, Pp3, Pp4, Pp5, Pp6, Pq1, Pq2, Pq3, Pq4, Pq5, Pq6, Pps, Pqs, Qm, Pm,
D1,D2,D3,D4,D5,D6,C1,C2,C3,C4,C5,C6,L1,L2,L3,L4,L5,L6, topplatf,baseplatf,legs, normal,T1,T2,T3,T4,T5,T6,T7,T8,Pext, Pp7,Cadd1,Cadd2,Cadd3:
Pp1 := sphere(Ps[1],Radius):
Pp2 := sphere(Ps[2],Radius):
Pp3 := sphere(Ps[3],Radius):
Pp4 := sphere(Ps[4],Radius):
Pp5 := sphere(Ps[5],Radius):
Pp6 := sphere(Ps[6],Radius):
Pps := Pp1, Pp2, Pp3, Pp4, Pp5, Pp6:
Pm := (1/6)*(Ps[1] + Ps[2] + Ps[3] + Ps[4] + Ps[5] + Ps[6]):
normal := convert(CrossProduct(Vector(Ps[1] - Ps[2]), Vector(Ps[1] - Ps[3])),list);
Pext := Pm + (3/sqrt(normal[1]^2 + normal[2]^2 + normal[3]^2))*normal;
Pp7 := sphere(Pext, Radius-0.05, color = gray, style = patchnogrid):
Cadd1 := line(Ps[2],Pext, color = gray, thickness = 2);
Cadd2 := line(Ps[4],Pext, color = gray, thickness = 2);
Cadd3 := line(Ps[6],Pext, color = gray, thickness = 2);
T1 := polygon([Ps[1],Ps[2],Pm], color = red):
T2 := polygon([Ps[1],Ps[6],Pm], color = red):
T3 := polygon([Ps[5],Ps[4],Pm], color = red):
T4 := polygon([Ps[2],Ps[5],Pm], color = red):
T5 := polygon([Ps[4],Ps[6],Pm], color = red):
topplatf:= Cadd1, Cadd2, Cadd3, Pp7, T1,T2, T3, T4, T5,Pps;
L1 :=line(Ps[1],Qs[1],color = legcol, thickness = 3):
L2 :=line(Ps[2],Qs[2],color = legcol, thickness = 3):
L3 :=line(Ps[3],Qs[3],color = legcol, thickness = 3):
L4 :=line(Ps[4],Qs[4],color = legcol, thickness = 3):
L5 :=line(Ps[5],Qs[5],color = legcol, thickness = 3):
L6 :=line(Ps[6],Qs[6],color = legcol, thickness = 3):
legs:= L1,L2,L3,L4,L5,L6;
return(topplatf,legs);
end proc:Ps:=[seq( Ls[2][3*i+1..3*(i+1)], i = 0..5)];Qs:=[seq( Ls[1][3*i+1..3*(i+1)], i = 0..5)];GLsplatfrm:=[seq(display(baseplatf,plotPlatform([seq( Ls[j][3*i+1..3*(i+1)] - Qcenter, i = 0..5)], Qsnew,blue,red,black)), j = 3..3)]:
display(platfrm[1]);Loading the trace curvestrpnt1:=readdata(cat(prefix,"trace.m2"), 3):
trcurve1:=pointplot3d(trpnt1, color = green, style = line, thickness = 4):
conline1:=pointplot3d([trpnt1[85],trpnt1[1]], color = green, style = line, thickness = 4):
trpnt2:=readdata(cat(prefix,"trace2.m2"), 3):
trcurve2:=pointplot3d(trpnt2, color = cyan, style = line, thickness = 4):
conline2:=pointplot3d([trpnt2[82],trpnt2[1]], color = cyan, style = line, thickness = 4):
extr0 := readdata(cat(prefix,"exttrace.m2"), 3):
extr := [seq(extr0[i] - Qcenter, i = 1.. nops(extr0))]:
trcurve3:= pointplot3d(extr, color = green, style = line, thickness = 4):
conline3:= pointplot3d([extr[nops(extr)],extr[1]], color = green, style = line, thickness = 4):
extr20 := readdata(cat(prefix,"exttrace2.m2"), 3):
extr2 := [seq(extr20[i] - Qcenter, i = 1.. nops(extr20))]:
trcurve4:= pointplot3d(extr2, color = cyan, style = line, thickness = 4):
conline4:= pointplot3d([extr2[nops(extr2)],extr2[1]], color = cyan, style = line, thickness = 4):
Setting up the plot optionsplotsetup(gif, plotoutput="/Users/floriangeiss/plot4/test"||(1000)||".gif",
plotoptions=`color,portrait,noborder,height=1024,width=1024`);display(plane,platfrm[2],trcurve4, conline4, orientation = [0,80], scaling = constrained, axes = boxed);plotsetup(bmp, plotoutput="/Users/floriangeiss/plot8/test"||(1)||".bmp",
plotoptions=`color,portrait,noborder,height=2048,width=2048`);
display(platfrm[1],trcurve3,conline3, trcurve4,conline4, orientation = [-28,85],view = [-6..6,-6..6,-6..6], scaling = constrained);nops(platfrm);k:=1;
for i from 1 to 1 do
plotsetup(gif, plotoutput="/Users/floriangeiss/plotpaper/shot"||(3000+i)||".gif",
plotoptions=`color,portrait,noborder,height=1024,width=1024`);
display(platfrm[i],trcurve3,conline3, trcurve4, conline4, orientation = [20,80],view = [-6..6,-6..6,-6..6], scaling = constrained);
plotsetup(inline);
end do;
i := 70; k := 1;
plotsetup(gif, plotoutput="/Users/floriangeiss/plot8/test"||(1002)||".gif",
plotoptions=`color,portrait,noborder,height=1024,width=1024`);
display(platfrm[i],trcurve3, conline3,trcurve4,conline4, orientation = [k*90 + 90/nops(Ls)*i,80],view = [-6..6,-6..6,-6..6], scaling = constrained);
plotsetup(inline);for i from 83 to 90 do
plotsetup(gif, plotoutput="/Users/floriangeiss/plot2/comp2"||(1090+i)||".gif",
plotoptions=`color,portrait,noborder`);
display(platfrm[82],trcurve1,trcurve2,conline1,conline2, orientation = [180 + i/2,45]);
plotsetup(inline);
end do;k:=8;
for i from 1 to 83 do
plotsetup(gif, plotoutput="/Users/floriangeiss/plot3/comp"||(1000+(k-1)*90 + i)||".gif",
plotoptions=`color,portrait,noborder`);
display(platfrm[i],trcurve1,conline1, trcurve2, conline2, orientation = [(k - 1)*45 + i/2,45]);
plotsetup(inline);
end do;
for i from 84 to 90 do
plotsetup(gif, plotoutput="/Users/floriangeiss/plot3/comp"||(1000+(k-1)*90 + i)||".gif",
plotoptions=`color,portrait,noborder`);
display(platfrm[82],trcurve1,conline1, trcurve2, conline2, orientation = [(k - 1)*45 + i/2,45]);
plotsetup(inline);
end do;
i1 := 1; i2 := 2; i3 := 6;
nrm:=convert(CrossProduct(Vector(Qsnew[i1] - Qsnew[i2]) , Vector(Qsnew[i1] - Qsnew[i3])),list):
eqpl:=nrm[1]*x + nrm[2]*y + nrm[3]*z - (nrm[1]*Qsnew[i1][1] + nrm[2]*Qsnew[i1][2] +nrm[3]*Qsnew[i1][3]);
subs({x = Qsnew[i1][1], y = Qsnew[i1][2], z = Qsnew[i1][3]}, eqpl);
subs({x = Qsnew[i2][1], y = Qsnew[i2][2], z = Qsnew[i2][3]}, eqpl);
subs({x = Qsnew[i3][1], y = Qsnew[i3][2], z = Qsnew[i3][3]}, eqpl);
plane:=implicitplot3d( eqpl, x= -6..3.5, y = -2.1..2.8, z =-4..1.1, color = white):
plotsetup(gif, plotoutput="/Users/floriangeiss/plot4/test"||(1000)||".gif",
plotoptions=`color,portrait,noborder,height=1024,width=1024`);
display(platfrm[2],trcurve3, conline3, orientation = [90,0], view = [-6..6,-6..6,-6..6], scaling = constrained);Qsnewimplicitplot3d( eqpl, x= -5..4, y = -2..3, z =-3..2, color = gray, axes = boxed, orientation = [160,45]);nrm[1]/nrm[3];nrmplane := implicitplot3d( z + 4, x = -6..4, y = -2.1..3, z = -4..4, color = gray, scaling = constrained): display(plane,L1,L2,L3,L4,L5,L6,Q1,Q2,Q3,Q4,Q5,Q6, orientation = [20,80]);plane := implicitplot3d( z + 4, x = -6..4, y = -2.1..3, z = -4..4, color = gray, scaling = constrained):
L1:= line(Qsnew[1], [Qsnew[1][1], Qsnew[1][2], -4], color = blue, thickness = 3):
L2:= line(Qsnew[2], [Qsnew[2][1], Qsnew[2][2], -4], color = blue, thickness = 3):
L3:= line(Qsnew[3], [Qsnew[3][1], Qsnew[3][2], -4], color = blue, thickness = 3):
L4:= line(Qsnew[4], [Qsnew[4][1], Qsnew[4][2], -4], color = blue, thickness = 3):
L5:= line(Qsnew[5], [Qsnew[5][1], Qsnew[5][2], -4], color = blue, thickness = 3):
L6:= line(Qsnew[6], [Qsnew[6][1], Qsnew[6][2], -4], color = blue, thickness = 3):
Q1:= sphere(Qsnew[1],0.1):
Q2:= sphere(Qsnew[2],0.1):
Q3:= sphere(Qsnew[3],0.1):
Q4:= sphere(Qsnew[4],0.1):
Q5:= sphere(Qsnew[5],0.1):
Q6:= sphere(Qsnew[6],0.1):
baseplatf := plane,L1,L2,L3,L4,L5,L6,Q1,Q2,Q3,Q4,Q5,Q6:Qs[4];Ps[4];with(plottools):Pts:=[seq( GLs[1][3*i+1..3*(i+1)] - Qcenter, i = 0..5)];
F := 1/6* ( Pts[1] + Pts[2] + Pts[3] + Pts[4] + Pts[5] + Pts[6]);T1 := polygon([Pts[1],Pts[2],F], color = red):
T2 := polygon([Pts[1],Pts[6],F], color = green):
T3 := polygon([Pts[5],Pts[4],F], color = blue):
T4 := polygon([Pts[2],Pts[5],F], color = gray):
T5 := polygon([Pts[4],Pts[6],F], color = yellow):display(T1,T2,T2,T3,T4,T5);