BeginPackage["FourBar`"] Coupler::usage = "Coupler[r1, r2, r3, r4, r5, alpha, t2, sgn, options] draws the coupler curve traced by the coupler point P. The coupler point has polar coordinates r5, alpha(in degrees) wrt point B of the mechanism(see usage notes). r1 is the length of the fixed link, r2 the length of the input link, r3 is the length of the coupler, and r4 the length of the output link. sgn is either 1 or -1 and determines whether the linkage is in its open or crossed configuration. The figure drawn includes a picture of the linkage at the position with input link angle equal to t2(in degrees)." FourBar::usage = "FourBar[r1, r2, r3, r4, sgn, options] plots the angular velocities and accelerations of the links 3, 4 when the input link is driven at a constant angular velocity. r1 is the length of the fixed link, r2 the length of the input link, r3 is the length of the coupler, and r4 the length of the output link. sgn is either 1 or -1 and determines whether the linkage is in its open or crossed configuration. The velocity and acceleration of a coupler point may also be plotted by setting the PointMotion option." Angles::usage = "Angles[r1, r2, r3, r4][theta2] returns the angles of the links 2, 3, 4 in degrees for the position when the angle of the input link 2 is theta2 in degrees. The first set is for sgn = -1 and the second set for sgn = 1. See also usage message for Coupler." InputRange::usage = "InputRange[r1, r2, r3, r4] returns the possible range of values of the input link angle in degrees. See also usage message for Coupler." PSFile::usage = "PSFile is an option for Coupler and FourBar. In Coupler, its default value is None, indicating that no PostScript file is to be saved; if specified as a string fname, Coupler will save the coupler curve and the linkage picture to a PostScript file fname.ps. In FourBar, the default value is {None, None}, meaning no PostScript files are to be saved; if PSFile is specified as {\"fname\", None}, the angular velocity plot is saved in a PostScript file fname.ps1 and the the angular acceleration plot is saved in fname.ps2; if PSFile is given as {None, \"filename\"}, then the velocity and acceleration plots for a specified point of the coupler link is saved to the PostScript file filename.ps. PSFile may also be {\"fname1\", \"fname2\"} to save both angular and linear velocity and acceleration plots." PointMotion::usage = "PointMotion is an option for FourBar with default setting None. If PointMotion is given as {r, theta}, theta in degrees, defining a point on the coupler wrt point B, then FourBar calculates its velocity and acceleration over a cycle of operation of the linkage." Coupler::input = "Input angle is not in `1`." Coupler::nossm = "Linkage has incompatible dimensions." Coupler::limits = "Error in input angle limit calculation." Coupler::invopt = "`1` is not a valid option for Coupler." FourBar::invopt = "`1` is not a valid option for FourBar." Options[FourBar] = {PSFile -> {None, None}, PointMotion -> None} Options[Coupler] = {PSFile -> None} Begin["`Private`"] Coupler[r1_, r2_, r3_, r4_, r5_, alpha_, t2_, sgn:1 | -1, opts___Rule] := Block[{treg, interval, psfile, angles, pts, alfa, rcp, rca, rmax, G, graph}, treg = crlim[r1, r2, r3, r4]; If[ treg === None, Message[Coupler::limits]; Return[] ]; interval = Interval @@ Transpose[treg]; interval = IntervalUnion[ Interval[ First[IntervalIntersection[interval, Interval[{pi, 2 pi}]]] - 2 pi ], interval]; If[ Not[ IntervalMemberQ[interval, t2 degrees] ], Message[Coupler::input, interval/degrees]; Return[] ]; If[ r1 >= r2 + r3 + r4, Message[Coupler::nossm]; Return[] ]; psfile = PSFile /. {opts} /. Options[Coupler]; If[ Not[ MatchQ[psfile, None | _String] ], Message[Coupler::invopt, opts]; Return[] ]; angles = pos[r1, r2, r3, r4, sgn][ t2 degrees ]; pts = Append[ {r2, r3, r4} Transpose[{Cos[angles], Sin[angles]}], {-r1, 0} ]; alfa = alpha degrees + Part[angles, 2]; rcp = r5 {Cos[alfa], Sin[alfa]}; rca = Plus @@ Take[pts, 2]; rmax = Max[r1, r2, r3, r4]; G = {points[{0, 0}, pts], points[First[pts], {rcp}], psp[{0, 0}, rmax], psp[{r1, 0}, rmax], Text["B", Scaled[-0.02 pts[[2]]/r3, First[pts]], {0, 0}], Text["A", {-3 rmax/100, 0}, {1, 0}], Text["C", Scaled[0.02 pts[[2]]/r3, rca], {0, 0}], Text["D", {3 rmax/100, 0} + { r1, 0}, {-1, 0}], Text["P", Scaled[0.02 rcp/r5 + {0.005, 0.005}, rcp + First[pts]], {0, 0}] }; graph = {Graphics[ Rectangle[ Scaled[{0, 0}], Scaled[{0.9, 0.9}], Show[coupler[r1, r2, r3, r4, r5, alpha, sgn], Graphics[{Thickness[0.001], G}], AspectRatio -> Automatic, DisplayFunction -> Identity, PlotRange -> All]] ], Graphics[ Rectangle[ Scaled[{0, 0.9}], Scaled[{1, 1}], Graphics[ Text[StringJoin["Coupler curve (Linkage at input angle = ", filtp[t2, 2], ")"], Scaled[{0.5, 0.95}], {0, 0}], PlotRange -> All]] ]}; If[ psfile === None, Show[graph, DisplayFunction -> $DisplayFunction, PlotRange -> All], Display[StringJoin["!psfix > ", psfile, ".ps"], Show[graph, DisplayFunction -> Identity, PlotRange -> All]] ] ] /; And @@ numberQ[{r1, r2, r3, r4, r5, alpha, t2}] Angles[lengths__][t2_?numberQ] := Block[ {x2, treg}, x2 = t2 degrees; treg = crlim[lengths]; If[ treg === None, Message[Coupler::limits]; Return[] ]; interval = Interval @@ Transpose[treg]; interval = IntervalUnion[ Interval[ First[IntervalIntersection[interval, Interval[{pi, 2 pi}]]] - 2 pi ], interval]; If[ Not[ IntervalMemberQ[interval, x2 ] ], Message[Coupler::input, interval/degrees]; Return[] ]; If[ First[{lengths}] >= Plus @@ Rest[{lengths}], Message[Coupler::nossm]; Return[] ]; {pos[lengths, -1][x2], pos[lengths, 1][x2]}/degrees ] /; And @@ numberQ[{lengths}] InputRange[lengths__] := Block[ {treg}, treg = crlim[lengths]; If[ treg === None, Message[Coupler::limits]; Return[] ]; If[ First[{lengths}] >= Plus @@ Rest[{lengths}], Message[Coupler::nossm]; Return[] ]; (Interval @@ Transpose[treg])/degrees ] /; And @@ numberQ[{lengths}] pos[r1_, r2_, r3_, r4_, sgn_][t2_] := Block[ {x2, C1, C2, C3, t4}, x2 = N[t2]; C1 = (r4 - r1)( r4 - r1 + 2 r2 Cos[x2]) + r2^2 - r3^2; C2 = 4 r2 r4 Sin[x2]; C3 = C1 - 4 r4 (r2 Cos[x2] - r1); t4 = If[ C3 == 0, pi, 2 ArcTan[(- C2 + sgn (C2^2 - 4 C1 C3)^(1/2))/(2 C3)] ]; {x2, Arg[ r1 - r2 Exp[I x2] - r4 Exp[I t4] ], t4} ] pos[x__, T2_List] := Transpose[Map[ pos[x], T2 ]] coupler[r1_, r2_, r3_, r4_, r5_, alpha_, sgn_] := Block[ {indata}, indata = {0, alpha degrees} + Drop[ data[r1, r2, r3, r4, sgn], -1]; Map[ListPlot[ #1, PlotJoined -> True, AspectRatio -> Automatic, DisplayFunction -> Identity, Axes -> {False, True}, PlotRange -> All, PlotStyle -> Thickness[0.002] ] &, Partition[Transpose[Apply[Plus, {{r2, r5} Cos[indata], {r2, r5} Sin[indata]}, 1]], Length[First[indata]]/2] ] ] data[r1_, r2_, r3_, r4_, sgn_] := data[r1, r2, r3, r4, sgn] = Block[{treg, delta}, treg = crlim[r1, r2, r3, r4] + 0.0001 {{1, 1}, {-1, -1}}; delta = (Last[treg] - First[treg])/101; pos[r1, r2, r3, r4, sgn, Join[ NestList[Plus[#1, First[delta]]&, Part[treg, 1, 1], 101], NestList[Plus[#1, Last[delta]]&, Part[treg, 1, 2], 101]]] ] crlim[r1_, r2_, r3_, r4_]:= Block[ {x = r1^2 + r2^2, L, U}, {L, U} = N[ Map[ArcCos, {x -(r3 + r4)^2, x - (r3 - r4)^2}/(2 r1 r2) ], 30 ]; If[ MatchQ[U, _Complex], U = 0]; If[ MatchQ[L, _Complex], L = pi]; {{U, 2 pi - L}, {L, 2 pi - U}} ] /; And @@ Positive[{r1, r2, r3, r4}] crlim[___] := None $DefaultFont = {"Times-Italic", 12} psp[center_, rmax_] := Block[{r = rmax/100}, {Circle[center, r], Circle[center, 2r, {0, pi}], Line[{ center + {2r, 0}, center + {2r, -2r} }], Line[{ center - {2r, 0}, center - {2r, 2r}}], points[center -{3r, 2r}, {{6r, 0}, {0, -r/2}, {-6r, 0}, {0, r/2}}] } ] points[origin_, vlist_] := Line[FoldList[ Plus, origin, vlist ]] pi = N[Pi] degrees = N[Degree] filtp[v_Integer, n_] := ToString[v] filtp[v_, n_] := Block[{x, pos, len}, x = ToString[N[v]]; pos = First @@ StringPosition[x, "."]; len = StringLength[x]; If[pos == len, StringTake[x, len - 1], StringTake[x, Min[n + pos, len] ]] ] opts[title_String] := Sequence @@ {PlotJoined -> True, DisplayFunction -> Identity, PlotStyle -> Thickness[0.001], Frame -> True, FrameLabel -> {"Input angle(rad)", " ", StringJoin["Angular ", title, " Ratios"], " "}}; opts2[title_String] := Sequence @@ {PlotJoined -> True, AspectRatio -> Automatic, DisplayFunction -> Identity, Frame -> True, FrameLabel -> {" ", " ", title, " "}, PlotRange -> All, PlotStyle -> Thickness[0.002]}; auxtext[text_, pts_] := Text[FontForm[ToString[InputForm[text]], {"Helvetica", 10}], Scaled[{0, -0.02}, pts], {0, 0}] auxvec[{x_, y_}] := {y, -x} SetAttributes[numberQ, Listable] numberQ[x_] := NumberQ[N[x]] FourBar[r1_, r2_, r3_, r4_, sgn:1 | -1, options___Rule] := Block[ {pdata, a13, a21, a32, plength, index, r2r3, r2r4, r4r3, r3r4, w3, w4, loc, graph, psfile, r5, t5, w3sq, w4sq}, pdata = data[r1, r2, r3, r4, sgn]; {a13, a21, a32} = pdata - RotateRight[pdata]; plength = Length[a13]/2; index = Ceiling[{1, 4, 2, 5, 1, 4, 2, 5} plength/3]; r2r3 = r2/r3; r2r4 = r2/r4; r4r3 = r4/r3; r3r4 = r3/r4; psfile = PSFile /. {options} /. Options[FourBar]; r5 = PointMotion /. {options} /. Options[FourBar]; If[ Not [ Or @@ (MatchQ[psfile, #] & /@ Flatten[ Outer[List, {None, _String}, {None, _String}], 1])], Message[FourBar::invopt, StringJoin["PSFile -> " , ToString[InputForm[psfile]]]]]; If[ r5 =!= None && Not[ MatchQ[ r5, {_?Positive, _?numberQ}] ], Message[FourBar::invopt, StringJoin["PointMotion -> " , ToString[InputForm[r5]]]]; Return[]]; {w3, w4} = Inner[ Divide, {r2r3 Sin[a13], r2r4 Sin[a21]}, Sin[a32], List ]; w3sq = w3^2; w4sq = w4^2; w4 = Join[ {w3, w4}, Inner[ Divide, Inner[Times, {w3sq, -w4sq}, Cos[a32], List] + {r2r3 Cos[a13] + r4r3 w4sq, - r2r4 Cos[a21] - w3sq r3r4}, Sin[a32], List ] ]; a21 = Part[w4, 3]; loc = MapThread[auxtext, {{3, 3, 4, 4, 3, 3, 4, 4}, Transpose[{Part[First[pdata], index], Flatten[MapThread[Part, {w4, Partition[index, 2]}]]}]}]; w4 = Join @@ (Transpose[{Partition[ First[pdata], plength ], #}] &) /@ (Partition[#, plength] & /@ w4); graph = {{Map[(ListPlot[Transpose[#], opts["Velocity"]] &), Take[w4, 4]], Graphics[Take[loc, 4], DisplayFunction -> Identity]}, {Map[(ListPlot[Transpose[#], opts["Acceleration"]] &), Take[w4, -4]], Graphics[Take[loc, -4], DisplayFunction -> Identity]}}; If[ First[psfile] === None, Map[Show[ #, DisplayFunction -> $DisplayFunction] &, graph], Display[StringJoin["!psfix > ", First[psfile], ".ps1"], Show[First[graph], DisplayFunction -> Identity]]; Display[StringJoin["!psfix > ", First[psfile], ".ps2"], Show[Last[graph], DisplayFunction -> Identity]] ]; If[r5 =!= None, {r5, t5} = r5 {1/r2, degrees}; a32 = r5 {-1, 1} Through[{Sin, Cos}[t5 + Part[pdata, 2]] ]; a13 = -Through[{Cos, Sin}[First[pdata]] ]; graph = {{Show[Map[ListPlot[ #1, opts2["Polar plot of coupler point velocity"]] &, Partition[Transpose[auxvec[a13] + Inner[Times, a32, w3, List]], plength] ], DisplayFunction -> Identity]}, {Show[Map[ListPlot[ #1, opts2["Polar plot of coupler point acceleration"]] &, Partition[Transpose[a13 - Inner[Times, auxvec[a32], w3sq, List] + Inner[Times, a32, a21, List]], plength] ], DisplayFunction -> Identity]}}; If[Last[psfile] === None, Show[ GraphicsArray[graph], DisplayFunction -> $DisplayFunction ], Display[StringJoin["!psfix > ", Last[psfile], ".ps"], Show[GraphicsArray[graph], DisplayFunction -> Identity]] ] ]; Null ] /; And @@ Join[ numberQ[{r1, r2, r3, r4}], Positive[{r1, r2, r3, r4}] ] End[] EndPackage[]