function result = pyramid_unit_o48_3d ( func ) %*****************************************************************************80 % %% PYRAMID_UNIT_O48_3D approximates an integral inside the unit pyramid in 3D. % % Discussion: % % An 48 point degree 7 formula, Stroud CN:C2:7-1, is used. % % The (X,Y,Z) integration region can be represented as: % % - ( 1 - Z ) <= X <= 1 - Z % - ( 1 - Z ) <= Y <= 1 - Z % 0 <= Z <= 1. % % When Z is zero, the integration region is a square lying in the (X,Y) % plane, centered at (0,0,0) with "radius" 1. As Z increases to 1, the % radius of the square diminishes, and when Z reaches 1, the square has % contracted to the single point (0,0,1). % % Licensing: % % This code is distributed under the GNU LGPL license. % % Modified: % % 01 April 2008 % % Author: % % John Burkardt % % Reference: % % Arthur Stroud, % Approximate Calculation of Multiple Integrals, % Prentice Hall, 1971, % ISBN: 0130438936, % LC: QA311.S85. % % Parameters: % % Input, external FUNC, the name of the user supplied function which % evaluates F(X,Y,Z), of the form % function value = func ( x, y, z ) % % Output, real RESULT, the approximate integral of the function. % order = 48; w = [ ... 2.01241939442682455E-002, ... 2.01241939442682455E-002, ... 2.01241939442682455E-002, ... 2.01241939442682455E-002, ... 2.60351137043010779E-002, ... 2.60351137043010779E-002, ... 2.60351137043010779E-002, ... 2.60351137043010779E-002, ... 1.24557795239745531E-002, ... 1.24557795239745531E-002, ... 1.24557795239745531E-002, ... 1.24557795239745531E-002, ... 1.87873998794808156E-003, ... 1.87873998794808156E-003, ... 1.87873998794808156E-003, ... 1.87873998794808156E-003, ... 4.32957927807745280E-002, ... 4.32957927807745280E-002, ... 4.32957927807745280E-002, ... 4.32957927807745280E-002, ... 1.97463249834127288E-002, ... 1.97463249834127288E-002, ... 1.97463249834127288E-002, ... 1.97463249834127288E-002, ... 5.60127223523590526E-002, ... 5.60127223523590526E-002, ... 5.60127223523590526E-002, ... 5.60127223523590526E-002, ... 2.55462562927473852E-002, ... 2.55462562927473852E-002, ... 2.55462562927473852E-002, ... 2.55462562927473852E-002, ... 2.67977366291788643E-002, ... 2.67977366291788643E-002, ... 2.67977366291788643E-002, ... 2.67977366291788643E-002, ... 1.22218992265373354E-002, ... 1.22218992265373354E-002, ... 1.22218992265373354E-002, ... 1.22218992265373354E-002, ... 4.04197740453215038E-003, ... 4.04197740453215038E-003, ... 4.04197740453215038E-003, ... 4.04197740453215038E-003, ... 1.84346316995826843E-003, ... 1.84346316995826843E-003, ... 1.84346316995826843E-003, ... 1.84346316995826843E-003 ]'; x = [ ... 0.88091731624450909, ... -0.88091731624450909, ... 0.0000000000000000, ... 0.0000000000000000, ... 0.70491874112648223, ... -0.70491874112648223, ... 0.0000000000000000, ... 0.0000000000000000, ... 0.44712732143189760, ... -0.44712732143189760, ... 0.0000000000000000, ... 0.0000000000000000, ... 0.18900486065123448, ... -0.18900486065123448, ... 0.0000000000000000, ... 0.0000000000000000, ... 0.36209733410322176, ... -0.36209733410322176, ... -0.36209733410322176, ... 0.36209733410322176, ... 0.76688932060387538, ... -0.76688932060387538, ... -0.76688932060387538, ... 0.76688932060387538, ... 0.28975386476618070, ... -0.28975386476618070, ... -0.28975386476618070, ... 0.28975386476618070, ... 0.61367241226233160, ... -0.61367241226233160, ... -0.61367241226233160, ... 0.61367241226233160, ... 0.18378979287798017, ... -0.18378979287798017, ... -0.18378979287798017, ... 0.18378979287798017, ... 0.38925011625173161, ... -0.38925011625173161, ... -0.38925011625173161, ... 0.38925011625173161, ... 7.76896479525748113E-02, ... -7.76896479525748113E-02, ... -7.76896479525748113E-02, ... 7.76896479525748113E-02, ... 0.16453962988669860, ... -0.16453962988669860, ... -0.16453962988669860, ... 0.16453962988669860 ]'; y = [ ... 0.0000000000000000, ... 0.0000000000000000, ... 0.88091731624450909, ... -0.88091731624450909, ... 0.0000000000000000, ... 0.0000000000000000, ... 0.70491874112648223, ... -0.70491874112648223, ... 0.0000000000000000, ... 0.0000000000000000, ... 0.44712732143189760, ... -0.44712732143189760, ... 0.0000000000000000, ... 0.0000000000000000, ... 0.18900486065123448, ... -0.18900486065123448, ... 0.36209733410322176, ... 0.36209733410322176, ... -0.36209733410322176, ... -0.36209733410322176, ... 0.76688932060387538, ... 0.76688932060387538, ... -0.76688932060387538, ... -0.76688932060387538, ... 0.28975386476618070, ... 0.28975386476618070, ... -0.28975386476618070, ... -0.28975386476618070, ... 0.61367241226233160, ... 0.61367241226233160, ... -0.61367241226233160, ... -0.61367241226233160, ... 0.18378979287798017, ... 0.18378979287798017, ... -0.18378979287798017, ... -0.18378979287798017, ... 0.38925011625173161, ... 0.38925011625173161, ... -0.38925011625173161, ... -0.38925011625173161, ... 7.76896479525748113E-02, ... 7.76896479525748113E-02, ... -7.76896479525748113E-02, ... -7.76896479525748113E-02, ... 0.16453962988669860, ... 0.16453962988669860, ... -0.16453962988669860, ... -0.16453962988669860 ]'; z = [ ... 4.85005494469969989E-02, ... 4.85005494469969989E-02, ... 4.85005494469969989E-02, ... 4.85005494469969989E-02, ... 0.23860073755186201, ... 0.23860073755186201, ... 0.23860073755186201, ... 0.23860073755186201, ... 0.51704729510436798, ... 0.51704729510436798, ... 0.51704729510436798, ... 0.51704729510436798, ... 0.79585141789677305, ... 0.79585141789677305, ... 0.79585141789677305, ... 0.79585141789677305, ... 4.85005494469969989E-02, ... 4.85005494469969989E-02, ... 4.85005494469969989E-02, ... 4.85005494469969989E-02, ... 4.85005494469969989E-02, ... 4.85005494469969989E-02, ... 4.85005494469969989E-02, ... 4.85005494469969989E-02, ... 0.23860073755186201, ... 0.23860073755186201, ... 0.23860073755186201, ... 0.23860073755186201, ... 0.23860073755186201, ... 0.23860073755186201, ... 0.23860073755186201, ... 0.23860073755186201, ... 0.51704729510436798, ... 0.51704729510436798, ... 0.51704729510436798, ... 0.51704729510436798, ... 0.51704729510436798, ... 0.51704729510436798, ... 0.51704729510436798, ... 0.51704729510436798, ... 0.79585141789677305, ... 0.79585141789677305, ... 0.79585141789677305, ... 0.79585141789677305, ... 0.79585141789677305, ... 0.79585141789677305, ... 0.79585141789677305, ... 0.79585141789677305 ]'; % % Quadrature. % quad = 0.0; for i = 1 : order quad = quad + w(i) * feval ( func, x(i), y(i), z(i) ); end % % Volume. % volume = pyramid_unit_volume_3d ( ); % % Result. % result = quad * volume; return end