Add user documentation and example notebook

.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
+

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:

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: