WECFullCode.m --- Script reproducing the numerical results of the paper
================================================================

Purpose
    Generates all numerical figures of the paper:
      ResponseCurves3.png      displacement response curves
      peakPeriod1.png          peak time intervals
      z1ExPhaseCurve1.png      phase plane (z1, v1)
      z2ExPhaseCurve1.png      phase plane (z2, v2)
      z1z2ExPhaseCurve.png     phase plane (z1, z2)
      v1v2ExPhaseCurve.png     phase plane (v1, v2)
      DifferentICCurves1.png   different initial data (z1, v1)
      DifferentICCurves2.png   different initial data (z2, v2)
      DiffICSol.png            differences of solutions for different initial data

Requirements
    MATLAB R2019a or later (uses readmatrix).
    No external toolbox is required.

How to run
    1. Place WECFullCode.m in the same directory as the data workbook:
         WEC_Data.xlsx
       This workbook contains one sheet per hydrodynamic dataset, each with the
       header Frequency_Hz in column A and the coefficient in column B:
         AddedMass          added mass [kg]
         RadiationDamping   radiation damping [N*s/m]
         ExcitedForce       wave excitation force per unit wave height [N/m]
    2. Run WECFullCode.m from that directory.

Note on output paths
    The script writes all figures into a subfolder named "figures" of the current
    working directory (it creates the folder if it does not exist). It uses only
    relative paths, so it can be run from any location without modification. To
    place the figures elsewhere, simply edit the variable outdir.

Parameter sources
    The physical parameters are taken from the cited reference:
    "Non-linear dynamic simulations of two-body wave energy converters via
     identification of viscous drag coefficients of different shapes of the
     submerged body based on numerical wave tank CFD simulation"
    The PTO damping coefficient c, the damping exponent alpha and the incident
    wave height H are user-defined values.

Additional remarks
    1. The script also contains a power-frequency sweep section
       ("%% Power-frequency curve"); its output is not included in the paper.
       The conclusions of the paper concern only the existence and uniqueness of
       the periodic solution and global attractivity.
    2. In that section the instantaneous absorbed power is computed as P = F_pto * v_rel:
         PTO damping force      F_pto  = c * |v_rel|^alpha * v_rel      (consistent with myODE)
         Instantaneous power    P_inst = c * |v_rel|^(alpha+2)
       Note: if the damping force were written as c*|v_rel|^alpha*sign(v_rel), the
       power would instead be c*|v_rel|^(alpha+1); the two conventions differ by a
       factor |v_rel| and must not be mixed. The unit of c used in this paper,
       kg*m^(-1/2)*s^(-1/2), corresponds to the former convention.
