Chapter 5: Functions & Modules - Building Reusable Scientific Code

def calculate_mean(values):
    """
    Calculate the arithmetic mean of a list of values.
    
    This is our first function - notice the structure!
    """
    # Note: assert is for development/debugging
    # In production, use explicit validation:
    assert len(values) > 0, "Cannot calculate mean of empty list"
    # Better production version would be:
    # if len(values) == 0:
    #     raise ValueError("Cannot calculate mean of empty list")
    
    total = sum(values)
    count = len(values)
    mean = total / count
    return mean

# Using the function
measurements = [23.5, 24.1, 23.8, 24.3, 23.9]
avg = calculate_mean(measurements)
print(f"Mean temperature: {avg:.2f}°C")

# The assert would raise an error with empty list:
# empty_mean = calculate_mean([])  
# AssertionError: Cannot calculate mean of empty list

def report_statistics(data, name="Dataset"):
    """Report basic statistics without returning values."""
    mean = sum(data) / len(data)
    minimum = min(data)
    maximum = max(data)
    
    print(f"Statistics for {name}:")
    print(f"  Mean: {mean:.3f}")
    print(f"  Range: [{minimum:.3f}, {maximum:.3f}]")
    # No return statement - returns None implicitly

# Report some calculations
temperatures = [20.1, 21.5, 19.8, 22.3, 20.9]
report_statistics(temperatures, "Temperature (°C)")

def analyze_data(values):
    """
    Calculate multiple statistics at once.
    
    Returns:
        mean, std_dev, min_val, max_val
    """
    n = len(values)
    mean = sum(values) / n
    
    # Calculate standard deviation
    variance = sum((x - mean)**2 for x in values) / n
    std_dev = variance ** 0.5
    
    return mean, std_dev, min(values), max(values)

# Analyze experimental data
data = [9.8, 10.1, 9.9, 10.2, 9.7, 10.0]
mean, std, min_val, max_val = analyze_data(data)

print(f"Mean: {mean:.3f}")
print(f"Std Dev: {std:.3f}")
print(f"Range: [{min_val}, {max_val}]")

def kinetic_energy_cgs(mass_g, velocity_cms):
    """
    Calculate kinetic energy in ergs.
    
    Parameters:
        mass_g: mass in grams
        velocity_cms: velocity in cm/s
    
    Returns:
        energy in ergs (g⋅cm²/s²)
    """
    energy_ergs = 0.5 * mass_g * velocity_cms**2
    return energy_ergs

# Example: electron moving at 1% speed of light
electron_mass = 9.109e-28  # grams
c_light = 2.998e10  # speed of light in cm/s
electron_velocity = 0.01 * c_light  # 1% of c
ke = kinetic_energy_cgs(electron_mass, electron_velocity)

print(f"Electron kinetic energy: {ke:.2e} ergs")
print(f"In eV: {ke/1.602e-12:.2f} eV")  # 1 eV = 1.602e-12 ergs

"""
DESIGN: Function to validate numerical calculation

PURPOSE: Ensure numerical results are physically reasonable
INPUT: value, expected_range
OUTPUT: boolean (True if valid)
PROCESS:
    - Check if value is finite
    - Check if value is within range
"""

def validate_result(value, min_val=None, max_val=None):
    """
    Validate a numerical result is reasonable.
    
    Simple validation - we'll expand this pattern later!
    """
    import math
    
    # Check if finite
    if not math.isfinite(value):
        return False
    
    # Check bounds if provided
    if min_val is not None and value < min_val:
        return False
    if max_val is not None and value > max_val:
        return False
    
    return True

# Test validation
energy = kinetic_energy_cgs(1.0, 100.0)
print(f"Energy = {energy} ergs")
print(f"Valid (positive)? {validate_result(energy, min_val=0)}")
print(f"Valid (in range)? {validate_result(energy, 0, 10000)}")

def escape_velocity(mass_g, radius_cm):
    """
    Calculate escape velocity from a spherical body.
    
    Parameters:
        mass_g: mass of the body in grams
        radius_cm: radius of the body in cm
    
    Returns:
        escape velocity in cm/s
    """
    G = 6.674e-8  # gravitational constant in cm³/(g⋅s²)
    v_escape = (2 * G * mass_g / radius_cm) ** 0.5
    return v_escape

# Earth parameters
earth_mass = 5.972e27    # grams
earth_radius = 6.371e8   # cm

# Correct calculation
v_correct = escape_velocity(earth_mass, earth_radius)

# WRONG: Arguments reversed
v_wrong = escape_velocity(earth_radius, earth_mass)

print(f"Correct escape velocity: {v_correct:.2e} cm/s = {v_correct/1e5:.1f} km/s")
print(f"Wrong (args reversed): {v_wrong:.2e} cm/s")
print(f"The wrong value is {v_correct/v_wrong:.0f}x too small!")

# Keywords make intent clear
v_safe = escape_velocity(radius_cm=earth_radius, mass_g=earth_mass)
print(f"\nWith keywords: {v_safe:.2e} cm/s")

def simulate_decay(initial_atoms, half_life_years=5730, time_years=0):
    """
    Simulate radioactive decay.
    
    Parameters:
        initial_atoms: starting number of atoms
        half_life_years: half-life in years (default: C-14 = 5730)
        time_years: elapsed time in years (default: 0)
    """
    import math
    remaining = initial_atoms * (0.5 ** (time_years / half_life_years))
    return remaining

# Different usage patterns
n0 = 1000000  # One million atoms

# Use all defaults (returns initial count)
print(f"At t=0: {simulate_decay(n0):.0f} atoms")

# Override just time (uses default half-life for C-14)
print(f"After 5730 years: {simulate_decay(n0, time_years=5730):.0f} atoms")

# Override all parameters (U-238 example)
print(f"U-238 after 1 billion years: {simulate_decay(n0, 4.5e9, 1e9):.0f} atoms")

# THE TRAP - DON'T DO THIS!
def add_measurement_buggy(value, data=[]):
    """BUGGY VERSION - default list is shared!"""
    data.append(value)
    return data

# Watch the disaster unfold
day1 = add_measurement_buggy(23.5)
print(f"Day 1: {day1}")

day2 = add_measurement_buggy(24.1)  # Surprise!
print(f"Day 2: {day2}")  # Contains BOTH days!

print(f"Same list? {day1 is day2}")  # True - it's the same object!

def add_measurement_fixed(value, data=None):
    """CORRECT VERSION - new list each time if needed."""
    if data is None:
        data = []  # Create new list when not provided
    data.append(value)
    return data

# Now it works correctly
day1 = add_measurement_fixed(23.5)
day2 = add_measurement_fixed(24.1)
print(f"Day 1: {day1}")
print(f"Day 2: {day2}")  # Separate lists!

# Can still provide existing list
combined = []
add_measurement_fixed(23.5, combined)
add_measurement_fixed(24.1, combined)
print(f"Combined: {combined}")

def calculate_statistics(*values, method='all'):
    """
    Calculate statistics for any number of values.
    
    *args collects positional arguments into a tuple.
    """
    if not values:
        return None
    
    n = len(values)
    mean = sum(values) / n
    
    if method == 'mean':
        return mean
    
    # Calculate variance and std dev
    variance = sum((x - mean)**2 for x in values) / n
    std_dev = variance ** 0.5
    
    if method == 'all':
        return {
            'n': n,
            'mean': mean,
            'std': std_dev,
            'min': min(values),
            'max': max(values)
        }

# Works with any number of arguments!
print(calculate_statistics(1, 2, 3, 4, 5))
print(f"Mean only: {calculate_statistics(2.5, 3.7, 1.2, method='mean'):.2f}")

def run_experiment(name, **parameters):
    """
    Run experiment with flexible parameters.
    
    **kwargs collects keyword arguments into a dictionary.
    """
    print(f"=== Experiment: {name} ===")
    
    # Default parameters in CGS units
    defaults = {
        'temperature_k': 293.15,       # Kelvin
        'pressure_dynes_cm2': 1.01325e6,  # 1 atm in dynes/cm²
        'trials': 10
    }
    
    # Update with user parameters
    config = defaults.copy()
    config.update(parameters)
    
    print("Configuration:")
    for key, value in config.items():
        print(f"  {key}: {value}")
    
    return f"Completed {config['trials']} trials"

# Simple call
result = run_experiment("Test A")

# Complex call with many options
result = run_experiment(
    "Test B",
    temperature_k=350,
    pressure_dynes_cm2=2e6,
    trials=50,
    catalyst="Platinum"
)

# WRONG - SyntaxError!
# def bad_function(x=10, y):  # Can't have required after optional
#     pass

# CORRECT - Required parameters first, then optional
def good_function(x, y=10):
    """Required 'x' must come before optional 'y'."""
    return x + y

# Usage
print(good_function(5))      # x=5, y=10 (default)
print(good_function(5, 20))  # x=5, y=20 (override)

# Example 1: Mixing positional-only and keyword-only (Python 3.8+)
def divide(numerator, denominator, /, *, precision=2):
    """
    numerator, denominator: must be positional
    precision: must be keyword
    """
    result = numerator / denominator
    return round(result, precision)

print(divide(10, 3))                    # Uses default precision
print(divide(10, 3, precision=4))       # Specifies precision
# divide(numerator=10, denominator=3)   # ERROR! Must be positional

# Example 2: Flexible argument collection
def process(*data_points, normalize=True, **options):
    """
    data_points: any number of positional args
    normalize: keyword-only parameter
    options: any extra keyword arguments
    """
    print(f"Processing {len(data_points)} points")
    print(f"Normalize: {normalize}")
    print(f"Options: {options}")

process(1, 2, 3, 4, normalize=False, threshold=0.5, verbose=True)

# Without keyword-only: unclear what True and False mean
def process_data_unclear(data, normalize, remove_outliers):
    """What do the boolean parameters mean at the call site?"""
    pass

# Call is confusing:
# process_data_unclear(measurements, True, False)  # What's True? What's False?

# With keyword-only: self-documenting
def process_data_clear(data, *, normalize=True, remove_outliers=False):
    """The * forces callers to name the boolean parameters."""
    if normalize:
        # Normalize the data...
        pass
    if remove_outliers:
        # Remove outliers...
        pass

# Call is clear:
# process_data_clear(measurements, normalize=True, remove_outliers=False)

# Pattern 1: Simple function with defaults
def calculate_interest(principal, rate=0.05, years=1):
    """Most common pattern - some required, some optional."""
    return principal * (1 + rate) ** years

# Pattern 2: Collecting variable arguments
def average(*numbers):
    """Accept any number of arguments."""
    if not numbers:
        return 0
    return sum(numbers) / len(numbers)

print(f"Average of 2, 4, 6, 8: {average(2, 4, 6, 8)}")

# Pattern 3: Force clarity with keyword-only
def create_particle(mass, velocity, *, charge=0, spin=0.5):
    """Mass and velocity are required, charge and spin must be named."""
    return {
        'mass': mass,
        'velocity': velocity, 
        'charge': charge,
        'spin': spin
    }

# Clear what each parameter means:
electron = create_particle(9.109e-28, 1e6, charge=-1.602e-19, spin=0.5)

# Pattern 4: Flexible keyword arguments
def plot(x, y, **options):
    """Accept x, y, and any number of plot options."""
    print(f"Plotting {len(x)} points")
    for key, value in options.items():
        print(f"  {key}: {value}")

plot([1, 2, 3], [2, 4, 6], color='red', linewidth=2, marker='o')

# The / marker (Python 3.8+): everything before is positional-only
def divide(dividend, divisor, /):
    """Both arguments MUST be passed positionally."""
    return dividend / divisor

print(divide(10, 2))           # Works
# print(divide(dividend=10, divisor=2))  # ERROR! Can't use keywords

# The * marker: everything after is keyword-only  
def configure(name, *, debug=False, verbose=False):
    """name can be positional or keyword, others must be keyword."""
    print(f"Configuring {name}")
    if debug:
        print("  Debug mode enabled")
    if verbose:
        print("  Verbose output enabled")

configure("MyApp", debug=True)  # Must name debug
# configure("MyApp", True)  # ERROR! Can't pass debug positionally

# Demonstrating LEGB with scientific context
speed_of_light = 2.998e10  # Global scope (cm/s in CGS)

def calculate_energy():
    # Enclosing scope
    rest_mass = 9.109e-28  # electron mass in grams
    
    def relativistic_factor(velocity_cms):
        # Local scope
        c = speed_of_light  # Accesses global (cm/s)
        beta = velocity_cms / c
        
        # Local calculation
        import math
        if beta >= 1:
            return float('inf')  # Cannot exceed speed of light
        gamma = 1 / math.sqrt(1 - beta**2)
        
        print(f"  Inside relativistic: β = {beta:.3f}, γ = {gamma:.3f}")
        return gamma
    
    # Use inner function
    v = 0.5 * speed_of_light  # Half light speed (cm/s)
    gamma = relativistic_factor(v)
    
    # Calculate relativistic energy (E = γmc²)
    energy_ergs = rest_mass * speed_of_light**2 * gamma
    print(f"Inside calculate: E = {energy_ergs:.2e} ergs")
    return energy_ergs

result = calculate_energy()
print(f"Global scope: c = {speed_of_light:.2e} cm/s")
print(f"Final energy: {result:.2e} ergs = {result/1.602e-12:.2f} eV")

counter = 0  # Global

def increment_wrong():
    """This will crash - Python sees assignment and assumes local!"""
    # counter += 1  # UnboundLocalError!
    print("Uncomment the line above to see the error")

def increment_with_global():
    """Explicitly use global variable."""
    global counter
    counter += 1
    return counter

def increment_better(current_count):
    """Best approach - avoid global state entirely."""
    return current_count + 1

# Test the better approach
count = 0
count = increment_better(count)
count = increment_better(count)
print(f"Count: {count}")

# Global approach (generally avoid)
increment_with_global()
increment_with_global()
print(f"Global counter: {counter}")

def create_integrator(method='rectangle'):
    """
    Factory function that creates specialized integrators.
    
    The returned function 'remembers' the method.
    """
    def integrate(func, a, b, n=100):
        """Numerically integrate func from a to b."""
        dx = (b - a) / n
        total = 0
        
        if method == 'rectangle':
            # Left rectangle rule
            for i in range(n):
                x = a + i * dx
                total += func(x) * dx
                
        elif method == 'midpoint':
            # Midpoint rule (more accurate)
            for i in range(n):
                x = a + (i + 0.5) * dx
                total += func(x) * dx
        
        return total
    
    # Return the customized function
    return integrate

# Create specialized integrators
rect_integrate = create_integrator('rectangle')
mid_integrate = create_integrator('midpoint')

# Test with x² from 0 to 1 (exact answer = 1/3)
def f(x):
    return x**2

rect_result = rect_integrate(f, 0, 1, 100)
mid_result = mid_integrate(f, 0, 1, 100)
exact = 1/3

print(f"Rectangle rule: {rect_result:.6f}")
print(f"Midpoint rule:  {mid_result:.6f}")
print(f"Exact:          {exact:.6f}")
print(f"Rectangle error: {abs(rect_result-exact):.6f}")
print(f"Midpoint error:  {abs(mid_result-exact):.6f}")

# Traditional function
def square(x):
    return x**2

# Equivalent lambda
square_lambda = lambda x: x**2

print(f"Traditional: {square(5)}")
print(f"Lambda: {square_lambda(5)}")

# Lambdas shine for sorting and filtering
data_points = [
    {'x': 1.5, 'y': 2.3, 'error': 0.1},
    {'x': 2.1, 'y': 4.7, 'error': 0.05},
    {'x': 0.8, 'y': 1.2, 'error': 0.2},
]

# Sort by x value
by_x = sorted(data_points, key=lambda p: p['x'])
print("\nSorted by x:")
for p in by_x:
    print(f"  x={p['x']}, y={p['y']}")

# Sort by relative error (error/y)
by_rel_error = sorted(data_points, key=lambda p: p['error']/p['y'])
print("\nSorted by relative error:")
for p in by_rel_error:
    rel_err = p['error']/p['y'] * 100
    print(f"  x={p['x']}, relative error={rel_err:.1f}%")

from functools import reduce

# Sample measurements in CGS units
measurements_cms = [981, 979, 983, 9900, 980, 982, 978]  # cm/s² (g values)

# FILTER: Remove outliers
mean_estimate = 980  # cm/s²
tolerance = 10  # cm/s²

good_data = list(filter(
    lambda x: abs(x - mean_estimate) < tolerance,
    measurements_cms
))

print(f"Original: {measurements_cms}")
print(f"Filtered: {good_data}")
print(f"Removed {len(measurements_cms) - len(good_data)} outliers")

# MAP: Apply calibration factor
calibration = 1.002  # 0.2% correction
calibrated = list(map(lambda x: x * calibration, good_data))
print(f"Calibrated: {[f'{x:.1f}' for x in calibrated]}")

# REDUCE: Calculate mean
mean = reduce(lambda a, b: a + b, calibrated) / len(calibrated)
print(f"Final mean: {mean:.1f} cm/s²")

# Note: These functional patterns are essential for JAX
# Pure functions enable automatic differentiation
# Map operations parallelize on GPUs
# Immutable data prevents side effects

def apply_operator(data, operator_func):
    """
    Apply any operation to data.
    
    operator_func determines the transformation.
    """
    return [operator_func(x) for x in data]

# Define different operations
def grams_to_kilograms(g):
    """Convert grams to kilograms."""
    return g / 1000.0

def dynes_to_newtons(dynes):
    """Convert dynes to Newtons (1 N = 10⁵ dynes)."""
    return dynes / 1e5

def normalize(x, mean=0, std=1):
    """Standardize value."""
    return (x - mean) / std

# Use different operations
masses_g = [100, 200, 300, 400]  # grams
masses_kg = apply_operator(masses_g, grams_to_kilograms)
print(f"Grams: {masses_g}")
print(f"Kilograms: {masses_kg}")

forces_dynes = [1e5, 2e5, 3e5, 4e5]  # dynes
forces_n = apply_operator(forces_dynes, dynes_to_newtons)
print(f"\nDynes: {forces_dynes}")
print(f"Newtons: {forces_n}")

# Create custom normalizer with closure
values = [2, 4, 6, 8, 10]
mean = sum(values) / len(values)
std = (sum((x - mean)**2 for x in values) / len(values)) ** 0.5

normalizer = lambda x: normalize(x, mean, std)
normalized = apply_operator(values, normalizer)
print(f"\nOriginal: {values}")
print(f"Normalized: {[f'{x:.2f}' for x in normalized]}")

# Method 1: Show the module contents (you would save this in a file)
module_code = '''
"""
statistics_tools.py - Basic statistical calculations.

This module provides fundamental statistical functions.
"""

import math

def mean(data):
    """Calculate arithmetic mean."""
    return sum(data) / len(data)

def variance(data):
    """Calculate population variance."""
    m = mean(data)
    return sum((x - m)**2 for x in data) / len(data)

def std_dev(data):
    """Calculate population standard deviation."""
    return math.sqrt(variance(data))

def std_error(data):
    """Calculate standard error of the mean."""
    return std_dev(data) / math.sqrt(len(data))

# Module-level code runs on import
print("Loading statistics_tools module...")

# Test code that only runs when script is executed directly
if __name__ == "__main__":
    test_data = [1, 2, 3, 4, 5]
    print(f"Test mean: {mean(test_data)}")
    print(f"Test std dev: {std_dev(test_data):.3f}")
'''

# Method 2: Programmatically create the file
with open('statistics_tools.py', 'w') as f:
    f.write(module_code)

print("Module file 'statistics_tools.py' has been created in your current directory")

# Different import methods

# Method 1: Import entire module
import statistics_tools

data = [10, 12, 11, 13, 10, 11, 12]
m = statistics_tools.mean(data)
s = statistics_tools.std_dev(data)
print(f"Method 1 - Mean: {m:.2f}, Std: {s:.2f}")

# Method 2: Import specific functions
from statistics_tools import mean, std_error

se = std_error(data)
print(f"Method 2 - Standard error: {se:.3f}")

# Method 3: Import with alias
import statistics_tools as stats

var = stats.variance(data)
print(f"Method 3 - Variance: {var:.3f}")

# Create physics_cgs.py module with all CGS units
physics_module = '''
"""
physics_cgs.py - Physics calculations in CGS units.

All calculations use CGS (centimeter-gram-second) units.
"""

# Physical constants in CGS
SPEED_OF_LIGHT = 2.998e10      # cm/s
PLANCK_CONSTANT = 6.626e-27    # erg⋅s  
BOLTZMANN = 1.381e-16          # erg/K
ELECTRON_MASS = 9.109e-28      # g
PROTON_MASS = 1.673e-24        # g
GRAVITATIONAL_CONSTANT = 6.674e-8  # cm³/(g⋅s²)

def kinetic_energy(mass_g, velocity_cms):
    """KE = ½mv² in ergs."""
    return 0.5 * mass_g * velocity_cms**2

def momentum(mass_g, velocity_cms):
    """p = mv in g⋅cm/s."""
    return mass_g * velocity_cms

def de_broglie_wavelength(mass_g, velocity_cms):
    """λ = h/p in cm."""
    p = momentum(mass_g, velocity_cms)
    return PLANCK_CONSTANT / p if p != 0 else float('inf')

def thermal_velocity(temp_k, mass_g):
    """RMS thermal velocity in cm/s."""
    import math
    return math.sqrt(3 * BOLTZMANN * temp_k / mass_g)

def photon_energy(wavelength_cm):
    """E = hc/λ in ergs."""
    return PLANCK_CONSTANT * SPEED_OF_LIGHT / wavelength_cm

def gravitational_force(m1_g, m2_g, distance_cm):
    """F = Gm₁m₂/r² in dynes."""
    if distance_cm == 0:
        return float('inf')
    return GRAVITATIONAL_CONSTANT * m1_g * m2_g / distance_cm**2

class Particle:
    """Simple particle with physics methods."""
    
    def __init__(self, mass_g, velocity_cms):
        self.mass = mass_g
        self.velocity = velocity_cms
    
    def kinetic_energy(self):
        return kinetic_energy(self.mass, self.velocity)
    
    def wavelength(self):
        return de_broglie_wavelength(self.mass, self.velocity)
'''

with open('physics_cgs.py', 'w') as f:
    f.write(physics_module)

# Now use it
import physics_cgs

# Electron at 1% light speed (all CGS units)
v_electron = 0.01 * physics_cgs.SPEED_OF_LIGHT  # cm/s
ke = physics_cgs.kinetic_energy(physics_cgs.ELECTRON_MASS, v_electron)
wavelength = physics_cgs.de_broglie_wavelength(physics_cgs.ELECTRON_MASS, v_electron)

print(f"Electron at 1% c:")
print(f"  Velocity = {v_electron:.2e} cm/s")
print(f"  KE = {ke:.2e} ergs")
print(f"  de Broglie λ = {wavelength:.2e} cm")

# Test gravitational force (Earth-Moon in CGS)
earth_mass_g = 5.972e27  # grams
moon_mass_g = 7.342e25   # grams  
earth_moon_distance_cm = 3.844e10  # cm

force = physics_cgs.gravitational_force(earth_mass_g, moon_mass_g, earth_moon_distance_cm)
print(f"\nEarth-Moon gravitational force: {force:.2e} dynes")

# GOOD: Clear, explicit imports
import numpy as np
import matplotlib.pyplot as plt
from math import sin, cos, pi

# BAD: Wildcard imports pollute namespace
# from numpy import *  # Adds 600+ names!
# from math import *   # Conflicts with numpy!

# Example of namespace collision
import math
import cmath  # Complex math module (introduced in Chapter 2)

# Clear which log we're using
real_log = math.log(10)      # Natural log of real number
complex_log = cmath.log(-1)   # Can handle complex numbers

print(f"math.log(10) = {real_log:.3f}")
print(f"cmath.log(-1) = {complex_log}")  # Returns complex number

# If we had used 'from math import *':
# log(10)  # Which log? math or cmath? Unclear!

# Level 1: Minimal (but essential)
def quadratic_formula(a, b, c):
    """Solve ax² + bx + c = 0."""
    discriminant = b**2 - 4*a*c
    if discriminant < 0:
        return None, None
    sqrt_disc = discriminant**0.5
    return (-b + sqrt_disc)/(2*a), (-b - sqrt_disc)/(2*a)

# Level 2: Professional documentation
def integrate_trapezoid(x_values, y_values):
    """
    Integrate using the trapezoidal rule.
    
    Parameters
    ----------
    x_values : list or array
        X coordinates of data points (must be sorted)
    y_values : list or array
        Y coordinates of data points
    
    Returns
    -------
    float
        Approximate integral of y with respect to x
    
    Examples
    --------
    >>> x = [0, 1, 2, 3]
    >>> y = [0, 1, 4, 9]  # y = x²
    >>> integrate_trapezoid(x, y)
    8.0  # Exact answer is 9
    
    Notes
    -----
    Uses the composite trapezoidal rule:
    ∫y dx ≈ Σ(y[i] + y[i+1])/2 * (x[i+1] - x[i])
    
    More accurate for smooth functions with
    many points.
    """
    if len(x_values) != len(y_values):
        raise ValueError("x and y must have same length")
    
    integral = 0
    for i in range(len(x_values) - 1):
        dx = x_values[i + 1] - x_values[i]
        avg_y = (y_values[i] + y_values[i + 1]) / 2
        integral += avg_y * dx
    
    return integral

# Test the function
x = [0, 0.5, 1.0, 1.5, 2.0]
y = [0, 0.25, 1.0, 2.25, 4.0]  # y = x²
result = integrate_trapezoid(x, y)
print(f"Integral of x² from 0 to 2: {result:.3f}")
print(f"Exact answer: {8/3:.3f}")
print(f"Error: {abs(result - 8/3):.3f}")

def quadratic_formula(a, b, c):
    """
    Solve quadratic equation ax² + bx + c = 0.
    
    Returns:
        tuple: (root1, root2) or (None, None) if no real roots
    """
    import math
    
    # Calculate discriminant
    discriminant = b**2 - 4*a*c
    
    # Check if real roots exist
    if discriminant < 0:
        return None, None
    elif discriminant == 0:
        # One repeated root
        root = -b / (2*a)
        return root, root
    else:
        # Two distinct roots
        sqrt_disc = math.sqrt(discriminant)
        root1 = (-b + sqrt_disc) / (2*a)
        root2 = (-b - sqrt_disc) / (2*a)
        return root1, root2

def test_quadratic():
    """Test quadratic formula solver."""
    
    # Test 1: Simple case (x² - 5x + 6 = 0)
    # Roots should be 2 and 3
    r1, r2 = quadratic_formula(1, -5, 6)
    assert abs(r1 - 3) < 1e-10 or abs(r1 - 2) < 1e-10
    assert abs(r2 - 3) < 1e-10 or abs(r2 - 2) < 1e-10
    print("✓ Test 1: Simple roots")
    
    # Test 2: No real roots (x² + 1 = 0)
    r1, r2 = quadratic_formula(1, 0, 1)
    assert r1 is None and r2 is None
    print("✓ Test 2: No real roots")
    
    # Test 3: Single root (x² - 2x + 1 = 0)
    # Root should be 1 (twice)
    r1, r2 = quadratic_formula(1, -2, 1)
    assert abs(r1 - 1) < 1e-10
    assert abs(r2 - 1) < 1e-10
    print("✓ Test 3: Repeated root")
    
    print("All tests passed!")

# Run tests
test_quadratic()

def kinetic_energy_typed(mass_g: float, velocity_cms: float) -> float:
    """
    Calculate kinetic energy with type hints.
    
    The ': float' annotations indicate expected types.
    The '-> float' indicates the return type.
    These are documentation only - Python doesn't enforce them!
    """
    return 0.5 * mass_g * velocity_cms**2

# Type hints help IDEs provide better autocomplete and catch potential errors
energy: float = kinetic_energy_typed(9.109e-28, 3e9)
print(f"Energy: {energy:.2e} ergs")

# Python still allows "wrong" types - hints aren't enforced!
# This works even though we pass integers instead of floats:
energy2 = kinetic_energy_typed(1, 100)  # Still works!
print(f"Integer inputs still work: {energy2} ergs")

# For more complex types, import from typing module
from typing import List, Tuple, Optional

def analyze_data_typed(values: List[float]) -> Tuple[float, float, Optional[float]]:
    """
    Analyze data with complex type hints.
    
    List[float] means a list of floats
    Tuple[float, float, Optional[float]] means it returns
    a tuple with two floats and maybe a third float (or None)
    """
    if not values:
        return 0.0, 0.0, None
    
    n = len(values)
    mean = sum(values) / n
    variance = sum((x - mean)**2 for x in values) / n
    std_dev = variance ** 0.5
    
    sorted_vals = sorted(values)
    median = sorted_vals[n // 2] if n % 2 == 1 else (sorted_vals[n//2-1] + sorted_vals[n//2]) / 2
    
    return mean, std_dev, median

# The IDE now knows exactly what types are returned
mean_val, std_val, median_val = analyze_data_typed([1.0, 2.0, 3.0, 4.0, 5.0])
print(f"Mean: {mean_val:.1f}, Std: {std_val:.2f}, Median: {median_val:.1f}")

import time

def simple_add(a, b):
    """Minimal work - overhead visible."""
    return a + b

def complex_calc(x):
    """Heavy computation - overhead negligible."""
    result = 0
    for i in range(100):
        result += (x + i) ** 0.5
    return result / 100

# Measure overhead
n_calls = 100000

# Simple function
start = time.perf_counter()
for _ in range(n_calls):
    simple_add(1.0, 2.0)
simple_time = time.perf_counter() - start

# Inline equivalent
start = time.perf_counter()
for _ in range(n_calls):
    result = 1.0 + 2.0
inline_time = time.perf_counter() - start

# Complex function
start = time.perf_counter()
for _ in range(n_calls // 100):  # Fewer calls (it's slower)
    complex_calc(1.0)
complex_time = time.perf_counter() - start

print("Function Call Analysis:")
print(f"Simple function: {simple_time*1000:.1f} ms")
print(f"Inline addition: {inline_time*1000:.1f} ms")
print(f"Overhead factor: {simple_time/inline_time:.1f}x")
print(f"\nComplex function: {complex_time*1000:.1f} ms")
print("\nLesson: Overhead only matters for trivial operations!")

from functools import lru_cache

# Fibonacci without memoization (exponentially slow!)
def fib_slow(n):
    if n < 2:
        return n
    return fib_slow(n-1) + fib_slow(n-2)

# Fibonacci with memoization (linear time!)
@lru_cache(maxsize=128)
def fib_fast(n):
    if n < 2:
        return n
    return fib_fast(n-1) + fib_fast(n-2)

# Add astrophysics-specific memoization example
@lru_cache(maxsize=128)
def planck_function(wavelength_nm, temp_k):
    """
    Cached Planck function for blackbody radiation.
    Expensive calculation cached for repeated calls.
    Common in spectral fitting where same temps recur.
    """
    import math
    h = 6.626e-27  # Planck constant in erg⋅s
    c = 2.998e10   # Speed of light in cm/s
    k = 1.381e-16  # Boltzmann constant in erg/K
    
    wavelength_cm = wavelength_nm * 1e-7
    
    # Planck function calculation
    exp_term = (h * c) / (wavelength_cm * k * temp_k)
    if exp_term > 700:  # Prevent overflow
        return 0
    
    numerator = 2 * h * c**2 / wavelength_cm**5
    denominator = math.exp(exp_term) - 1
    return numerator / denominator

# Compare performance
import time

n = 30

start = time.perf_counter()
result_slow = fib_slow(n)
time_slow = time.perf_counter() - start

start = time.perf_counter()
result_fast = fib_fast(n)
time_fast = time.perf_counter() - start

print(f"Fibonacci({n}) = {result_fast}")
print(f"Without memoization: {time_slow*1000:.1f} ms")
print(f"With memoization: {time_fast*1000:.3f} ms")
print(f"Speedup: {time_slow/time_fast:.0f}x")

# Check cache statistics
print(f"\nCache info: {fib_fast.cache_info()}")

# Test Planck function caching
print(f"\nPlanck function at 500nm, 5778K: {planck_function(500, 5778):.2e} erg/(s⋅cm²⋅cm⋅sr)")
print(f"Cache info: {planck_function.cache_info()}")