# ============================================================================
# DIDACTIC WORKSHEET: DERIVATION OF THE ROTATION MATRIX ABOUT THE X‑AXIS
#
# This worksheet guides you through:
#   1. 2D rotation in the yz‑plane (point on a circle)
#   2. Extension to 3D: rotation of a vector about the x‑axis
#   3. Rotation of a cube about the x‑axis with a dynamic matrix display
#
# Run each section separately to see the theory and the animations.
# ============================================================================

restart;
with(plots):        # for plots and animations
with(plottools):    # for lines, polygons, spheres
with(LinearAlgebra):# for matrix multiplication

# ----------------------------------------------------------------------------
# PART 1: 2D ROTATION IN THE YZ‑PLANE
# ----------------------------------------------------------------------------
print("============================================================================");
print("PART 1: 2D rotation in the yz‑plane");
print("We consider a point P = (y, z) in the yz‑plane.");
print("We want to rotate it counter‑clockwise by an angle θ.");
print("");

print("Step 1 – Polar coordinates:");
print("   y = r·cos(φ),   z = r·sin(φ)");
print("where r is the distance from the origin and φ is the angle with the positive y‑axis.");
print("");

print("Step 2 – Rotation by θ adds θ to the angle:");
print("   y' = r·cos(φ+θ),   z' = r·sin(φ+θ)");
print("");

print("Step 3 – Use the addition formulas:");
print("   cos(φ+θ) = cosφ·cosθ - sinφ·sinθ");
print("   sin(φ+θ) = sinφ·cosθ + cosφ·sinθ");
print("");

print("Step 4 – Multiply by r and substitute y = r·cosφ, z = r·sinφ:");
print("   y' = y·cosθ - z·sinθ");
print("   z' = y·sinθ + z·cosθ");
print("");

print("Step 5 – Write as a matrix multiplication:");
print("   [ y' ]   [ cosθ  -sinθ ] [ y ]");
print("   [ z' ] = [ sinθ   cosθ ] [ z ]");
print("The 2×2 matrix is the 2D rotation matrix R_2D(θ).");
print("============================================================================");
print("");

print("Animation: A point starting at (1,0) rotates through a full circle.");
print("Below the plot you see the numerical matrix R_2D(θ).");
print("");

N := 32;  # number of frames
frames_2d := []:
for i from 0 to N do
    theta := i * (2*Pi/N);
    R2 := Matrix([[cos(theta), -sin(theta)], [sin(theta), cos(theta)]]);
    rotated := convert(R2 . Vector([1,0]), list);
    
    # Graphical elements
    p := pointplot([rotated], symbol=solidcircle, symbolsize=20, color="red");
    l := line([0,0], rotated, color="blue", linestyle=dot);
    circ := plot([cos(t), sin(t), t=0..2*Pi], color="gray", thickness=1);
    y_axis := plot([[-1.5,0],[1.5,0]], color="black", thickness=1);
    z_axis := plot([[0,-1.5],[0,1.5]], color="black", thickness=1);
    y_label := textplot([1.6,0.1, "y"], font=["Times",12]);
    z_label := textplot([0.1,1.6, "z"], font=["Times",12]);
    matrix_text := sprintf("R_2D = [%7.3f  %7.3f]\n        [%7.3f  %7.3f]",
                           cos(theta), -sin(theta), sin(theta), cos(theta));
    mat_plot := textplot([-1.2, -1.3, matrix_text], font=["Courier",10]);
    
    frames_2d := [op(frames_2d),
                  display(p, l, circ, y_axis, z_axis, y_label, z_label, mat_plot,
                          scaling=constrained, view=[-1.5..1.5,-1.5..1.5],
                          title=cat("θ = ", sprintf("%.2f", theta), " rad"))];
end do:
anim_2d := display(frames_2d, insequence=true):
print(anim_2d);
print("Click the play button. The point rotates and the matrix updates.");
print("");

# ----------------------------------------------------------------------------
# PART 2: EXTENSION TO 3D – ROTATION ABOUT THE X‑AXIS
# ----------------------------------------------------------------------------
print("============================================================================");
print("PART 2: From 2D to 3D – rotation about the x‑axis");
print("In 3D we have coordinates (x, y, z).");
print("When rotating about the x‑axis:");
print("  - The x‑coordinate remains unchanged (the axis is the x‑axis).");
print("  - The y and z coordinates rotate exactly as in the yz‑plane.");
print("Therefore:");
print("   x' = x");
print("   y' = y·cosθ - z·sinθ");
print("   z' = y·sinθ + z·cosθ");
print("");

print("Step 6 – Write as a 3×3 matrix:");
print("   [ x' ]   [ 1    0     0 ] [ x ]");
print("   [ y' ] = [ 0  cosθ -sinθ ] [ y ]");
print("   [ z' ]   [ 0  sinθ  cosθ ] [ z ]");
print("This is the rotation matrix R_x(θ).");
print("============================================================================");
print("");

print("Animation: A red vector initially at (0,1,0) rotates about the x‑axis.");
print("The vector stays in the yz‑plane; its x‑coordinate remains 0.");
print("");

Rx := theta -> Matrix([[1,0,0],[0,cos(theta),-sin(theta)],[0,sin(theta),cos(theta)]]):
frames_vec := []:
N := 32;
for i from 0 to N do
    theta := i * (2*Pi/N);
    v0 := Vector([0,1,0]);
    v_rot := Rx(theta) . v0;
    tip := convert(v_rot, list);
    # Draw vector as a line + a small sphere at the tip
    line_vec := line([0,0,0], tip, color="red", thickness=3);
    tip_sphere := sphere(tip, 0.08, color="red");
    # Coordinate axes as lines
    x_axis := line([-1.5,0,0], [3,0,0], color="black", thickness=1);
    y_axis := line([0,-1.5,0], [0,3,0], color="black", thickness=1);
    z_axis := line([0,0,-1.5], [0,0,3], color="black", thickness=1);
    axis_text := textplot3d([[3.2,0,0,"x"],[0,3.2,0,"y"],[0,0,3.2,"z"]], font=["Times",12]);
    frames_vec := [op(frames_vec),
                   display(line_vec, tip_sphere, x_axis, y_axis, z_axis, axis_text,
                           scaling=constrained, view=[-1.5..3,-1.5..3,-1.5..3],
                           orientation=[70,60], title=cat("θ = ", sprintf("%.2f", theta), " rad"))];
end do:
anim_vec := display(frames_vec, insequence=true):
print(anim_vec);
print("Click the play button. The red vector rotates in the yz‑plane.");
print("");

# ----------------------------------------------------------------------------
# PART 3: ROTATING A CUBE ABOUT THE X‑AXIS
# ----------------------------------------------------------------------------
print("============================================================================");
print("PART 3: Rotating a cube about the x‑axis");
print("We apply the matrix R_x(θ) to every vertex of a cube.");
print("The cube has side length 2 and its centre at the origin.");
print("A dynamic caption below the plot shows the current numerical matrix.");
print("============================================================================");
print("");

# Cube definition
cube_vertices := [
    [-1,-1,-1], [ 1,-1,-1], [ 1, 1,-1], [-1, 1,-1],  # bottom face (z = -1)
    [-1,-1, 1], [ 1,-1, 1], [ 1, 1, 1], [-1, 1, 1]   # top face (z =  1)
];
edges := [
    [1,2],[2,3],[3,4],[4,1],  # bottom
    [5,6],[6,7],[7,8],[8,5],  # top
    [1,5],[2,6],[3,7],[4,8]   # vertical
];

draw_cube := proc(vertices)
    local edge, lines, faces;
    lines := seq( line(vertices[edge[1]], vertices[edge[2]], color="Black", thickness=2), edge in edges );
    faces := [
        polygon([vertices[1],vertices[2],vertices[3],vertices[4]], color="Red", transparency=0.5),
        polygon([vertices[5],vertices[6],vertices[7],vertices[8]], color="Green", transparency=0.5),
        polygon([vertices[1],vertices[2],vertices[6],vertices[5]], color="Blue", transparency=0.5),
        polygon([vertices[2],vertices[3],vertices[7],vertices[6]], color="Yellow", transparency=0.5),
        polygon([vertices[3],vertices[4],vertices[8],vertices[7]], color="Cyan", transparency=0.5),
        polygon([vertices[4],vertices[1],vertices[5],vertices[8]], color="Magenta", transparency=0.5)
    ];
    return display(lines, faces);
end:

frames_cube := []:
N := 32;
for i from 0 to N do
    theta := i * (2*Pi/N);
    rotated_vertices := map(v -> convert(Rx(theta) . Vector(v), list), cube_vertices);
    M := Rx(theta);
    matrix_caption := cat(
        "R_x(θ) = \n",
        sprintf("[%7.3f  %7.3f  %7.3f]", M[1,1], M[1,2], M[1,3]), "\n",
        sprintf("[%7.3f  %7.3f  %7.3f]", M[2,1], M[2,2], M[2,3]), "\n",
        sprintf("[%7.3f  %7.3f  %7.3f]", M[3,1], M[3,2], M[3,3])
    );
    frames_cube := [op(frames_cube),
                    display(draw_cube(rotated_vertices),
                            title = cat("θ = ", sprintf("%.2f", theta), " rad"),
                            orientation = [45,45,0],
                            lightmodel = light4,
                            scaling = constrained,
                            axes = normal,
                            labels = ["x","y","z"],
                            view = [-2.5..2.5, -2.5..2.5, -2.5..2.5],
                            caption = matrix_caption
                           )]:
end do:
anim_cube := display(frames_cube, insequence=true):
print(anim_cube);
print("Click the play button. The cube rotates about the x‑axis.");
print("The caption shows the matrix R_x(θ) with current numerical values.");
print("============================================================================");
print("");

print("END OF WORKSHEET.");
print("You have now seen step by step how the rotation matrix about the x‑axis is derived,");
print("from a 2D point rotation, through a 3D vector, to a full cube.");