% Nonlinear sound equations in 1D. PDE is represented as 2N coupled ODEs
% on an N-element spatial grid with periodic boundaries.
function dqdt = nonlinear(q, t, Delta, c)
% q     = state vector. Should be a row vector, the first half being
%         the dimensionless velocity v/c, and the second half the 
%         dimensionless density (rho-<rho>)/<rho>.
% t     = time. Not used. Provided for compatibility with ODE solvers.
% Delta = grid spacing.
% c     = wave speed.
%
% To set up initial conditions: velocity = density gives a wave traveling
% in the +x direction, while velocity = -density gives a wave traveling in the
% -x direction.

N = numel(q)/2; % Number of grid points
gamma = 5/3; % Monatomic gas

% We implement periodic boundary conditions by wrapping around both 
% ends, so each of our two variables has N+2 elements.
velocity = q(1:N);
density = q((N+1):(2*N));
vel = [ velocity(N), velocity, velocity(1) ]; % with wrapped ends
den = [  density(N),  density,  density(1) ]; % with wrapped ends

% Calculate spatial derivatives:
dx_vel = ( den(3:(N+2)) - den(1:N) ) / (2*Delta);
dx_den = ( vel(3:(N+2)) - vel(1:N) ) / (2*Delta);
   % Note how the periodic BCs work --- derivs have N elements.

% PDEs:
dt_vel = -c * ( velocity.*dx_vel + dx_den.*(1+density).^gamma );
dt_den = -c * ( dx_vel.*(1+density) + velocity.*dx_den ) ;

dqdt = [dt_vel, dt_den];

