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!
No products are associated with this question.
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)
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!