function carriergreenspan()

%% determine wave

sigma = 0:.1:20;
lambda = 0:.1:30;

A = 0.5;
g = 9.81;
h0 = 5;
a = 1/25;
T = 32;

omega = 1;
l0 = T^2*g*a;

phi = nan(length(sigma),length(lambda));
phi_l = nan(size(phi));
phi_s = nan(size(phi));

for is = 1:length(sigma)
    for il = 1:length(lambda)
        s = sigma(is);
        l = lambda(il);
        
        x = -pi:.1:pi;
        y = exp(1i*omega*s*sin(x));
        J0 = 1/(2*pi)*trapz(x,y);
        
        y = exp(1i*(x+omega*s*sin(x)));
        Jm1 = 1/(2*pi)*trapz(x,y);

        phi(is,il) = A*J0*cos(omega*l);
        phi_l(is,il) = -omega*A*J0*sin(omega*l);
        phi_s(is,il) = omega*A*Jm1*cos(omega*l);
    end
end

[L S] = meshgrid(sigma, lambda);

v = 1./S.*phi_s; v(1,:) = 0;
x = real(phi_l./4-S.^2./16-v.^2./2);
eta = real(phi_l./4-v.^2./2);
t = real(L./2-v);

%% plot

close all;

figure;

%for i = 1:size(eta,2)
for it = 0:.1:max(max(t))
    
    neta = nan(size(t,1));
    nx = nan(size(t,1));
    for j = 1:size(t,1)
        neta(j) = interp1(t(j,:), eta(j,:), it);
        nx(j) = interp1(t(j,:), x(j,:), it);
    end
    
    subplot(211);
    plot(100*[-1 1], 100*[-1 1], 'k'); hold on;
    %plot(x(:,i), eta(:,i), '-xr'); hold off;
    plot(nx, neta, '-xr'); hold off;
    set(gca, 'XLim', [-1.2 .3], 'YLim', [-.3 .3]);
    xlabel('x'); ylabel('\eta');
    grid on;
    
    %tt = sqrt(l0./a./g).*[min(t(:,i)) mean(t(:,i)) max(t(:,i)) max(t(:,i))-min(t(:,i))];
    %title(['t = ' num2str(tt)]);
    title(['t = ' num2str(sqrt(l0./a./g)*it)]);
    
    subplot(212);
    plot(1380+100*[-1 1], h0+100*[-1 1]*a, 'k'); hold on;
    %plot(1380+l0*x(:,i), h0+l0*a*eta(:,i), '-xr'); hold off;
    plot(1380+l0*nx, h0+l0*a*neta, '-xr'); hold off;
    set(gca, 'XLim', [500 1500], 'YLim', h0+[-2.5 2.5]);
    xlabel('x^*'); ylabel('\eta^*');
    grid on;
    
    drawnow;
    pause(.1);
end
