diff --git a/examples/environmental_contours_example.mlx b/examples/environmental_contours_example.mlx
index 47c84caf..6bcfa6fd 100644
Binary files a/examples/environmental_contours_example.mlx and b/examples/environmental_contours_example.mlx differ
diff --git a/mhkit/wave/resource/environmental_contours.m b/mhkit/wave/resource/environmental_contours.m
index 779c3118..5c7fd86e 100644
--- a/mhkit/wave/resource/environmental_contours.m
+++ b/mhkit/wave/resource/environmental_contours.m
@@ -34,24 +34,45 @@
 %         Desired return period (years) for calculation of environmental
 %         contour, can be a scalar or a vector.
 %
-%     PCA: Structure (optional)
-% 	      principal component analysis dictionary from previous function
-%         call. When supplied the function will skip the PCA calculation
-%         for the passe x1, and x2.
-%         to call: environmental_contour(x1,x2,dt,period,"PCA",PCA)
+%     method: string or list
+%         Copula method to apply. Options include ['PCA','gaussian',
+%         'gumbel', 'clayton', 'rosenblatt', 'nonparametric_gaussian',
+%         'nonparametric_clayton', 'nonparametric_gumbel', 'bivariate_KDE'
+%         'bivariate_KDE_log']
 %
-%     bin_size : double (optional)
-%         Data points in each bin
-%         to call: environmental_contour(x1,x2,dt,period,"bin_size",bin_size)
-%
-%     nb_steps : int (optional)
-%         Discretization of the circle in the normal space used for
-%         IFORM calculation.
-%         to call: environmental_contour(x1,x2,dt,period,"nb_steps",nb_steps)
-%
-%     return_PCA: boolean
-% 	      Default False, if True will retun the PCA dictionary
-%         to call: environmental_contour(x1,x2,dt,period,"return_PCA",return_PCA)
+%     **OPTIONS
+%        min_bin_count: int
+%            Passed to _copula_parameters to sets the minimum number of
+%            bins allowed. Default = 40.
+%        initial_bin_max_val: int, double
+%            Passed to _copula_parameters to set the max value of the
+%            first bin. Default = 1.
+%        bin_val_size: int, double
+%            Passed to _copula_parameters to set the size of each bin
+%            after the initial bin.  Default 0.25.
+%        nb_steps: int
+%            Discretization of the circle in the normal space is used for
+%            copula component calculation. Default nb_steps=1000.
+%        bandwidth:
+%            Must specify bandwidth for bivariate KDE method.
+%            Default = None.
+%        Ndata_bivariate_KDE: int
+%            Must specify bivariate KDE method. Defines the contoured
+%            space from which samples are taken. Default = 100.
+%        max_x1: double
+%            Defines the max value of x1 to discretize the KDE space
+%        max_x2: double
+%            Defines the max value of x2 to discretize the KDE space
+%        PCA: dict
+%            If provided, the principal component analysis (PCA) on x1,
+%            x2 is skipped. The PCA will be the same for a given x1, x2
+%            therefore this step may be skipped if multiple calls to
+%            environmental contours are made for the same x1, x2 pair.
+%            The PCA dict may be obtained by setting return_fit=True when
+%            calling the PCA method.
+%        return_fit: boolean
+%            Will return fitting parameters used for each method passed.
+%            Default False.
 %
 %
 % Returns
@@ -68,9 +89,16 @@
     dt
     period
     method
-    options.PCA = py.None;
-    options.bin_size = 250;
+    options.min_bin_count = 40;
+    options.initial_bin_max_val = 1;
+    options.bin_val_size = 0.25;
     options.nb_steps = 1000;
+    options.bandwidth = py.None;
+    options.Ndata_bivariate_KDE = 100;
+    options.max_x1 = py.None;
+    options.max_x2 = py.None;
+    options.PCA = py.None;
+    options.PCA_bin_size = 250;
     options.return_fit = py.False;
 end
 
@@ -90,8 +118,12 @@
 end
 
 data = py.mhkit.wave.contours.environmental_contours(py.numpy.array(x1),py.numpy.array(x2),...
-    int32(dt),period_py,method,pyargs('PCA',options.PCA,'bin_size',int32(options.bin_size),...
-    'nb_steps',int32(options.nb_steps),'return_fit',options.return_fit));
+    int32(dt),period_py,method,pyargs('min_bin_count',int32(options.min_bin_count),...
+    'initial_bin_max_val',options.initial_bin_max_val,'bin_val_size',options.bin_val_size,...
+    'nb_steps',int32(options.nb_steps),'bandwidth',options.bandwidth,...
+    'Ndata_bivariate_KDE',options.Ndata_bivariate_KDE,'max_x1',options.max_x1,...
+    'max_x2',options.max_x2,'PCA',options.PCA,'PCA_bin_size',int32(options.PCA_bin_size),...
+    'return_fit',options.return_fit));
 
 data_struct = struct(data);