function xp=radialpart129to37(x)
% xp=RADIALPART129TO37(x)
%
% INPUT:
%
% x         The three-dimensional array of size MxNx129
%
% OUTPUT:
%
% xp        The three-dimensional array of size MxNx37 which also
%           includes the Jacobian already
%
% Assuming 129 equally spaced radii (edges), numerically integrate to 37
% layers, according to hardwired layers provided in here.

% Written by Ignace Loris (igloris@vub.ac.be) on 19/06/2009
% Last modified by fjsimons-at-alum.mit.edu, 07/14/2009

% Error check on hardwired dimensions
nx3=size(x,3);
if nx3~=129
  error('Input should have size 129 in third dimension!')
end

% Prepare array
xp=zeros(size(x,1),size(x,2),37);

% Hardwire the radii, in km, and produce the jacobian, in km^2
% The Earth's total radius
R_E=6371;
% The Earth's core radius (the core-mantle-boundary)
R_CMB=3481.4;
% The thickness of the layers
dr=(R_E-R_CMB)/(nx3-1);
% The location of the interfaces themselves
r=linspace(R_CMB,R_E,nx3);
% Produce the Jacobian and include the sampling interval
Jdr=r.^2*dr;
% Multiply this into the input matrix already
prefactor=repmat(Jdr,[size(x,2) 1 size(x,1)]);
% Special provision if x is effecively unidimensionsl
if size(x,1)~=1 && size(x,2)~=1
  prefactor=shiftdim(prefactor,2);
else
  prefactor=reshape(prefactor,size(x));
end
x=x.*prefactor;

% Perform the numerical integration:
% Use trapezoid rule for 22.575 km layers:
xp(:,:,37)=(x(:,:,129)+x(:,:,128))/2;
xp(:,:,36)=(x(:,:,128)+x(:,:,127))/2;

% Use Simpson's rule for 45.15 km layers:
xp(:,:,35)=(x(:,:,125)+4*x(:,:,126)+x(:,:,127))/3;
xp(:,:,31)=(x(:,:,111)+4*x(:,:,112)+x(:,:,113))/3;
xp(:,:,30)=(x(:,:,109)+4*x(:,:,110)+x(:,:,111))/3;
xp(:,:,27)=(x(:,:, 99)+4*x(:,:,100)+x(:,:,101))/3;
xp(:,:,26)=(x(:,:, 97)+4*x(:,:, 98)+x(:,:, 99))/3;
xp(:,:, 2)=(x(:,:,  3)+4*x(:,:,  4)+x(:,:,  5))/3;
xp(:,:, 1)=(x(:,:,  1)+4*x(:,:,  2)+x(:,:,  3))/3;

% These layers are the thickest: 90.3 km
% Use extended Simpson's rule, with weights [1 4 2 4 1]/3;
for i=2:24
  xp(:,:,i+1)=(x(:,:,4*i-3)+4*x(:,:,4*i-2)+2*x(:,:,4*i-1)+4*x(:,:,4*i)+x(:,:,4*i+1))/3;
end
for i=26:27
  xp(:,:,i+2)=(x(:,:,4*i-3)+4*x(:,:,4*i-2)+2*x(:,:,4*i-1)+4*x(:,:,4*i)+x(:,:,4*i+1))/3;
end
for i=29:31
  xp(:,:,i+3)=(x(:,:,4*i-3)+4*x(:,:,4*i-2)+2*x(:,:,4*i-1)+4*x(:,:,4*i)+x(:,:,4*i+1))/3;
end