MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

# "Distance" from repeatedly changing origin in matrix

Asked by Christian jensen on 10 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!

## 1 Comment

Sean de Wolski on 10 May 2012

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.

## Products

No products are associated with this question.

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!