I've been bootstrapping myself to the very alien world of Mathematica and there came my first WTF moment:
RRP[P_, O_, Phi_, L_, M_] :=
If[d <= L, O + l2*v0, {}] //. {
v0 -> {Cos[Phi], Sin[Phi]},
OP -> P - O,
l1 -> Dot[OP, v0],
v1 -> l1*v0,
d -> Norm[OP - v1],
l2 -> l1 + M Sqrt[L^2 - d^2]
};
What I wanted to achieve is to recursively apply these rules(which I felt similar to a context-free grammar) to form the desired expression. There was no cyclic dependency in the rules, nor were alternatives, so I kinda expected it to expand correctly. But it didn't!
Edit: I'm trying to calculate the position of the internal joint in a planar RRP linkage. Expressed in procedural language it's probably like this(pseudocode):
/* P: position of external R joint
* O: reference point of P joint
* Phi: direction of P joint
* L: length of the R-R bar
* M: +1/-1
*/
function RRP(P, O, Phi, L, M){
OP=P-O;
v0=[cos(Phi), sin(Phi)];
l1=dot(OP, v0);
v1=l1*v0;
d=norm(OP-v1);
l2=l1+M*sqrt(L^2-d^2);
if(d<=L) return O+l2*v0;
else return [];
}
Since I haven't figured out how to even debug here, I tried to output intermediate variables(or non-terminals in CFG terminology):
RRP[P_, O_, Phi_, L_, M_] :=
{d, OP, v1, OP - v1} //. {
v0 -> {Cos[Phi], Sin[Phi]},
OP -> P - O,
l1 -> Dot[OP, v0],
v1 -> l1*v0,
d -> Norm[OP - v1],
l2 -> l1 + M Sqrt[L^2 - d^2]
};
RRP[{0, 1}, {5, 0}, 0, Sqrt[2], 1]
And it gave me this:
{{6, Sqrt[26]}, {-5, 1}, {-5, 0}, {{0, -5}, {6, 1}}}
OP
and v1
were as expected. But WTF was OP-v1
? A Matrix?
I figure I can't do much with OP
so I tried manually substituting v1
in the d->
rule:
RRP[P_, O_, Phi_, L_, M_] :=
{d, OP, v1, OP - l1*v0} //. {
v0 -> {Cos[Phi], Sin[Phi]},
OP -> P - O,
l1 -> Dot[OP, v0],
v1 -> l1*v0,
d -> Norm[OP - l1*v0],
l2 -> l1 + M Sqrt[L^2 - d^2]
};
RRP[{0, 1}, {5, 0}, 0, Sqrt[2], 1]
Now it gives me the correct answer:
{1, {-5, 1}, {-5, 0}, {0, 1}}
So now I'm VERY confused about this,
I don't even know if I'm following the proper/idiomatic way of constructing complex functions, thanks to the arcane Help that came with Mathematica... Could anyone please either point out the problem, or kindly enlighten me on how I could construct a equivalent function in "The Mathematica Way"?
EDIT: After some lookup I figured out that the evaluation order for OP - v1
was to blame. So I added Hold[]
around the culprit and //ReleaseHold
the whole thing after //.
finishes building the expression. It works:
RRP[P_, O_, \[Phi]_, L_, M_] := (
If[d <= L, O + l2*v0, {}] //. {
v0 -> {Cos[\[Phi]], Sin[\[Phi]]},
OP -> P - O,
l1 -> Dot[OP, v0],
v1 -> l1*v0,
d -> Norm[Hold[OP - v1]],
l2 -> l1 + M Sqrt[L^2 - d^2]
}) // ReleaseHold;
RRP[{0, 1}, {5, 0}, 0, Sqrt[2], 1]
But I am still confused:
Could I know where I need Hold[]
when I write the code, not as a workaround after I found out something wrong with it?