Lab 04: File I/O and Python Functions#

1. Demonstration#

1.1 Example 1: Quadratic Equation Solver#

Consider a function roots(a, b, c) that returns the roots of quadratic equation:

\[ ax^2 + bx + c = 0 \]

Using the quadratic formula:

\[ x = \frac{-b \pm \sqrt{b^2 - 4ac}}{2a} \]

We’ll assume the equation has two real roots (discriminant ≥ 0).

import math

def roots(a, b, c):

    
    return root1, root2
# Example 1: x^2 - 5x + 6 = 0 (roots: 3, 2)
print("Example 1: x² - 5x + 6 = 0")
Example 1: x² - 5x + 6 = 0
# Example 2: 2x^2 + 4x - 6 = 0 (roots: 1, -3)
print("Example 2: 2x² + 4x - 6 = 0")
Example 2: 2x² + 4x - 6 = 0
# Example 3: x^2 - 4 = 0 (roots: 2, -2)
print("Example 3: x² - 4 = 0")
Example 3: x² - 4 = 0

1.2 Example 2: Gravitational Force Calculator with Filtering#

Define a function to compute the gravitational force:

\[ F = G\frac{m_1 m_2}{r^2} \]

where:

  • G = 6.674 × 10⁻¹¹ N·m²/kg² (gravitational constant)

  • m₁, m₂ = masses in kg

  • r = distance between centers in meters

  • F = gravitational force in Newtons

We’ll loop over several masses and print results only if F < 1 N.

1.3 Example 3: Force Calculation with File I/O#

Define a function to compute the force F = ma that takes two arguments (m and a) and return the value of F in Newtons. Loop over several masses from a list. Print the results on the screen and save them into the file “masses_file.dat”.


2. Practice Problems#

Now it’s your turn! Complete the following problems using functions, loops, and conditional statements.

Problem 1: Electrostatic Force Calculator with Filtering#

Create a function electrostatic_force(q1, q2, r) that calculates the electrostatic force between two charges using Coulomb’s Law:

\[ F = k_e \frac{|q_1 q_2|}{r^2} \]

where:

  • \(k_e\) = 8.99 × 10⁹ N·m²/C² (Coulomb’s constant)

  • q₁, q₂ = charges in Coulombs (C)

  • r = distance between charges in meters

  • F = electrostatic force in Newtons

Tasks:

  1. Define the function to return the electrostatic force

  2. Test with the following charge pairs (all with q₁ = 1.0 × 10⁻⁶ C):

    • q₂ values (C): [1.0e-6, 2.0e-6, 5.0e-6, 1.0e-5, 5.0e-5, 1.0e-4]

    • r values (m): [0.01, 0.02, 0.05, 0.1, 0.2, 0.5]

  3. Loop over all combinations (6 test cases total)

  4. Print results ONLY if F > 100 N

  5. Count and report how many cases meet this criterion

  6. Create a formatted table showing q₁, q₂, r, and F for filtered results

Format: Similar to the gravitational force example, show a status column indicating whether each case was printed or filtered.

# Test parameters
q1 = 1.0e-6  # C
q2_values = [1.0e-6, 2.0e-6, 5.0e-6, 1.0e-5, 5.0e-5, 1.0e-4]
r_values = [0.01, 0.02, 0.05, 0.1, 0.2, 0.5]

Problem 2: Equations of State - Comparing Gas Models#

Using the following equations of state, define Python functions returning the pressure of a gas (p), given its volume (V), amount (n) and temperature (T). Use SI units.

Ideal Gas Equation:

\[ p = \frac{nRT}{V} \]

Van der Waals Equation:

\[ p = \frac{nRT}{V - nb} - \frac{n^2a}{V^2} \]

Dieterici Equation:

\[ p = \frac{nRT \exp\left(-\frac{na}{RTV}\right)}{V - nb} \]

Tasks:

  1. Define three functions:

    • p_ideal(n, V, T, R)

    • p_vdw(n, V, T, R, a, b) (Van der Waals)

    • p_dieterici(n, V, T, R, a, b) (Dieterici)

  2. Compare the pressure for 0.5 mol of CO₂ at three temperatures confined to a volume of 10 L:

    • T = 273.15 K, 303.15 K, and 333.15 K (use a list)

  3. Use the following constants:

    • R = 8.3144626 J/(K·mol)

    • Van der Waals parameters: a = 3.640 L²·bar/mol², b = 0.04267 L/mol

    • Dieterici parameters: a = 4.692 L²·bar/mol², b = 0.04639 L/mol

  4. Important unit conversions:

    • Convert volume from L to m³: V_m3 = V_L / 1000

    • Convert a from L²·bar/mol² to Pa·m⁶/mol²: multiply by 100

    • Convert b from L/mol to m³/mol: divide by 1000

  5. Use a for loop to iterate over temperatures and print results for each T

  6. For each temperature, print:

    • Temperature in K

    • Ideal gas pressure in Pa

    • Van der Waals pressure in Pa

    • Dieterici pressure in Pa

Hint: Use math.exp() for the exponential function in Dieterici equation.

# Constants
R = 8.3144626  # J/(K·mol)
n = 0.5        # mol CO2
V_L = 10       # Volume in L

# Unit conversions
V_m3 = V_L / 1000  # L to m³

# Van der Waals parameters for CO2 (converted to SI)
a_vdw = 3.640 / 1000     # L²·bar/mol² to Pa·m⁶/mol²
b_vdw = 0.04267 / 1000   # L/mol to m³/mol

# Dieterici parameters for CO2 (converted to SI)
a_diet = 4.692 / 1000    # L²·bar/mol² to Pa·m⁶/mol²
b_diet = 0.04639 / 1000  # L/mol to m³/mol

# Temperatures to test
temperatures = [273.15, 303.15, 333.15]

Problem 3: Van der Waals Equation with File Input#

This part was done in the previous lab. Using the following van der Waals equation of states, define the Python function returning the pressure of a gas (p), given its volume (V), amount (n) and temperature (T). Use 1.0 mol of CO₂ at different temperatures from the file (“temperatures.dat”) confined to a volume of 0.5 L. Use R = 0.082057 [L atm K⁻¹ mol⁻¹]. Take the van der Waals parameters a = 3.640 L² bar/mol² and b = 0.04267 L/mol.

Van der Waals equation:

\[ p = \frac{nRT}{V - nb} - \frac{n^2a}{V^2} \]

First, create the file “temperatures.dat” containing the following:

# Temperature, K
50
70
120
150
200
250
# Create the temperatures.dat file
with open('temperatures.dat', 'w') as f:
    f.write("# Temperature, K\n")
    f.write("50\n")
    f.write("70\n")
    f.write("120\n")
    f.write("150\n")
    f.write("200\n")
    f.write("250\n")

print("File 'temperatures.dat' created successfully")
File 'temperatures.dat' created successfully
# Constants
R = 0.082057   # L·atm/(K·mol)
n = 1.0        # mol CO2
V = 0.5        # L
a = 3.640      # L²·bar/mol²
b = 0.04267    # L/mol