NFYK22003U Computational Atmosphere and Ocean Dynamics

It is very much a work in progress! Please do not share it yet, as there should be a better version soon.
Course Information
Schedule 2025
Sibling: NFYB21003U Fluid Mechanics
Next: Lecture 1 – Introduction and Transport Processes
"""
        vonNeuman_elephant.py
            "With four parameters I can fit an elephant,
               and with five I can make him wiggle his trunk."
        
        Original Version
            Author: Piotr A. Zolnierczuk
            Retrieved on 14 September 2011 from
            http://www.johndcook.com/blog/2011/06/21/how-to-fit-an-elephant/
        Modified to wiggle trunk:
            2 October 2011 by David Bailey (http://www.physics.utoronto.ca/~dbailey)
        
        Based on the paper:
            "Drawing an elephant with four complex parameters", by
            Jurgen Mayer, Khaled Khairy, and Jonathon Howard,
            Am. J. Phys. 78, 648 (2010), DOI:10.1119/1.3254017
        
            The paper does not specify how the wiggle parameter controls the
            trunk, so a guess was made.
        
            Attributed to von Neumann by Enrico Fermi, as quoted by
              Freeman Dyson in "A meeting with Enrico Fermi" in
              Nature 427 (22 January 2004) p. 297 
        """
import matplotlib
matplotlib.use('Agg')  # Use non-interactive backend

import matplotlib.pyplot as plt
import numpy as np
import os

# Elephant parameters
p = [50 - 30j, 18 + 8j, 12 - 10j, -14 - 60j, 40 + 20j]

def fourier(t, C):
    f = np.zeros(t.shape)
    for k in range(len(C)):
        f += C.real[k]*np.cos(k*t) + C.imag[k]*np.sin(k*t)
    return f

def elephant(t, p):
    npar = 6
    Cx = np.zeros((npar,), dtype='complex')
    Cy = np.zeros((npar,), dtype='complex')
    Cx[1] = p[0].real * 1j
    Cy[1] = p[3].imag + p[0].imag * 1j
    Cx[2] = p[1].real * 1j
    Cy[2] = p[1].imag * 1j
    Cx[3] = p[2].real
    Cy[3] = p[2].imag * 1j
    Cx[5] = p[3].real
    x = np.append(fourier(t, Cy), [p[4].imag])
    y = -np.append(fourier(t, Cx), [-p[4].imag])
    return x, y

# Create output directory
os.makedirs("frames", exist_ok=True)

# Plot and save 10 frames of the wiggling trunk
for i in range(10):
    t_body = np.linspace(0.4 + 1.3*np.pi, 2*np.pi + 0.9*np.pi, 1000)
    x_body, y_body = elephant(t_body, p)

    t_trunk = np.linspace(2*np.pi + 0.9*np.pi, 0.4 + 3.3*np.pi, 1000)
    x_trunk, y_trunk = elephant(t_trunk, p)

    for ii in range(len(y_trunk) - 1):
        y_trunk[ii] -= np.sin((x_trunk[ii] - x_trunk[0]) * np.pi / len(y_trunk)) * np.sin(i) * p[4].real

    fig, ax = plt.subplots(figsize=(6, 6))
    ax.plot(x_body, y_body, 'r.', label='Body')
    ax.plot(x_trunk, y_trunk, 'r.', label='Trunk')
    ax.set_aspect('equal')
    ax.axis('off')
    plt.tight_layout()
    plt.savefig(f"frames/elephant_frame_{i:03d}.png", dpi=150)
    plt.close()

print("Saved 10 frames in ./frames")