SphericalSpherical

model SphericalSpherical "Spherical - spherical joint aggregation (1 constraint, no potential states) with an optional point mass in the middle"
    import Modelica.Mechanics.MultiBody.Types;

    extends Interfaces.PartialTwoFrames;

    parameter Boolean animation = true "= true, if animation shall be enabled";
    parameter Boolean showMass = true "= true, if mass shall be shown (provided animation = true and m > 0)";
    parameter Boolean computeRodLength = false "= true, if rodLength shall be computed during initialization (see info)";
    parameter SI.Length rodLength(min = Modelica.Constants.eps, fixed = not computeRodLength, start = 1) "Distance between the origins of frame_a and frame_b (if computeRodLength=true, guess value)";
    parameter SI.Mass m(min = 0) = 0 "Mass of rod (= point mass located in middle of rod)";
    parameter SI.Diameter sphereDiameter = world.defaultJointLength "Diameter of spheres representing the spherical joints"
        annotation (Dialog(
            tab = "Animation",
            group = "if animation = true",
            enable = animation));
    input Types.Color sphereColor = Modelica.Mechanics.MultiBody.Types.Defaults.JointColor "Color of spheres representing the spherical joints"
        annotation (Dialog(
            colorSelector = true,
            tab = "Animation",
            group = "if animation = true",
            enable = animation));
    parameter SI.Diameter rodDiameter = sphereDiameter / Types.Defaults.JointRodDiameterFraction "Diameter of rod connecting the two spherical joint"
        annotation (Dialog(
            tab = "Animation",
            group = "if animation = true",
            enable = animation));
    input Types.Color rodColor = Modelica.Mechanics.MultiBody.Types.Defaults.RodColor "Color of rod connecting the two spherical joints"
        annotation (Dialog(
            colorSelector = true,
            tab = "Animation",
            group = "if animation = true",
            enable = animation));
    parameter SI.Diameter massDiameter = sphereDiameter "Diameter of sphere representing the mass point"
        annotation (Dialog(
            tab = "Animation",
            group = "if animation = true and showMass = true and m > 0",
            enable = animation and showMass and 0 < m));
    input Types.Color massColor = Modelica.Mechanics.MultiBody.Types.Defaults.BodyColor "Color of sphere representing the mass point"
        annotation (Dialog(
            colorSelector = true,
            tab = "Animation",
            group = "if animation = true and showMass = true and m > 0",
            enable = animation and showMass and 0 < m));
    input Types.SpecularCoefficient specularCoefficient = world.defaultSpecularCoefficient "Reflection of ambient light (= 0: light is completely absorbed)"
        annotation (Dialog(
            tab = "Animation",
            group = "if animation = true",
            enable = animation));
    parameter Boolean kinematicConstraint = true "= false, if no constraint shall be defined, due to analytically solving a kinematic loop (\"false\" should not be used by user, but only by MultiBody.Joints.Assemblies joints)"
        annotation (Dialog(tab = "Advanced"));
    Real constraintResidue = rRod_0 * rRod_0 - rodLength * rodLength "Constraint equation of joint in residue form: Either length constraint (= default) or equation to compute rod force (for analytic solution of loops in combination with Internal.RevoluteWithLengthConstraint/PrismaticWithLengthConstraint)"
        annotation (Dialog(
            tab = "Advanced",
            enable = not kinematicConstraint));
    parameter Boolean checkTotalPower = false "= true, if total power flowing into this component shall be determined (must be zero)"
        annotation (Dialog(tab = "Advanced"));
    SI.Force f_rod "Constraint force in direction of the rod (positive on frame_a, when directed from frame_a to frame_b)";
    SI.Position rRod_0[3] "Position vector from frame_a to frame_b resolved in world frame";
    SI.Position rRod_a[3] "Position vector from frame_a to frame_b resolved in frame_a";
    Real eRod_a[3](each final unit = "1") "Unit vector in direction from frame_a to frame_b, resolved in frame_a";
    SI.Position r_CM_0[3] "Dummy if m==0, or position vector from world frame to mid-point of rod, resolved in world frame";
    SI.Velocity v_CM_0[3] "First derivative of r_CM_0";
    SI.Force f_CM_a[3] "Dummy if m==0, or inertial force acting at mid-point of rod due to mass point acceleration, resolved in frame_a";
    SI.Force f_CM_e[3] "Dummy if m==0, or projection of f_CM_a onto eRod_a, resolved in frame_a";
    SI.Force f_b_a1[3] "Force acting at frame_b, but without force in rod, resolved in frame_a";
    SI.Power totalPower "Total power flowing into this element, if checkTotalPower=true (otherwise dummy)";
protected
    Visualizers.Advanced.Shape shape_rod(shapeType = "cylinder", color = rodColor, specularCoefficient = specularCoefficient, length = rodLength, width = rodDiameter, height = rodDiameter, lengthDirection = eRod_a, widthDirection = {0, 1, 0}, r = frame_a.r_0, R = frame_a.R) if world.enableAnimation and animation;
    Visualizers.Advanced.Shape shape_a(shapeType = "sphere", color = sphereColor, specularCoefficient = specularCoefficient, length = sphereDiameter, width = sphereDiameter, height = sphereDiameter, lengthDirection = eRod_a, widthDirection = {0, 1, 0}, r_shape = -eRod_a * (0.5 * sphereDiameter), r = frame_a.r_0, R = frame_a.R) if world.enableAnimation and animation;
    Visualizers.Advanced.Shape shape_b(shapeType = "sphere", color = sphereColor, specularCoefficient = specularCoefficient, length = sphereDiameter, width = sphereDiameter, height = sphereDiameter, lengthDirection = eRod_a, widthDirection = {0, 1, 0}, r_shape = eRod_a * (rodLength - 0.5 * sphereDiameter), r = frame_a.r_0, R = frame_a.R) if world.enableAnimation and animation;
    Visualizers.Advanced.Shape shape_mass(shapeType = "sphere", color = massColor, specularCoefficient = specularCoefficient, length = massDiameter, width = massDiameter, height = massDiameter, lengthDirection = eRod_a, widthDirection = {0, 1, 0}, r_shape = eRod_a * (0.5 * rodLength - 0.5 * sphereDiameter), r = frame_a.r_0, R = frame_a.R) if world.enableAnimation and animation and showMass and 0 < m;
equation
    if checkTotalPower then 
        totalPower = frame_a.f * Frames.resolve2(frame_a.R, der(frame_a.r_0)) + frame_b.f * Frames.resolve2(frame_b.R, der(frame_b.r_0)) + (-m) * (der(v_CM_0) - world.gravityAcceleration(r_CM_0)) * v_CM_0 + frame_a.t * Frames.angularVelocity2(frame_a.R) + frame_b.t * Frames.angularVelocity2(frame_b.R);
    else 
        totalPower = 0;
    end if;
    if kinematicConstraint then 
        rRod_0 = transpose(frame_b.R.T) * (frame_b.R.T * frame_b.r_0) - transpose(frame_a.R.T) * (frame_a.R.T * frame_a.r_0);
    else 
        rRod_0 = frame_b.r_0 - frame_a.r_0;
    end if;
    if 0 < m then 
        r_CM_0 = frame_a.r_0 + 0.5 * rRod_0;
        v_CM_0 = der(r_CM_0);
        f_CM_a = m * Frames.resolve2(frame_a.R, der(v_CM_0) - world.gravityAcceleration(r_CM_0));
        f_CM_e = f_CM_a * eRod_a * eRod_a;
        frame_a.f = 0.5 * (f_CM_a - f_CM_e) + f_rod * eRod_a;
        f_b_a1 = 0.5 * (f_CM_a + f_CM_e);
        frame_b.f = Frames.resolveRelative(f_b_a1 - f_rod * eRod_a, frame_a.R, frame_b.R);
    else 
        r_CM_0 = zeros(3);
        v_CM_0 = zeros(3);
        f_CM_a = zeros(3);
        f_CM_e = zeros(3);
        f_b_a1 = zeros(3);
        frame_a.f = f_rod * eRod_a;
        frame_b.f = -Frames.resolveRelative(frame_a.f, frame_a.R, frame_b.R);
    end if;
    eRod_a = rRod_a / rodLength;
    rRod_a = Frames.resolve2(frame_a.R, rRod_0);
    constraintResidue = 0;
    frame_a.t = zeros(3);
    frame_b.t = zeros(3);

    annotation (
        Icon(
            coordinateSystem(
                preserveAspectRatio = true,
                extent = {
                    {-100, -100}, 
                    {100, 100}}),
            graphics = {
                Ellipse(
                    extent = {
                        {-95, -40}, 
                        {-15, 40}},
                    fillPattern = FillPattern.Sphere,
                    fillColor = {192, 192, 192}), 
                Ellipse(
                    extent = {
                        {-84, -30}, 
                        {-24, 30}},
                    fillColor = {255, 255, 255},
                    fillPattern = FillPattern.Solid), 
                Ellipse(
                    extent = {
                        {15, -40}, 
                        {95, 40}},
                    fillPattern = FillPattern.Sphere,
                    fillColor = {192, 192, 192}), 
                Ellipse(
                    extent = {
                        {25, -29}, 
                        {85, 30}},
                    lineColor = {128, 128, 128},
                    fillColor = {255, 255, 255},
                    fillPattern = FillPattern.Solid), 
                Text(
                    extent = {
                        {-150, 90}, 
                        {150, 50}},
                    textString = "%name",
                    lineColor = {0, 0, 255}), 
                Rectangle(
                    extent = {
                        {-40, 40}, 
                        {41, -41}},
                    lineColor = {255, 255, 255},
                    fillColor = {255, 255, 255},
                    fillPattern = FillPattern.Solid), 
                Rectangle(
                    extent = {
                        {-51, 6}, 
                        {48, -4}},
                    fillPattern = FillPattern.HorizontalCylinder,
                    fillColor = {192, 192, 192}), 
                Ellipse(
                    extent = {
                        {-68, 15}, 
                        {-39, -13}},
                    fillPattern = FillPattern.Sphere,
                    fillColor = {192, 192, 192}), 
                Ellipse(
                    extent = {
                        {39, 14}, 
                        {68, -14}},
                    fillPattern = FillPattern.Sphere,
                    fillColor = {192, 192, 192}), 
                Text(
                    extent = {
                        {-150, -60}, 
                        {150, -90}},
                    textString = "%rodLength")}),
        Diagram(
            coordinateSystem(
                preserveAspectRatio = true,
                extent = {
                    {-100, -100}, 
                    {100, 100}}),
            graphics = {
                Ellipse(
                    extent = {
                        {-98, -40}, 
                        {-18, 40}},
                    fillColor = {160, 160, 164},
                    fillPattern = FillPattern.Solid), 
                Ellipse(
                    extent = {
                        {-88, -30}, 
                        {-28, 30}},
                    lineColor = {160, 160, 164},
                    fillColor = {255, 255, 255},
                    fillPattern = FillPattern.Solid), 
                Ellipse(
                    extent = {
                        {18, -40}, 
                        {98, 40}},
                    lineColor = {160, 160, 164},
                    fillColor = {160, 160, 164},
                    fillPattern = FillPattern.Solid), 
                Ellipse(
                    extent = {
                        {29, -30}, 
                        {89, 29}},
                    lineColor = {192, 192, 192},
                    fillColor = {255, 255, 255},
                    fillPattern = FillPattern.Solid), 
                Line(
                    points = {
                        {-56, -60}, 
                        {46, -60}},
                    color = {0, 0, 255}), 
                Polygon(
                    points = {
                        {56, -60}, 
                        {41, -54}, 
                        {41, -66}, 
                        {56, -60}},
                    lineColor = {0, 0, 255},
                    fillColor = {0, 0, 255},
                    fillPattern = FillPattern.Solid), 
                Text(
                    extent = {
                        {-37, -63}, 
                        {33, -79}},
                    textString = "rodLength",
                    lineColor = {0, 0, 255}), 
                Rectangle(
                    extent = {
                        {-40, 41}, 
                        {40, -40}},
                    lineColor = {255, 255, 255},
                    fillColor = {255, 255, 255},
                    fillPattern = FillPattern.Solid), 
                Rectangle(
                    extent = {
                        {-51, 6}, 
                        {48, -4}},
                    fillPattern = FillPattern.HorizontalCylinder,
                    fillColor = {192, 192, 192}), 
                Ellipse(
                    extent = {
                        {-71, 15}, 
                        {-42, -13}},
                    fillPattern = FillPattern.Solid), 
                Ellipse(
                    extent = {
                        {42, 14}, 
                        {71, -14}},
                    fillPattern = FillPattern.Solid), 
                Line(
                    points = {
                        {-56, -71}, 
                        {-56, 1}},
                    color = {0, 0, 255}), 
                Line(
                    points = {
                        {56, -72}, 
                        {56, 0}},
                    color = {0, 0, 255}), 
                Polygon(
                    points = {
                        {11, 1}, 
                        {-1, 4}, 
                        {-1, -2}, 
                        {11, 1}},
                    lineColor = {0, 0, 255}), 
                Line(
                    points = {
                        {-56, 1}, 
                        {-1, 1}},
                    color = {0, 0, 255}), 
                Text(
                    extent = {
                        {-32, -4}, 
                        {4, -29}},
                    textString = "eRod_a")}),
        Documentation(info = "<html>\n<p>\nJoint that has a spherical joint on each of its two ends.\nThe rod connecting the two spherical joints is approximated by a\npoint mass that is located in the middle of the rod. When the mass\nis set to zero (default), special code for a massless body is generated.\nIn the following default animation figure, the two spherical joints are\nrepresented by two red spheres, the connecting rod by a grey cylinder\nand the point mass in the middle of the rod by a light blue sphere:\n</p>\n\n<p>\n<img src=\"modelica://Modelica/Resources/Images/Mechanics/MultiBody/Joints/SphericalSpherical.png\" alt=\"model Joints.SphericalSpherical\">\n</p>\n\n<p>\nThis joint introduces <strong>one constraint</strong> defining that the distance between\nthe origin of frame_a and the origin of frame_b is constant (= rodLength).\nIt is highly recommended to use this joint in loops\nwhenever possible, because this enhances the efficiency\nconsiderably due to smaller systems of non-linear algebraic\nequations.\n</p>\n<p>\nIt is sometimes desirable to <strong>compute</strong> the <strong>rodLength</strong>\nof the connecting rod during initialization. For this, parameter\n<strong>computeLength</strong> has to be set to <strong>true</strong> and instead <strong>one</strong> other,\neasier to determine, position variable in the same loop\nneeds to have a fixed attribute of <strong>true</strong>. For example,\nif a loop consists of one Revolute joint, one Prismatic joint and\na SphericalSpherical joint, one may fix the start values of the revolute\njoint angle and of the relative distance of the prismatic joint\nin order to compute the rodLength of the rod.\n</p>\n<p>\nIt is not possible to connect other components, such as a body with mass\nproperties or a special visual shape object to the rod connecting\nthe two spherical joints. If this is needed, use instead joint Joints.<strong>UniversalSpherical</strong>\nthat has this property.\n</p>\n</html>"));
end SphericalSpherical;