PTC_Twiss 2x2 coupling matrix is incorrect

[ldeniau]

First, the alpha functions should not be scaled by onedp here: https://github.com/MethodicalAcceleratorDesign/MAD-X/blob/88211ac4a0022d29b24a64374592fb72cadc72cf/src/madx_ptc_twiss.f90#L2179

Second, the algebra is wrong here: https://github.com/MethodicalAcceleratorDesign/MAD-X/blob/88211ac4a0022d29b24a64374592fb72cadc72cf/src/madx_ptc_twiss.f90#L2215

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).