function f = partition(z);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function f = partition(z)
% This function evaluates the partition function for the
% harmonic oscillator potential at the points specified
% in the rows of z.  The number of columns of z determines
% myd, the number of particles. The partition function is
% evaluated relative to a particle at the position specified
% in the global variable mya.  The variables myk and
% mydelta are parameters in the function evaluation.
% partition.m Dianne P. O'Leary   09/2004
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global myd mya myk mydelta

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Evaluate the harmonic oscillator potential.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

a = mya;
V =  .5*myk*(z.*z);
Va = .5*myk*a^2;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute first part of exponent, determined by the
% distance between particles.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

dz = sum((z(:,2:myd)-z(:,1:myd-1)).^2,2) + ...
         (a-z(:,1)).^2 + (z(:,myd)-a).^2;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute second part of exponent, determined by the
% harmonic oscillator potential.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

vpart = 2*sum(V,2) + 2*Va;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute the argument for the exponentials.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

eg = - .5/mydelta * dz - .5*mydelta * vpart;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute the function.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

f = (1/sqrt(2*pi*mydelta))^(myd+1) *exp(eg);