diff --git a/.gitignore b/.gitignore
index f07f617..f7dd7f8 100644
--- a/.gitignore
+++ b/.gitignore
@@ -20,41 +20,37 @@ parts/
 sdist/
 var/
 wheels/
-pip-wheel-metadata/
-share/python-wheels/
 *.egg-info/
 .installed.cfg
 *.egg
-MANIFEST
 
 # PyInstaller
 *.manifest
 *.spec
 
+# Installer logs
+pip-log.txt
+pip-delete-this-directory.txt
+
 # Unit test / coverage reports
 htmlcov/
 .tox/
-.nox/
 .coverage
 .coverage.*
 .cache
 nosetests.xml
 coverage.xml
 *.cover
-*.py,cover
 .hypothesis/
 .pytest_cache/
 
-# Environments
-.env
-.venv
-env/
-venv/
-ENV/
-env.bak/
-venv.bak/
+# Jupyter Notebook
+.ipynb_checkpoints
+
+# pyenv
+.python-version
 
-# IDE
+# IDEs
 .vscode/
 .idea/
 *.swp
@@ -65,6 +61,4 @@ venv.bak/
 .DS_Store
 Thumbs.db
 
-# Temporary files
-/tmp/
-*.tmp
+
diff --git a/README.md b/README.md
index 5093580..ca20ac7 100644
--- a/README.md
+++ b/README.md
@@ -1,95 +1,202 @@
-# pyDiwasp
+# pyDIWASP
 
-[![CI](https://github.com/SBFRF/pyDIWASP/actions/workflows/ci.yml/badge.svg)](https://github.com/SBFRF/pyDIWASP/actions/workflows/ci.yml)
+Python conversion of the DIWASP package (DIrectional WAve SPectrum analysis Version 1.4)
 
-conversion of diwasp package (DIWASP: DIrectional WAve SPectrum analysis Version 1.4) for python
-converted from https://github.com/metocean/diwasp
+pyDIWASP is a Python implementation of the DIWASP toolbox for estimating directional wave spectra from wave measurement data. It provides functions to analyze wave measurements from arrays of instruments (e.g., pressure sensors, current meters, wave gauges) and compute the directional distribution of wave energy.
 
-I would LOVE help making this into better package of the original diwasp tool. Please check issues for needed functionality adds.
+Converted from: https://github.com/metocean/diwasp
+
+## Features
+
+- **Directional Wave Analysis**: Estimate directional wave spectra from instrument array data
+- **Multiple Estimation Methods**: Supports IMLM and EMEP methods
+- **Flexible Input**: Works with various instrument types (pressure, velocity, elevation sensors)
+- **Visualization**: Built-in plotting functions for 3D and polar spectral plots
+- **Wave Statistics**: Calculate significant wave height, peak period, and dominant direction
 
 ## Installation
 
-### From source
+### Basic Installation
+
 ```bash
+# Clone the repository
 git clone https://github.com/SBFRF/pyDIWASP.git
 cd pyDIWASP
-pip install -r requirements.txt
+
+# Install dependencies
+pip install numpy scipy matplotlib
 ```
 
 ### Requirements
-- Python 3.8+
-- NumPy >= 1.20.0, < 2.0
-- SciPy >= 1.7.0, < 2.0
-- Matplotlib >= 3.3.0, < 4.0
 
-## Testing
+- Python 3.6+
+- NumPy
+- SciPy
+- Matplotlib
+
+## Quick Start
+
+```python
+import numpy as np
+from dirspec import dirspec
+
+# Define instrument data structure
+ID = {
+    'layout': np.array([[0, 10, 20], [0, 0, 0], [0, 0, 0]]),  # x, y, z positions
+    'datatypes': ['pres', 'pres', 'pres'],  # instrument types
+    'depth': 10.0,  # water depth in meters
+    'fs': 2.0,  # sampling frequency in Hz
+    'data': wave_data  # nsamples x ninstruments array
+}
+
+# Define spectral matrix structure
+SM = {
+    'freqs': np.linspace(0.05, 0.5, 50),  # frequency bins in Hz
+    'dirs': np.linspace(-np.pi, np.pi, 36)  # direction bins in radians
+}
+
+# Define estimation parameters
+EP = {
+    'method': 'IMLM',  # estimation method
+    'iter': 100  # number of iterations
+}
+
+# Compute directional spectrum
+SMout, EPout = dirspec(ID, SM, EP)
+```
 
-The package includes a comprehensive test suite that documents the existing capabilities:
+For more detailed examples, see the [example notebook](examples/pyDIWASP_example.ipynb).
 
-```bash
-# Install test dependencies
-pip install pytest pytest-cov
+## Main Functions
 
-# Run all tests
-pytest tests/ -v
+### dirspec.py
+Main function for directional wave spectrum analysis.
 
-# Run tests with coverage
-pytest tests/ -v --cov=. --cov-report=term
+**Usage:**
+```python
+SMout, EPout = dirspec(ID, SM, EP, Options)
 ```
 
-**Test Coverage:**
-- 25 tests covering core functionality
-- Core functions: wavenumber calculations, wave height calculations, data validation
-- API functions: spectrum info, interpolation, file I/O
-- Integration tests: full directional analysis workflow
-
-## Toolbox contents:
-### Main functions:
-- dirspec.m           Main function for directional wave analysis
-- readspec.m          Reads in DIWASP format spectrum files
-- writespec.m         Writes DIWASP format spectrum files
-- plotspec.m          Plots DIWASP spectrums
-- testspec.m          Testing function for the estimation methods
-- makespec.m          Makes a fake spectrum and generates fake data for testing dirspec.m
-- infospec.m          Returns information about a directional spectrum
-- data_structures.m   is a help file describing the new Version 1.1 data structures
-
-## Private functions (some can be used as stand alone functions):
-### The transfer functions
-- /private/elev.m
-- /private/pres.m
-- /private/velx.m
-- /private/vely.m
-- /private/velz.m
-- /private/slpx.m
-- /private/slpy.m
-- /private/vels.m
-- /private/accs.m
-
-All transfer functions are now implemented in pyDIWASP!
-
-### The estimation functions
-- /private/DFTM.m
-- /private/EMLM.m
-- /private/IMLM.m
-- /private/EMEP.m
-- /private/BDM.m
-
-All estimation methods are now implemented in pyDIWASP!
-
-### Miscellaneous functions
-- /private/smoothspec.m
-- /private/wavenumber.m
-- /private/makerandomsea.m
-- /private/makewavedata.m
-- /private/Hsig.m
-- /private/gsamp.m
-- /private/check_data.m
-  
-
-carying original license agreement and copyright
-
-## License agreement
+**Parameters:**
+- `ID`: Instrument data structure (dict)
+- `SM`: Spectral matrix structure (dict)
+- `EP`: Estimation parameters structure (dict)
+- `Options`: Optional parameters (list of key-value pairs)
+
+**Returns:**
+- `SMout`: Output spectral matrix with computed directional spectrum
+- `EPout`: Estimation parameters with actual values used
+
+### infospec.py
+Calculates and displays information about a directional spectrum.
+
+**Returns:** Significant wave height (Hsig), peak period (Tp), direction of peak period (DTp), and dominant direction (Dp).
+
+### plotspec.py
+Plots the spectral matrix in 3D or polar form.
+
+**Plot types:**
+1. 3D surface plot
+2. Polar plot
+3. 3D surface plot (compass bearings)
+4. Polar plot (compass bearings)
+
+### writespec.py
+Writes directional spectrum data to file in DIWASP format.
+
+### interpspec.py
+Interpolates a spectrum onto a different frequency/direction grid.
+
+## Private Functions
+
+The `private/` directory contains internal functions used by the main analysis routines:
+
+### Transfer Functions
+Calculate instrument response to waves:
+- `elev.py` - Surface elevation transfer function
+- `pres.py` - Pressure sensor transfer function
+- `velx.py`, `vely.py` - Horizontal velocity transfer functions
+- `wavenumber.py` - Dispersion relation solver
+
+### Estimation Methods
+Directional spectrum estimation algorithms:
+- `IMLM.py` - Iterated Maximum Likelihood Method (default)
+- `EMEP.py` - Extended Maximum Entropy Principle
+
+### Utility Functions
+- `smoothspec.py` - Spectral smoothing
+- `hsig.py` - Significant wave height calculation
+- `check_data.py` - Data structure validation
+- `diwasp_csd.py` - Cross-spectral density estimation
+- `spectobasis.py` - Spectral basis conversion
+
+## Data Structures
+
+### Instrument Data Structure (ID)
+Dictionary with the following fields:
+- `layout`: 3×N array of instrument positions [x; y; z] in meters
+- `datatypes`: List of instrument types ('pres', 'elev', 'velx', 'vely', etc.)
+- `depth`: Water depth in meters
+- `fs`: Sampling frequency in Hz
+- `data`: N_samples × N_instruments array of measurements
+
+### Spectral Matrix Structure (SM)
+Dictionary with the following fields:
+- `freqs`: Array of frequency values
+- `dirs`: Array of direction values (radians or degrees)
+- `S`: Frequency × Direction array of spectral density (output only)
+- `funit`: Frequency unit ('Hz' or 'rad/s')
+- `dunit`: Direction unit ('rad' or 'deg')
+- `xaxisdir`: Reference axis direction (default: 90 = East)
+
+### Estimation Parameters Structure (EP)
+Dictionary with the following fields:
+- `method`: Estimation method ('IMLM' or 'EMEP')
+- `nfft`: FFT length (auto-calculated if not specified)
+- `dres`: Directional resolution (default: 180)
+- `iter`: Number of iterations for iterative methods (default: 100)
+- `smooth`: Spectral smoothing ('ON' or 'OFF', default: 'ON')
+
+## Estimation Methods
+
+### IMLM (Iterated Maximum Likelihood Method)
+- **Default method**
+- Iteratively refines directional spectrum estimate
+- Good balance of accuracy and computational efficiency
+- Recommended for most applications
+
+### EMEP (Extended Maximum Entropy Principle)
+- Based on maximum entropy principle
+- Works well for narrow directional spreads
+- Suitable for swell-dominated conditions
+
+**Note:** The original DIWASP MATLAB toolbox includes additional methods (EMLM, DFTM, BDM) that are not yet implemented in this Python version. Contributions to add these methods are welcome!
+
+## Examples
+
+See the [examples directory](examples/) for Jupyter notebooks demonstrating:
+- Basic directional spectrum estimation
+- Working with different instrument configurations
+- Comparing estimation methods
+- Visualization options
+
+## Contributing
+
+Contributions are welcome! Please feel free to submit pull requests or open issues for:
+- Bug fixes
+- Documentation improvements
+- New features
+- Additional examples
+
+## References
+
+All implemented calculation algorithms are described in:
+
+Hashimoto, N. (1997). "Analysis of the directional wave spectrum from field data". 
+In: *Advances in Coastal Engineering Vol. 3*. Ed: Liu, P.L-F. 
+Pub: World Scientific, Singapore.
+
+## Original Copyright and License
 DIWASP, is free software; you can redistribute it and/or modify it under the terms of the 
 GNU General Public License as published by the Free Software Foundation. 
 However, the DIWASP license includes the following addendum concerning its usage:
diff --git a/dirspec.py b/dirspec.py
index da4d393..feb302e 100644
--- a/dirspec.py
+++ b/dirspec.py
@@ -91,7 +91,7 @@ def dirspec(ID, SM, EP, Options_=None):
     ndat, szd = np.shape(ID['data'])
 
     #get resolution of FFT - if not specified, calculate a sensible value
-    if 'nfft' not in EP:
+    if 'nfft' not in EP or not EP['nfft']:
         nfft = int(2 ** (8 + np.round(np.log2(ID['fs']))))
         EP['nfft'] = nfft
     else:
diff --git a/examples/pyDIWASP_example.ipynb b/examples/pyDIWASP_example.ipynb
new file mode 100644
index 0000000..4b372e3
--- /dev/null
+++ b/examples/pyDIWASP_example.ipynb
@@ -0,0 +1,609 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# pyDIWASP Example: Directional Wave Spectrum Analysis\n",
+    "\n",
+    "This notebook demonstrates how to use pyDIWASP to analyze directional wave spectra from instrument array data.\n",
+    "\n",
+    "## Overview\n",
+    "\n",
+    "pyDIWASP is a Python implementation of the DIWASP (DIrectional WAve SPectrum) toolbox. It estimates the directional distribution of wave energy from measurements by an array of wave instruments.\n",
+    "\n",
+    "### What is a Directional Wave Spectrum?\n",
+    "\n",
+    "A directional wave spectrum describes how wave energy is distributed across:\n",
+    "- **Frequencies** (or wave periods)\n",
+    "- **Directions** (where waves are coming from)\n",
+    "\n",
+    "This is more informative than a 1D frequency spectrum because it tells us not just how much energy is at each frequency, but also which direction that energy is traveling."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Setup\n",
+    "\n",
+    "First, import the necessary libraries:\n",
+    "\n",
+    "**Note on imports:** This notebook assumes you're running it from the `examples/` directory. If you encounter import errors, you can either:\n",
+    "1. Install pyDIWASP in development mode: `pip install -e .` (from repo root)\n",
+    "2. Ensure your working directory is the `examples/` folder when starting Jupyter\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import sys\n",
+    "import os\n",
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "\n",
+    "# Add pyDIWASP to Python path\n",
+    "# This assumes you're running the notebook from the examples/ directory\n",
+    "# The parent directory (..) contains the pyDIWASP modules\n",
+    "if os.getcwd().endswith('examples'):\n",
+    "    sys.path.insert(0, '..')\n",
+    "elif os.path.basename(os.getcwd()) != 'pyDIWASP':\n",
+    "    print('Warning: Notebook should be run from examples/ directory or repo root')\n",
+    "\n",
+    "from dirspec import dirspec\n",
+    "from infospec import infospec\n",
+    "from plotspec import plotspec\n",
+    "\n",
+    "# Set random seed for reproducibility\n",
+    "np.random.seed(42)\n",
+    "\n",
+    "print(\"Libraries imported successfully!\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Example 1: Basic Directional Spectrum Analysis\n",
+    "\n",
+    "In this example, we'll create synthetic wave data representing a wave field with:\n",
+    "- A dominant wave direction from the east (90 degrees)\n",
+    "- Peak period around 10 seconds (0.1 Hz)\n",
+    "- Wave height around 2 meters\n",
+    "\n",
+    "We'll use a simple 3-sensor array to measure pressure at different locations."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Step 1: Generate Synthetic Wave Data\n",
+    "\n",
+    "We'll create a simple wave field using superposition of sinusoidal components:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Simulation parameters\n",
+    "duration = 1024  # seconds\n",
+    "fs = 2.0  # sampling frequency (Hz)\n",
+    "depth = 10.0  # water depth (m)\n",
+    "\n",
+    "# Time array\n",
+    "t = np.arange(0, duration, 1/fs)\n",
+    "nsamples = len(t)\n",
+    "\n",
+    "# Wave parameters\n",
+    "peak_freq = 0.1  # Hz (10 second period)\n",
+    "peak_dir = np.pi/2  # radians (from East)\n",
+    "wave_amplitude = 1.0  # meters\n",
+    "\n",
+    "# Instrument array layout (3 pressure sensors in an L-shape)\n",
+    "# Layout is 3xN array: [x positions; y positions; z positions]\n",
+    "layout = np.array([\n",
+    "    [0, 10, 0],      # x positions (meters)\n",
+    "    [0, 0, 10],      # y positions (meters)\n",
+    "    [-depth, -depth, -depth]  # z positions (meters, negative = below surface)\n",
+    "])\n",
+    "\n",
+    "ninst = layout.shape[1]\n",
+    "\n",
+    "print(f\"Simulation setup:\")\n",
+    "print(f\"  Duration: {duration} seconds\")\n",
+    "print(f\"  Sampling rate: {fs} Hz\")\n",
+    "print(f\"  Number of samples: {nsamples}\")\n",
+    "print(f\"  Number of instruments: {ninst}\")\n",
+    "print(f\"  Water depth: {depth} m\")\n",
+    "print(f\"\\nInstrument positions (x, y, z):\")\n",
+    "for i in range(ninst):\n",
+    "    print(f\"  Sensor {i+1}: ({layout[0,i]:.1f}, {layout[1,i]:.1f}, {layout[2,i]:.1f}) m\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Generate synthetic wave data\n",
+    "# Simple model: superposition of wave components\n",
+    "\n",
+    "# Create several wave components with different frequencies and directions\n",
+    "frequencies = np.array([0.08, 0.10, 0.12, 0.15])  # Hz\n",
+    "directions = np.array([np.pi/2, np.pi/2, np.pi/2 + np.pi/6, np.pi/3])  # radians\n",
+    "amplitudes = np.array([0.5, 1.0, 0.3, 0.2])  # meters\n",
+    "\n",
+    "# Calculate wavenumbers using dispersion relation: omega^2 = g*k*tanh(k*h)\n",
+    "g = 9.81  # gravity (m/s^2)\n",
+    "wavenumbers = []\n",
+    "for f in frequencies:\n",
+    "    omega = 2 * np.pi * f\n",
+    "    # Solve iteratively for k\n",
+    "    k = omega**2 / g  # deep water approximation as initial guess\n",
+    "    for _ in range(10):\n",
+    "        k = omega**2 / (g * np.tanh(k * depth))\n",
+    "    wavenumbers.append(k)\n",
+    "wavenumbers = np.array(wavenumbers)\n",
+    "\n",
+    "# Initialize data array\n",
+    "data = np.zeros((nsamples, ninst))\n",
+    "\n",
+    "# Generate pressure signal at each instrument location\n",
+    "for i in range(len(frequencies)):\n",
+    "    f = frequencies[i]\n",
+    "    theta = directions[i]\n",
+    "    A = amplitudes[i]\n",
+    "    k = wavenumbers[i]\n",
+    "    omega = 2 * np.pi * f\n",
+    "    \n",
+    "    # Phase at each instrument location\n",
+    "    for j in range(ninst):\n",
+    "        x = layout[0, j]\n",
+    "        y = layout[1, j]\n",
+    "        z = layout[2, j]\n",
+    "        \n",
+    "        # Spatial phase\n",
+    "        phase_xy = k * (x * np.cos(theta) + y * np.sin(theta))\n",
+    "        \n",
+    "        # Pressure response (includes depth attenuation)\n",
+    "        pressure_response = np.cosh(k * (z + depth)) / np.cosh(k * depth)\n",
+    "        \n",
+    "        # Add wave component\n",
+    "        random_phase = np.random.uniform(0, 2*np.pi)\n",
+    "        data[:, j] += A * pressure_response * np.cos(omega * t - phase_xy + random_phase)\n",
+    "\n",
+    "# Add some noise\n",
+    "noise_level = 0.05  # 5% noise\n",
+    "data += noise_level * np.random.randn(nsamples, ninst)\n",
+    "\n",
+    "print(f\"\\nGenerated wave data shape: {data.shape}\")\n",
+    "print(f\"Data range: [{data.min():.2f}, {data.max():.2f}] (pressure units)\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Visualize the Generated Data"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Plot time series from each instrument\n",
+    "fig, axes = plt.subplots(ninst, 1, figsize=(12, 6), sharex=True)\n",
+    "for i in range(ninst):\n",
+    "    axes[i].plot(t[:500], data[:500, i])\n",
+    "    axes[i].set_ylabel(f'Sensor {i+1}\\n(pressure)')\n",
+    "    axes[i].grid(True, alpha=0.3)\n",
+    "\n",
+    "axes[-1].set_xlabel('Time (seconds)')\n",
+    "fig.suptitle('Synthetic Wave Measurements (first 250 seconds)', fontsize=14)\n",
+    "plt.tight_layout()\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Step 2: Set Up the Instrument Data Structure\n",
+    "\n",
+    "The Instrument Data (ID) structure contains all information about the measurement array:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Create instrument data structure\n",
+    "ID = {\n",
+    "    'layout': layout,\n",
+    "    'datatypes': ['pres', 'pres', 'pres'],  # all pressure sensors\n",
+    "    'depth': depth,\n",
+    "    'fs': fs,\n",
+    "    'data': data\n",
+    "}\n",
+    "\n",
+    "print(\"Instrument Data (ID) structure created:\")\n",
+    "print(f\"  Layout shape: {ID['layout'].shape}\")\n",
+    "print(f\"  Data types: {ID['datatypes']}\")\n",
+    "print(f\"  Depth: {ID['depth']} m\")\n",
+    "print(f\"  Sampling frequency: {ID['fs']} Hz\")\n",
+    "print(f\"  Data shape: {ID['data'].shape}\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Step 3: Define the Spectral Matrix Structure\n",
+    "\n",
+    "The Spectral Matrix (SM) structure defines the frequency and directional resolution of the output:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Define output spectral matrix\n",
+    "SM = {\n",
+    "    'freqs': np.linspace(0.05, 0.3, 30),  # frequency bins (Hz)\n",
+    "    'dirs': np.linspace(-np.pi, np.pi, 36)  # direction bins (radians)\n",
+    "}\n",
+    "\n",
+    "print(\"Spectral Matrix (SM) structure created:\")\n",
+    "print(f\"  Frequency range: {SM['freqs'][0]:.3f} - {SM['freqs'][-1]:.3f} Hz\")\n",
+    "print(f\"  Number of frequency bins: {len(SM['freqs'])}\")\n",
+    "print(f\"  Direction range: {SM['dirs'][0]:.2f} - {SM['dirs'][-1]:.2f} radians\")\n",
+    "print(f\"  Number of direction bins: {len(SM['dirs'])}\")\n",
+    "print(f\"  Directional resolution: {np.degrees(SM['dirs'][1] - SM['dirs'][0]):.1f} degrees\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Step 4: Set Estimation Parameters\n",
+    "\n",
+    "The Estimation Parameters (EP) structure controls the analysis method:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Define estimation parameters\n",
+    "EP = {\n",
+    "    'method': 'IMLM',  # Iterated Maximum Likelihood Method\n",
+    "    'iter': 100,       # number of iterations\n",
+    "    'smooth': 'ON'     # enable spectral smoothing\n",
+    "}\n",
+    "\n",
+    "print(\"Estimation Parameters (EP) structure created:\")\n",
+    "print(f\"  Method: {EP['method']}\")\n",
+    "print(f\"  Iterations: {EP['iter']}\")\n",
+    "print(f\"  Smoothing: {EP['smooth']}\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Step 5: Compute the Directional Spectrum\n",
+    "\n",
+    "Now we can run the main analysis function:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Compute directional spectrum\n",
+    "print(\"Computing directional spectrum...\")\n",
+    "print(\"=\" * 50)\n",
+    "\n",
+    "# Set plot type to 0 (no automatic plotting) so we can customize\n",
+    "options = ['PLOTTYPE', 0]\n",
+    "\n",
+    "SMout, EPout = dirspec(ID, SM, EP, options)\n",
+    "\n",
+    "print(\"=\" * 50)\n",
+    "print(\"\\nDirectional spectrum computed successfully!\")\n",
+    "print(f\"Output spectrum shape: {SMout['S'].shape}\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Step 6: Analyze Results\n",
+    "\n",
+    "Use `infospec` to get key wave statistics:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Get wave information\n",
+    "Hsig, Tp, DTp, Dp = infospec(SMout)\n",
+    "\n",
+    "print(\"\\n\" + \"=\"*50)\n",
+    "print(\"WAVE ANALYSIS RESULTS\")\n",
+    "print(\"=\"*50)\n",
+    "print(f\"Significant Wave Height (Hsig): {Hsig:.2f} m\")\n",
+    "print(f\"Peak Period (Tp): {Tp:.2f} s\")\n",
+    "print(f\"Direction at Peak (DTp): {DTp:.2f} rad = {np.degrees(DTp):.1f}\u00b0\")\n",
+    "print(f\"Dominant Direction (Dp): {Dp:.2f} rad = {np.degrees(Dp):.1f}\u00b0\")\n",
+    "print(\"=\"*50)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Step 7: Visualize the Directional Spectrum\n",
+    "\n",
+    "Create custom visualizations of the results:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Create a 1D frequency spectrum (integrated over directions)\n",
+    "freq_spectrum = np.sum(np.real(SMout['S']), axis=1) * (SMout['dirs'][1] - SMout['dirs'][0])\n",
+    "\n",
+    "# Create figure with subplots\n",
+    "fig, axes = plt.subplots(1, 2, figsize=(14, 5))\n",
+    "\n",
+    "# Plot 1D frequency spectrum\n",
+    "axes[0].plot(SMout['freqs'], freq_spectrum, 'b-', linewidth=2)\n",
+    "axes[0].axvline(1/Tp, color='r', linestyle='--', label=f'Peak Period = {Tp:.1f}s')\n",
+    "axes[0].set_xlabel('Frequency (Hz)', fontsize=12)\n",
+    "axes[0].set_ylabel('Spectral Density (m\u00b2/Hz)', fontsize=12)\n",
+    "axes[0].set_title('Frequency Spectrum (Directionally Integrated)', fontsize=13)\n",
+    "axes[0].grid(True, alpha=0.3)\n",
+    "axes[0].legend()\n",
+    "\n",
+    "# Plot 1D directional distribution (integrated over frequencies)\n",
+    "dir_spectrum = np.sum(np.real(SMout['S']), axis=0) * (SMout['freqs'][1] - SMout['freqs'][0])\n",
+    "axes[1].plot(np.degrees(SMout['dirs']), dir_spectrum, 'g-', linewidth=2)\n",
+    "axes[1].axvline(np.degrees(Dp), color='r', linestyle='--', label=f'Dominant Dir = {np.degrees(Dp):.1f}\u00b0')\n",
+    "axes[1].set_xlabel('Direction (degrees)', fontsize=12)\n",
+    "axes[1].set_ylabel('Spectral Density (m\u00b2/degree)', fontsize=12)\n",
+    "axes[1].set_title('Directional Distribution (Frequency Integrated)', fontsize=13)\n",
+    "axes[1].grid(True, alpha=0.3)\n",
+    "axes[1].legend()\n",
+    "\n",
+    "plt.tight_layout()\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Create 2D contour plot of directional spectrum\n",
+    "fig, ax = plt.subplots(figsize=(10, 8))\n",
+    "\n",
+    "# Convert directions to degrees for better readability\n",
+    "dirs_deg = np.degrees(SMout['dirs'])\n",
+    "\n",
+    "# Create meshgrid\n",
+    "F, D = np.meshgrid(SMout['freqs'], dirs_deg)\n",
+    "\n",
+    "# Plot filled contours\n",
+    "levels = 20\n",
+    "contourf = ax.contourf(F, D, np.real(SMout['S']).T, levels=levels, cmap='viridis')\n",
+    "contour = ax.contour(F, D, np.real(SMout['S']).T, levels=levels, colors='k', alpha=0.3, linewidths=0.5)\n",
+    "\n",
+    "# Add colorbar\n",
+    "cbar = plt.colorbar(contourf, ax=ax, label='Spectral Density (m\u00b2s/deg)')\n",
+    "\n",
+    "# Mark peak\n",
+    "ax.plot(1/Tp, np.degrees(DTp), 'r*', markersize=15, label=f'Peak: T={Tp:.1f}s, \u03b8={np.degrees(DTp):.1f}\u00b0')\n",
+    "\n",
+    "ax.set_xlabel('Frequency (Hz)', fontsize=12)\n",
+    "ax.set_ylabel('Direction (degrees)', fontsize=12)\n",
+    "ax.set_title('2D Directional Wave Spectrum', fontsize=14, fontweight='bold')\n",
+    "ax.legend(loc='upper right')\n",
+    "ax.grid(True, alpha=0.3)\n",
+    "\n",
+    "plt.tight_layout()\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Step 8: Use Built-in Plotting Functions\n",
+    "\n",
+    "pyDIWASP includes built-in plotting functions for standard visualizations:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# 3D surface plot\n",
+    "plotspec(SMout, ptype=1)\n",
+    "plt.title('3D Surface Plot - Cartesian Directions', fontsize=14, fontweight='bold')\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Polar contour plot\n",
+    "plotspec(SMout, ptype=2)\n",
+    "plt.title('Polar Contour Plot - Cartesian Directions', fontsize=14, fontweight='bold')\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Example 2: Comparing Different Estimation Methods\n",
+    "\n",
+    "pyDIWASP supports several estimation methods. Let's compare two common ones:\n",
+    "- **IMLM**: Iterated Maximum Likelihood Method (default)\n",
+    "- **EMEP**: Extended Maximum Entropy Principle"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Compute spectrum using EMEP method\n",
+    "EP_emep = {\n",
+    "    'method': 'EMEP',\n",
+    "    'iter': 100,\n",
+    "    'smooth': 'ON'\n",
+    "}\n",
+    "\n",
+    "print(\"Computing directional spectrum using EMEP method...\")\n",
+    "print(\"=\" * 50)\n",
+    "SMout_emep, _ = dirspec(ID, SM, EP_emep, ['PLOTTYPE', 0])\n",
+    "print(\"=\" * 50)\n",
+    "print(\"Complete!\\n\")\n",
+    "\n",
+    "# Get statistics for EMEP\n",
+    "Hsig_emep, Tp_emep, DTp_emep, Dp_emep = infospec(SMout_emep)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Compare the two methods\n",
+    "fig, axes = plt.subplots(1, 2, figsize=(14, 6))\n",
+    "\n",
+    "# IMLM method\n",
+    "F, D = np.meshgrid(SMout['freqs'], np.degrees(SMout['dirs']))\n",
+    "c1 = axes[0].contourf(F, D, np.real(SMout['S']).T, levels=15, cmap='viridis')\n",
+    "axes[0].plot(1/Tp, np.degrees(DTp), 'r*', markersize=12, label=f'Peak')\n",
+    "axes[0].set_xlabel('Frequency (Hz)')\n",
+    "axes[0].set_ylabel('Direction (degrees)')\n",
+    "axes[0].set_title(f'IMLM Method\\nHsig={Hsig:.2f}m, Tp={Tp:.1f}s', fontweight='bold')\n",
+    "plt.colorbar(c1, ax=axes[0], label='Density')\n",
+    "axes[0].legend()\n",
+    "axes[0].grid(True, alpha=0.3)\n",
+    "\n",
+    "# EMEP method\n",
+    "c2 = axes[1].contourf(F, D, np.real(SMout_emep['S']).T, levels=15, cmap='viridis')\n",
+    "axes[1].plot(1/Tp_emep, np.degrees(DTp_emep), 'r*', markersize=12, label=f'Peak')\n",
+    "axes[1].set_xlabel('Frequency (Hz)')\n",
+    "axes[1].set_ylabel('Direction (degrees)')\n",
+    "axes[1].set_title(f'EMEP Method\\nHsig={Hsig_emep:.2f}m, Tp={Tp_emep:.1f}s', fontweight='bold')\n",
+    "plt.colorbar(c2, ax=axes[1], label='Density')\n",
+    "axes[1].legend()\n",
+    "axes[1].grid(True, alpha=0.3)\n",
+    "\n",
+    "plt.tight_layout()\n",
+    "plt.show()\n",
+    "\n",
+    "print(\"\\nComparison of Results:\")\n",
+    "print(\"=\" * 50)\n",
+    "print(f\"{'Parameter':<20} {'IMLM':<15} {'EMEP':<15}\")\n",
+    "print(\"=\" * 50)\n",
+    "print(f\"{'Hsig (m)':<20} {Hsig:<15.2f} {Hsig_emep:<15.2f}\")\n",
+    "print(f\"{'Tp (s)':<20} {Tp:<15.2f} {Tp_emep:<15.2f}\")\n",
+    "print(f\"{'DTp (degrees)':<20} {np.degrees(DTp):<15.1f} {np.degrees(DTp_emep):<15.1f}\")\n",
+    "print(f\"{'Dp (degrees)':<20} {np.degrees(Dp):<15.1f} {np.degrees(Dp_emep):<15.1f}\")\n",
+    "print(\"=\" * 50)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Summary\n",
+    "\n",
+    "This notebook demonstrated:\n",
+    "\n",
+    "1. **Setting up input data structures** (ID, SM, EP)\n",
+    "2. **Running directional spectrum analysis** using `dirspec()`\n",
+    "3. **Extracting wave statistics** using `infospec()`\n",
+    "4. **Visualizing results** with both custom plots and built-in functions\n",
+    "5. **Comparing different estimation methods** (IMLM vs EMEP)\n",
+    "\n",
+    "### Key Takeaways\n",
+    "\n",
+    "- pyDIWASP requires three main input structures: ID (instrument data), SM (spectral matrix), and EP (estimation parameters)\n",
+    "- The IMLM method is the default and works well for most applications\n",
+    "- Different estimation methods may give slightly different results\n",
+    "- The output includes both the full 2D spectrum and integrated statistics\n",
+    "\n",
+    "### Next Steps\n",
+    "\n",
+    "To use pyDIWASP with your own data:\n",
+    "1. Load your wave measurement data\n",
+    "2. Set up the ID structure with your instrument configuration\n",
+    "3. Define appropriate frequency and direction ranges in SM\n",
+    "4. Choose an estimation method in EP\n",
+    "5. Run `dirspec()` and analyze the results\n",
+    "\n",
+    "For more information, see the [pyDIWASP README](../README.md)."
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.8.0"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}
\ No newline at end of file