Implement this in your Matlab
clear
clc
% Loading the surface
getd = @(p)path(p,path);
getd('toolbox_signal/');
getd('toolbox_general/');
getd('toolbox_graph/');
getd('toolbox_wavelet_meshes/');
% getting the surface
% [vertex,faces] = read_mesh('bunny');
% [vertex,faces] = read_mesh('elephant-50kv');
[vertex,faces] = read_mesh('cow');
N = size(vertex,2);
% underlying geography
tau = ones(N,1);
A_bar = ones(N,1);
U_bar = ones(N,1);
% parameters
sigma = 9;
alpha = 0;
beta = 0;
% Calculating bilateral trade costs
% Now we calculate the distances
distance = NaN(N,N); % this is the matrix of all the bilateral distances
h = waitbar(0,'Calculating bilateral trade costs');
for n=1:N
waitbar(n/N,h)
pstarts = n;
options = struct('W',tau,'start_points',pstarts);
[D,S,Q] = perform_fast_marching_mesh(vertex, faces, pstarts, options);
distance(n,:) = D';
end
close(h)
% trade costs are the exponent
distance = log(1.5) * (distance ./ mean(mean(distance))); % making the average trade cost about 1.5
T = exp(distance);
% Determining equilibrium distribution of population
% creating the amalgamation variables
gamma_1 = ((sigma-1)/(2*sigma-1))*(1-(sigma-1)*(alpha-beta))-beta*(sigma-1);
gamma_2 = (sigma-1)-((sigma-1)^2)/(2*sigma-1);
gamma_3 = ((sigma-1)^2)/(2*sigma-1);
gamma_4 = ((sigma-1)/(2*sigma-1))*(1-(sigma-1)*(alpha-beta))+alpha*(sigma-1);
% iterating to find an equilibrium
diff = 1;
tol = 0.0001;
L = ones(N,1);
h = waitbar(0,'Determining equilibrium distribution of economic activity');
while diff>tol
waitbar(1-diff,h)
L_new = ((U_bar.^gamma_2).*(A_bar.^gamma_3).*...
mean((T.^(1-sigma)).*(ones(size(A_bar))*((A_bar.^gamma_2).*(U_bar.^gamma_3).*(L.^gamma_4))'),2)).^(1/gamma_1);
L_new = L_new ./ mean(L_new(L_new>0)); % normalizing to have mean of one
diff = norm(L - L_new);
L = L_new;
end
close(h)
w = ((L.^(1-(sigma-1)*(alpha-beta))).*(A_bar.^(1-sigma)).*(U_bar.^(sigma-1))).^(1/(1-2*sigma));
w = w ./ mean(w(L_new>0)); % normalizing to have mean of one
w(isnan(w)) = 1; % normalizing wages to 1 where no one lives
% Making the figure
clear options
options = struct('face_vertex_color',L);
figure(1)
clf
h = plot_mesh(vertex,faces, options);
hold on
rotate(h,[0 1 1],50)
colormap(jet)
shading flat
lighting none
material dull
colorbar
view([25, 90])
hold off
print -dpng cow_popdistn_051713
