Experiment with a scalar formulation

```julia
using Test
using Multibody
using ModelingToolkit
import ModelingToolkitStandardLibrary.Mechanical.Rotational
using OrdinaryDiffEq
using LinearAlgebra
using JuliaSimCompiler

t = Multibody.t
D = Differential(t)
world = Multibody.world
W(args...; kwargs...) = Multibody.world

@mtkmodel TestUSR begin
    @components begin
        world = W()
        j1 = JointUSR(positive_branch=true, use_arrays=false)
        fixed = FixedTranslation(r=[1,0,0])
        b1 = Body(isroot=false, neg_w=true)
        p1 = Prismatic(state_priority=100)
    end
    @equations begin
        connect(world.frame_b, j1.frame_a, fixed.frame_a)
        connect(fixed.frame_b, p1.frame_a)
        connect(p1.frame_b, j1.frame_b)
        connect(j1.frame_im, b1.frame_a)
    end
end

@nAmed model = TestUSR()
model = complete(model)
ss = structural_simplify(IRSystem(model))
prob = ODEProblem(ss, [model.b1.a_0[1]=>0.0, D(D(model.p1.s))=>0.0], (0.0, 1.0))
sol = solve(prob, FBDF(autodiff=true))
```
works

src/fancy_joints.jl

@@ -443,12 +443,7 @@ end
         0 .~ collect(frame_a.f .+ resolve1(Rrel, frame_b.f))
         0 .~ collect(frame_a.tau .+ resolve1(Rrel, frame_b.tau))
 
-        if use_arrays
-            angle ~ compute_angle2(length_constraint, e_array, r_a_array, r_b_array, positive_branch)[1]
-            # angle ~ Symbolics.term(compute_angle2, length_constraint, e_array, r_a_array, r_b_array, positive_branch, type=Real)
-        else
-            angle ~ compute_angle2(length_constraint, e, r_a, r_b, positive_branch)[1]
-        end
+        angle ~ _compute_angle2(length_constraint, e, r_a, r_b, positive_branch)
     ]
 
     sys = ODESystem(eqs, t; name=:nothing, systems)#, parameter_dependencies = [positive_branch => select_branch(length_constraint, e, phi_offset + phi_guess, r_a, r_b)])  # JuliaSimCompiler ignores parameter dependencies, the user has to provide it instead
@@ -739,27 +734,33 @@ function compute_angle(L::Real, e::AbstractVector, r_a::AbstractVector, r_b::Abs
     atan(ksin1, kcos1)
 end
 
-function compute_angle2(L::Real, e::AbstractVector, r_a::AbstractVector, r_b::AbstractVector, positive_branch)
+function _compute_angle2(L::Real, e::AbstractVector, r_a::AbstractVector, r_b::AbstractVector, positive_branch)
+    @show e,r_a
     e_r_a = e'r_a
     e_r_b = e'r_b
     A = -2*(r_b'r_a - e_r_b'e_r_a)
     B = 2*r_b'cross(e, r_a)
     C = r_a'r_a + r_b'r_b - L^2 - 2*e_r_b'e_r_a
     k1 = A^2 + B^2
     k1a = k1 - C^2
-    k1a > 1e-10 || error("Singular position of loop (either no or two analytic solutions; the mechanism has lost one-degree-of freedom in this position). Try to use another joint-assembly component. In most cases it is best to let joints outside of the JointXXX component be revolute and _not_ prismatic joints. If this also leads to singular positions, it could be that this kinematic loop cannot be solved analytically. In this case you have to build up the loop with basic joints (_no_ aggregation JointXXX components) and rely on dynamic state selection, i.e., during simulation the states will be dynamically selected in such a way that in no position a degree of freedom is lost.")
+    #k1a > 1e-10 || error("Singular position of loop (either no or two analytic solutions; the mechanism has lost one-degree-of freedom in this position). Try to use another joint-assembly component. In most cases it is best to let joints outside of the JointXXX component be revolute and _not_ prismatic joints. If this also leads to singular positions, it could be that this kinematic loop cannot be solved analytically. In this case you have to build up the loop with basic joints (_no_ aggregation JointXXX components) and rely on dynamic state selection, i.e., during simulation the states will be dynamically selected in such a way that in no position a degree of freedom is lost.")
     k1b = max(k1a, 1.0e-12)
     k2 = sqrt(k1b)
     kcos1 = -A*C + B*k2*ifelse(positive_branch == true, 1, -1)
     ksin1 = -B*C + A*k2*ifelse(positive_branch == true, -1, 1)
-    [atan(ksin1, kcos1)]
+    atan(ksin1, kcos1)
 end
 
-@register_array_symbolic compute_angle2(L::Real, e::AbstractVector, r_a::AbstractVector, r_b::AbstractVector, positive_branch::Bool) begin
-    size = (1, )
-    eltype = Float64
+function compute_angle2(L::Real, e::AbstractVector, r_a::AbstractVector, r_b::AbstractVector, positive_branch::Bool)
+    _compute_angle2(L, e, r_a, r_b, positive_branch)
 end
 
+@register_symbolic compute_angle2(L::Real, e::AbstractVector, r_a::AbstractVector, r_b::AbstractVector, positive_branch::Bool)
+#@register_array_symbolic compute_angle2(L::Real, e::AbstractVector, r_a::AbstractVector, r_b::AbstractVector, positive_branch::Bool) begin
+#    size = (1, )
+#    eltype = Float64
+#end
+
 # @register_symbolic compute_angle(L::Real, e::AbstractVector, r_a::AbstractVector, r_b::AbstractVector, positive_branch::Bool)::Real
 
 # @register_symbolic compute_angle(L::Num, e1::Num, e1::Num, e2::Num, r_a1::Num, r_a2::Num, r_a3::Num, r_b1::Num, r_b2::Num, r_b3::Num, positive_branch)
@@ -769,4 +770,4 @@ end
 #     r_a = SA[r_a1, r_a2, r_a3]
 #     r_b = SA[r_b1, r_b2, r_b3]
 #     compute_angle(L, e, r_a, r_b, positive_branch)
-# end
+# end