Your d, g are arrays of same size as t....You can use mesh grid on d/g... Are you looking for some thing like below?
clc;
S = 21;
t = linspace(-pi,pi,S);
knownPoints = [cos(t);sin(t)]; %(x,y)
f = exp(1i*t);
d = zeros (1,S);
calculationPoint = [10;9];
for idx=1:S
tmp1 = calculationPoint(1,1)-knownPoints(1,idx);
tmp2 = calculationPoint(2,1)-knownPoints(2,idx);
d(idx) = sqrt(tmp1^2+tmp2^2);
end
g = d.^2 + 1./d.^3;
[T,G] = meshgrid(t,g) ;
z = 1/S*sum(g.*f); % z function is of parameter calculation point