% Generates a solution of the Einstein's equation for the case studied
% in the Belinski-Zakharov transform
%
% The grid size is [Time]x[Space] with steps [dt]x[dz]
% PlotType can be either
%     '3D' - to obtain a numerical solution of the equation
%     'Exact' - to obtain the one-soliton solution using the
%     Belinski-Zakharov transform
%     'Error' - to see the difference between the numerical and the
%     analytic solutions
% The variable [Periodic] is a boolean that denotes whether we use periodic
% boundary conditions. Note that one can change the initial conditions
% inside the function.

function SimulateOneSoliton(Time, Space, dt, dz, PlotType, Periodic)
    fig = figure;
    Dimensions = 2; % Number of spatial + time dimensions

    t0 = 10;
    z0 = 80;
    
    GridSizeZ = round(Space / dz);
    GridSizeT = round(Time / dt);
    
%    if (GridSizeT > GridSizeZ/2)
%        sprintf('Warning: Time parameter is larger than the simulated triangle height. Setting GridSizeT=GridSizeZ.')
%        GridSizeT = GridSizeZ / 2;
%        Time = GridSizeT * dt;
%    end
    
    z = (z0 + dz):dz:(z0 + Space);
    t = (t0 + dt):dt:(t0 + Time);
    
    epsilon = 0.01;  % Our small parameter
    g = zeros(GridSizeT, GridSizeZ, Dimensions, Dimensions);  % The metric tensor, i.e. g_{ij}
    gExact = zeros(GridSizeT, GridSizeZ, Dimensions, Dimensions);  % The exact solution for the metric tensor, i.e. g_{ij}
    f = ones(GridSizeT, GridSizeZ); % The g_33 and -g_00 term in the metric tensor.
    eta = zeros(GridSizeT, GridSizeZ, Dimensions, Dimensions); % The zero order term of the metric tensor
    Alpha = zeros(GridSizeT, GridSizeZ); % Our alpha parameter (alpha^2 = det(g))
    A = zeros(Dimensions, Dimensions);
    B = zeros(Dimensions, Dimensions);
    gzp = zeros(Dimensions, Dimensions); % The metric tensor at the next spatial element
    gzm = zeros(Dimensions, Dimensions); % The metric tensor at the previous spatial element
    gCurrent = zeros(Dimensions, Dimensions); % The current metric tensor
    gInverse = zeros(Dimensions, Dimensions); % The inverse of the current metric tensor
    gtp = zeros(Dimensions, Dimensions); % The metric tensor at the next time moment
    gtm = zeros(Dimensions, Dimensions); % The metric tensor at the previous time moment
    s1 = 1; s2 = 0; C1 = 1; C2 = 0; % The constants of integration for the exact solution
    MetricIsSingular = false;
    MaxError = 0; % A measure of the maximum deviation of g_11 from the exact solution gExact_11
    
    if (dz < dt)
        dt = dz;
        sprintf('Warning: Initial dz was smaller than dz. Setting dt=dz')
    end
    
    function ComputeMetric(tCoordinate)
        for zCoordinate=1:GridSizeZ
            if (g(tCoordinate,zCoordinate,1,1) * g(tCoordinate,zCoordinate,2,2) - g(tCoordinate,zCoordinate,1,2) * g(tCoordinate,zCoordinate,2,1) <= 0)
                MetricIsSingular = true;
                break;
            end
            
            Alpha(tCoordinate,zCoordinate) = sqrt(g(tCoordinate,zCoordinate,1,1) * g(tCoordinate,zCoordinate,2,2) - g(tCoordinate,zCoordinate,1,2) * g(tCoordinate,zCoordinate,2,1));
        end
    end

    function ComputeAB(tCoordinate,zCoordinate)
        DgDt(:,:) = (g(tCoordinate+1,zCoordinate,:,:) - g(tCoordinate-1,zCoordinate,:,:))/(2*dt);

        zNext = zCoordinate + 1;
        if (zNext > GridSizeZ)
            zNext = zNext - GridSizeZ;
        end

        zPrevious = zCoordinate - 1;
        if (zPrevious < 1)
            zPrevious = zPrevious + GridSizeZ;
        end

        DgDz(:,:) = (g(tCoordinate,zNext,:,:) - g(tCoordinate,zPrevious,:,:))/(2*dz);

        invg = zeros(Dimensions, Dimensions);
        a = g(tCoordinate,zCoordinate,1,1);
        b = g(tCoordinate,zCoordinate,1,2);
        c = g(tCoordinate,zCoordinate,2,1);
        d = g(tCoordinate,zCoordinate,2,2);
        gdet = (Alpha(tCoordinate,zCoordinate))^2;

        invg(1,1) =  d/gdet; invg(1,2) = -b/gdet;
        invg(2,1) = -c/gdet; invg(2,2) =  a/gdet;

        A = -Alpha(tCoordinate,zCoordinate) * (DgDz + DgDt) * invg;
        B = Alpha(tCoordinate,zCoordinate) * (DgDz - DgDt) * invg;
    end

    % We first find the exact solution, just to compare it with ours later
    % on.
    for i=1:GridSizeT % Corresponds to t
        for j=1:GridSizeZ % Corresponds to z
            if (z(j) > t(i))
                r = log(2*z(j)^2/t(i)^2 - 1 - 2*(z(j)^2/t(i)^2 * (z(j)^2/t(i)^2 - 1))^0.5);
            else
                r = 0;
            end

            gExact(i,j,1,1) = cosh(s1 * r + C2)/cosh((s1-1/2) * r + C2) * t(i)^(2*s1);
            gExact(i,j,2,2) = cosh(s2 * r - C2)/cosh((1/2-s2) * r + C2) * t(i)^(2*s2);
            gExact(i,j,1,2) = -sinh(r/2)/cosh((s1-1/2)*r + C2)*t(i);
            gExact(i,j,2,1) = -sinh(r/2)/cosh((s1-1/2)*r + C2)*t(i);
        end
    end

    % Next we fix the initial conditions that correspond to a simple
    % soliton solution
    f(1, :) = 1;
    f(2, :) = 1;

    for j=1:GridSizeZ % Corresponds to z
        g(1,j,1,1) = gExact(1,j,1,1);
        g(1,j,2,2) = gExact(1,j,2,2);
        g(1,j,1,2) = gExact(1,j,1,2);
        g(1,j,2,1) = gExact(1,j,2,1);
        
        g(2,j,1,1) = gExact(2,j,1,1);
        g(2,j,2,2) = gExact(2,j,2,2);
        g(2,j,1,2) = gExact(2,j,1,2);
        g(2,j,2,1) = gExact(2,j,2,1);
    end
    
    ComputeMetric(1);

    if (strcmp(PlotType, '2D'))
        subplot(2,2,1);
        g11plot = plot(z, g(1,:,1,1));
        axis([-pi pi 0 2]);
        xlabel('z (Space)');
        ylabel('g_{11}');

        subplot(2,2,2);
        g12plot = plot(z, g(1,:,1,2));
        axis([-pi pi 0 2]);
        xlabel('z (Space)');
        ylabel('g_{12}');

        subplot(2,2,3);
        g21plot = plot(z, g(1,:,2,1));
        axis([-pi pi 0 2]);
        xlabel('z (Space)');
        ylabel('g_{21}');

        subplot(2,2,4);
        g22plot = plot(z, g(1,:,2,2));
        axis([-pi pi 0 2]);
        xlabel('z (Space)');
        ylabel('g_{22}');

        set(fig,'userdata',dt);
    end
    
    if (~strcmp(PlotType, 'Exact'))
        for i=2:(GridSizeT-1) % Corresponds to the time t
            if (strcmp(PlotType,'2D'))
                dt = get(fig,'userdata');
            end

            if (~strcmp(PlotType, 'Exact'))
                ComputeMetric(i);
            end

            if MetricIsSingular
                t(i)-t0
                i = GridSizeT;
                break;
            end

            for j=1:GridSizeZ % Corresponds to z
                % We make sure jp (j plus) and jm (j minus) will be within our
                % grid by making the z-axis periodic
                jp = j + 1;
                if (jp > GridSizeZ)
                    if Periodic
                        jp = jp - GridSizeZ;
                    else
                        jp = j;
                    end
                end

                jm = j - 1;
                if (jm < 1)
                    if Periodic
                        jm = jm + GridSizeZ;
                    else
                        jm = j;
                    end
                end

                gzp = squeeze(g(i,jp,:,:));
                gzm = squeeze(g(i,jm,:,:));
                gCurrent = squeeze(g(i,j,:,:));
                gInverse = inv(gCurrent);
                gtm = squeeze(g(i-1,j,:,:));

                gtp = 2 * gCurrent - gtm;
                gtp = gtp + (dt/dz)^2 * (gzp - 2 * gCurrent + gzm);
                gtp = gtp + (dt/(2*dz))^2 * (Alpha(i,jp) - Alpha(i,jm)) * (gzp - gzm) / Alpha(i,j);
                gtp = gtp - ((Alpha(i,j) - Alpha(i-1,j)) * (gCurrent - gtm)) / Alpha(i,j);
                gtp = gtp + (dt/(2*dz))^2 * (gzp - gzm) * gInverse * (gzp - gzm);
                gtp = gtp - (gCurrent - gtm) * gInverse * (gCurrent - gtm);
                g(i+1,j,:,:) = gtp;
                
                if ((abs(g(i+1,j,1,1) - gExact(i+1,j,1,1)) > MaxError) && (i+1 <= j) && (i+1 <= (GridSizeZ - j + 1)))
                    MaxError = abs(g(i+1,j,1,1) - gExact(i+1,j,1,1));
                end
            end

            if (strcmp(PlotType, '2D'))
                set(g11plot,'ydata',g(i,:,1,1));
                set(g12plot,'ydata',g(i,:,1,2));
                set(g21plot,'ydata',g(i,:,2,1));
                set(g22plot,'ydata',g(i,:,2,2));

                drawnow
            end
        end
    end

    if (~strcmp(PlotType, 'Exact'))
        ComputeMetric(GridSizeT);
    end
        
    if (~strcmp(PlotType, '2D'))
        subplot(2,2,1);
        if (strcmp(PlotType, '3D'))
            surf(z,t,g(:,:,1,1), 'FaceColor', 'interp', 'EdgeColor', 'none', 'FaceLighting', 'phong');
        elseif (strcmp(PlotType, 'Exact'))
            surf(z,t,gExact(:,:,1,1), 'FaceColor', 'interp', 'EdgeColor', 'none', 'FaceLighting', 'phong');
        elseif (strcmp(PlotType, 'Error'))
            surf(z,t,(g(:,:,1,1)-gExact(:,:,1,1))./gExact(:,:,1,1), 'FaceColor', 'interp', 'EdgeColor', 'none', 'FaceLighting', 'phong');
        end
        xlabel('z (Space)');
        ylabel('t (Time)');
        zlabel('g_{11}');
        colorbar;
        subplot(2,2,2);
        if (strcmp(PlotType, '3D'))
            surf(z,t,g(:,:,1,2), 'FaceColor', 'interp', 'EdgeColor', 'none', 'FaceLighting', 'phong');
        elseif (strcmp(PlotType, 'Exact'))
            surf(z,t,gExact(:,:,1,2), 'FaceColor', 'interp', 'EdgeColor', 'none', 'FaceLighting', 'phong');
        elseif (strcmp(PlotType, 'Error'))
            surf(z,t,(g(:,:,1,2)-gExact(:,:,1,2))./gExact(:,:,1,2), 'FaceColor', 'interp', 'EdgeColor', 'none', 'FaceLighting', 'phong');
        end
        xlabel('z (Space)');
        ylabel('t (Time)');
        zlabel('g_{12}');
        colorbar;
        subplot(2,2,3);
        if (strcmp(PlotType, '3D'))
            surf(z,t,g(:,:,2,1), 'FaceColor', 'interp', 'EdgeColor', 'none', 'FaceLighting', 'phong');
        elseif (strcmp(PlotType, 'Exact'))
            surf(z,t,gExact(:,:,2,1), 'FaceColor', 'interp', 'EdgeColor', 'none', 'FaceLighting', 'phong');
        elseif (strcmp(PlotType, 'Error'))
            surf(z,t,(g(:,:,2,1)-gExact(:,:,2,1))./gExact(:,:,2,1), 'FaceColor', 'interp', 'EdgeColor', 'none', 'FaceLighting', 'phong');
        end
        xlabel('z (Space)');
        ylabel('t (Time)');
        zlabel('g_{21}');
        colorbar;
        subplot(2,2,4);
        if (strcmp(PlotType, '3D'))
            surf(z,t,g(:,:,2,2), 'FaceColor', 'interp', 'EdgeColor', 'none', 'FaceLighting', 'phong');
        elseif (strcmp(PlotType, 'Exact'))
            surf(z,t,gExact(:,:,2,2), 'FaceColor', 'interp', 'EdgeColor', 'none', 'FaceLighting', 'phong');
        elseif (strcmp(PlotType, 'Error'))
            surf(z,t,(g(:,:,2,2)-gExact(:,:,2,2))./gExact(:,:,2,2), 'FaceColor', 'interp', 'EdgeColor', 'none', 'FaceLighting', 'phong');
        end
        xlabel('z (Space)');
        ylabel('t (Time)');
        zlabel('g_{22}');
        colorbar;
    end
    
    if (strcmp(PlotType, 'Error'))
        MaxError
    end
end
