CoCalc News

# Essential physics imports
import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate, optimize
import sympy as sp

# Welcome to computational physics!
print("Welcome to Computational Physics with CoCalc!")

# Physical constants (in SI units)
c = 299792458          # Speed of light (m/s)
h = 6.62607015e-34     # Planck constant (J⋅s)
hbar = h / (2 * np.pi) # Reduced Planck constant
k_B = 1.380649e-23     # Boltzmann constant (J/K)
e = 1.602176634e-19    # Elementary charge (C)
m_e = 9.1093837015e-31 # Electron mass (kg)
m_p = 1.67262192369e-27 # Proton mass (kg)

print(f"Speed of light: {c:.0e} m/s")
print(f"Planck constant: {h:.3e} J⋅s")
print(f"Electron mass: {m_e:.3e} kg")

# Classical mechanics: projectile motion
def projectile_motion():
    """
    Simulate projectile motion with air resistance
    """
    
    # Parameters
    g = 9.81          # Gravitational acceleration (m/s²)
    v0 = 50           # Initial velocity (m/s)
    angle = 45        # Launch angle (degrees)
    m = 1.0           # Mass (kg)
    b = 0.1           # Air resistance coefficient
    
    # Convert angle to radians
    theta = np.radians(angle)
    
    # Initial conditions
    vx0 = v0 * np.cos(theta)
    vy0 = v0 * np.sin(theta)
    
    # Differential equation: F = ma = -mg - bv
    def equations_of_motion(t, state):
        x, y, vx, vy = state
        
        # Air resistance force
        v_magnitude = np.sqrt(vx**2 + vy**2)
        if v_magnitude > 0:
            ax = -b * vx * v_magnitude / m
            ay = -g - b * vy * v_magnitude / m
        else:
            ax = 0
            ay = -g
        
        return [vx, vy, ax, ay]
    
    # Solve the differential equation
    t_span = (0, 10)
    initial_state = [0, 0, vx0, vy0]
    
    # Event function to stop when projectile hits ground
    def hit_ground(t, state):
        return state[1]  # y-coordinate
    hit_ground.terminal = True
    hit_ground.direction = -1
    
    solution = integrate.solve_ivp(
        equations_of_motion, t_span, initial_state,
        events=hit_ground, dense_output=True, rtol=1e-8
    )
    
    # Extract trajectory
    t_flight = solution.t
    trajectory = solution.y
    
    # Plot trajectory
    plt.figure(figsize=(10, 6))
    plt.plot(trajectory[0], trajectory[1], 'b-', linewidth=2, label='With air resistance')
    
    # Compare with no air resistance (analytical solution)
    t_no_air = np.linspace(0, t_flight[-1], 100)
    x_no_air = vx0 * t_no_air
    y_no_air = vy0 * t_no_air - 0.5 * g * t_no_air**2
    
    plt.plot(x_no_air, y_no_air, 'r--', linewidth=2, label='No air resistance')
    
    plt.xlabel('Horizontal Distance (m)')
    plt.ylabel('Height (m)')
    plt.title(f'Projectile Motion (v₀={v0} m/s, θ={angle}°)')
    plt.legend()
    plt.grid(True, alpha=0.3)
    plt.show()
    
    # Calculate range
    range_with_air = trajectory[0][-1]
    range_no_air = v0**2 * np.sin(2*theta) / g
    
    print(f"Range with air resistance: {range_with_air:.1f} m")
    print(f"Range without air resistance: {range_no_air:.1f} m")
    print(f"Air resistance reduces range by: {(1 - range_with_air/range_no_air)*100:.1f}%")
    
    return t_flight, trajectory

# Execute projectile motion simulation
time, trajectory = projectile_motion()

# Simple harmonic motion and damped oscillations
def harmonic_oscillator():
    """
    Study simple harmonic motion and damping effects
    """
    
    # Parameters
    omega_0 = 2.0     # Natural frequency (rad/s)
    gamma = 0.1       # Damping coefficient
    F0 = 1.0          # Driving force amplitude
    omega_d = 1.8     # Driving frequency
    
    # Equation of motion: m*x'' + gamma*x' + k*x = F*cos(omega_d*t)
    # Using omega_0^2 = k/m, we get: x'' + gamma*x' + omega_0^2*x = (F0/m)*cos(omega_d*t)
    
    def oscillator_equation(t, state):
        x, v = state
        
        # Driving force
        driving_force = F0 * np.cos(omega_d * t)
        
        # Acceleration
        a = -2*gamma*v - omega_0**2*x + driving_force
        
        return [v, a]
    
    # Time array
    t = np.linspace(0, 20, 1000)
    
    # Solve for different initial conditions
    solutions = {}
    initial_conditions = [
        (1.0, 0.0, "x₀=1, v₀=0"),
        (0.0, 2.0, "x₀=0, v₀=2"),
        (0.5, 1.0, "x₀=0.5, v₀=1")
    ]
    
    plt.figure(figsize=(12, 8))
    
    for i, (x0, v0, label) in enumerate(initial_conditions):
        solution = integrate.solve_ivp(
            oscillator_equation, (0, 20), [x0, v0],
            t_eval=t, rtol=1e-8
        )
        
        solutions[label] = solution
        
        plt.subplot(2, 2, i+1)
        plt.plot(solution.t, solution.y[0], 'b-', linewidth=1.5)
        plt.xlabel('Time (s)')
        plt.ylabel('Position (m)')
        plt.title(f'Damped Oscillator: {label}')
        plt.grid(True, alpha=0.3)
    
    # Phase space plot
    plt.subplot(2, 2, 4)
    for label, solution in solutions.items():
        plt.plot(solution.y[0], solution.y[1], linewidth=1.5, label=label)
    
    plt.xlabel('Position (m)')
    plt.ylabel('Velocity (m/s)')
    plt.title('Phase Space')
    plt.legend()
    plt.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.show()
    
    # Analyze frequency response
    frequencies = np.linspace(0.1, 4.0, 100)
    amplitude_response = []
    
    for omega in frequencies:
        # Steady-state amplitude for driven oscillator
        amplitude = F0 / np.sqrt((omega_0**2 - omega**2)**2 + (2*gamma*omega)**2)
        amplitude_response.append(amplitude)
    
    plt.figure(figsize=(10, 6))
    plt.plot(frequencies, amplitude_response, 'r-', linewidth=2)
    plt.axvline(omega_0, color='b', linestyle='--', alpha=0.7, label=f'ω₀ = {omega_0}')
    plt.xlabel('Driving Frequency (rad/s)')
    plt.ylabel('Steady-State Amplitude')
    plt.title('Frequency Response of Damped Harmonic Oscillator')
    plt.legend()
    plt.grid(True, alpha=0.3)
    plt.show()
    
    return solutions, amplitude_response

# Execute oscillator analysis
oscillator_solutions, frequency_response = harmonic_oscillator()

# Electromagnetic field visualization
def electromagnetic_fields():
    """
    Visualize electric and magnetic fields
    """
    
    # Electric field of point charges
    def electric_field_point_charges():
        """Visualize electric field of multiple point charges"""
        
        # Define point charges: (x, y, charge)
        charges = [
            (2, 2, 1.0),    # Positive charge
            (-2, -2, -1.0), # Negative charge
            (0, 3, 0.5),    # Smaller positive charge
        ]
        
        # Create grid
        x = np.linspace(-4, 4, 20)
        y = np.linspace(-4, 4, 20)
        X, Y = np.meshgrid(x, y)
        
        # Calculate electric field
        Ex = np.zeros_like(X)
        Ey = np.zeros_like(Y)
        V = np.zeros_like(X)  # Electric potential
        
        k = 8.99e9  # Coulomb constant (N⋅m²/C²)
        
        for charge_x, charge_y, q in charges:
            # Distance from each grid point to the charge
            dx = X - charge_x
            dy = Y - charge_y
            r = np.sqrt(dx**2 + dy**2)
            
            # Avoid division by zero
            r = np.where(r < 0.1, 0.1, r)
            
            # Electric field components
            Ex += k * q * dx / r**3
            Ey += k * q * dy / r**3
            
            # Electric potential
            V += k * q / r
        
        # Plot electric field and equipotential lines
        plt.figure(figsize=(12, 5))
        
        # Electric field vectors
        plt.subplot(1, 2, 1)
        plt.quiver(X, Y, Ex, Ey, np.sqrt(Ex**2 + Ey**2), cmap='viridis', alpha=0.7)
        
        # Mark charge locations
        for charge_x, charge_y, q in charges:
            color = 'red' if q > 0 else 'blue'
            size = abs(q) * 100
            plt.scatter(charge_x, charge_y, c=color, s=size, edgecolors='black')
        
        plt.xlabel('x (m)')
        plt.ylabel('y (m)')
        plt.title('Electric Field')
        plt.axis('equal')
        plt.colorbar(label='Field Magnitude')
        
        # Equipotential lines
        plt.subplot(1, 2, 2)
        contour = plt.contour(X, Y, V, levels=20, colors='blue', alpha=0.6)
        plt.clabel(contour, inline=True, fontsize=8)
        
        # Mark charge locations
        for charge_x, charge_y, q in charges:
            color = 'red' if q > 0 else 'blue'
            size = abs(q) * 100
            plt.scatter(charge_x, charge_y, c=color, s=size, edgecolors='black')
        
        plt.xlabel('x (m)')
        plt.ylabel('y (m)')
        plt.title('Electric Potential')
        plt.axis('equal')
        
        plt.tight_layout()
        plt.show()
        
        return X, Y, Ex, Ey, V
    
    # Magnetic field of current loop
    def magnetic_field_current_loop():
        """Magnetic field of a circular current loop"""
        
        # Parameters
        I = 1.0      # Current (A)
        R = 1.0      # Loop radius (m)
        mu_0 = 4*np.pi*1e-7  # Permeability of free space
        
        # Calculate field along axis
        z = np.linspace(-3, 3, 100)
        Bz = mu_0 * I * R**2 / (2 * (R**2 + z**2)**(3/2))
        
        plt.figure(figsize=(10, 6))
        
        plt.subplot(1, 2, 1)
        plt.plot(z, Bz*1e6, 'b-', linewidth=2)
        plt.xlabel('Distance along axis (m)')
        plt.ylabel('Magnetic Field (μT)')
        plt.title(f'Magnetic Field of Current Loop (I={I}A, R={R}m)')
        plt.grid(True, alpha=0.3)
        
        # Field lines visualization (simplified)
        plt.subplot(1, 2, 2)
        theta = np.linspace(0, 2*np.pi, 100)
        loop_x = R * np.cos(theta)
        loop_y = R * np.sin(theta)
        
        # Draw current loop
        plt.plot(loop_x, loop_y, 'r-', linewidth=3, label='Current Loop')
        
        # Simplified field line representation
        r_field = np.linspace(0.2, 3, 10)
        for r in r_field:
            # Approximate field lines
            phi = np.linspace(0, 2*np.pi, 50)
            field_x = r * np.cos(phi)
            field_y = 0.5 * r * np.sin(phi)
            plt.plot(field_x, field_y, 'b-', alpha=0.5, linewidth=1)
            plt.plot(field_x, -field_y, 'b-', alpha=0.5, linewidth=1)
        
        plt.xlabel('x (m)')
        plt.ylabel('y (m)')
        plt.title('Magnetic Field Lines')
        plt.axis('equal')
        plt.legend()
        
        plt.tight_layout()
        plt.show()
        
        return z, Bz
    
    # Execute electromagnetic field calculations
    electric_results = electric_field_point_charges()
    magnetic_results = magnetic_field_current_loop()
    
    return electric_results, magnetic_results

# Execute electromagnetic field analysis
em_results = electromagnetic_fields()

# Thermodynamics and kinetic theory
def thermal_physics():
    """
    Explore thermodynamics and statistical mechanics
    """
    
    # Ideal gas law and kinetic theory
    def kinetic_theory_simulation():
        """Simulate kinetic theory of gases"""
        
        # Parameters
        N = 1000         # Number of particles
        T = 300          # Temperature (K)
        m = 4.65e-26     # Mass of N2 molecule (kg)
        
        # Maxwell-Boltzmann distribution
        def maxwell_boltzmann(v, T, m):
            """Maxwell-Boltzmann speed distribution"""
            k_B = 1.38e-23
            factor = 4 * np.pi * (m / (2 * np.pi * k_B * T))**(3/2)
            return factor * v**2 * np.exp(-m * v**2 / (2 * k_B * T))
        
        # Generate random velocities according to Maxwell-Boltzmann
        # (Simplified: using normal distribution for each component)
        sigma = np.sqrt(k_B * T / m)  # Standard deviation for each component
        
        vx = np.random.normal(0, sigma, N)
        vy = np.random.normal(0, sigma, N)
        vz = np.random.normal(0, sigma, N)
        
        # Calculate speeds
        speeds = np.sqrt(vx**2 + vy**2 + vz**2)
        
        # Plot speed distribution
        plt.figure(figsize=(12, 8))
        
        plt.subplot(2, 2, 1)
        plt.hist(speeds, bins=50, density=True, alpha=0.7, label='Simulation')
        
        # Theoretical Maxwell-Boltzmann distribution
        v_theory = np.linspace(0, np.max(speeds), 200)
        mb_theory = maxwell_boltzmann(v_theory, T, m)
        plt.plot(v_theory, mb_theory, 'r-', linewidth=2, label='Theory')
        
        plt.xlabel('Speed (m/s)')
        plt.ylabel('Probability Density')
        plt.title(f'Maxwell-Boltzmann Distribution (T={T}K)')
        plt.legend()
        plt.grid(True, alpha=0.3)
        
        # Calculate characteristic speeds
        v_avg = np.mean(speeds)
        v_rms = np.sqrt(np.mean(speeds**2))
        v_mp = np.sqrt(2 * k_B * T / m)  # Most probable speed
        
        plt.axvline(v_avg, color='blue', linestyle='--', label=f'Average: {v_avg:.0f} m/s')
        plt.axvline(v_rms, color='green', linestyle='--', label=f'RMS: {v_rms:.0f} m/s')
        plt.axvline(v_mp, color='red', linestyle='--', label=f'Most probable: {v_mp:.0f} m/s')
        plt.legend()
        
        # Energy distribution
        plt.subplot(2, 2, 2)
        kinetic_energies = 0.5 * m * speeds**2
        plt.hist(kinetic_energies / k_B, bins=50, density=True, alpha=0.7)
        plt.xlabel('Kinetic Energy / k_B (K)')
        plt.ylabel('Probability Density')
        plt.title('Kinetic Energy Distribution')
        plt.grid(True, alpha=0.3)
        
        # Temperature dependence
        plt.subplot(2, 2, 3)
        temperatures = [200, 300, 400, 500]
        v_range = np.linspace(0, 1500, 200)
        
        for T_temp in temperatures:
            mb_dist = maxwell_boltzmann(v_range, T_temp, m)
            plt.plot(v_range, mb_dist, linewidth=2, label=f'T = {T_temp}K')
        
        plt.xlabel('Speed (m/s)')
        plt.ylabel('Probability Density')
        plt.title('Temperature Dependence')
        plt.legend()
        plt.grid(True, alpha=0.3)
        
        # Pressure calculation from kinetic theory
        plt.subplot(2, 2, 4)
        volumes = np.linspace(0.01, 0.1, 100)  # m³
        n_moles = 1.0  # mol
        R = 8.314  # J/(mol⋅K)
        
        # Ideal gas law: PV = nRT
        pressures_ideal = n_moles * R * T / volumes
        
        plt.plot(volumes*1000, pressures_ideal/1000, 'b-', linewidth=2, label='Ideal Gas Law')
        plt.xlabel('Volume (L)')
        plt.ylabel('Pressure (kPa)')
        plt.title(f'P-V Diagram (T={T}K)')
        plt.grid(True, alpha=0.3)
        plt.legend()
        
        plt.tight_layout()
        plt.show()
        
        print(f"Theoretical values at T = {T}K:")
        print(f"Average speed: {np.sqrt(8*k_B*T/(np.pi*m)):.0f} m/s")
        print(f"RMS speed: {np.sqrt(3*k_B*T/m):.0f} m/s")
        print(f"Most probable speed: {np.sqrt(2*k_B*T/m):.0f} m/s")
        
        return speeds, kinetic_energies
    
    # Heat transfer
    def heat_transfer_simulation():
        """Simulate heat conduction"""
        
        # One-dimensional heat equation: ∂T/∂t = α ∂²T/∂x²
        # where α is thermal diffusivity
        
        # Parameters
        L = 1.0          # Length (m)
        alpha = 1e-4     # Thermal diffusivity (m²/s)
        dx = 0.01        # Spatial step
        dt = 0.1         # Time step
        t_max = 3600     # Maximum time (s)
        
        # Grid
        x = np.arange(0, L + dx, dx)
        nx = len(x)
        
        # Check stability condition
        stability = alpha * dt / dx**2
        print(f"Stability parameter: {stability:.3f} (should be ≤ 0.5)")
        
        # Initial condition: hot at one end, cold at the other
        T = np.zeros(nx)
        T[0] = 100    # Hot end (°C)
        T[-1] = 0     # Cold end (°C)
        
        # Time stepping
        times = np.arange(0, t_max + dt, dt)
        T_history = []
        
        for t in times[::int(len(times)/10)]:  # Store every 10th time step
            T_history.append(T.copy())
            
            # Update temperature using finite difference
            T_new = T.copy()
            for i in range(1, nx-1):
                T_new[i] = T[i] + alpha * dt / dx**2 * (T[i+1] - 2*T[i] + T[i-1])
            
            T = T_new
        
        # Plot temperature evolution
        plt.figure(figsize=(10, 6))
        
        plt.subplot(1, 2, 1)
        for i, T_snapshot in enumerate(T_history):
            time_label = f't = {i * t_max / 10 / 60:.0f} min'
            plt.plot(x, T_snapshot, linewidth=2, label=time_label)
        
        plt.xlabel('Position (m)')
        plt.ylabel('Temperature (°C)')
        plt.title('Heat Conduction Over Time')
        plt.legend()
        plt.grid(True, alpha=0.3)
        
        # Steady-state solution (linear)
        plt.subplot(1, 2, 2)
        T_steady = 100 * (1 - x / L)  # Linear temperature profile
        plt.plot(x, T_history[-1], 'b-', linewidth=2, label='Numerical solution')
        plt.plot(x, T_steady, 'r--', linewidth=2, label='Analytical steady-state')
        plt.xlabel('Position (m)')
        plt.ylabel('Temperature (°C)')
        plt.title('Final Temperature Profile')
        plt.legend()
        plt.grid(True, alpha=0.3)
        
        plt.tight_layout()
        plt.show()
        
        return x, T_history
    
    # Execute thermal physics simulations
    speeds, energies = kinetic_theory_simulation()
    x_heat, T_heat = heat_transfer_simulation()
    
    return speeds, energies, x_heat, T_heat

# Execute thermal physics analysis
thermal_results = thermal_physics()

# Dimensional analysis tools
def dimensional_analysis():
    """
    Explore dimensional analysis and scaling laws
    """
    
    # Fundamental dimensions
    dimensions = {
        'length': 'L',
        'mass': 'M', 
        'time': 'T',
        'electric_current': 'I',
        'temperature': 'Θ',
        'amount': 'N',
        'luminous_intensity': 'J'
    }
    
    # Common physical quantities and their dimensions
    quantities = {
        'velocity': 'L T⁻¹',
        'acceleration': 'L T⁻²',
        'force': 'M L T⁻²',
        'energy': 'M L² T⁻²',
        'power': 'M L² T⁻³',
        'pressure': 'M L⁻¹ T⁻²',
        'electric_field': 'M L T⁻³ I⁻¹',
        'magnetic_field': 'M T⁻² I⁻¹'
    }
    
    print("DIMENSIONAL ANALYSIS EXAMPLES")
    print("="*40)
    
    # Example 1: Pendulum period
    print("Example 1: Simple Pendulum Period")
    print("Variables: length L, gravity g, mass m")
    print("Dimensions: [L] = L, [g] = L T⁻², [m] = M")
    print("Proposed form: T = L^a × g^b × m^c")
    print("Dimension equation: T = L^a × (L T⁻²)^b × M^c")
    print("                   T = L^(a+b) × T^(-2b) × M^c")
    print("Matching dimensions:")
    print("  Time: 1 = -2b  →  b = -1/2")
    print("  Length: 0 = a + b  →  a = 1/2")
    print("  Mass: 0 = c  →  c = 0")
    print("Result: T ∝ √(L/g)")
    
    # Numerical verification
    L_values = np.array([0.5, 1.0, 1.5, 2.0])  # lengths in meters
    g = 9.81
    T_predicted = 2 * np.pi * np.sqrt(L_values / g)
    
    plt.figure(figsize=(10, 6))
    plt.subplot(1, 2, 1)
    plt.plot(L_values, T_predicted, 'bo-', linewidth=2, markersize=8)
    plt.xlabel('Length (m)')
    plt.ylabel('Period (s)')
    plt.title('Pendulum Period vs Length')
    plt.grid(True, alpha=0.3)
    
    # Example 2: Drag force scaling
    print("\nExample 2: Drag Force at High Speed")
    print("Variables: velocity v, density ρ, area A, coefficient C_d")
    print("Expected form: F_drag = ½ × C_d × ρ × A × v²")
    
    # Show velocity scaling
    velocities = np.linspace(10, 100, 20)
    rho = 1.225  # Air density (kg/m³)
    A = 2.0      # Cross-sectional area (m²)
    C_d = 0.3    # Drag coefficient
    
    F_drag = 0.5 * C_d * rho * A * velocities**2
    
    plt.subplot(1, 2, 2)
    plt.plot(velocities, F_drag, 'r-', linewidth=2)
    plt.xlabel('Velocity (m/s)')
    plt.ylabel('Drag Force (N)')
    plt.title('Drag Force vs Velocity (∝ v²)')
    plt.grid(True, alpha=0.3)
    
    plt.tight_layout()
    plt.show()
    
    return quantities, T_predicted, F_drag

# Execute dimensional analysis
dim_analysis = dimensional_analysis()

# Error analysis and uncertainty propagation
def error_analysis():
    """
    Demonstrate error analysis and uncertainty propagation
    """
    
    print("ERROR ANALYSIS IN PHYSICS")
    print("="*30)
    
    # Example: Measuring gravitational acceleration
    # Using pendulum: g = 4π²L/T²
    
    # Measurements with uncertainties
    L = 1.000          # Length (m)
    dL = 0.001         # Uncertainty in length (m)
    T = 2.006          # Period (s)
    dT = 0.010         # Uncertainty in period (s)
    
    # Calculate g
    g = 4 * np.pi**2 * L / T**2
    
    # Uncertainty propagation: g = 4π²L/T²
    # Relative uncertainty: δg/g = δL/L + 2δT/T
    relative_uncertainty = dL/L + 2*dT/T
    dg = g * relative_uncertainty
    
    print(f"Measurement results:")
    print(f"Length: L = {L:.3f} ± {dL:.3f} m")
    print(f"Period: T = {T:.3f} ± {dT:.3f} s")
    print(f"Calculated g = {g:.2f} ± {dg:.2f} m/s²")
    print(f"Accepted value: 9.81 m/s²")
    print(f"Percent error: {abs(g - 9.81)/9.81 * 100:.1f}%")
    
    # Monte Carlo error propagation
    def monte_carlo_errors(n_samples=10000):
        """Use Monte Carlo method for error propagation"""
        
        # Generate random samples
        L_samples = np.random.normal(L, dL, n_samples)
        T_samples = np.random.normal(T, dT, n_samples)
        
        # Calculate g for each sample
        g_samples = 4 * np.pi**2 * L_samples / T_samples**2
        
        # Statistics
        g_mean = np.mean(g_samples)
        g_std = np.std(g_samples)
        
        print(f"\nMonte Carlo results ({n_samples} samples):")
        print(f"g = {g_mean:.2f} ± {g_std:.2f} m/s²")
        
        # Plot distribution
        plt.figure(figsize=(10, 6))
        
        plt.subplot(1, 2, 1)
        plt.hist(g_samples, bins=50, density=True, alpha=0.7, edgecolor='black')
        plt.axvline(g_mean, color='red', linestyle='--', linewidth=2, label=f'Mean: {g_mean:.2f}')
        plt.axvline(9.81, color='green', linestyle='--', linewidth=2, label='True value: 9.81')
        plt.xlabel('g (m/s²)')
        plt.ylabel('Probability Density')
        plt.title('Distribution of g Values')
        plt.legend()
        plt.grid(True, alpha=0.3)
        
        # Error vs number of measurements
        plt.subplot(1, 2, 2)
        n_measurements = np.logspace(1, 4, 20).astype(int)
        standard_errors = []
        
        for n in n_measurements:
            # Standard error decreases as 1/√n
            standard_error = g_std / np.sqrt(n)
            standard_errors.append(standard_error)
        
        plt.loglog(n_measurements, standard_errors, 'bo-', linewidth=2)
        plt.loglog(n_measurements, 1/np.sqrt(n_measurements), 'r--', 
                  linewidth=2, label='∝ 1/√n')
        plt.xlabel('Number of Measurements')
        plt.ylabel('Standard Error (m/s²)')
        plt.title('Error Reduction with More Measurements')
        plt.legend()
        plt.grid(True, alpha=0.3)
        
        plt.tight_layout()
        plt.show()
        
        return g_samples, g_mean, g_std
    
    # Execute Monte Carlo analysis
    g_samples, g_mean, g_std = monte_carlo_errors()
    
    return g, dg, g_samples

# Execute error analysis
error_results = error_analysis()

# Welcome to mathematical computing!
# Let's start with basic symbolic mathematics

# Working with symbolic variables
x, y, z = var('x y z')

# Symbolic expressions and algebra
expr = (x + 1)^3
expanded = expand(expr)
factored = factor(x^3 + 3*x^2 + 3*x + 1)

print("Original expression:", expr)
print("Expanded form:", expanded)
print("Factored form:", factored)

# Basic calculus
f = x^2 + 3*x + 2
derivative = diff(f, x)
integral = integrate(f, x)

print("\nFunction: f(x) =", f)
print("Derivative: f'(x) =", derivative)
print("Integral: ∫f(x)dx =", integral)

Original expression: (x + 1)^3
Expanded form: x^3 + 3*x^2 + 3*x + 1
Factored form: (x + 1)^3
Function: f(x) = x^2 + 3*x + 2
Derivative: f'(x) = 2*x + 3
Integral: ∫f(x)dx = 1/3*x^3 + 3/2*x^2 + 2*x

# Number theory explorations
print("=== NUMBER THEORY ADVENTURES ===")

# Prime numbers and factorization
n = 2023
print("\nIs", n, "prime?", is_prime(n))
print("Prime factorization of", n, ":", factor(n))

# Finding primes
primes_up_to_100 = [p for p in range(2, 101) if is_prime(p)]
print("\nFirst few primes:", primes_up_to_100[:10])
print("Number of primes up to 100:", len(primes_up_to_100))

# The Euclidean algorithm
a, b = 48, 18
gcd_result = gcd(a, b)
print("\ngcd(", a, ",", b, ") =", gcd_result)

# Extended Euclidean algorithm
g, u, v = xgcd(a, b)
print("Extended:", u, "·", a, " + ", v, "·", b, " = ", g)

=== NUMBER THEORY ADVENTURES ===

Is 2023 prime? False
Prime factorization of 2023 : 7 * 17^2

First few primes: [2, 3, 5, 7, 11, 13, 17, 19, 23, 29]
Number of primes up to 100: 25

gcd( 48 , 18 ) = 6
Extended: -1 · 48  +  3 · 18  =  6

# Group theory basics
print("=== ALGEBRAIC STRUCTURES ===")

# Modular arithmetic
print("\nModular arithmetic:")
for i in range(12):
    print(i, "mod 5 =", i % 5)

# Working with permutations
# Create a simple permutation group
S3 = SymmetricGroup(3)
print("\nSymmetric group S3 has order:", S3.order())
print("Elements of S3:", list(S3))

# Group operations
g1 = S3([2, 1, 3])  # Swap first two elements
g2 = S3([1, 3, 2])  # Swap last two elements
product = g1 * g2
print("\nPermutation multiplication:")
print(g1, "*", g2, "=", product)

=== ALGEBRAIC STRUCTURES ===

Modular arithmetic:
0 mod 5 = 0
1 mod 5 = 1
2 mod 5 = 2
3 mod 5 = 3
4 mod 5 = 4
5 mod 5 = 0
6 mod 5 = 1
7 mod 5 = 2
8 mod 5 = 3
9 mod 5 = 4
10 mod 5 = 0
11 mod 5 = 1

Symmetric group S3 has order: 6
Elements of S3: [(), (2,3), (1,2), (1,2,3), (1,3,2), (1,3)]

Permutation multiplication:
(1,2) * (2,3) = (1,2,3)

# Calculus explorations
print("=== CALCULUS AND ANALYSIS ===")

# Limits
x = var('x')
limit_expr = sin(x)/x
limit_result = limit(limit_expr, x=0)
print("lim(x→0) sin(x)/x =", limit_result)

# Series expansions
f = e^x
taylor_series = f.taylor(x, 0, 5)  # Taylor series around x=0, degree 5
print("\nTaylor series of e^x:", taylor_series)

# Definite integrals
integral_result = integrate(x^2 * e^(-x^2), x, -oo, oo)
print("\n∫_{-∞}^{∞} x²e^(-x²) dx =", integral_result)

# Plotting functions
p1 = plot(sin(x), (x, -2*pi, 2*pi), color='blue', legend_label='sin(x)')
p2 = plot(cos(x), (x, -2*pi, 2*pi), color='red', legend_label='cos(x)')
combined_plot = p1 + p2
combined_plot.show()