function [ u, v, p, b1, b2 ] = mhd_solution ( x, y )

%*****************************************************************************80
%
%% mhd_solution(): Hartmann exact solution of a steady 2D MHD equation.
%
%  Discussion:
%
%    We consider the domain XxY = [0,L]x[-1,+1].
%    The horizontal velocity u(x,y) is constant along any fixed y value.
%    The vertical velocity v(x,y) is zero everywhere.
%    The magnetic field B satisfies B(x,y) = (b(y),1).
%    The boundary walls y=-1 and y=+1 have no slip boundary conditions:
%      u(x,-1) = u(x,+1) = 0
%    The inlet x=0 and outlet x=L have Neumann boundary conditions:
%      -p n + 1/Re (Del u,n) = - pd n
%    where n is the unit normal vector and the pressure pd is given.
%    The walls are perfectly insulating.
%    On the boundary, we impose the tangential part of the magnetic field:
%      B x n = Bd x n
%    where the external transverse magnetic field Bd is given
%    as Bd(x,y)=(0,1).
%
%      pd(x,y) = - G x - 1/2 S b(y)^2 + p0
%      u''(x,y) + Re S b'(y) = - G Re
%      v(x,y) = 0
%      b''(x,y) + R m u'(x,y) = 0
%
%    The system to be solved is
%
%      -1/Re del^2 u + del p + u dot del u + s B cross curl B = f
%      div u = 0
%      1/Rm curl ( curl ( B ) ) - curl ( u cross B ) = h
%      div B = 0
%
%    with homogeneous boundary conditions
%
%      1/Rm curl B cross n = 0
%      B dot n = 0
%      u = 0
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    23 August 2024
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Jean-Frederic Gerbeau, Claude Le Bris, Tony Lelievre,
%    Mathematical methods for the magnetohydrodynamics
%    of liquid metals,
%    Oxford University Press, 2006,
%    ISBN 978-0-19-856665-6
%
%  Input:
%
%    real X(:), Y(:): the spatial locations.
%
%  Output:
%
%    real U(:), V(:), P(:), B1(:), B2(:): the solution components
%    (velocity, pressure, magnetic field).
%
  [ m, n ] = size ( x );

  [ G, Ha, L, p0, Re, Rm, S ] = mhd_parameters ( );

  u = G * Re / Ha / tanh ( Ha ) ...
    * ( 1.0 - cosh ( y * Ha ) / cosh ( Ha ) );
  v = zeros ( m, n );
  b = G / S * ( sinh ( y * Ha ) / sinh ( Ha ) - y );
  b1 = b;
  b2 = ones ( m, n );
  p = - G * x - 0.5 * S * b.^2;

  return
end

