function [par_fit,fit_func,f_fine] = fit_JONSWAP(f,E,par_0)
% This function finds the Pierson-Moskowitz spectrum that minimizes the
% mean squared error between the observed data (f,E) and the model. This is
% accompished by finding that parameter fp that minimizes the erro.
    par_fit = fminsearch(@compSSE,par_0);
    df = f(2)-f(1);
    f_fine = min(f):df/10:max(f);
    fit_func = JS_spec(f_fine,par_fit);
    function mse = compSSE(par)
%         if sum(par>0)==5
%             E_fit = JS_spec(f,par);
%             mse = sum((E-E_fit).^2);
%         else
%             mse = Inf;
%         end
            E_fit = JS_spec(f,par);
            mse = sum((E-E_fit).^2);
%         if sum(par>0)<5 || par(3)<0.03
%             mse = 1.1*mse;
%         end

    end
end

function E_J = JS_spec(f,par) % JONSWAP spectral density
    alpha   = par(1);
    fp      = par(2);
    sigma_a = par(3);
    sigma_b = par(4);
    gamma   = par(5);
    g = 9.8;
        
    top = f>fp;  % one where f> fp zero otherwise;
    bot = ~top;  % one where f<=fp zero otherwise;
    sigma = sigma_a*bot + sigma_b*top;
    u   = f/fp;
    PM  = alpha*(g/(2*pi)^2)^2*f.^-5.*exp(-1.25*u.^-4);
    G   = gamma.^(exp(-0.5*((u-1)./sigma).^2));
    E_J = PM.*G;
end

