First of all please excuse the long wall of test coming up. I need to solve the following set of equations which arise from a coordinate transformation:
Given v_n, v_m, x_p, z_p, y_p, s, d, solve the system below for v_yp, v_zp:
I: v_n = v_yp * (1/(s*d)) + v_xp * y_p * (1/(s*d^2))
II: v_m = v_zp * (1/(s*d)) + v_xp * z_p * (1/(s*d^2))
III: v_xp = -1/Jp_11 * (Jp_12 * v_yp + Jp_13 * v_zp)
The terms Jp_ij are the entries of a matrix Jp and can be written as:
Jp_11 = sin(theta) * cos(phi) * cos(b_0) + cos(theta) * sin(b_0)
Jp_12 = sin(theta) * sin(phi)
Jp_13 = -sin(theta) * cos(phi) * sin(b_0) + cos(theta) * cos(b_0)
b_0 is given as well. Now I want to solve I,II,III for v_yp and v_zp. I tried the following commands in SymPy:
In [156]: EQ1 = v_n - (v_yp * (1/(s*d)) + yp * v_xp * (1/(s*d**2)))
In [157]: EQ2 = v_m - (v_zp * (1/(s*d)) + zp * v_xp * (1/(s*d**2)))
In [158]: EQ3 = v_xp + 1/Jp.row(0).col(0).tolist()[0][0] * (Jp.row(0).col(1).tolist()[0][0] * v_yp + Jp.row(0).col(2).tolist()[0][0] * v_zp)
In [159]: EQ1
Out[159]: v_n - v_yp/(d*s) - v_xp*yp/(d**2*s)
In [160]: EQ2
Out[160]: v_m - v_zp/(d*s) - v_xp*zp/(d**2*s)
In [161]: EQ3
Out[161]: v_xp + (v_yp*sin(phi)*sin(theta) + v_zp*(-sin(b0)*sin(theta)*cos(phi) + cos(b0)*cos(theta)))/(sin(b0)*cos(theta) + sin(theta)*cos(b0)*cos(phi))
In [162]: sympy.linsolve([EQ1,EQ2,EQ3], (v_xp, v_yp, v_zp))
Out[162]: {(0, d*s*v_n, d*s*v_m)}
In [163]: sympy.linsolve([EQ1,EQ2,EQ3], (v_yp, v_zp))
Out[163]: EmptySet()
In [164]: sympy.linsolve([EQ1,EQ2], (v_yp, v_zp))
Out[164]: {(d*s*v_n - v_xp*yp/d, d*s*v_m - v_xp*zp/d)}
In [165]: sympy.solve([EQ1,EQ2,EQ3], (v_yp, v_zp))
Out[165]: []
Since the third equation already sets v_xp and v_yp,v_zp in relation I'd like to solve the system in terms of v_yp and v_zp only. Line 163 does not give me the answer I want and line 162's output is unexpected: v_xp is a velocity in direction to the observer of the Sun and should be zero only when looking at the disk center (where by assumption the radial velocity, equalling v_xp (only at that point), was set to v_r = 0).
I also tried plugging in v_xp (III) manually into I and II.
In [175]: v_xp = - (v_yp*sympy.sin(phi)*sympy.sin(theta) + v_zp*(-sympy.sin(b0)*sympy.sin(theta)*sympy.cos(phi) + sympy.cos(b0)*sympy.cos(theta)))/(sympy.sin(b0)*sympy.cos(theta) + sympy.sin(theta)*sympy.cos(b0)*sympy.cos(phi))
In [176]: v_xp
Out[176]: (-v_yp*sin(phi)*sin(theta) - v_zp*(-sin(b0)*sin(theta)*cos(phi) + cos(b0)*cos(theta)))/(sin(b0)*cos(theta) + sin(theta)*cos(b0)*cos(phi))
In [178]: EQ4 = v_n - v_yp/(s*d) - v_xp*yp/(s*d**2)
In [179]: EQ4
Out[179]: v_n - v_yp/(d*s) - yp*(-v_yp*sin(phi)*sin(theta) - v_zp*(-sin(b0)*sin(theta)*cos(phi) + cos(b0)*cos(theta)))/(d**2*s*(sin(b0)*cos(theta) + sin(theta)*cos(b0)*cos(phi)))
In [181]: EQ5 = v_m - v_zp/(s*d) - v_xp*zp/(s*d**2)
In [182]: EQ5
Out[182]: v_m - v_zp/(d*s) - zp*(-v_yp*sin(phi)*sin(theta) - v_zp*(-sin(b0)*sin(theta)*cos(phi) + cos(b0)*cos(theta)))/(d**2*s*(sin(b0)*cos(theta) + sin(theta)*cos(b0)*cos(phi)))
In [183]: sympy.linsolve([EQ4,EQ5], (v_yp, v_zp))
Out[183]: {(d*s*v_n, d*s*v_m)}
Again we have the result that would correspond to v_xp = 0.
However, when doing the calculations manually, the expressions for v_yp, v_zp become some complicated terms. If needed I can post them later. I was doing the symbolic computation because I wanted to check the result and see whether it could still be simplified.
Solving the equations does not work as expected. Why?
PUSH: I tried the same task in MATLAB: Here is the code and the result. Note that I replaced Jp_ij by Jij for simplicity.
Code:
syms xp yp zp v_xp v_yp v_zp s d v_n v_m b0 phi theta J11 J12 J13;
EQ1 = v_xp == (-v_yp*sin(phi)*sin(theta) - v_zp*(-sin(b0)*sin(theta)*cos(phi) + cos(b0)*cos(theta)))/(sin(b0)*cos(theta) + sin(theta)*cos(b0)*cos(phi));
EQ2 = v_n == v_yp/(d*s) + v_xp*yp/(s*d^2);
EQ3 = v_m == v_zp/(s*d) + v_xp*zp/(s*d^2);
[sol_vyp,sol_vzp] = solve([EQ1,EQ2,EQ3],[v_yp,v_zp]);
sol_vyp, sol_vzp
[sol_vyp] = solve([EQ2],[v_yp]);
sol_vyp
EQ4 = v_n == v_yp/(d*s) + yp*(-v_yp*sin(phi)*sin(theta) - v_zp*(-sin(b0)*sin(theta)*cos(phi) + cos(b0)*cos(theta)))/(d^2*s*(sin(b0)*cos(theta) + sin(theta)*cos(b0)*cos(phi)));
EQ5 = v_m == v_zp/(d*s) + zp*(-v_yp*sin(phi)*sin(theta) - v_zp*(-sin(b0)*sin(theta)*cos(phi) + cos(b0)*cos(theta)))/(d^2*s*(sin(b0)*cos(theta) + sin(theta)*cos(b0)*cos(phi)));
[sol_vyp] = solve([EQ4],[v_yp]);
sol_vyp
[sol_vyp,sol_vzp] = solve([EQ4,EQ5],[v_yp,v_zp]);
sol_vyp, sol_vzp
EQ6 = v_n == v_yp/(d*s) + yp*(-v_yp*J12/J11 - v_zp*J13/J11)/(s*d^2);
EQ7 = v_m == v_zp/(d*s) + zp*(-v_yp*J12/J11 - v_zp*J13/J11)/(s*d^2);
[sol_vyp] = solve([EQ6],[v_yp]);
sol_vyp
[sol_vyp,sol_vzp] = solve([EQ6,EQ7],[v_yp,v_zp]);
sol_vyp, sol_vzp
Result:
>> test_transform
sol_vyp =
Empty sym: 0-by-1
sol_vzp =
Empty sym: 0-by-1
sol_vyp =
d*s*(v_n - (v_xp*yp)/(d^2*s))
sol_vyp =
(v_n + (v_zp*yp*(cos(b0)*cos(theta) - cos(phi)*sin(b0)*sin(theta)))/(d^2*s*(sin(b0)*cos(theta) + cos(b0)*cos(phi)*sin(theta))))/(1/(d*s) - (yp*sin(phi)*sin(theta))/(d^2*s*(sin(b0)*cos(theta) + cos(b0)*cos(phi)*sin(theta))))
sol_vyp =
(d*s*(d*v_n*sin(b0)*cos(theta) + v_m*yp*cos(b0)*cos(theta) - v_n*zp*cos(b0)*cos(theta) + d*v_n*cos(b0)*cos(phi)*sin(theta) - v_m*yp*cos(phi)*sin(b0)*sin(theta) + v_n*zp*cos(phi)*sin(b0)*sin(theta)))/(d*sin(b0)*cos(theta) - zp*cos(b0)*cos(theta) - yp*sin(phi)*sin(theta) + d*cos(b0)*cos(phi)*sin(theta) + zp*cos(phi)*sin(b0)*sin(theta))
sol_vzp =
(d*s*(v_n*zp*sin(phi)*sin(theta) - v_m*yp*sin(phi)*sin(theta) + d*v_m*sin(b0)*cos(theta) + d*v_m*cos(b0)*cos(phi)*sin(theta)))/(d*sin(b0)*cos(theta) - zp*cos(b0)*cos(theta) - yp*sin(phi)*sin(theta) + d*cos(b0)*cos(phi)*sin(theta) + zp*cos(phi)*sin(b0)*sin(theta))
sol_vyp =
(v_n + (J13*v_zp*yp)/(J11*d^2*s))/(1/(d*s) - (J12*yp)/(J11*d^2*s))
sol_vyp =
-(d*s*(J11*d*v_n + J13*v_m*yp - J13*v_n*zp))/(J12*yp - J11*d + J13*zp)
sol_vzp =
-(d*s*(J11*d*v_m - J12*v_m*yp + J12*v_n*zp))/(J12*yp - J11*d + J13*zp)
With exception of the first line, the results are as expected. I'd need some background info why this does not work in python though, please. As a last check, I also tried it the following way in python, keeping the matrix symbols:
In [220]: J11,J12,J13 = sympy.symbols('J11 J12 J13')
In [222]: EQ6 = v_n - v_yp/(d*s) - yp*(-v_yp*J12/J11 - v_zp*J13/J11)/(s*d**2)
In [223]: EQ7 = v_m - v_zp/(d*s) - zp*(-v_yp*J12/J11 - v_zp*J13/J11)/(s*d**2)
In [228]: sympy.linsolve([EQ6, EQ7], [v_yp, v_zp])
Out[228]: {(d*s*v_n, d*s*v_m)}
Once more, I don't know why SymPy does not give me the answer I expect.
PUSH: sympy.linsolve can solve relatively simple systems. Somehow at a certain level of complexity it cracks down, while sympy.solve keeps going.
In [236]: a, b, c, d, e, f = sympy.symbols('a b c d e f')
In [246]: sympy.linsolve([v_n - a*v_yp - yp*(b*v_yp)], v_yp)
Out[246]: {(v_n/(a + b*yp),)}
In [247]: sympy.linsolve([v_n - a*v_yp - yp*(b*v_yp + c * v_zp)], v_yp)
Out[247]: {(v_n/a,)}
In [245]: sympy.solve([v_n - a*v_yp - yp*(b*v_yp + c * v_zp)], v_yp)
Out[245]: {v_yp: (-c*v_zp*yp + v_n)/(a + b*yp)}
So with solve() I can retrieve the answers in the correct form.
In [255]: sympy.linsolve([EQ6], [v_yp])
Out[255]: {(d*s*v_n,)}
In [256]: sympy.solve([EQ6], [v_yp])
Out[256]: {v_yp: (J11*d**2*s*v_n + J13*v_zp*yp)/(J11*d - J12*yp)}
In [257]: sympy.solve([EQ6, EQ7], [v_yp, v_zp])
Out[257]:
{v_yp: d*s*(J11*d*v_n + J13*v_m*yp - J13*v_n*zp)/(J11*d - J12*yp - J13*zp),
v_zp: d*s*(J11*d*v_m - J12*v_m*yp + J12*v_n*zp)/(J11*d - J12*yp - J13*zp)}
In [260]: sympy.solve([EQ4, EQ5], [v_yp, v_zp])
Out[260]:
{v_yp: d*s*(d*v_n*sin(b0)*cos(theta) + d*v_n*sin(theta)*cos(b0)*cos(phi) - v_m*yp*sin(b0)*sin(theta)*cos(phi) + v_m*yp*cos(b0)*cos(theta) + v_n*zp*sin(b0)*sin(theta)*cos(phi) - v_n*zp*cos(b0)*cos(theta))/(d*sin(b0)*cos(theta) + d*sin(theta)*cos(b0)*cos(phi) - yp*sin(phi)*sin(theta) + zp*sin(b0)*sin(theta)*cos(phi) - zp*cos(b0)*cos(theta)),
v_zp: d*s*(d*v_m*sin(b0)*cos(theta) + d*v_m*sin(theta)*cos(b0)*cos(phi) - v_m*yp*sin(phi)*sin(theta) + v_n*zp*sin(phi)*sin(theta))/(d*sin(b0)*cos(theta) + d*sin(theta)*cos(b0)*cos(phi) - yp*sin(phi)*sin(theta) + zp*sin(b0)*sin(theta)*cos(phi) - zp*cos(b0)*cos(theta))}
