First, the alpha functions should not be scaled by onedp here:
|
ax=kx*tw%alfa(1,1) * onedp -tw%alfa(1,2) * onedp /kx; |
Second, the algebra is wrong here:
|
cosvp = (ax*ay-bx*cy)/(ay*ay+cy*cy) |
cosvp = (ax*ay-bx*cy)/(ay*ay+cy*cy) ! ok
sinvp = (ax*cy+ay*bx)/(ay*ay+cy*cy) ! ok
cosvm = (ax*ay+by*cx)/(ax*ax+cx*cx) ! bug: (ax*ay + cx*by)/(ay*ay + by*by)
sinvm = (ax*by-ay*cx)/(ax*ax+cx*cx) ! bug: (ay*cx - ax*cy)/(ay*ay + by*by)
cosv2 = sqrt((1 + cosvp*cosvm - sinvp*sinvm)/2) ! bug: ± sqrt((1 + cosvp*cosvm + sinvp*sinvm)/2)
sinv2 = -sqrt((1 - cosvp*cosvm + sinvp*sinvm)/2) ! bug: ± sqrt((1 - cosvp*cosvm - sinvp*sinvm)/2)
cosv1 = -sqrt((1 + cosvp*cosvm + sinvp*sinvm)/2) ! bug: ± sqrt((1 + cosvp*cosvm - sinvp*sinvm)/2)
sinv1 = -sqrt((1 - cosvp*cosvm - sinvp*sinvm)/2) ! bug: ± sqrt((1 - cosvp*cosvm + sinvp*sinvm)/2)
Third, the disambiguation of the signs for the four square roots is complicated, it would be better to let the complex functions do the job (usually by going back and forth to polar coordinates).