function f = fxyz1 ( n, x, y, z )

%*****************************************************************************80
%
%% FXYZ1 is the first 3D example.
%
%  Discussion:
%
%    This example is known as the 3D Shepp-Logan phantom.
%
%    It should be plotted on [-1,+1] x [-1,+1.5] x [-1.5,+1.5].
%
%    Seventeen objects are modeled by ellipses of various gray levels,
%    including:
%
%     1: Outer skull
%     2: Inner skull
%     3: Left eye
%     4: Right eye
%     5: Nose
%     6: Mouth
%     7: Left ear
%     8: Right ear
%     9: Left small tumor
%    10: Center small tumor
%    11: Right small tumor
%    12: Old f
%    13: Old g
%    14: Old e
%    15: Right ventricle
%    16: Left ventricle
%    17: Blood clot
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    21 September 2014
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Larry Shepp,
%    Computerized tomography and nuclear magnetic resonance,
%    Journal of Computer Assisted Tomography,
%    Volume 4, Number 1, February 1980, pages 94-107.
%
%  Parameters:
%
%    Input, integer N, the number of points.
%
%    Input, real X(N), Y(N), Z(N), the arguments.
%
%    Output, real F(N), the function values.
%
  x0 = [  0.0000,  0.0000,  0.2583, -0.2583,  0.0000, ...
          0.0000,  0.7076, -0.7076, -0.0800,  0.0000, ...
          0.0600,  0.0000,  0.0000,  0.0000,  0.2200, ...
         -0.2200,  0.5600 ]';
  y0 = [  0.0000, -0.0184,  0.7534,  0.7534,  1.1398, ...
          0.0000, -0.1378, -0.1378, -0.6050, -0.6050, ...
         -0.6050,  0.1000, -0.1000,  0.3500,  0.0000, ...
          0.0000, -0.4000 ]';
  z0 = [  0.0000, -0.0185,  0.0000,  0.0000, -0.1957, ...
         -0.7620, -0.1905, -0.1905,  0.3810,  0.3810, ...
          0.3810,  0.3810,  0.1270,  0.3810,  0.3810, ...
          0.3810,  0.3810 ]';
  a1 = [  0.7233,  0.7008,  0.1270,  0.1270,  0.1270, ...
          0.4575,  0.0635,  0.0635,  0.0460,  0.0230, ...
          0.0230,  0.0460,  0.2100,  0.1100,  0.1600, ...
          0.1600,  0.0300 ]';
  a2 = [  0.9644,  0.9246,  0.1270,  0.1270,  0.3400, ...
          0.6099,  0.3175,  0.3175,  0.0230,  0.0230, ...
          0.0460,  0.0460,  0.2581,  0.2500,  0.3100, ...
          0.4100,  0.2000 ]';
  a3 = [  1.2700,  1.2241,  0.1270,  0.1270,  0.1700, ...
          0.5080,  0.3175,  0.3175,  0.0230,  0.0460, ...
          0.0230,  0.0460,  0.2581,  0.2300,  0.2540, ...
          0.3810,  0.2000 ]';
  g =  [  2.0000, -0.9800, -1.0000, -1.0000,  1.5000, ...
         -1.0000,  1.0000,  1.0000,  0.0100,  0.0100, ...
          0.0100,  0.0100,  0.0100,  0.0100, -0.0200, ...
         -0.0200,  0.0300 ]';
  v11 = [ 1.0000,  1.0000,  1.0000,  1.0000,  1.0000, ...
          1.0000,  0.9903, -0.9903,  1.0000,  1.0000, ...
          1.0000,  1.0000,  1.0000,  1.0000,  0.9511, ...
         -0.9511,  0.9192 ]';
  v12 = [ 0.0000,  0.0000,  0.0000,  0.0000,  0.0000, ...
          0.0000, -0.1085, -0.1085,  0.0000,  0.0000, ...
          0.0000,  0.0000,  0.0000,  0.0000, -0.3090, ...
         -0.3090, -0.3381 ]';
  v13 = [ 0.0000,  0.0000,  0.0000,  0.0000,  0.0000, ...
          0.0000, -0.0865, -0.0865,  0.0000,  0.0000, ...
          0.0000,  0.0000,  0.0000,  0.0000,  0.0000, ...
          0.0000,  0.2020 ]';
  v21 = [ 0.0000,  0.0000,  0.0000,  0.0000,  0.0000, ...
          0.0000,  0.1089, -0.1089,  0.0000,  0.0000, ...
          0.0000,  0.0000,  0.0000,  0.0000,  0.3090, ...
         -0.3090,  0.3452 ]';
  v22 = [ 1.0000,  1.0000,  1.0000,  1.0000,  0.5446, ...
          1.0000,  0.9941,  0.9941,  1.0000,  1.0000, ...
          1.0000,  1.0000,  1.0000,  1.0000,  0.9511, ...
          0.9511,  0.9385 ]';
  v23 = [ 0.0000,  0.0000,  0.0000,  0.0000, -0.8387, ...
          0.0000,  0.0000,  0.0000,  0.0000,  0.0000, ...
          0.0000,  0.0000,  0.0000,  0.0000,  0.0000, ...
          0.0000,  0.0000 ]';
  v31 = [ 0.0000,  0.0000,  0.0000,  0.0000,  0.0000, ...
          0.0000,  0.0860, -0.0860,  0.0000,  0.0000, ...
          0.0000,  0.0000,  0.0000,  0.0000,  0.0000, ...
          0.0000,  0.1896 ]';
  v32 = [ 0.0000,  0.0000,  0.0000,  0.0000,  0.8387, ...
          0.0000, -0.0094, -0.0094,  0.0000,  0.0000, ...
          0.0000,  0.0000,  0.0000,  0.0000,  0.0000, ...
          0.0000, -0.0697 ]';
  v33 = [ 1.0000,  1.0000,  1.0000,  1.0000,  0.5446, ...
          1.0000,  0.9963,  0.9963,  1.0000,  1.0000, ...
          1.0000,  1.0000,  1.0000,  1.0000,  1.0000, ...
          1.0000, -0.9794 ]';

  f = zeros ( size ( x ) );

  for e = 1 : 17

    c = ( ( ( x - x0(e) ) * v11(e) ...
          + ( y - y0(e) ) * v12(e) ...
          + ( z - z0(e) ) * v13(e) ) / a1(e) ) .^2 ...
      + ( ( ( x - x0(e) ) * v21(e) ...
          + ( y - y0(e) ) * v22(e) ...
          + ( z - z0(e) ) * v23(e) ) / a2(e) ) .^2 ...
      + ( ( ( x - x0(e) ) * v31(e) ...
          + ( y - y0(e) ) * v32(e) ...
          + ( z - z0(e) ) * v33(e) ) / a3(e) ) .^2;

    i = find ( c <= 1.0 );

    f(i) = f(i) + g(e);

  end

  return
end

