Chapter 3: Control Flow & Logic

Learning Objectives¶

By the end of this chapter, you will be able to:

• (1) Design algorithms using structured pseudocode before writing any Python code, progressing through three levels of refinement.

• (2) Implement conditional statements (if/elif/else) with proper guard clauses and edge case handling.

• (3) Choose appropriate loop structures (for vs while) based on problem requirements and iteration patterns.

• (4) Master all six comparison operators (>, <, >=, <=, ==, !=) and three logical operators (and, or, not).

• (5) Handle floating-point comparisons safely using math.isclose() instead of direct equality checks.

• (6) Debug logic errors systematically using IPython’s %debug magic command and strategic assert statements.

• (7) Write efficient list comprehensions for simple transformations while recognizing when explicit loops are clearer.

• (8) Recognize and apply five universal algorithmic patterns: accumulation, filtering, mapping, searching, and convergence.

• (9) Build defensive code with guard clauses that validate inputs and fail fast with clear error messages.

Prerequisites Check¶

Chapter Overview¶

Programming is fundamentally about teaching computers to make decisions and repeat tasks. When you write an if statement or a loop, you’re translating human logic into instructions a machine can follow. But here’s the critical insight that separates computational thinkers from mere coders: the logic must be designed before it’s implemented. This chapter transforms you from someone who writes code to someone who designs algorithms.

We’ll start with the lost art of pseudocode — not as a bureaucratic exercise, but as the difference between code that works by accident and code that works by design. You’ll learn to recognize universal patterns that appear across all of computational physics: iteration, accumulation, filtering, mapping, and reduction. These patterns will appear in every project you build, from N-body simulations to neural networks. Whether you’re folding light curves to find exoplanet periods or iterating until your stellar model converges, these patterns form the backbone of computational astronomy.

The control flow structures we explore here are where your numerical calculations from Chapter 2 become dynamic algorithms. Every convergence test, every adaptive timestep, every Monte Carlo acceptance criterion depends on mastering these concepts deeply, not just syntactically. By chapter’s end, you’ll see code not as a sequence of commands, but as a carefully orchestrated flow of decisions and iterations that solve real scientific problems. You’ll write algorithms that could process data from the James Webb Space Telescope or control the adaptive optics on the next generation of ground-based observatories.

3.1 Algorithmic Thinking: The Lost Art of Pseudocode¶

Most students jump straight from problem to code, then wonder why they spend hours debugging. Professional computational scientists spend more time thinking than typing. Pseudocode is how we think precisely about algorithms without getting distracted by syntax. Think of it as your algorithm’s blueprint — you wouldn’t build a telescope without optical designs, so why write code without algorithmic designs?

Why Pseudocode Matters in Scientific Computing¶

Consider this scenario: You need to implement adaptive timestepping for an orbital integrator. Without pseudocode, you’ll likely write code, run it, watch orbits spiral incorrectly, debug for hours, and maybe get it working through trial and error. With pseudocode, you’ll identify edge cases, boundary conditions, and logical flaws before writing a single line of Python.

# WITHOUT PSEUDOCODE (typical student approach):
# "I'll figure it out as I code..."
def integrate_naive(state, t_end):
    dt = 0.01
    while state.time < t_end:
        new_state = step(state, dt)
        error = estimate_error(state, new_state)
        if error > tolerance:
            dt = dt * 0.5  # Seems reasonable?
        state = new_state
    return state
# Wait, this doesn't work... infinite loop when error is bad!
# Also, dt never increases... hours of debugging ahead

Now let’s see how pseudocode reveals problems immediately! This is exactly how professional astronomers design algorithms for everything from orbit calculations to galaxy simulations. You’re about to learn the same systematic approach used at NASA’s Jet Propulsion Laboratory for spacecraft navigation and at the Space Telescope Science Institute for Hubble’s scheduling algorithms.

The Three Levels of Pseudocode Refinement¶

Professional algorithm development happens in stages, each revealing different issues. Don’t worry if this feels strange at first — every programmer has felt that way! But once you embrace pseudocode, you’ll save countless hours of debugging. Let’s build this skill together:

Level 1: Conceptual Overview (The Big Picture)

WHILE simulation not done:       # WHILE means "repeat as long as condition is true"
    Take a step
    Check IF step was good        # IF means "only do this when condition is true"
    Adjust timestep

This level helps you understand the overall flow. The WHILE construct creates a loop that continues until some condition becomes false. The IF construct makes a decision based on a condition. Already, we can ask critical questions: What defines “done”? What makes a step “good”? How much should we adjust? These questions matter!

Level 2: Structural Detail (The Flow)

FUNCTION adaptive_integrate(initial_state, end_time):  # FUNCTION groups reusable code
    state ← initial_state                              # ← means "assign value to variable"
    dt ← estimate_initial_timestep(state)
    
    WHILE time < end_time:                             # Loop continues while time hasn't reached end
        DO:                                             # DO-UNTIL creates a loop that runs at least once
            trial_step = integrate(state, dt)
            error = compute_error(trial_step)
        UNTIL error < tolerance OR dt < dt_min         # OR means "either condition can be true"
        
        state = trial_step
        dt = adjust_timestep(error, dt)
    
    RETURN state                                        # RETURN sends value back to caller

Now we see the retry logic and minimum timestep safeguard. The DO-UNTIL construct ensures we attempt at least one integration step. The OR operator means either condition being true will exit the inner loop. FUNCTION defines a reusable block of code that can be called with arguments and RETURN a result.

Level 3: Implementation-Ready (Stage 1: Core Logic)

FUNCTION adaptive_integrate(initial_state, end_time, tolerance):
    state ← initial_state
    dt ← estimate_initial_timestep(state)
    
    WHILE state.time < end_time:
        trial_state ← rk4_step(state, dt)
        error ← estimate_error(state, trial_state)
        
        IF error < tolerance:                          # Decision point
            state ← trial_state
            dt ← min(dt * 1.5, dt_max)                # Can grow
        ELSE:                                          # ELSE handles "otherwise" case
            dt ← max(dt * 0.5, dt_min)                # Must shrink

Level 3: Implementation-Ready (Stage 2: Add Safety)

FUNCTION adaptive_integrate(initial_state, end_time, tolerance):
    state ← initial_state
    dt ← estimate_initial_timestep(state)
    dt_min ← 1e-10 * (end_time - initial_state.time)
    dt_max ← 0.1 * (end_time - initial_state.time)
    
    WHILE state.time < end_time:
        step_accepted ← False                          # Boolean flag (True/False)
        attempts ← 0
        
        WHILE NOT step_accepted AND attempts < MAX_ATTEMPTS:  # NOT inverts, AND requires both
            trial_state ← rk4_step(state, dt)
            error ← estimate_error(state, trial_state)
            
            IF error < tolerance:
                step_accepted ← True
                state ← trial_state
            ELSE:
                dt ← max(dt * 0.5, dt_min)
            attempts ← attempts + 1

Each refinement level reveals new issues and solutions. The NOT operator inverts a boolean value (True becomes False, False becomes True). The AND operator requires both conditions to be true. This is computational thinking in action!

3.2 Boolean Logic in Scientific Computing¶

Every decision in your code ultimately reduces to true or false. But in scientific computing, these decisions often involve floating-point numbers, where equality is treacherous and precision is limited. Let’s master this fundamental building block that underlies everything from data quality checks to convergence criteria!

The Complete Set of Comparison Operators¶

Python provides six comparison operators that return boolean values (True or False):

# All comparison operators in Python
temperature = 5778  # Kelvin (Sun's surface)

# The six fundamental comparisons
print(f"Greater than: {temperature > 6000}")           # False
print(f"Less than: {temperature < 6000}")              # True
print(f"Greater or equal: {temperature >= 5778}")      # True
print(f"Less or equal: {temperature <= 5778}")         # True
print(f"Equal to: {temperature == 5778}")              # True
print(f"Not equal to: {temperature != 6000}")          # True

# Chaining comparisons (Python's elegant feature!)
print(f"\nMain sequence star? {3000 < temperature < 50000}")  # True
# This is equivalent to: (3000 < temperature) AND (temperature < 50000)

The != operator (not equal) is particularly useful for filtering out sentinel values or checking if something has changed. The ability to chain comparisons like 3000 < temperature < 50000 is a Python feature that makes code more readable and matches mathematical notation.

The Three Logical Operators: AND, OR, NOT¶

Python’s logical operators combine or modify boolean values:

# Demonstrating all three logical operators
is_bright = True
is_variable = False

# AND: Both conditions must be true
print(f"Bright AND variable: {is_bright and is_variable}")  # False

# OR: At least one condition must be true  
print(f"Bright OR variable: {is_bright or is_variable}")     # True

# NOT: Inverts the boolean value
print(f"NOT bright: {not is_bright}")                       # False
print(f"NOT variable: {not is_variable}")                   # True

# Complex combinations (real telescope scheduling logic!)
observable = True
weather_good = True
calibrated = False

can_observe = observable and weather_good and (not calibrated or calibrated)
print(f"\nCan observe? {can_observe}")

# Truth table demonstration
print("\nTruth Table for AND:")
for a in [True, False]:
    for b in [True, False]:
        print(f"  {a:5} AND {b:5} = {a and b}")

Special Comparison Operators: is, in¶

Python has two special operators that are incredibly useful in scientific programming:

import math

# The 'is' operator checks identity (same object in memory)
a = [1, 2, 3]
b = [1, 2, 3]
c = a

print(f"a == b: {a == b}")  # True - same values
print(f"a is b: {a is b}")  # False - different objects
print(f"a is c: {a is c}")  # True - same object

# Special case: None should always use 'is'
result = None
print(f"Checking None: {result is None}")  # Preferred
# Don't use: result == None

# The 'in' operator checks membership
stellar_types = ['O', 'B', 'A', 'F', 'G', 'K', 'M']
print(f"\nIs 'G' a stellar type? {'G' in stellar_types}")  # True
print(f"Is 'X' a stellar type? {'X' in stellar_types}")    # False

# Works with strings too!
filename = "observations_2024.fits"
print(f"Is FITS file? {'fits' in filename}")  # True

The Walrus Operator: Assignment Expressions (Python 3.8+)¶

Python 3.8 introduced the walrus operator (:=) which allows assignment within expressions:

# Traditional approach - two steps
data = get_observations()  # Assume this function exists
if len(data) > 100:
    print(f"Large dataset: {len(data)} observations")
    # Notice we call len(data) twice!

# With walrus operator - assign and test in one line
if (n := len(data)) > 100:
    print(f"Large dataset: {n} observations")
    # Now n contains the length, no need to recalculate!

# Real astronomical example (simulated)
def check_observation_quality(observations):
    """Check if we have enough high-quality observations."""
    # Without walrus operator:
    good_obs = [obs for obs in observations if obs['snr'] > 5]
    if len(good_obs) >= 10:
        print(f"Found {len(good_obs)} good observations")
        return good_obs
    
    # With walrus operator (Python 3.8+):
    if (good_count := len(good_obs)) >= 10:
        print(f"Found {good_count} good observations")
        return good_obs
    
    return None

# Useful in while loops too
# while (line := file.readline()):  # Read and check in one step
#     process(line)

print("Note: Walrus operator requires Python 3.8+")
print("It's useful but not essential - all code can be written without it")

The Floating-Point Equality Trap¶

Never use == with floating-point numbers! Even tiny rounding errors break equality:

# The floating-point trap strikes!
calculated_temp = 5778.0000000001
expected_temp = 5778.0

print(f"Calculated == Expected? {calculated_temp == expected_temp}")  # False!
print(f"Tiny difference: {calculated_temp - expected_temp:.2e}")

# The solution: tolerance-based comparison
def safe_equal(a, b, rel_tol=1e-9, abs_tol=1e-12):
    """
    Safe floating-point comparison for scientific computing.
    Used in actual telescope pointing systems!
    """
    # Handle exact equality (includes infinities)
    if a == b:
        return True
    
    # Handle NaN (NaN != NaN by IEEE standard)
    if math.isnan(a) or math.isnan(b):
        return False
    
    # Handle infinity cases
    if math.isinf(a) or math.isinf(b):
        return a == b
    
    # Normal comparison with tolerance
    # Note: math.isclose() does this internally, but understanding it matters!
    return abs(a - b) <= max(rel_tol * max(abs(a), abs(b)), abs_tol)

# Test our safe comparison
print(f"\n0.1 + 0.2 == 0.3? {0.1 + 0.2 == 0.3}")  # False!
print(f"Safe equal? {safe_equal(0.1 + 0.2, 0.3)}")  # True!
print(f"math.isclose? {math.isclose(0.1 + 0.2, 0.3)}")  # True!

Short-Circuit Evaluation: Order Matters!¶

Python’s and and or operators use short-circuit evaluation — they stop evaluating as soon as the result is determined:

# Short-circuit evaluation prevents errors and saves computation
stellar_data = []  # Empty for demonstration

# WRONG - will crash if data is empty!
# if stellar_data[0] > 0 and len(stellar_data) > 0:  # IndexError!

# CORRECT - checks length first (short-circuits if empty)
if len(stellar_data) > 0 and stellar_data[0] > 0:
    print("First star has positive measurement")
else:
    print("No data or first measurement not positive")

# Even more Pythonic (empty list evaluates to False)
if stellar_data and stellar_data[0] > 0:
    print("First star has positive measurement")

# OR also short-circuits
def expensive_check():
    print("  Running expensive calculation...")
    return True

# This won't call expensive_check() because True or anything is True
print("\nShort-circuit OR demonstration:")
result = True or expensive_check()  # expensive_check never runs!
print(f"Result: {result}")

# But this will call it
result = False or expensive_check()  # expensive_check must run
print(f"Result: {result}")

def check_collision_risk(satellite1, satellite2):
    """
    Simplified collision risk logic similar to conjunction assessment.
    Real systems use complex probability calculations, but follow 
    similar efficiency principles.
    
    Based on standard conjunction assessment practices where:
    - Initial screening uses simple distance checks
    - Detailed analysis only for close approaches
    - Probability calculations only when necessary
    """
    
    # Check cheap calculations first
    if distance > screening_threshold:  # ~10-50 km typical
        return "No risk"  # 99% filtered here
    
    # Only if close, calculate relative velocity
    if closing_velocity < 0:
        return "Moving apart"
    
    # Only if approaching, compute collision probability
    if collision_probability > 1e-4:  # Standard threshold
        return "COLLISION RISK!"
    
    return "Monitor"

3.3 Conditional Statements: Teaching Computers to Decide¶

Conditional statements are where your code makes decisions. In scientific computing, these decisions often involve numerical thresholds, convergence criteria, and boundary conditions. Let’s build your intuition for writing robust conditionals that could run on spacecraft or control telescopes!

The if Statement: Your First Decision Maker¶

The if statement is the simplest conditional — it executes code only when a condition is true:

# Basic if statement
magnitude = 4.5

if magnitude < 6.0:
    print(f"Star is visible to naked eye (mag {magnitude})")

# Nothing happens if condition is false
magnitude = 8.0
if magnitude < 6.0:
    print("This won't print")
    
# Multiple statements in if block
stellar_mass = 10.0  # Solar masses

if stellar_mass > 8:
    print("Massive star detected!")
    print(f"Mass: {stellar_mass} Msun")
    print("Will end as supernova")
    remnant = "neutron star or black hole"

The if-else Statement: Binary Decisions¶

The else clause provides an alternative when the condition is false:

# Binary decision with if-else
redshift = 0.8

if redshift < 0.1:
    classification = "nearby galaxy"
else:
    classification = "distant galaxy"
    
print(f"z = {redshift}: {classification}")

# You can have multiple statements in each block
observation_snr = 3.5  # Signal-to-noise ratio

if observation_snr >= 5.0:
    print("High quality detection")
    process_immediately = True
    confidence = "high"
else:
    print("Low SNR - needs verification")
    process_immediately = False
    confidence = "low"

The elif Statement: Multiple Choices¶

The elif (else if) statement allows multiple conditions to be checked in sequence:

def classify_stellar_remnant(mass_solar):
    """
    Determine stellar remnant type based on initial mass.
    Demonstrates guard clauses and elif chains.
    Based on Chandrasekhar limit and stellar evolution theory.
    """
    # Guard clauses - validate input FIRST
    if mass_solar <= 0:
        raise ValueError(f"Mass must be positive: {mass_solar}")
    
    if not math.isfinite(mass_solar):
        raise ValueError(f"Mass must be finite: {mass_solar}")
    
    # Main classification logic with elif chain
    if mass_solar < 0.08:
        remnant = "brown dwarf (failed star)"
    elif mass_solar < 8:
        remnant = "white dwarf"
    elif mass_solar < 25:
        remnant = "neutron star"
    else:  # Final catch-all
        remnant = "black hole"
    
    # Add uncertainty near boundaries (real issue in astronomy!)
    boundaries = [0.08, 8, 25]
    min_distance = min(abs(mass_solar - b) for b in boundaries)
    
    if min_distance < 0.5:
        remnant += " (near boundary - uncertain)"
    
    return remnant

# Test our classifier
print(classify_stellar_remnant(1.0))   # Our Sun's fate
print(classify_stellar_remnant(7.8))   # Near boundary!
print(classify_stellar_remnant(30))    # Massive star fate

Guard Clauses: Fail Fast, Fail Clear¶

Guard clauses handle special cases immediately, preventing deep nesting and making code clearer. This pattern is essential for scientific code where invalid inputs can cause subtle bugs hours into a simulation!

def calculate_orbital_period(a, M, validate=True):
    """
    Kepler's third law with comprehensive validation.
    This defensive style could have saved Mars Climate Orbiter!
    
    Parameters:
        a: Semi-major axis [AU]
        M: Central mass [solar masses]
    """
    # Guard clauses handle problems immediately
    if validate:
        if a <= 0:
            raise ValueError(f"Semi-major axis must be positive: {a} AU")
        if M <= 0:
            raise ValueError(f"Mass must be positive: {M} Msun")
        
        # Check for orbit inside Schwarzschild radius!
        rs_au = 2.95e-8 * M  # Schwarzschild radius in AU
        if a < rs_au:
            raise ValueError(f"Orbit inside black hole event horizon: a={a} AU, Rs={rs_au} AU")
    
    # Main calculation - only runs if guards pass
    G_au_msun = 39.478  # G in AU³/Msun/year²
    period_years = math.sqrt(a**3 / M)  # Simplified Kepler's third law
    
    # Sanity check result
    if validate and period_years > 13.8e9:
        import warnings
        warnings.warn(f"Period exceeds age of universe: {period_years:.2e} years")
    
    return period_years

# Test with real systems
print(f"Earth: {calculate_orbital_period(1.0, 1.0):.2f} years")
print(f"Mercury: {calculate_orbital_period(0.387, 1.0):.2f} years")
print(f"Proxima Centauri b: {calculate_orbital_period(0.0485, 0.122):.3f} years")

The Ternary Operator: Compact Conditionals¶

Python’s ternary operator provides a compact way to write simple if-else statements:

# Ternary operator: value_if_true if condition else value_if_false
magnitude = 3.5
visibility = "visible" if magnitude < 6.0 else "not visible"
print(f"Star with magnitude {magnitude} is {visibility}")

# Useful for setting defaults based on conditions
exposure_time = 30  # seconds
quality = "good" if exposure_time > 10 else "poor"

# Can be nested but don't overdo it!
stellar_class = "G"
temperature = 5778 if stellar_class == "G" else (7500 if stellar_class == "A" else 3500)
print(f"Class {stellar_class} star: ~{temperature}K")

def process_thrust_data(force, units):
    """This guard clause would have saved $327 million!"""
    
    # Validate units BEFORE processing
    valid_units = {'N': 1.0, 'lbf': 4.448222}  # Conversion factors
    
    if units not in valid_units:
        raise ValueError(f"Unknown units: {units}. Use 'N' or 'lbf'")
    
    # Convert to standard units (Newtons)
    force_newtons = force * valid_units[units]
    
    # Additional sanity check
    if force_newtons > 1000:  # Typical max thruster force
        warnings.warn(f"Unusually high thrust: {force_newtons} N")
    
    return force_newtons

# The actual error: ground software sent lbf, flight expected N
# Result: 4.45× error accumulated over months!

3.4 Loops: The Heart of Scientific Computation¶

Now that you’ve mastered making decisions with conditionals, let’s make your code repeat tasks efficiently! Loops are where your programs gain superpowers — they’re the difference between analyzing one star and analyzing millions. Every N-body simulation, every light curve analysis, every Monte Carlo calculation depends on loops. The patterns you learn here will appear in every algorithm you write for the rest of your career.

The for Loop: Iterating Over Sequences¶

The for loop iterates over any sequence (list, tuple, string, range):

# Basic for loop
stellar_types = ['O', 'B', 'A', 'F', 'G', 'K', 'M']

for spectral_class in stellar_types:
    print(f"Class {spectral_class} star")

# Using range() for counting
print("\nCounting with range:")
for i in range(5):  # 0, 1, 2, 3, 4 (not 5!)
    print(f"Observation {i}")

# Range with start, stop, step
print("\nEvery 2nd hour from 20:00 to 02:00:")
for hour in range(20, 26, 2):  # 20, 22, 24
    print(f"{hour:02d}:00")

The Accumulator Pattern in Astronomy¶

The accumulator pattern is fundamental to scientific computing:

# Calculating center of mass for a star cluster
star_masses = [1.2, 0.8, 2.1, 0.5, 1.5]  # Solar masses
star_positions = [0.1, 0.3, 0.5, 0.7, 0.9]  # Parsecs from origin

total_mass = 0
weighted_position = 0

for mass, position in zip(star_masses, star_positions):
    total_mass += mass
    weighted_position += mass * position

center_of_mass = weighted_position / total_mass
print(f"Cluster center of mass: {center_of_mass:.3f} pc")
print(f"Total cluster mass: {total_mass:.1f} Msun")

# This same pattern calculates:
# - Barycenter of binary systems (how we detect exoplanets!)
# - Average stellar metallicity in galaxies
# - Integrated luminosity functions
# - Weighted mean magnitudes

Common for Loop Patterns¶

Python provides several useful functions for loop patterns:

# enumerate() gives you index and value
magnitudes = [10.2, 10.1, 9.5, 10.3, 8.2, 10.2]

print("Finding bright events with enumerate:")
for i, mag in enumerate(magnitudes):
    if mag < 9.0:
        print(f"  Alert! Index {i}: magnitude {mag}")

# zip() for parallel iteration
times = [0, 1, 2, 3, 4]  # seconds
positions = [0, 4.9, 19.6, 44.1, 78.4]  # meters

print("\nParallel iteration with zip:")
for t, x in zip(times, positions):
    if t > 0:  # Avoid division by zero
        velocity = x / t
        print(f"  t={t}s: v={velocity:.1f} m/s")

The while Loop: Conditional Iteration¶

The while loop continues as long as a condition remains true:

# Basic while loop
count = 0
while count < 3:
    print(f"Iteration {count}")
    count += 1  # Don't forget to update!

# While loop for convergence
print("\nConvergence example:")
value = 100.0
target = 1.0
iteration = 0

while abs(value - target) > 0.01 and iteration < 100:  # Safety limit!
    value = value * 0.9 + target * 0.1  # Gradual approach
    iteration += 1
    if iteration <= 3 or iteration % 10 == 0:  # Print selectively
        print(f"  Iter {iteration}: value = {value:.3f}")

print(f"Converged to {value:.3f} after {iteration} iterations")

Loop Control: break, continue, and else¶

Python provides additional loop control statements:

# break: Exit loop early
print("Using break to find first detection:")
observations = [0.1, 0.3, 0.2, 5.8, 0.4, 6.2]

for obs in observations:
    if obs > 1.0:
        print(f"First significant detection: {obs}")
        break  # Stop searching
        
# continue: Skip to next iteration
print("\nUsing continue to skip bad data:")
measurements = [1.2, -999, 2.3, -999, 3.4]

for value in measurements:
    if value == -999:  # Sentinel value
        continue  # Skip this iteration
    print(f"Processing: {value}")

# else clause: Runs if loop completes without break
print("\nLoop else clause:")
search_list = [1, 2, 3, 4, 5]
target = 7

for item in search_list:
    if item == target:
        print("Found!")
        break
else:  # This runs because we didn't break
    print(f"Target {target} not found in list")

The pass Statement: Placeholder¶

The pass statement does nothing — useful as a placeholder:

# pass as placeholder
for i in range(3):
    if i == 1:
        pass  # TODO: Add processing here later
    else:
        print(f"Processing {i}")

# Often used in exception handling
try:
    risky_operation = 1 / 1  # No error this time
except ZeroDivisionError:
    pass  # Silently ignore this specific error

print("Continued after pass")

Nested Loops: Processing 2D Data¶

Loops can be nested to process multi-dimensional data:

# Nested loops for 2D grid (like CCD image)
print("Processing 3x3 pixel grid:")
for row in range(3):
    for col in range(3):
        pixel_value = row * 3 + col  # Simulated pixel value
        print(f"({row},{col})={pixel_value}", end="  ")
    print()  # New line after each row

# More realistic: Finding peaks in 2D data
data_2d = [
    [1, 2, 1],
    [2, 9, 2],  # Peak at (1,1)
    [1, 2, 1]
]

print("\nFinding peaks in 2D array:")
for i in range(len(data_2d)):
    for j in range(len(data_2d[i])):
        if data_2d[i][j] > 5:
            print(f"Peak at ({i},{j}): value={data_2d[i][j]}")

3.5 List Comprehensions: Elegant and Efficient¶

Now that you’ve mastered loops, let’s evolve them into something even more powerful! List comprehensions are Python’s gift to scientific programmers. They transform verbose loops into concise, readable, and faster expressions.

From Loop to Comprehension¶

# Traditional loop approach
squares_loop = []
for x in range(10):
    if x % 2 == 0:  # Even numbers only
        squares_loop.append(x**2)
print(f"Loop result: {squares_loop}")

# List comprehension - same result, clearer intent!
squares_comp = [x**2 for x in range(10) if x % 2 == 0]
print(f"Comprehension: {squares_comp}")

# The anatomy of a list comprehension:
# [expression for item in sequence if condition]
#      ↓           ↓         ↓           ↓
#  Transform   Variable  Source    Filter (optional)

Real Astronomical Applications¶

# Filter and transform magnitude data from a survey
magnitudes = [12.3, 15.1, 13.7, 16.2, 14.5, 17.1, 11.8, 18.5, 13.2]

# Get fluxes for observable stars (mag < 16, typical small telescope limit)
observable_fluxes = [10**(-0.4 * mag) 
                     for mag in magnitudes 
                     if mag < 16.0]

print(f"Observable star count: {len(observable_fluxes)}/{len(magnitudes)}")
print(f"Brightest flux: {max(observable_fluxes):.2e}")
print(f"Faintest flux: {min(observable_fluxes):.2e}")

# Dictionary comprehension (bonus!)
star_dict = {f"star_{i}": mag 
             for i, mag in enumerate(magnitudes) 
             if mag < 15}
print(f"\nBright stars dictionary: {star_dict}")

When NOT to Use Comprehensions¶

# BAD: Too complex - unreadable!
# result = [process(x) if condition(x) else alternative(y) 
#           for x, y in zip(list1, list2) 
#           if validate(x) and check(y)]

# GOOD: Clear loop for complex logic
def classify_galaxies(redshifts, luminosities):
    """When logic is complex, loops are clearer!"""
    classifications = []
    for z, L in zip(redshifts, luminosities):
        if z < 0.5 and L > 1e10:
            classifications.append("nearby bright")
        elif z < 0.5:
            classifications.append("nearby faint")
        elif L > 1e10:
            classifications.append("distant bright")
        else:
            classifications.append("distant faint")
    return classifications

# Clear and maintainable!
z_values = [0.1, 0.8, 0.3, 1.2]
L_values = [5e9, 2e10, 8e10, 3e9]
galaxy_types = classify_galaxies(z_values, L_values)
print(f"Galaxy classifications: {galaxy_types}")

3.6 Advanced Control Flow Patterns¶

Now let’s explore powerful patterns that appear throughout scientific computing. These aren’t just code tricks — they’re fundamental algorithmic building blocks!

Welford’s Algorithm: Numerically Stable Statistics¶

def running_statistics(data_stream):
    """
    Calculate mean and variance in single pass.
    Uses Welford's algorithm (1962) for numerical stability.
    Essential for processing streaming telescope data!
    """
    import math

    n = 0
    mean = 0.0
    M2 = 0.0
    
    for value in data_stream:
        n += 1
        delta = value - mean
        mean += delta / n
        delta2 = value - mean
        M2 += delta * delta2
    
    if n < 2:
        return mean, float('nan')
    
    variance = M2 / (n - 1)
    return mean, variance

# Test with problematic data (large baseline with small variations)
photometry = [1e8, 1e8 + 1, 1e8 + 2, 1e8 - 1, 1e8 + 0.5]  
mean, var = running_statistics(photometry)
print(f"Stable algorithm: mean={mean:.1f}, std={math.sqrt(var):.2f}")

# The naive approach would lose precision!
naive_mean = sum(photometry) / len(photometry)
naive_var = sum(x**2 for x in photometry) / len(photometry) - naive_mean**2
print(f"Naive (problematic): mean={naive_mean:.1f}, std={math.sqrt(abs(naive_var)):.2f}")

Why is Welford’s algorithm numerically stable?

The naive approach (sum all values, then divide) accumulates large sums that can lose precision. For values like [1e8, 1e8+1, 1e8+2], the sum becomes ~3e8, and the tiny variations (1, 2) get lost in floating-point representation.

Welford’s algorithm maintains a running mean and updates it incrementally with small deltas. Instead of computing (1e8 + 1e8 + 1e8)/3, it computes:

• mean = 1e8

• mean += (1e8 - 1e8)/2 = 0

• mean += (1e8 - 1e8)/3 = 0 This keeps all arithmetic operations on similar scales, preserving precision.

The variance calculation similarly avoids subtracting large nearly-equal numbers (catastrophic cancellation) by accumulating squared deviations incrementally.

3.7 Debugging Control Flow¶

Logic errors are the hardest bugs because the code runs without crashing but produces wrong results. Let’s build your debugging arsenal!

Strategic Print Debugging¶

def debug_convergence(initial, target, rate, max_iter=20):
    """
    Example of strategic debug output.
    Shows exactly where and why algorithms succeed or fail.
    """
    
    current = initial
    history = []
    
    for iteration in range(max_iter):
        old = current
        current = current * (1 - rate) + target * rate
        change = current - old
        history.append(current)
        
        # Strategic output - not everything!
        if iteration < 3 or iteration % 5 == 0:
            print(f"Iter {iteration:2d}: {old:.4f} → {current:.4f} (Δ={change:+.5f})")
        
        # Convergence check
        if abs(current - target) < 1e-6:
            print(f"✓ CONVERGED at iteration {iteration}")
            return current
        
        # Detect problems early
        if len(history) > 3:
            if (history[-1] > history[-2] < history[-3]):
                print(f"⚠ OSCILLATION detected")
                
    print(f"✗ FAILED after {max_iter} iterations")
    return current

# Test the algorithm
print("Testing convergence:")
result = debug_convergence(0, 100, 0.1)

Using Assertions for Validation¶

The assert statement helps catch bugs during development:

def calculate_magnitude_average(magnitudes):
    """Calculate average magnitude with assertions for debugging."""
    
    # Assertions document and enforce assumptions
    assert len(magnitudes) > 0, "Need at least one magnitude"
    assert all(isinstance(m, (int, float)) for m in magnitudes), "All must be numbers"
    assert all(0 < m < 30 for m in magnitudes), "Magnitudes must be reasonable"
    
    # Safe to proceed after assertions
    return sum(magnitudes) / len(magnitudes)

# Test with good data
good_mags = [10.2, 10.5, 10.3]
avg = calculate_magnitude_average(good_mags)
print(f"Average magnitude: {avg:.2f}")

# WARNING: Assertions can be disabled in production with python -O
# This means they should NEVER be used for actual validation!

def kepler_planet_search(star_id, light_curve):
    """Simplified Kepler planet detection algorithm"""
    
    # Guard clause - data quality check
    if len(light_curve) < 1000:
        return None
    
    # Remove outliers (cosmic rays, etc.)
    cleaned = [point for point in light_curve 
               if abs(point - median) < 5 * sigma]
    
    # Search for periodic dips
    best_period = None
    best_depth = 0
    
    for trial_period in range(1, 365):  # Days
        folded = fold_light_curve(cleaned, trial_period)
        depth = measure_transit_depth(folded)
        
        if depth > best_depth and depth > 3 * noise_level:
            best_period = trial_period
            best_depth = depth
    
    # Validate as planet (not eclipsing binary)
    if best_period:
        if is_v_shaped(folded):  # Binary check
            return None
        if depth > 0.5:  # Too deep
            return None
            
        return {'period': best_period, 'depth': best_depth}
    
    return None

def kepler_planet_search(star_id, light_curve):
    """Simplified Kepler planet detection algorithm"""
    
    # Guard clause - data quality check
    if len(light_curve) < 1000:
        return None
    
    # Remove outliers (cosmic rays, etc.)
    cleaned = [point for point in light_curve 
               if abs(point - median) < 5 * sigma]
    
    # Search for periodic dips
    best_period = None
    best_depth = 0
    
    for trial_period in range(1, 365):  # Days
        folded = fold_light_curve(cleaned, trial_period)
        depth = measure_transit_depth(folded)
        
        if depth > best_depth and depth > 3 * noise_level:
            best_period = trial_period
            best_depth = depth
    
    # Validate as planet (not eclipsing binary)
    if best_period:
        if is_v_shaped(folded):  # Binary check
            return None
        if depth > 0.5:  # Too deep
            return None
            
        return {'period': best_period, 'depth': best_depth}
    
    return None

This ran on 150,000+ stars for 4 years! Your code uses the same patterns that revealed the universe is full of planets! ([])

:::{admonition} 🛠️ Debug This! The Telescope Priority Bug
:class: challenge

A telescope scheduling system has a subtle bug in its priority logic. Can you find and fix it?

```python
def assign_telescope_priority(observation):
    """
    Assign priority based on object type and time sensitivity.
    Higher numbers = higher priority.
    
    THIS CODE HAS A BUG - find it!
    """
    magnitude = observation['magnitude']
    obj_type = observation['type']
    time_critical = observation['time_critical']
    
    # Assign base priority by object type
    if obj_type == 'supernova':
        priority = 100
    elif obj_type == 'variable_star' and magnitude < 12:
        priority = 70
    elif obj_type == 'variable_star':  # This has a problem!
        priority = 50
    elif obj_type == 'asteroid' and time_critical:
        priority = 80
    elif obj_type == 'galaxy':
        priority = 30
    else:
        priority = 10
    
    # Boost for bright objects
    if magnitude < 10:
        priority += 20
    
    return priority

# Test case that reveals the bug:
obs = {'type': 'variable_star', 'magnitude': 8, 'time_critical': True}
print(f"Priority: {assign_telescope_priority(obs)}")
# Expected: 70 + 20 = 90 (bright variable star)
# Actually gets: 90 (seems right... or is it?)

obs2 = {'type': 'variable_star', 'magnitude': 11.5, 'time_critical': True}
print(f"Priority: {assign_telescope_priority(obs2)}")
# Expected: 70 (bright variable star condition)
# Actually gets: 70 (correct!)

obs3 = {'type': 'variable_star', 'magnitude': 14, 'time_critical': True}
print(f"Priority: {assign_telescope_priority(obs3)}")
# Expected: 50 (dim variable star)
# What does it actually get?
```

def assign_telescope_priority_fixed(observation):
    """Fixed version with clearer logic."""
    magnitude = observation['magnitude']
    obj_type = observation['type']
    time_critical = observation['time_critical']
    
    # Assign base priority by object type
    if obj_type == 'supernova':
        priority = 100
    elif obj_type == 'asteroid' and time_critical:
        priority = 80
    elif obj_type == 'variable_star':
        # Nested logic is clearer for subcategories
        if magnitude < 12:
            priority = 70
        else:
            priority = 50
    elif obj_type == 'galaxy':
        priority = 30
    else:
        priority = 10
    
    # Boost for bright objects
    if magnitude < 10:
        priority += 20
    
    return priority

Main Takeaways¶

What an incredible journey you’ve just completed! You’ve transformed from someone who writes code line by line to someone who designs algorithms systematically. This transformation mirrors the evolution every computational scientist goes through, from tentative beginner to confident algorithm designer.

You started by learning to think in pseudocode, a skill that gives you the power to design before you code. Those three levels of refinement you practiced are your blueprint for success. Every hour you invest in pseudocode saves many hours of debugging. When you design your next algorithm for analyzing galaxy spectra or simulating stellar evolution, you’ll catch logical flaws on paper instead of after hours of computation.

The complete set of comparison and logical operators you’ve mastered — from simple greater-than checks to complex boolean combinations with and, or, and not — gives you the full vocabulary for expressing any logical condition. You understand that == is dangerous with floats, that is checks identity not equality, and that in elegantly tests membership. These aren’t just syntax details; they’re the building blocks of every data validation, every convergence check, every quality filter you’ll ever write.

Your understanding of conditional statements goes beyond syntax to defensive programming philosophy. Those guard clauses you learned could literally prevent spacecraft crashes or save irreplaceable telescope time. The elif chains you practiced will classify astronomical objects, determine observing strategies, and control instrument settings. Every conditional you write is a decision that shapes how your code responds to the infinite variety of real data.

The loop patterns you’ve mastered are universal across computational physics. That accumulator pattern using Welford’s algorithm? It’s calculating photometric precision in the TESS pipeline right now. The convergence pattern with safety limits? It’s finding equilibrium in stellar models. The nested loops you practiced? They’re processing CCD images from every major observatory. Whether using for loops to iterate through catalogs, while loops to converge solutions, or list comprehensions to filter data, you now have the full toolkit.

Most importantly, you’ve learned that bugs aren’t failures — they’re learning opportunities. Every infinite loop teaches you about termination conditions. Every off-by-one error reinforces proper indexing. The debugging strategies you’ve developed, from strategic print statements to assertions, will serve you throughout your career. Even the experts at NASA and ESA use these same techniques.

Remember that every major computational achievement relies on these fundamentals. The control flow patterns you’ve learned detected gravitational waves at LIGO, discovered thousands of exoplanets with Kepler, and process images from JWST. You’re not just learning Python syntax — you’re joining a tradition of computational thinking that enables humanity’s greatest discoveries.

Definitions¶

Accumulator Pattern: An algorithmic pattern where values are iteratively combined into a running total or aggregate, fundamental to reductions and statistical calculations in astronomical data processing.

Adaptive Refinement: A universal pattern where parameters are adjusted based on quality metrics, with safeguards against infinite refinement, appearing in timestepping, mesh refinement, and optimization throughout computational astrophysics.

and: Logical operator that returns True only if both operands are true, using short-circuit evaluation.

assert: Statement that raises an AssertionError if a condition is false, used for debugging and documenting assumptions during development (not for production validation).

Boolean Logic: The system of true/false values and logical operations (and, or, not) that underlies all conditional execution in programs.

break: Statement that immediately exits the current loop, skipping any remaining iterations.

Conditional Statement: A control structure (if/elif/else) that executes different code blocks based on whether conditions evaluate to true or false.

continue: Statement that skips the rest of the current loop iteration and proceeds to the next iteration.

elif: “Else if” statement that checks an additional condition when the previous if or elif was false.

else: Clause that executes when all previous if/elif conditions were false, or when a loop completes without breaking.

for: Loop that iterates over elements in a sequence or iterable object.

Guard Clause: A conditional statement at the beginning of a function that handles special cases or invalid inputs immediately, preventing deep nesting.

if: Statement that executes code only when a specified condition is true.

in: Operator that tests membership in a sequence or collection.

is: Operator that tests object identity (same object in memory), not just equality of values.

List Comprehension: A concise Python syntax for creating lists: [expression for item in iterable if condition].

not: Logical operator that inverts a boolean value (True becomes False, False becomes True).

or: Logical operator that returns True if at least one operand is true, using short-circuit evaluation.

pass: Null statement that does nothing, used as a placeholder where syntax requires a statement.

Pseudocode: A human-readable description of an algorithm that focuses on logic and structure without syntactic details.

sentinel value: A special marker value (like -999, ‘END’, or None) that signals the end of data or a special condition, allowing loops to know when to stop processing

Short-circuit Evaluation: The behavior where logical operators stop evaluating as soon as the result is determined.

Walrus Operator (:=): Assignment expression operator (Python 3.8+) that assigns a value to a variable as part of an expression, allowing both assignment and testing in a single statement.

while: Loop that continues executing as long as a specified condition remains true.

Key Takeaways¶

✓ Pseudocode reveals logical flaws before they become bugs — always design before implementing

✓ Master all six comparison operators (>, <, >=, <=, ==, !=) and three logical operators (and, or, not)

✓ Never use == with floating-point numbers; always use tolerance-based comparisons like math.isclose()

✓ The is operator checks identity, not equality — use it for None checks

✓ The in operator elegantly tests membership in sequences or strings

✓ Guard clauses handle special cases first, making main logic clearer

✓ for loops iterate over sequences, while loops continue until a condition becomes false

✓ break exits loops early, continue skips to the next iteration, else runs if loop completes

✓ List comprehensions are faster than loops for simple transformations but become unreadable for complex logic

✓ Short-circuit evaluation in and/or prevents errors and improves performance

✓ The accumulator pattern is fundamental to scientific computing, appearing in all statistical calculations

✓ Always include maximum iteration limits in while loops to prevent infinite loops

✓ Welford’s algorithm solves numerical stability issues in streaming statistics

✓ Use assert statements to document and enforce assumptions during development

✓ Every major astronomical discovery relies on the control flow patterns you’ve learned

Quick Reference Tables¶

Python Module & Method Reference (Chapter 3 Additions)¶

New Built-in Functions¶

Logical Testing

• all(iterable) - Returns True if all elements are true

• any(iterable) - Returns True if any element is true

• isinstance(obj, type) - Check if object is of specified type

Loop Support

• enumerate(iterable, start=0) - Returns index-value pairs

• zip(*iterables) - Combines multiple iterables for parallel iteration

• range(start, stop, step) - Generate arithmetic progression

Control Flow Keywords¶

Conditionals

• if - Execute block if condition is true

• elif - Check additional condition if previous was false

• else - Execute if all previous conditions were false

Loops

• for - Iterate over sequence

• while - Loop while condition is true

• break - Exit loop immediately

• continue - Skip to next iteration

• else - Execute if loop completes without break

Other

• pass - Null operation placeholder

• assert - Raise AssertionError if condition is false

Operators¶

Comparison

• >, <, >=, <=, ==, != - Numerical comparisons

• is, is not - Identity comparisons

• in, not in - Membership testing

Logical

• and - Logical AND with short-circuit evaluation

• or - Logical OR with short-circuit evaluation

• not - Logical NOT (inversion)

New Math Module Functions¶

import math

• math.isclose(a, b, rel_tol=1e-9, abs_tol=0.0) - Safe floating-point comparison

• math.isfinite(x) - Check if neither infinite nor NaN

• math.isnan(x) - Check if value is NaN

• math.isinf(x) - Check if value is infinite

Debugging Support¶

IPython Magic Commands

• %debug - Enter debugger after exception

• %pdb - Automatic debugger on exceptions

Debugger Commands (when in pdb)

• p variable - Print variable value

• pp variable - Pretty-print variable

• l - List code around current line

• n - Next line

• s - Step into function

• c - Continue execution

• u/d - Move up/down call stack

• q - Quit debugger

Next Chapter Preview¶

You’ve conquered control flow — now get ready for the next level! Chapter 4 will reveal how to organize data efficiently using Python’s powerful data structures. You’ll discover when to use lists versus dictionaries versus sets, and more importantly, you’ll understand why these choices can make your algorithms run 100 times faster or 100 times slower.

Imagine trying to find a specific star in a catalog of millions. With a list, you’d check each star one by one — taking minutes or hours. With a dictionary, you’ll find it instantly — in microseconds! The data structures you’ll learn next are the difference between simulations that finish in minutes and ones that run for days.

The control flow patterns you’ve mastered here will operate on the data structures you’ll learn next. Your loops will iterate through dictionaries of astronomical objects. Your conditionals will filter sets of observations. Your comprehensions will transform lists of measurements into meaningful results. Together, control flow and data structures give you the power to handle the massive datasets of modern astronomy — from Gaia’s billion-star catalog to the petabytes of data from the Square Kilometre Array.

Get excited — Chapter 4 is where your code goes from processing dozens of data points to handling thousands efficiently!