par.a = 75e-3/2;    % radius [m]
par.h = 1.004e-3;   % thickness [m]
par.E = 73.2e9;     % Young's modulus [Pa]
par.nu = 0.155;     % Poisson's ratio
par.rho = 2202;     % density [kg/m^3]

%Calculate fundamental modes of the disk
[freqs, modes, shapes, x, y] = disk_frequencies(par, 10000, 1, 'shapes', 0.5e-3);

%Now we extract the force profile from the COMSOL model
model = run('faster.m');

%Form coordinate grid to match the one used for the mode shapes
x0 = -37.5 : .5 : 37.5;
y0 = -37.5 : .5 : 37.5;
z0 = 1 : 0.1 : 2;
[xg, yg, zg] = meshgrid(x0, y0, z0);
xx = [xg(:),yg(:),zg(:)]';

%Extract force profile in the upward direction F=alpha*grad(E^2)
fy = mphinterp(model,'d(es.Ex^2 + es.Ey^2 + es.Ez^2, z)', 'coord', xx);

%Grid the force data to be in the same form as the shapes
force = griddata(xx(1, :), xx(2, :), xx(3, :), fy, xg, yg, zg);

%{
figure()
subplot(121)
pcolor(x,y, force)
shading flat
axis equal
axis tight
subplot(122)
pcolor(x,y, shapes{3})
shading flat
axis equal
axis tight
%}


dx = xg(1, 2, 1) - xg(1, 1, 1);
dy = yg(2, 1, 1) - yg(1, 1, 1);
dz = zg(1, 1, 2) - zg(1, 1, 1);

%Numerically integrate the profile for each mode.
a = size(modes);
b = a(1);
product = zeros(1, b);
for n = 1 : a(1)
    product(1, n) = sum(shapes{n}(~isnan(shapes{n})) .* force(~isnan(shapes{n}))) * dx * dy * dz;
end