Chapter 6: OOP Fundamentals - Organizing Scientific Code

# Can you predict what this prints?
def modify(data):
    data['value'] = data['value'] * 2
    return data

measurement = {'value': 100, 'unit': 'K'}
result = modify(measurement)
print(f"Original: {measurement['value']}")  # What value?
print(f"Result: {result['value']}")         # What value?

# Approach 1: Functions with dictionaries (what you know)
def create_particle(mass, x, y, vx, vy):
    """Create a particle dictionary."""
    return {
        'mass': mass,  # grams
        'x': x, 'y': y,  # cm
        'vx': vx, 'vy': vy  # cm/s
    }

def kinetic_energy(particle):
    """Calculate kinetic energy in ergs."""
    v_squared = particle['vx']**2 + particle['vy']**2
    return 0.5 * particle['mass'] * v_squared

def update_position(particle, dt):
    """Update position based on velocity."""
    particle['x'] += particle['vx'] * dt
    particle['y'] += particle['vy'] * dt
    return particle

# Using functions
p1 = create_particle(1.0, 0, 0, 10, 5)
energy = kinetic_energy(p1)
p1 = update_position(p1, 0.1)
print(f"Energy: {energy:.1f} ergs")

# Approach 2: Object-Oriented (what you're learning)
class Particle:
    """A particle with position and velocity."""
    
    def __init__(self, mass, x, y, vx, vy):
        """Initialize particle state."""
        self.mass = mass  # grams
        self.x = x        # cm
        self.y = y        # cm
        self.vx = vx      # cm/s
        self.vy = vy      # cm/s
    
    def kinetic_energy(self):
        """Calculate kinetic energy in ergs."""
        v_squared = self.vx**2 + self.vy**2
        return 0.5 * self.mass * v_squared
    
    def update_position(self, dt):
        """Update position based on velocity."""
        self.x += self.vx * dt
        self.y += self.vy * dt

# Using objects
p2 = Particle(1.0, 0, 0, 10, 5)
energy = p2.kinetic_energy()
p2.update_position(0.1)
print(f"Energy: {energy:.1f} ergs")

# Define a class (blueprint)
class Measurement:
    """A scientific measurement with uncertainty."""
    
    def __init__(self, value, error):
        """Initialize measurement with value and error."""
        self.value = value
        self.error = error
    
    def relative_error(self):
        """Calculate relative error as percentage."""
        if self.value == 0:
            return float('inf')
        return abs(self.error / self.value) * 100

# Create objects (instances)
temp = Measurement(293.15, 0.1)  # Temperature in Kelvin
# Pressure in dyne/cm² (CGS unit for pressure)
# Note: CGS (centimeter-gram-second) is traditional in astronomy
pressure = Measurement(1.01325e6, 500)  # 1 atm = 1.01325e6 dyne/cm²

print(f"Temperature: {temp.value} ± {temp.error} K")
print(f"Pressure: {pressure.value:.2e} ± {pressure.error} dyne/cm²")
print(f"Pressure relative error: {pressure.relative_error():.3f}%")

class Counter:
    """Demonstrates how self works."""
    
    def __init__(self):
        """Initialize counter at zero."""
        self.count = 0  # Each object gets its own count
    
    def increment(self):
        """Increment this counter by one."""
        self.count += 1  # self refers to the specific object
    
    def get_count(self):
        """Return current count value."""
        return self.count

# Each object maintains independent state
c1 = Counter()
c2 = Counter()

c1.increment()
c1.increment()
c2.increment()

print(f"Counter 1: {c1.get_count()}")  # 2
print(f"Counter 2: {c2.get_count()}")  # 1

# What Python actually does when you call c1.increment():
# Counter.increment(c1)  # c1 becomes 'self' in the method

class Simulation:
    """Demonstrates instance vs class attributes."""
    
    # Class attribute - shared by all simulations
    speed_of_light = 2.998e10  # cm/s (CGS)
    total_runs = 0
    
    def __init__(self, name, particle_count):
        """Initialize a new simulation."""
        # Instance attributes - unique to each simulation
        self.name = name
        self.particle_count = particle_count
        self.time = 0.0
        self.box_size_cm = 1.0e-4  # 1 micron box
        # Increment shared counter
        Simulation.total_runs += 1
    
    def advance(self, dt):
        """Advance simulation by dt seconds."""
        self.time += dt

# Create simulations
sim1 = Simulation("Test A", 1000)
sim2 = Simulation("Test B", 5000)

print(f"Total simulations: {Simulation.total_runs}")
print(f"Sim1 particles: {sim1.particle_count}")
print(f"Sim2 particles: {sim2.particle_count}")
print(f"Speed of light: {Simulation.speed_of_light:.2e} cm/s")
print(f"Box size: {sim1.box_size_cm:.2e} cm")

# Note: In production code, always import at file top:
# import math
# Here we import locally for pedagogical clarity

import math  # Best practice: import at the top

# First: A simple class with basic methods
class Sample:
    """Scientific sample with basic operations."""
    
    def __init__(self, mass_g, volume_cm3):
        """Initialize with mass in grams, volume in cm³."""
        self.mass = mass_g
        self.volume = volume_cm3
    
    def density(self):
        """Calculate density in g/cm³."""
        return self.mass / self.volume
    
    def is_denser_than_water(self):
        """Check if denser than water (1 g/cm³)."""
        return self.density() > 1.0

# Using simple methods
iron = Sample(7.87, 1.0)  # Iron sample
print(f"Iron density: {iron.density()} g/cm³")
print(f"Sinks in water: {iron.is_denser_than_water()}")

class Vector2D:
    """A 2D vector for physics calculations."""
    
    def __init__(self, x, y):
        """Initialize vector components in cm."""
        self.x = x
        self.y = y
    
    def magnitude(self):
        """Calculate vector magnitude (length)."""
        return (self.x**2 + self.y**2)**0.5
    
    def normalize(self):
        """Scale vector to unit length (magnitude = 1)."""
        mag = self.magnitude()
        if mag > 0:
            self.x /= mag
            self.y /= mag
    
    def dot(self, other):
        """Calculate dot product with another vector."""
        return self.x * other.x + self.y * other.y
    
    def angle_with(self, other):
        """Calculate angle with another vector in radians."""
        dot_product = self.dot(other)
        mags = self.magnitude() * other.magnitude()
        if mags == 0:
            return 0
        cos_angle = dot_product / mags
        # Clamp to [-1, 1] to avoid numerical errors
        cos_angle = max(-1, min(1, cos_angle))
        return math.acos(cos_angle)

# Using vector methods
v1 = Vector2D(3, 4)
v2 = Vector2D(1, 0)

print(f"v1 magnitude: {v1.magnitude()} cm")
print(f"Dot product: {v1.dot(v2)} cm²")
print(f"Angle: {v1.angle_with(v2):.2f} radians")

v1.normalize()  # Makes magnitude = 1
print(f"After normalization: ({v1.x:.2f}, {v1.y:.2f})")
print(f"New magnitude: {v1.magnitude():.2f}")

class DataProcessor:
    """Demonstrates different method types."""
    
    version = "1.0"
    
    def __init__(self, data):
        """Initialize with data to process."""
        self.data = data
    
    # Instance method - needs self
    def process(self):
        """Process this object's data."""
        return sum(self.data) / len(self.data)
    
    # Class method - gets class, not instance
    @classmethod
    def from_file(cls, filename):
        """Create instance from file."""
        # Simulate file reading
        data = [1, 2, 3, 4, 5]
        return cls(data)  # Create new instance
    
    # Static method - doesn't need self or cls
    @staticmethod
    def validate_data(data):
        """Check if data is valid."""
        return len(data) > 0 and all(isinstance(x, (int, float)) for x in data)
    
    # Static method with physics example
    @staticmethod
    def celsius_to_kelvin(celsius):
        """Convert temperature - no instance needed."""
        return celsius + 273.15

# Using different method types
processor = DataProcessor([10, 20, 30])
print(f"Average: {processor.process()}")

# Class method creates new instance
# processor2 = DataProcessor.from_file("data.txt")

# Static methods work without instance
valid = DataProcessor.validate_data([1, 2, 3])
print(f"Data valid: {valid}")

# Physics static method
temp_k = DataProcessor.celsius_to_kelvin(25)
print(f"Temperature: {temp_k} K")

class Circle:
    """Circle with computed properties."""
    
    def __init__(self, radius):
        """Initialize with radius, using validation."""
        # This now uses the setter for validation!
        self.radius = radius
    
    @property
    def radius(self):
        """Get radius in cm."""
        return self._radius
    
    @radius.setter
    def radius(self, value):
        """Set radius with validation."""
        if value <= 0:
            raise ValueError(f"Radius must be positive, got {value}")
        self._radius = value
    
    @property
    def area(self):
        """Computed area in cm² using πr²."""
        return math.pi * self._radius**2
    
    @property
    def circumference(self):
        """Computed circumference in cm using 2πr."""
        return 2 * math.pi * self._radius

# Properties look like attributes but run code
circle = Circle(5)
print(f"Radius: {circle.radius} cm")
print(f"Area: {circle.area:.2f} cm²")
print(f"Circumference: {circle.circumference:.2f} cm")

# Changing radius automatically updates computed properties
circle.radius = 10
print(f"New area: {circle.area:.2f} cm²")

# Validation prevents invalid states (even in __init__ now!)
try:
    bad_circle = Circle(-5)  # This will raise an error
except ValueError as e:
    print(f"Error: {e}")

class Temperature:
    """Temperature with automatic unit conversion."""
    
    def __init__(self, kelvin):
        """Initialize with temperature in Kelvin."""
        self._kelvin = kelvin
    
    @property
    def kelvin(self):
        """Get temperature in Kelvin."""
        return self._kelvin
    
    @kelvin.setter
    def kelvin(self, value):
        """Set temperature with validation."""
        if value < 0:
            raise ValueError("Below absolute zero!")
        self._kelvin = value
    
    @property
    def celsius(self):
        """Get temperature in Celsius."""
        return self._kelvin - 273.15
    
    @celsius.setter
    def celsius(self, value):
        """Set temperature in Celsius."""
        # Use kelvin setter for validation
        self.kelvin = value + 273.15
    
    @property
    def fahrenheit(self):
        """Get temperature in Fahrenheit."""
        return self._kelvin * 9/5 - 459.67

# Same temperature, different units
temp = Temperature(300)
print(f"Water at {temp.kelvin:.1f} K")
print(f"  = {temp.celsius:.1f} °C")
print(f"  = {temp.fahrenheit:.1f} °F")

# Set in any unit, always consistent
temp.celsius = 100
print(f"Boiling: {temp.kelvin:.1f} K")

class Observation:
    """Astronomical observation with computed conditions."""
    
    def __init__(self, target, altitude_deg, moon_days_since_new=0):
        """Initialize observation with target and conditions."""
        self.target = target
        self.altitude_deg = altitude_deg  # Object altitude in degrees
        self.moon_days = moon_days_since_new
    
    @property
    def airmass(self):
        """Calculate airmass from altitude (simplified Pickering 2002)."""
        import math
        if self.altitude_deg <= 0:
            return float('inf')
        # Convert to zenith angle
        zenith_deg = 90 - self.altitude_deg
        zenith_rad = math.radians(zenith_deg)
        # Simplified plane-parallel approximation
        # Real calculation uses Pickering's formula
        return 1.0 / math.cos(zenith_rad) if zenith_deg < 85 else 38.0
    
    @property
    def moon_phase(self):
        """Calculate moon illumination fraction (0=new, 1=full)."""
        import math
        # Simplified: actual needs ephemeris
        # Moon cycle is ~29.53 days
        phase_angle = 2 * math.pi * self.moon_days / 29.53
        # Illumination varies as (1 + cos(phase))/2
        return (1 + math.cos(phase_angle)) / 2
    
    @property
    def sky_brightness_category(self):
        """Categorize sky brightness based on moon."""
        if self.moon_phase < 0.25:
            return "Dark"
        elif self.moon_phase < 0.75:
            return "Grey"
        else:
            return "Bright"

# Example observation
obs = Observation("M31", altitude_deg=60, moon_days_since_new=7)
print(f"Observing {obs.target}:")
print(f"  Airmass: {obs.airmass:.2f}")
print(f"  Moon illumination: {obs.moon_phase:.1%}")
print(f"  Sky conditions: {obs.sky_brightness_category}")

class Fraction:
    """A fraction with arithmetic operations."""
    
    def __init__(self, numerator, denominator):
        """Initialize and reduce fraction."""
        if denominator == 0:
            raise ValueError("Denominator cannot be zero")
        self.num = numerator
        self.den = denominator
        self._reduce()
    
    def _reduce(self):
        """Reduce to lowest terms using GCD (Euclidean algorithm)."""
        # Euclidean algorithm: repeatedly replace larger number 
        # with remainder when divided by smaller number
        a, b = abs(self.num), abs(self.den)
        while b:  # Continue until remainder is 0
            a, b = b, a % b  # GCD(a,b) = GCD(b, a mod b)
        self.num //= a
        self.den //= a
    
    def __str__(self):
        """Human-readable string for print()."""
        return f"{self.num}/{self.den}"
    
    def __repr__(self):
        """Unambiguous string for debugging."""
        return f"Fraction({self.num}, {self.den})"
    
    def __float__(self):
        """Convert to float."""
        return self.num / self.den
    
    def __add__(self, other):
        """Add two fractions: a/b + c/d = (ad+bc)/bd"""
        new_num = self.num * other.den + other.num * self.den
        new_den = self.den * other.den
        return Fraction(new_num, new_den)
    
    def __eq__(self, other):
        """Check equality: a/b == c/d if ad == bc"""
        return self.num * other.den == other.num * self.den

# Using special methods
f1 = Fraction(1, 2)
f2 = Fraction(1, 3)

print(f"f1 = {f1}")  # Calls __str__
print(f"f1 + f2 = {f1 + f2}")  # Calls __add__
print(f"f1 as float: {float(f1)}")  # Calls __float__
print(f"f1 == Fraction(2,4): {f1 == Fraction(2, 4)}")  # Calls __eq__

class DataSet:
    """Collection that acts like a built-in container."""
    
    def __init__(self, values=None):
        """Initialize with optional values."""
        self.values = values if values else []
    
    def __len__(self):
        """Support len(dataset)."""
        return len(self.values)
    
    def __getitem__(self, index):
        """Support dataset[index]."""
        return self.values[index]
    
    def __setitem__(self, index, value):
        """Support dataset[index] = value."""
        self.values[index] = value
    
    def __contains__(self, value):
        """Support 'value in dataset'."""
        return value in self.values
    
    def __iter__(self):
        """Support for loops."""
        return iter(self.values)
    
    def __bool__(self):
        """Support if dataset: (True if not empty)."""
        return len(self.values) > 0

# Acts like a built-in container
data = DataSet([10, 20, 30, 40, 50])

print(f"Length: {len(data)}")
print(f"First: {data[0]}")
print(f"Contains 30: {30 in data}")
print(f"Is non-empty: {bool(data)}")

data[1] = 25
for value in data:
    print(value, end=" ")

class RunningStatistics:
    """Maintains statistics as data arrives."""
    
    def __init__(self):
        """Initialize empty statistics."""
        self.count = 0
        self.sum = 0
        self.sum_sq = 0
    
    def add_value(self, x):
        """Add new value to statistics."""
        self.count += 1
        self.sum += x
        self.sum_sq += x**2
    
    @property
    def mean(self):
        """Current mean."""
        return self.sum / self.count if self.count > 0 else 0
    
    @property
    def variance(self):
        """Current variance."""
        if self.count < 2:
            return 0
        mean = self.mean
        return (self.sum_sq - self.count * mean**2) / (self.count - 1)

# Object maintains state between calls
stats = RunningStatistics()
for value in [1, 2, 3, 4, 5]:
    stats.add_value(value)
    print(f"After {value}: mean={stats.mean:.1f}, var={stats.variance:.1f}")

class Galaxy:
    """Model a galaxy with properties."""
    
    def __init__(self, name, distance_mpc, redshift):
        """Initialize galaxy with observed properties."""
        self.name = name
        self.distance_mpc = distance_mpc  # Megaparsecs
        self.redshift = redshift
    
    def recession_velocity(self):
        """Calculate recession velocity in km/s."""
        c = 3e5  # km/s
        return c * self.redshift  # Simple approximation
    
    def lookback_time_gyr(self):
        """Estimate lookback time in Gyr."""
        # Simplified: t ≈ z / H0
        H0 = 70  # km/s/Mpc
        return self.redshift * 1000 / H0

m31 = Galaxy("Andromeda", 0.78, -0.001)  # Blueshifted!
print(f"{m31.name}: v={m31.recession_velocity():.0f} km/s")

# No need for a class here
def celsius_to_kelvin(celsius):
    """Convert Celsius to Kelvin."""
    return celsius + 273.15

def calculate_orbital_period(semi_major_axis_au):
    """Kepler's third law: P² = a³ for period in years."""
    return semi_major_axis_au ** 1.5

# Simple functions are clearer than unnecessary classes
temp_k = celsius_to_kelvin(25)
period = calculate_orbital_period(1.0)  # Earth
print(f"Temperature: {temp_k} K")
print(f"Period: {period} years")

# These don't need to remember anything between calls
def mean(values):
    """Calculate mean."""
    return sum(values) / len(values)

def standard_deviation(values):
    """Calculate standard deviation."""
    m = mean(values)
    variance = sum((x - m)**2 for x in values) / len(values)
    return variance**0.5

data = [1, 2, 3, 4, 5]
print(f"Mean: {mean(data)}")
print(f"Std: {standard_deviation(data):.2f}")

class Instrument:
    """Scientific instrument for debugging demo."""
    
    def __init__(self, name, wavelength_nm):
        """Initialize instrument with name and wavelength."""
        self.name = name
        self.wavelength_nm = wavelength_nm
        self._calibrated = False
    
    def calibrate(self):
        """Calibrate instrument."""
        self._calibrated = True
        return f"{self.name} calibrated"

# Create instrument
spectrometer = Instrument("HARPS", 500)

# Introspection tools
print(f"Type: {type(spectrometer)}")
print(f"Class name: {spectrometer.__class__.__name__}")
print(f"Is Instrument?: {isinstance(spectrometer, Instrument)}")

# Check attributes
print(f"\nHas 'calibrate'?: {hasattr(spectrometer, 'calibrate')}")
print(f"Wavelength: {getattr(spectrometer, 'wavelength_nm', 'N/A')}")

# List all attributes (filtering out special ones)
attrs = [a for a in dir(spectrometer) if not a.startswith('_')]
print(f"\nPublic attributes: {attrs}")

# Instance __dict__ vs dir()
# __dict__ shows instance attributes only
# dir() shows everything the object can access
print(f"\nInstance __dict__: {vars(spectrometer)}")
print(f"dir() has {len(dir(spectrometer))} items (includes inherited)")
print(f"__dict__ has {len(vars(spectrometer))} items (instance only)")

from typing import Optional, List

class ObservationLog:
    """Demonstrates type hints with classes."""
    
    def __init__(self, observer: str, target: str) -> None:
        """Initialize observation log with type hints."""
        self.observer: str = observer
        self.target: str = target
        self.entries: List[str] = []
        self.completion_time: Optional[float] = None
    
    def add_entry(self, entry: str) -> None:
        """Add log entry."""
        self.entries.append(entry)
    
    def get_entry_count(self) -> int:
        """Return number of entries."""
        return len(self.entries)
    
    def mark_complete(self, time_hours: float) -> None:
        """Mark observation complete."""
        self.completion_time = time_hours

# Type hints help IDEs and documentation
log = ObservationLog("J. Smith", "M31")
log.add_entry("Target acquired")
log.add_entry("Exposure started")
log.mark_complete(2.5)

print(f"Observer: {log.observer}")
print(f"Entries: {log.get_entry_count()}")
print(f"Duration: {log.completion_time} hours")