"Distance" from repeatedly changing origin in matrix
Asked by Christian jensen on 10 May 2012
Latest activity Commented on by Christian jensen on 13 May 2012
Hi guys. I have a matrix of densities in x and y (a map) from which I want to generate a gravity anomaly map. To do this I use a function which takes in the densities and x,y coordinates from the density matrix and spits out the gravity anomaly at the desired point. For each point of gravity anomaly I need to run through the whole density matrix with that point at the origin, measuring "distance" (x and y) from the current origin. To make the whole gravity anomaly matrix then, I need to change the point of origin for every gravity anomaly measurement.
My big question is; how do I measure this distance? I find it extremely difficult to even explain this problem proporly.
Said in another way: I want to loop through a matrix and measure the x-x_i and y-y_i distance relative to a chosen origin (x,y).
I'm thinking elaborate loops, in loops, in loops, but I hope to avoid this. If someone has an idea, or can direct me to a thread I've overlooked, I would be glad to hear it. Thanks!
Products
No products are associated with this question.
2 Answers
Answer by Christian jensen on 10 May 2012
This is what I've managed to do so far. The result should look something like this, but this method is far too slow since I need to calculate 160x160 blocks. Not just 20x20 as in this example.
test = ones(20,20); % Density matrix test(3,5) = 2000; % One block with 2000 kg/m3 density test(19,12) = 1500; % Another block with 1500 kg/m3 density [l,w] = size(test); test = test(:); z = 10000; % Constant depth dist_x = zeros(w*l,1); % Preassigned matrix dist_y = zeros(w*l,1); % Preassigned matrix delta_g = zeros(w,l); % Preassigned matrix G = 6.67428e-11; % Gravitational constant V = a*a*4*a; % Cube/grid volume. Every cube is 1x1x4 km
for k = 1:l % x-axis point of origin
for h = 1:w % y-axis point of origin
x = h;
y = k;
count = 1;
% Pre-calculating distance vectors
for i = 1:l
for j = 1:w
dist_x(count) = abs(x-i)*1000;
dist_y(count) = abs(y-j)*1000;
count = count + 1;
end
end
temp = 0;
% Calculating gravity anomaly in every point with respect to
% current origin and summing in delta_g(h,k)
for i = 1:l*w
M = test(i)*V; % Mass
r = sqrt(dist_x(i)^2 + dist_y(i)^2 + z^2); % Distance from source
temp = temp + ((G*M*z/r^3)*(1 - (35/18)*(a/r)^4))*10^5;
end
delta_g(h,k) = temp;
end
end
figure
imagesc(delta_g)
1 Comment
Sean de Wolski on 10 May 2012
don't have time to look too much into this right now. But consider these few points:
M = test(i)*V. V doesn't change and you do this to each element in test. Thus, outside of the for-loop:
M = test*V;
Then just reference M(i) when necessary, you won't need test any more.
Next, the whole for i=1:lw, could be easily cut down. How many of those variables never change, yet you recalculate each time. For example, z^2, or 35/17, etc. None of these values ever change so there's no reason to recalculate them (l*w)^2 times.
finally, temp could be redone as couple of sums to get rid of this loop all together
r = sqrt(distx.^2+disty.^2+z^2)
delta_g = reshape(sum(G*M*the_rest_of_the_stuff,2),l,w); %i.e a matrix multiplication of M and the rest of stuff reshaped to your choice shape.
Answer by Christian jensen on 13 May 2012
Thank you very much for your inputs. I've taken every one of them into consideration. I found a way to completely get rid of loops:
function delta_g = grav_anom(a,rho,x,y,z) % input a = 1/2 of side cube side length G = 6.67428e-11; % Gravitaional constant V = (2*a)^3; % Cube is actually 1x1x4 km M = rho*V; % Mass grid = 160; [X,Y] = meshgrid(0:1:grid-1,0:1:grid-1); X = X(:); Y = Y(:); r = sqrt(((X-x)*1000).^2 + ((Y-y)*1000).^2 + (z*1000)^2); % Distance from source delta_g = ((G*M*z*1000./r.^3).*(1 - (35/18).*(a^4./r.^4))).*10^5; delta_g = reshape(delta_g,grid,grid);
Once again your help was much appriciated. Cheers!
2 Comments
1 Comment
Direct link to this comment:
http://www.mathworks.de/matlabcentral/answers/37968#comment_78617
Please provide us with small data example sets with what you expect (can use ones or magic squares - but something we can work with).
This looks like it will be a good bsxfun() problem.