function Ass6

  final_time = 5.0;
  num_points_i = 100;
  num_points_j = 100;
  num_cells_i = num_points_i - 1;
  num_cells_j = num_points_j - 1;
  CFL = 0.75;
  R_min = 1.0;
  R_max = 12.0;
  h = 1;
  ux  = 5;
  uy  = 0.0;
  g = 9.81;
  rho = 1000;
  U_free_stream = zeros(1,1,3);
  U_free_stream(1,1,1) = h;
  U_free_stream(1,1,2) = h*ux;
  U_free_stream(1,1,3) = h*uy;

  %Initialize Geometry
  delta_R     = (R_max-R_min)/(num_points_j-1);
  delta_theta = 2.0*pi/(num_points_i-1);
  x     = zeros(num_points_i,num_points_j);
  y     = zeros(num_points_i,num_points_j);
  U     = zeros(num_points_i-1,num_points_j-1, 3);
  areas = zeros(num_points_i-1,num_points_j-1);

  for i = 1:num_points_i
    theta = (i-1)*delta_theta;
    sin_theta = sin(theta);
    cos_theta = cos(theta);
    for j = 1:num_points_j
      R = R_min + (j-1)*delta_R;
      x(i,j) = R*cos_theta;
      y(i,j) = R*sin_theta;
    end
  end

  %Initialize Cell-Centred Data
  for i = 1:num_cells_i
    for j = 1:num_cells_j
        U(i,j,:) = U_free_stream;
        areas(i,j) = compute_cell_area(x(i  , j  ),y(i  , j ), ...
                                       x(i+1, j  ),y(i+1, j ), ...
                                       x(i+1, j+1),y(i+1, j+1), ...
                                       x(i  , j+1),y(i  , j+1));
    end
  end

  %Initialize Time
  time = 0;
  iteration = 0;

  %Time march
  while time < final_time
    %Initialize
    dt = 1.0e6;
    dUdt = zeros(size(U));
    %Compute i-direction flux contributions to dUdt
    for j = 1:num_cells_j
      for i = 1:num_cells_i
        [nx ny edge_length] = unit_normal(x(i,j+1),y(i,j+1),x(i,j),y(i,j));
        if i == 1
          Ul = U(num_cells_i,j,:);
        else
          Ul = U(i-1,j,:);
        end
        Ul_rot = rotate(Ul,nx,ny);
        Ur = U(i,j,:);
        Ur_rot = rotate(Ur,nx,ny);
        [flux_rot fastest_wave] = Local_Lax_Friedrichs(Ul_rot,Ur_rot,g);
        flux = rotate(flux_rot,nx,-ny);
        if i == 1
          dUdt(num_cells_i,j,:) = dUdt(num_cells_i,j,:) - flux*edge_length/areas(num_cells_i,j);
        else
          dUdt(i-1,j,:) = dUdt(i-1,j,:) - flux*edge_length/areas(i-1,j);
        end
        dUdt(i,j,:) = dUdt(i,j,:) + flux*edge_length/areas(i,j);
        cell_width = areas(i,j)/edge_length;
        dt = min(dt, CFL*cell_width/fastest_wave);
      end
    end

    %Compute j-direction Flux Contributions to dUdt
    for j = 1:num_cells_j+1
      for i = 1:num_cells_i
        [nx ny edge_length] = unit_normal(x(i,j),y(i,j),x(i+1,j),y(i+1,j));
        if j <= num_cells_j
          Ur = U(i,j,:);
        else
          Ur = U_free_stream;
        end
        Ur_rot = rotate(Ur,nx,ny);
        if j == 1
          Ul_rot = reflect(Ur_rot);
        else
          Ul = U(i,j-1,:);
          Ul_rot = rotate(Ul,nx,ny);
        end
        [flux_rot fastest_wave] = Local_Lax_Friedrichs(Ul_rot,Ur_rot,g);
        flux = rotate(flux_rot,nx,-ny);
        if j > 1
          dUdt(i,j-1,:) = dUdt(i,j-1,:) - flux*edge_length/areas(i,j-1);
        end
        if j <= num_cells_j
          dUdt(i,j,:) = dUdt(i,j,:) + flux*edge_length/areas(i,j);
          cell_width = areas(i,j)/edge_length;
          dt = min(dt, CFL*cell_width/fastest_wave);
        end
      end
    end

    %Update
    fit = 0;
    fii = 0;
    fjt = 0;
    fij = 0;
    ft = ((fit)^2 + (fjt)^2)^(0.5)
    U = U + dt*dUdt;
    iteration = iteration + 1;
    time = time + dt;
    for i = 1:num_cells_i
      %Unit Normals and Edge Length for Range of i and j = 1
      [nx ny edge_length] = unit_normal(x(i,1),y(i,1),x(i+1,1),y(i+1,1));
      %Force Calculations
      fii = -0.5*nx*rho*g*edge_length*(U(i,1,1))^2;
      fij = -0.5*ny*rho*g*edge_length*(U(i,1,1))^2;
      %Adjusting Total Forces
      fit = fit+fii
      fjt = fjt+fij
    end
    %interpolate_and_plot1(x,y,h_field(U),1);
    %interpolate_and_plot2(x,y,ux_field(U),2);
    %Output Update Every Time Step
    disp(sprintf('timestep %i, time = %0.5e, dt = %0.5e, norm(dUdt[1]) = %f', ...
                 iteration, time, dt, norm(dUdt(:,:,1))));
  end %While Time March
end

%Local Lax-Friedrichs Flux function
function [Flux fastest_wave] = Local_Lax_Friedrichs(Ul_rot,Ur_rot,g)
  fastest_wave = max([max_lambda_x(Ul_rot,g) max_lambda_x(Ur_rot,g)]);
  Fl = Flux_x(Ul_rot,g);
  Fr = Flux_x(Ur_rot,g);
  Flux = 0.5*(Fl+Fr) - 0.5*fastest_wave*(Ur_rot-Ul_rot);
end

%P.D.E. Specific Functions
function Flux_x = Flux_x(U,g)
  uy = U(3)/U(1);
  ux = U(2)/U(1);
  Flux_x = zeros(1,1,3);
  Flux_x(1,1,1) = U(2);
  Flux_x(1,1,2) = U(1)*ux^2+0.5*g*U(1)^2;
  Flux_x(1,1,3) = U(1)*ux*uy;
end

function Urot = rotate(U,nx,ny)
  Urot = U;
  Urot(2) = U(2)*nx+U(3)*ny;
  Urot(3) = -U(2)*ny+U(3)*nx;
end

function Ureflected = reflect(U)
  Ureflected = U;
  Ureflected(2) = -1.0*Ureflected(2);
end

function l = max_lambda_x(U,g)
  h = U(1);
  ux = U(2)/U(1);
  l1 = abs(ux-(g*h)^(0.5));
  l2 = abs(ux+(g*h)^(0.5));
  l = max([l1,l2]);
end

%Geometry Functions
function l = length(x1,y1,x2,y2)
  l = ((x2-x1)^2+(y2-y1)^2)^(0.5);
end

function [nx ny edge_length] = unit_normal(x1,y1,x2,y2)
  edge_length = length(x1,y1,x2,y2);
  nx =  (y2-y1)/edge_length;
  ny = -(x2-x1)/edge_length;
end

function A = compute_cell_area(x1,y1,x2,y2,x3,y3,x4,y4)
  A = 0.5*abs(x1*y2-x2*y1+x2*y3-x3*y2+x3*y4-x4*y3+x4*y1-x1*y4);
end

%Output-Related Functions
function h = h_field(U)
  h = U(:,:,1);
end

function ux = ux_field(U)
  ux = U(:,:,2)./U(:,:,1);
end

function uy = uy_field(U)
  uy = U(:,:,3)./U(:,:,1);
end

function interpolated = interpolate(x,y,W)
  %This Function Does a Crude Interpolation of Data From Cell
    %Centres to Nodes for Plotting.
  [num_i,num_j] = size(x);
  interpolated = zeros(size(x));
  %Interior Nodes are All the Same
  for i = 2:num_i-1
    for j = 2:num_j-1
      interpolated(i,j) = 0.25*( W(i,j) + W(i-1,j) + W(i,j-1) + W(i-1,j-1));
    end
  end
  %j_min and j_max
  for i = 2:num_i-1
      interpolated(i,1)     = 0.5*( W(i,1) + W(i-1,1));
      interpolated(i,num_j) = 0.5*( W(i,num_j-1) + W(i-1,num_j-1));
  end
  %i_min and i_max
  for j = 2:num_j-1
      interpolated(1,j)     = 0.5*( W(1,j) + W(1,j-1));
      interpolated(num_i,j) = 0.5*( W(num_i-1,j) + W(num_i-1,j-1));
  end
  %Four Corners
  interpolated(1,1) = W(1,1);
  interpolated(1,num_j) = W(1,num_j-1);
  interpolated(num_i,1) = W(num_i-1,1);
  interpolated(num_i,num_j) = W(num_i-1,num_j-1);

end

function interpolate_and_plot1(x,y,W,fig_num)
  W_interp = interpolate(x,y,W);
  figure(fig_num)
  contourf(x,y,W_interp);
  axis('square');
  colorbar;
  drawnow;
end

function interpolate_and_plot2(x,y,W,fig_num)
  W_interp = interpolate(x,y,W);
  figure(fig_num)
  contourf(x,y,W_interp);
  axis('square');
  colorbar;
  drawnow;

end