diff --git a/.github/copilot-instructions.md b/.github/copilot-instructions.md
index c0f8e4b..ad055d9 100644
--- a/.github/copilot-instructions.md
+++ b/.github/copilot-instructions.md
@@ -19,18 +19,70 @@ This is a **Julia linear algebra computing project** using DrWatson for reproduc
 # Main module uses @reexport for clean interface
 using Reexport
 @reexport using GeometryBasics, Plots, LinearAlgebra, RationalRoots, Symbolics
+
+# Configure plotting for both interactive and headless environments
+if haskey(ENV, "CI") || get(ENV, "GKSwstype", "") == "100"
+    # CI or headless environment - use headless mode
+    ENV["GKSwstype"] = "100"
+    gr(show=false)
+else
+    # Interactive environment - normal plotting
+    gr()
+end
+
 # Comprehensive exports for all functions
+# Pure computational functions (no plotting dependencies)
+export calculate_param_line
+# Integrated plotting functions (computation + visualization)
 export distance_2_points, center_of_gravity, barycentric_coord, plot_param_line
 # ... matrix functions, line functions, etc.
 ```
 
+### Environment-Aware Module Loading
+```julia
+# The module automatically detects CI vs interactive environments
+if haskey(ENV, "CI") || get(ENV, "GKSwstype", "") == "100"
+    ENV["GKSwstype"] = "100"  # Headless plotting
+    gr(show=false)
+end
+```
+
+## Linear_Algebra CI Testing Approach
+
+The superior CI testing strategy (aligned with Math_Foundations) consists of three components:
+
+### 1. Module-Level Headless Detection
+Configure plotting environment in the main module (`Linear_Algebra.jl`) at load time:
+```julia
+# Automatic CI detection and headless configuration
+if haskey(ENV, "CI") || get(ENV, "GKSwstype", "") == "100"
+    ENV["GKSwstype"] = "100"  # Force headless mode
+    gr(show=false)           # Disable plot display
+end
+```
+
+### 2. Manual GKS Configuration in Tests
+Set `ENV["GKSwstype"] = "100"` in test files before loading the module:
+```julia
+# In test files - Configure headless mode before loading module
+ENV["GKSwstype"] = "100"  # Force headless plotting for CI
+using DrWatson, Test
+@quickactivate "Linear_Algebra"
+using Linear_Algebra
+```
+
+### 3. Separated Computational/Plotting Logic with Robust Testing
+- **Pure computational functions** (`calculate_*`): Test mathematical logic directly, no try-catch
+- **Plotting functions** (`plot_*`): Test with try-catch fallback for CI compatibility
+- **Integration testing**: Verify both computation and visualization work together
+
 ### Test Setup (Uses @quickactivate)
 ```julia
 # Tests use DrWatson @quickactivate pattern
 using DrWatson, Test
 @quickactivate "Linear_Algebra"
+# Load the Linear_Algebra package
 using Linear_Algebra
-using GeometryBasics, LinearAlgebra
 ```
 
 ### CI-Compatible Plotting Pattern
@@ -77,18 +129,47 @@ end
 17. Maintain consistency with mathematical notation
 
 ### Testing Patterns
-18. **Comprehensive Coverage**: 49 tests covering ~100% of functions
-19. **CI-Safe**: Plotting tests skip in headless environments
+18. **Comprehensive Coverage**: Test coverage includes all mathematical functions
+19. **CI-Safe**: Plotting tests work in both local and headless environments
 20. **Edge Cases**: Test mathematical edge cases (orthogonal vectors, zero angles, etc.)
 21. **Type Testing**: Verify return types (Point2f, AbstractVector, matrices)
 22. **Numerical Precision**: Use `atol=1e-10` for floating-point comparisons
-23. **Error Handling**: Use try-catch for functions that might fail in CI
+23. **CI-Compatible Testing Pattern**: Separate computational logic from plotting, test math directly without try-catch, only use try-catch for visualization:
+```julia
+# Test computational logic directly (NO try-catch - mathematical errors should fail)
+@testset "Pure Computational Tests" begin
+    points = calculate_param_line(p, q, 3)
+    @test length(points) == 3
+    @test typeof(points) == Vector{Point2f}
+    # Test mathematical correctness without plotting dependencies
+end
+
+# Test integration (plotting + computation) with CI-safe fallback
+@testset "Integration Tests" begin
+    try
+        # Test the plotting function (includes computation + visualization)
+        result = plot_param_line(p, q, 3)
+        @test typeof(result) == Vector{Point2f}
+        @test length(result) == 3
+    catch e
+        # Only catch plotting-related errors, not computational errors
+        if contains(string(e), "display") || contains(string(e), "GKS") || isa(e, ArgumentError)
+            @test hasmethod(plot_param_line, (Point2f, Point2f, Int64))
+        else
+            # Re-throw computational errors - these should fail the test
+            rethrow(e)
+        end
+    end
+end
+```
 
 ### Code Organization
 24. **Two-File Structure**: Basic operations in `linear_algebra_basic.jl`, matrices in `linear_algebra_transform.jl`
 25. **Consistent Naming**: Functions end with descriptive suffixes (`_matrix`, `_line`, `_coord`)
 26. **Symbolic Variants**: Provide `_symbolic` versions for algebraic manipulation
 27. **Export Everything**: All public functions exported from main module
+28. **Testing Structure**: Modular test files (`test_linear_algebra_basic.jl`, `test_linear_algebra_transform.jl`)
+29. **Follow Math_Foundations Pattern**: Separate computational logic from plotting, use three-tier testing approach
 
 ## Dependencies & Performance
 
@@ -121,10 +202,18 @@ CI=true julia --project=. test/runtests.jl
 julia --project=. docs/make.jl
 ```
 
+### Julia Compilation Considerations
+- **Be Patient with First Runs**: Julia often needs to precompile packages and rebuild project cache on first run. when running a Julia command in the CLI for the first time, it may take a while to precompile the packages and build the project cache, so you won't see the results of running the command for a while.
+- **Typical First Run**: May take 15-30 seconds for precompilation before tests actually start
+- **Example Expected Output**: `Precompiling DrWatson... 3 dependencies successfully precompiled in 17 seconds`
+- **Subsequent Runs**: Much faster once cache is built
+- **Don't Cancel Early**: Allow time for compilation phase to complete
+- **IMPORTANT**: This applies to ALL Julia commands including CI testing with `CI=true julia --project=. test/runtests.jl`
+
 ### CI Considerations
 - Tests automatically detect CI environment via ENV variables
-- Plotting tests skip gracefully in headless mode
-- 49 tests in local mode, 47 tests in CI mode (plotting tests reduced)
+- Plotting tests skip gracefully in headless mode  
+- 68 tests pass in both local and CI modes (plotting tests with fallbacks)
 - Test execution time: ~15-16 seconds
 
 ## Git Best Practices
@@ -202,6 +291,10 @@ distance_2_points(p::Point, q::Point) -> Float64
 center_of_gravity(p::Point, q::Point, t) -> Point
 vector_angle_cos(p::Vector, q::Vector) -> Float64
 orthproj(v::Vector, w::Vector) -> Vector
+# Pure computational functions (no plotting dependencies)
+calculate_param_line(p::Point, q::Point, n::Int64) -> Vector{Point2f}
+# Integrated plotting functions (computation + visualization)
+plot_param_line(p::Point, q::Point, n::Int64) -> Vector{Point2f}
 ```
 
 ### Matrix Transformations
diff --git a/src/Linear_Algebra.jl b/src/Linear_Algebra.jl
index 9b8395c..943b584 100644
--- a/src/Linear_Algebra.jl
+++ b/src/Linear_Algebra.jl
@@ -1,10 +1,21 @@
 module Linear_Algebra
 using Reexport
 @reexport using  GeometryBasics, Plots, LinearAlgebra, RationalRoots, Symbolics
-# Set GR as the default plotting backend
-gr()
+
+# Configure plotting for both interactive and headless environments
+if haskey(ENV, "CI") || get(ENV, "GKSwstype", "") == "100"
+    # CI or headless environment - use headless mode
+    ENV["GKSwstype"] = "100"
+    gr(show=false)
+else
+    # Interactive environment - normal plotting
+    gr()
+end
 
 # Exports...
+# Pure computational functions (no plotting dependencies)
+export calculate_param_line
+# Integrated plotting functions (computation + visualization)
 export distance_2_points, center_of_gravity, barycentric_coord, plot_param_line
 export vector_angle_cos, is_orthogonal, polar_unit, orthproj, reflection, rotation
 export point_in_implicit_line, parametric_to_implicit_line, implicit_to_parametric_line, explicit_line
diff --git a/src/linear_algebra_basic.jl b/src/linear_algebra_basic.jl
index 5624261..78b7f01 100644
--- a/src/linear_algebra_basic.jl
+++ b/src/linear_algebra_basic.jl
@@ -14,6 +14,8 @@ Creates center of gravity of points `p` and `q`using parametric equation of a li
 function center_of_gravity(p::Point, q::Point, t)
     v = q - p
     r = p + (t * v)
+    # Ensure we return a Point2f type
+    return Point2f(r)
 end
 
 """
@@ -25,18 +27,33 @@ function barycentric_coord(p::Point, q::Point, r::Point)
 end
 
 """
-    plot_param_line(p::Point, q::Point, n::Int64) → [Point]
-Creates `n` points on a line defined by `p` and `q`, using the parametric equation of a line, then plot
+    calculate_param_line(p::Point, q::Point, n::Int64) → [Point]
+Creates `n` points on a line defined by `p` and `q`, using the parametric equation of a line.
+Pure computational function - no plotting dependencies.
 """
-
-function plot_param_line(p::Point, q::Point, n::Int64)
+function calculate_param_line(p::Point, q::Point, n::Int64)
     Ps = Point[]
     for i = 1:n
         t = i / n
         r = center_of_gravity(p, q, t)
-        s = scatter!(r, legend=false)
         push!(Ps, r)
     end
+    Ps
+end
+
+"""
+    plot_param_line(p::Point, q::Point, n::Int64) → [Point]
+Creates `n` points on a line defined by `p` and `q`, using the parametric equation of a line, then plot
+"""
+
+function plot_param_line(p::Point, q::Point, n::Int64)
+    # Use computational function for calculations
+    Ps = calculate_param_line(p, q, n)
+    
+    # Add plotting functionality
+    for point in Ps
+        s = scatter!(point, legend=false)
+    end
     s = plot!(Ps, legend=false)
     display(s)
     Ps
diff --git a/test/runtests.jl b/test/runtests.jl
index 85ced14..22a2caa 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,278 +1,16 @@
 using DrWatson, Test
 @quickactivate "Linear_Algebra"
+
+# Load the Linear_Algebra package
 using Linear_Algebra
-using GeometryBasics, LinearAlgebra
 
 # Run test suite
 println("Starting tests")
 ti = time()
 
 @testset "Linear_Algebra tests" begin
-    
-    @testset "Basic Linear Algebra Functions" begin
-        # Test distance_2_points
-        p1 = Point2f(0.0, 0.0)
-        p2 = Point2f(3.0, 4.0)
-        @test distance_2_points(p1, p2) ≈ 5.0
-        
-        # Test center_of_gravity
-        p = Point2f(0.0, 0.0)
-        q = Point2f(2.0, 2.0)
-        cog = center_of_gravity(p, q, 0.5)
-        @test cog ≈ Point2f(1.0, 1.0)
-        
-        # Test barycentric_coord
-        p = Point2f(0.0, 0.0)
-        q = Point2f(4.0, 0.0)
-        r = Point2f(1.0, 0.0)
-        t = barycentric_coord(p, q, r)
-        @test t ≈ 1.0/3.0
-        
-        # Test vector_angle_cos
-        v1 = [1.0, 0.0]
-        v2 = [0.0, 1.0]
-        @test vector_angle_cos(v1, v2) ≈ 0.0  # cos(90°) = 0
-        
-        v3 = [1.0, 0.0]
-        v4 = [1.0, 0.0]
-        @test vector_angle_cos(v3, v4) ≈ 1.0  # cos(0°) = 1
-        
-        # Test is_orthogonal
-        p1 = Point2f(1.0, 0.0)
-        p2 = Point2f(0.0, 1.0)
-        @test is_orthogonal(p1, p2) == true
-        
-        p3 = Point2f(1.0, 1.0)
-        p4 = Point2f(1.0, 0.0)
-        @test is_orthogonal(p3, p4) == false
-        
-        # Test polar_unit
-        v = [3.0, 4.0]
-        unit_v = polar_unit(v)
-        @test norm(unit_v) ≈ 1.0
-        @test unit_v ≈ [0.6, 0.8]
-        
-        # Test orthproj
-        v = [1.0, 0.0]  # project onto unit vector
-        w = [3.0, 4.0]  # vector to project
-        proj = orthproj(v, w)
-        # orthproj(v,w) = (dot(w,v)/norm(v)^2) * v = (3.0/1.0) * [1,0] = [3,0]
-        @test proj ≈ [3.0, 0.0]
-        
-        # Test reflection (vector function)
-        v = [1.0, 0.0]  # vector to reflect
-        w = [0.0, 1.0]  # reflection axis
-        refl = reflection(v, w)
-        @test length(refl) == 2  # should return a 2D vector
-        
-        # Test rotation (vector function)
-        θ = 90  # 90 degrees (the function expects degrees)
-        v = [1.0, 0.0]
-        rotated = rotation(θ, v)
-        @test rotated ≈ [0.0, 1.0] atol=1e-10
-        
-        # Test plot_param_line (skip in CI environments due to plotting)
-        p = Point2f(0.0, 0.0)
-        q = Point2f(1.0, 1.0)
-        
-        # Skip plotting tests in CI environments
-        if get(ENV, "CI", "false") == "true" || get(ENV, "GITHUB_ACTIONS", "false") == "true"
-            # In CI, just test that the function exists
-            @test hasmethod(plot_param_line, (typeof(p), typeof(q), Int64))
-        else
-            # Local testing - allow plotting but capture any display issues
-            try
-                points = plot_param_line(p, q, 3)
-                @test length(points) == 3
-                @test isa(points[1], Array) || isa(points[1], AbstractVector)
-                @test isa(points[end], Array) || isa(points[end], AbstractVector)
-            catch e
-                # If plotting fails for any reason, just test function existence
-                @test hasmethod(plot_param_line, (typeof(p), typeof(q), Int64))
-                println("Note: Plotting test skipped due to: $e")
-            end
-        end
-    end
-    
-    @testset "Transformation Matrix Functions" begin
-        # Test projection_matrix
-        x = [1.0, 0.0]
-        P = projection_matrix(x)
-        @test size(P) == (2, 2)
-        @test P * x ≈ x  # projecting onto itself should return itself
-        
-        # Test rotation_matrix  
-        d = 90  # 90 degrees (the function takes degrees, not radians)
-        R = rotation_matrix(d)
-        @test size(R) == (2, 2)
-        # Rotating [1,0] by 90° should give [0,1]
-        v = [1.0, 0.0]
-        rotated = R * v
-        @test rotated ≈ [0.0, 1.0] atol=1e-10
-        
-        # Test stretch_matrix
-        n = 2.0
-        S = stretch_matrix(n)
-        @test size(S) == (2, 2)
-        v = [1.0, 1.0]
-        stretched = S * v
-        @test stretched ≈ [2.0, 2.0]  # both components scaled by n
-        
-        # Test reflection_matrix
-        U = [1.0, 0.0]  # reflect across x-axis
-        R = reflection_matrix(U)
-        @test size(R) == (2, 2)
-        v = [1.0, 1.0]
-        reflected = R * v
-        @test reflected[1] ≈ 1.0   # x-component unchanged
-        @test reflected[2] ≈ -1.0  # y-component flipped
-        
-        # Test additional matrix functions
-        # Test projection_matrix_polar
-        n = 1.0
-        P_polar = projection_matrix_polar(n)
-        @test size(P_polar) == (2, 2)
-        
-        # Test projection_matrix_transpose
-        u = [1.0, 0.0]
-        P_transpose = projection_matrix_transpose(u)
-        @test size(P_transpose) == (2, 2)
-        
-        # Test rotation_matrix_ns (non-symbolic)
-        θ = π/4  # 45 degrees
-        R_ns = rotation_matrix_ns(θ)
-        @test size(R_ns) == (2, 2)
-        
-        # Test reflection_matrix_rational
-        U_rational = [1.0, 1.0]
-        R_rational = reflection_matrix_rational(U_rational)
-        @test size(R_rational) == (2, 2)
-        
-        # Test symbolic matrix functions (should return symbolic expressions)
-        P_sym = projection_matrix_symbolic()
-        @test size(P_sym) == (2, 2)
-        
-        P_sym_polar = projection_matrix_symbolic_polar()
-        @test size(P_sym_polar) == (2, 2)
-        
-        R_sym = rotation_matrix_symbolic()
-        @test size(R_sym) == (2, 2)
-        
-        S_sym = stretch_matrix_symbolic()
-        @test size(S_sym) == (2, 2)
-        
-        Refl_sym = reflection_matrix_symbolic()
-        @test size(Refl_sym) == (2, 2)
-    end
-    
-    @testset "Line Geometry Functions" begin
-        # Test parametric_to_implicit_line
-        p = Point2f(0.0, 0.0)
-        v = [1.0, 1.0]  # direction vector for y = x line
-        a, b, c = parametric_to_implicit_line(p, v)
-        # For line y = x, implicit form should be x - y = 0
-        @test abs(a) ≈ abs(b)  # coefficients should have same magnitude
-        
-        # Test distance_to_implicit_line
-        # Distance from point (1,0) to line x - y = 0 (y = x)
-        a, b, c = 1.0, -1.0, 0.0
-        r = Point2f(1.0, 0.0)
-        dist = distance_to_implicit_line(a, b, c, r)
-        expected_dist = 1.0 / sqrt(2.0)  # distance formula for point to line
-        @test dist ≈ expected_dist
-        
-        # Test explicit_line
-        p1 = Point2f(0.0, 0.0)
-        p2 = Point2f(1.0, 2.0)
-        m, b = explicit_line(p1, p2)
-        @test m ≈ 2.0  # slope
-        @test b ≈ 0.0  # y-intercept
-        
-        # Test point_in_implicit_line
-        p = Point2f(0.0, 0.0)
-        q = Point2f(1.0, 1.0)
-        x = Point2f(0.5, 0.5)  # point on line between p and q
-        result = point_in_implicit_line(p, q, x)
-        @test isa(result, Bool) || isa(result, Number)  # function returns boolean or number
-        
-        # Test implicit_to_parametric_line
-        a, b, c = 1.0, -1.0, 0.0  # line x - y = 0
-        p_param, v_param = implicit_to_parametric_line(a, b, c)
-        @test isa(p_param, AbstractVector) || isa(p_param, Point2f)
-        @test length(v_param) == 2
-        
-        # Test implicit_line_point_normal_form
-        a, b, c = 3.0, 4.0, 5.0
-        â, b̂, ĉ = implicit_line_point_normal_form(a, b, c)
-        # Should normalize the coefficients
-        norm_factor = sqrt(a^2 + b^2)
-        @test â ≈ a/norm_factor
-        @test b̂ ≈ b/norm_factor
-        
-        # Test distance_to_parametric_line
-        p = Point2f(0.0, 0.0)  # point on line
-        v = [1.0, 0.0]         # direction vector
-        r = Point2f(0.0, 1.0)  # test point
-        dist_param = distance_to_parametric_line(p, v, r)
-        @test dist_param ≈ 1.0  # distance from (0,1) to x-axis
-        
-        # Test foot_of_line functions
-        P = Point2f(0.0, 0.0)
-        v = [1.0, 0.0]
-        R = Point2f(1.0, 1.0)
-        
-        # Test foot_of_line with tuple return
-        foot, dist = foot_of_line(P, v, R)
-        @test isa(foot, AbstractVector) || isa(foot, Point2f)
-        @test isa(dist, Number)
-        
-        # Test foot_of_line with boolean parameter
-        result = foot_of_line(P, v, R, false)
-        @test isa(result, Point2f) || isa(result, Tuple)
-        
-        # Test foot_of_line with two points (skip due to undefined variable)
-        # A = Point2f(0.0, 0.0)
-        # B = Point2f(1.0, 0.0)
-        # foot_2pt = foot_of_line(A, B, false)
-        # @test isa(foot_2pt, Point2f) || isa(foot_2pt, Tuple)
-    end
-    
-    @testset "Advanced Line Functions" begin
-        # Test intersection functions (just check they execute without error)
-        v1 = [1.0, 0.0]  # horizontal line direction
-        w1 = [0.0, 1.0]  # vertical line direction
-        p1 = Point2f(0.0, 1.0)  # point on first line
-        q1 = Point2f(1.0, 0.0)  # point on second line
-        try
-            intersection = intersection_2_parametric_lines(v1, w1, p1, q1)
-            @test true  # if it doesn't error, that's good
-        catch
-            @test false  # if it errors, that's bad
-        end
-        
-        # Test intersection_2_implicit_lines
-        # Line 1: x + y = 1
-        # Line 2: x - y = 0
-        a₁, b₁, c₁ = 1.0, 1.0, -1.0
-        a₂, b₂, c₂ = 1.0, -1.0, 0.0
-        try
-            intersection_impl = intersection_2_implicit_lines(a₁, b₁, c₁, a₂, b₂, c₂)
-            @test true  # if it doesn't error, that's good
-        catch
-            @test false  # if it errors, that's bad
-        end
-        
-        # Test Linear_Algebra.rationalize function
-        x = π
-        try
-            rational_approx = Linear_Algebra.rationalize(x; sigdigits=4)
-            @test isa(rational_approx, Rational) || isa(rational_approx, Number)
-        catch
-            # If the function doesn't exist or has different signature, just pass
-            @test true
-        end
-    end
-    
+    include("test_linear_algebra_basic.jl")
+    include("test_linear_algebra_transform.jl")
 end
 
 ti = time() - ti
diff --git a/test/test_linear_algebra_basic.jl b/test/test_linear_algebra_basic.jl
new file mode 100644
index 0000000..98468d8
--- /dev/null
+++ b/test/test_linear_algebra_basic.jl
@@ -0,0 +1,235 @@
+using Test
+
+# Configure plotting for headless CI environments BEFORE loading Linear_Algebra
+ENV["GKSwstype"] = "100"  # Set GKS to use headless mode
+
+using Linear_Algebra
+using GeometryBasics, LinearAlgebra
+
+# Ensure plots directory exists for plotting tests
+if !isdir("plots")
+    mkdir("plots")
+end
+
+@testset "linear_algebra_basic.jl Tests" begin
+    
+    @testset "Basic Point and Vector Operations" begin
+        # Test distance_2_points
+        p1 = Point2f(0.0, 0.0)
+        p2 = Point2f(3.0, 4.0)
+        @test distance_2_points(p1, p2) ≈ 5.0
+        
+        # Test center_of_gravity
+        p = Point2f(0.0, 0.0)
+        q = Point2f(2.0, 2.0)
+        cog = center_of_gravity(p, q, 0.5)
+        @test cog ≈ Point2f(1.0, 1.0)
+        
+        # Test barycentric_coord
+        p = Point2f(0.0, 0.0)
+        q = Point2f(4.0, 0.0)
+        r = Point2f(1.0, 0.0)
+        t = barycentric_coord(p, q, r)
+        @test t ≈ 1.0/3.0
+        
+        # Test vector_angle_cos
+        v1 = [1.0, 0.0]
+        v2 = [0.0, 1.0]
+        @test vector_angle_cos(v1, v2) ≈ 0.0  # cos(90°) = 0
+        
+        v3 = [1.0, 0.0]
+        v4 = [1.0, 0.0]
+        @test vector_angle_cos(v3, v4) ≈ 1.0  # cos(0°) = 1
+        
+        # Test is_orthogonal
+        p1 = Point2f(1.0, 0.0)
+        p2 = Point2f(0.0, 1.0)
+        @test is_orthogonal(p1, p2) == true
+        
+        p3 = Point2f(1.0, 1.0)
+        p4 = Point2f(1.0, 0.0)
+        @test is_orthogonal(p3, p4) == false
+        
+        # Test polar_unit
+        v = [3.0, 4.0]
+        unit_v = polar_unit(v)
+        @test norm(unit_v) ≈ 1.0
+        @test unit_v ≈ [0.6, 0.8]
+        
+        # Test orthproj
+        v = [1.0, 0.0]  # project onto unit vector
+        w = [3.0, 4.0]  # vector to project
+        proj = orthproj(v, w)
+        # orthproj(v,w) = (dot(w,v)/norm(v)^2) * v = (3.0/1.0) * [1,0] = [3,0]
+        @test proj ≈ [3.0, 0.0]
+        
+        # Test reflection (vector function)
+        v = [1.0, 0.0]  # vector to reflect
+        w = [0.0, 1.0]  # reflection axis
+        refl = reflection(v, w)
+        @test length(refl) == 2  # should return a 2D vector
+        
+        # Test rotation (vector function)
+        θ = 90  # 90 degrees (the function expects degrees)
+        v = [1.0, 0.0]
+        rotated = rotation(θ, v)
+        @test rotated ≈ [0.0, 1.0] atol=1e-10
+    end
+    
+    # COMPUTATIONAL TESTS - These should NEVER use try-catch, ensuring math is always tested
+    @testset "Pure Computational Functions (Mathematics Only)" begin
+        @testset "calculate_param_line" begin
+            # Test parametric line generation
+            p = Point2f(0.0, 0.0)
+            q = Point2f(2.0, 2.0)
+            
+            # Test 3 points on line
+            points = calculate_param_line(p, q, 3)
+            @test length(points) == 3
+            @test points[1] ≈ Point2f(2.0/3.0, 2.0/3.0)  # t=1/3
+            @test points[2] ≈ Point2f(4.0/3.0, 4.0/3.0)  # t=2/3
+            @test points[3] ≈ Point2f(2.0, 2.0)          # t=1 (endpoint)
+            
+            # Test single point
+            single_point = calculate_param_line(p, q, 1)
+            @test length(single_point) == 1
+            @test single_point[1] ≈ Point2f(2.0, 2.0)  # t=1
+            
+            # Test return type
+            @test isa(points, Vector)
+            @test all(point -> isa(point, Point2f), points)
+        end
+    end
+    
+    # INTEGRATION TESTS - These test computation + plotting together with safe fallbacks
+    @testset "Integrated Plotting Functions (Computation + Visualization)" begin
+        @testset "plot_param_line integration" begin
+            p = Point2f(0.0, 0.0)
+            q = Point2f(1.0, 1.0)
+            
+            # Test that function runs and returns correct computational results
+            # even if plotting fails in CI environments
+            try
+                points = plot_param_line(p, q, 3)
+                # Mathematical correctness is still tested even if plotting fails
+                @test length(points) == 3
+                @test isa(points[1], Point2f)
+                @test isa(points[end], Point2f)
+                # Test that points lie on the expected line
+                @test points[1] ≈ Point2f(1.0/3.0, 1.0/3.0) atol=1e-10
+                @test points[2] ≈ Point2f(2.0/3.0, 2.0/3.0) atol=1e-10
+                @test points[3] ≈ Point2f(1.0, 1.0) atol=1e-10
+            catch e
+                # Only catch plotting-related errors, not computational errors
+                if contains(string(e), "display") || contains(string(e), "GKS") || isa(e, ArgumentError)
+                    @test hasmethod(plot_param_line, (Point2f, Point2f, Int64))
+                else
+                    # Re-throw computational errors - these should fail the test
+                    rethrow(e)
+                end
+            end
+        end
+    end
+    
+    @testset "Line Geometry Functions" begin
+        # Test parametric_to_implicit_line
+        p = Point2f(0.0, 0.0)
+        v = [1.0, 1.0]  # direction vector for y = x line
+        a, b, c = parametric_to_implicit_line(p, v)
+        # For line y = x, implicit form should be x - y = 0
+        @test abs(a) ≈ abs(b)  # coefficients should have same magnitude
+        
+        # Test distance_to_implicit_line
+        # Distance from point (1,0) to line x - y = 0 (y = x)
+        a, b, c = 1.0, -1.0, 0.0
+        r = Point2f(1.0, 0.0)
+        dist = distance_to_implicit_line(a, b, c, r)
+        expected_dist = 1.0 / sqrt(2.0)  # distance formula for point to line
+        @test dist ≈ expected_dist
+        
+        # Test explicit_line
+        p1 = Point2f(0.0, 0.0)
+        p2 = Point2f(1.0, 2.0)
+        m, b = explicit_line(p1, p2)
+        @test m ≈ 2.0  # slope
+        @test b ≈ 0.0  # y-intercept
+        
+        # Test point_in_implicit_line
+        p = Point2f(0.0, 0.0)
+        q = Point2f(1.0, 1.0)
+        x = Point2f(0.5, 0.5)  # point on line between p and q
+        result = point_in_implicit_line(p, q, x)
+        @test isa(result, Bool) || isa(result, Number)  # function returns boolean or number
+        
+        # Test implicit_to_parametric_line
+        a, b, c = 1.0, -1.0, 0.0  # line x - y = 0
+        p_param, v_param = implicit_to_parametric_line(a, b, c)
+        @test isa(p_param, AbstractVector) || isa(p_param, Point2f)
+        @test length(v_param) == 2
+        
+        # Test implicit_line_point_normal_form
+        a, b, c = 3.0, 4.0, 5.0
+        â, b̂, ĉ = implicit_line_point_normal_form(a, b, c)
+        # Should normalize the coefficients
+        norm_factor = sqrt(a^2 + b^2)
+        @test â ≈ a/norm_factor
+        @test b̂ ≈ b/norm_factor
+        
+        # Test distance_to_parametric_line
+        p = Point2f(0.0, 0.0)  # point on line
+        v = [1.0, 0.0]         # direction vector
+        r = Point2f(0.0, 1.0)  # test point
+        dist_param = distance_to_parametric_line(p, v, r)
+        @test dist_param ≈ 1.0  # distance from (0,1) to x-axis
+        
+        # Test foot_of_line functions
+        P = Point2f(0.0, 0.0)
+        v = [1.0, 0.0]
+        R = Point2f(1.0, 1.0)
+        
+        # Test foot_of_line with tuple return
+        foot, dist = foot_of_line(P, v, R)
+        @test isa(foot, AbstractVector) || isa(foot, Point2f)
+        @test isa(dist, Number)
+        
+        # Test foot_of_line with boolean parameter
+        result = foot_of_line(P, v, R, false)
+        @test isa(result, Point2f) || isa(result, Tuple)
+    end
+    
+    @testset "Advanced Line Functions" begin
+        # Test intersection functions (just check they execute without error)
+        v1 = [1.0, 0.0]  # horizontal line direction
+        w1 = [0.0, 1.0]  # vertical line direction
+        p1 = Point2f(0.0, 1.0)  # point on first line
+        q1 = Point2f(1.0, 0.0)  # point on second line
+        try
+            intersection = intersection_2_parametric_lines(v1, w1, p1, q1)
+            @test true  # if it doesn't error, that's good
+        catch
+            @test false  # if it errors, that's bad
+        end
+        
+        # Test intersection_2_implicit_lines
+        # Line 1: x + y = 1
+        # Line 2: x - y = 0
+        a₁, b₁, c₁ = 1.0, 1.0, -1.0
+        a₂, b₂, c₂ = 1.0, -1.0, 0.0
+        try
+            intersection_impl = intersection_2_implicit_lines(a₁, b₁, c₁, a₂, b₂, c₂)
+            @test true  # if it doesn't error, that's good
+        catch
+            @test false  # if it errors, that's bad
+        end
+        
+        # Test Linear_Algebra.rationalize function
+        x = π
+        try
+            rational_approx = Linear_Algebra.rationalize(x; sigdigits=4)
+            @test isa(rational_approx, Rational) || isa(rational_approx, Number)
+        catch
+            # If the function doesn't exist or has different signature, just pass
+            @test true
+        end
+    end
+end
diff --git a/test/test_linear_algebra_transform.jl b/test/test_linear_algebra_transform.jl
new file mode 100644
index 0000000..6a76725
--- /dev/null
+++ b/test/test_linear_algebra_transform.jl
@@ -0,0 +1,135 @@
+using Test
+
+# Configure plotting for headless CI environments BEFORE loading Linear_Algebra
+ENV["GKSwstype"] = "100"  # Set GKS to use headless mode
+
+using Linear_Algebra
+using GeometryBasics, LinearAlgebra
+
+@testset "linear_algebra_transform.jl Tests" begin
+    
+    @testset "Transformation Matrix Functions" begin
+        # Test projection_matrix
+        x = [1.0, 0.0]
+        P = projection_matrix(x)
+        @test size(P) == (2, 2)
+        @test P * x ≈ x  # projecting onto itself should return itself
+        
+        # Test rotation_matrix  
+        d = 90  # 90 degrees (the function takes degrees, not radians)
+        R = rotation_matrix(d)
+        @test size(R) == (2, 2)
+        # Rotating [1,0] by 90° should give [0,1]
+        v = [1.0, 0.0]
+        rotated = R * v
+        @test rotated ≈ [0.0, 1.0] atol=1e-10
+        
+        # Test stretch_matrix
+        n = 2.0
+        S = stretch_matrix(n)
+        @test size(S) == (2, 2)
+        v = [1.0, 1.0]
+        stretched = S * v
+        @test stretched ≈ [2.0, 2.0]  # both components scaled by n
+        
+        # Test reflection_matrix
+        U = [1.0, 0.0]  # reflect across x-axis
+        R = reflection_matrix(U)
+        @test size(R) == (2, 2)
+        v = [1.0, 1.0]
+        reflected = R * v
+        @test reflected[1] ≈ 1.0   # x-component unchanged
+        @test reflected[2] ≈ -1.0  # y-component flipped
+        
+        # Test additional matrix functions
+        # Test projection_matrix_polar
+        n = 1.0
+        P_polar = projection_matrix_polar(n)
+        @test size(P_polar) == (2, 2)
+        
+        # Test projection_matrix_transpose
+        u = [1.0, 0.0]
+        P_transpose = projection_matrix_transpose(u)
+        @test size(P_transpose) == (2, 2)
+        
+        # Test rotation_matrix_ns (non-symbolic)
+        θ = π/4  # 45 degrees
+        R_ns = rotation_matrix_ns(θ)
+        @test size(R_ns) == (2, 2)
+        
+        # Test reflection_matrix_rational
+        U_rational = [1.0, 1.0]
+        R_rational = reflection_matrix_rational(U_rational)
+        @test size(R_rational) == (2, 2)
+    end
+    
+    @testset "Symbolic Matrix Functions" begin
+        # Test symbolic matrix functions (should return symbolic expressions)
+        P_sym = projection_matrix_symbolic()
+        @test size(P_sym) == (2, 2)
+        
+        P_sym_polar = projection_matrix_symbolic_polar()
+        @test size(P_sym_polar) == (2, 2)
+        
+        R_sym = rotation_matrix_symbolic()
+        @test size(R_sym) == (2, 2)
+        
+        S_sym = stretch_matrix_symbolic()
+        @test size(S_sym) == (2, 2)
+        
+        Refl_sym = reflection_matrix_symbolic()
+        @test size(Refl_sym) == (2, 2)
+    end
+    
+    @testset "Matrix Mathematical Properties" begin
+        # Test mathematical properties of transformation matrices
+        
+        # Rotation matrix should be orthogonal (R^T * R = I)
+        d = 45  # 45 degrees
+        R = rotation_matrix(d)
+        I_approx = R' * R
+        @test I_approx ≈ [1.0 0.0; 0.0 1.0] atol=1e-10
+        
+        # Projection matrix should be idempotent (P^2 = P)
+        x = [1.0, 1.0]
+        P = projection_matrix(x)
+        P_squared = P * P
+        @test P_squared ≈ P atol=1e-10
+        
+        # Reflection matrix should be its own inverse (R^2 = I)
+        U = [1.0, 0.0]
+        Refl = reflection_matrix(U)
+        Refl_squared = Refl * Refl
+        @test Refl_squared ≈ [1.0 0.0; 0.0 1.0] atol=1e-10
+        
+        # Stretch matrix should scale uniformly
+        n = 3.0
+        S = stretch_matrix(n)
+        v = [1.0, 1.0]
+        scaled = S * v
+        @test norm(scaled) ≈ n * norm(v) atol=1e-10
+    end
+    
+    @testset "Edge Cases and Special Values" begin
+        # Test rotation by 0 degrees (should be identity)
+        R_zero = rotation_matrix(0)
+        @test R_zero ≈ [1.0 0.0; 0.0 1.0] atol=1e-10
+        
+        # Test rotation by 180 degrees
+        R_180 = rotation_matrix(180)
+        v = [1.0, 0.0]
+        rotated_180 = R_180 * v
+        @test rotated_180 ≈ [-1.0, 0.0] atol=1e-10
+        
+        # Test stretch by 1 (should be identity)
+        S_one = stretch_matrix(1.0)
+        @test S_one ≈ [1.0 0.0; 0.0 1.0] atol=1e-10
+        
+        # Test projection onto unit vector
+        unit_x = [1.0, 0.0]
+        P_unit = projection_matrix(unit_x)
+        test_vector = [3.0, 4.0]
+        projected = P_unit * test_vector
+        @test projected ≈ [3.0, 0.0] atol=1e-10  # Should project onto x-axis
+    end
+end