You are now following this question
- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.
simple nested loops error
3 views (last 30 days)
Show older comments
Hi,
I want to do modfication(nested loop) in the current code limts[r,c].
I have twelve files avg_mat(i = 12) to Process . But each file needs to use diffrent r and c values in the code. I need to mutiply r and c with cosine theta for every file.
where theta increases from 0deg to max deg and decreases to 0deg at the end with the increment (max/number of files). Please help.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
results = zeros(numl(avg_mat),1);
for i=1:numel(avg_mat)
writematrix(avg_mat{i}, ['avg_val_' num2str(i)]);
R = avg_mat{i}(:,1);
C = avg_mat{i}(:,2);
F = avg_mat{i}(:,5);
maxi = 90;
for t = [0:(maxi/numl(avg_mat)):maxi:-(maxi/numl(avg_mat)):0]
r(t) = linspace(min(R), max(R), 1000)*cos(t);%need modification based the file
c(t) = linspace(min(C), max(C), 1000)*cos(t);%need modification based on the file
[Rg, Cg] = meshgrid(r(t), c(t));
Fg(i,t) = griddata(R, C, F, Rg, Cg);
result(i,t) = trapz(c(t), trapz(r(t), Fg, 2));
end
end
Accepted Answer
Walter Roberson
on 27 Apr 2020
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
navg = numel(avg_mat);
r = cell(navg,1);
t = cell(navg,1);
maxi = 90;
t1 = linspace(0, maxi, navg);
t2 = fliplr(t1(1:end-1));
tvals = [t1, t2];
nt = length(tvals);
results = zeros(navg,nt);
Fg = cell(navg,nt);
for i = 1:navg
writematrix(avg_mat{i}, ['avg_val_' num2str(i)]);
R = avg_mat{i}(:,1);
C = avg_mat{i}(:,2);
F = avg_mat{i}(:,5);
for tidx = 1:nt
t = tvals(tidx);
r{t} = linspace(min(R), max(R), 1000)*cosd(t);
c{t} = linspace(min(C), max(C), 1000)*cosd(t);
[Rg, Cg] = meshgrid(r, c);
Fg{i,tidx} = griddata(R, C, F, Rg, Cg);
result(i,t) = trapz(c{t}, trapz(r{t}, Fg{i,tidx}, 2));
end
end
I suspect that you do not really need to record all those r and c and Fg values, but recording them is at least useful during the debugging phase.
8 Comments
MS
on 27 Apr 2020
Edited: MS
on 27 Apr 2020
Thanks.
I got this error. please help me to correct it
Array indices must be positive integers or logical values.
Error in average_folder_new (line 44)
r{t} = linspace(min(R), max(R), 1000)*cosd(t);
Walter Roberson
on 28 Apr 2020
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
navg = numel(avg_mat);
r = cell(navg,1);
t = cell(navg,1);
maxi = 90;
t1 = linspace(0, maxi, navg);
t2 = fliplr(t1(1:end-1));
tvals = [t1, t2];
nt = length(tvals);
results = zeros(navg,nt);
Fg = cell(navg,nt);
for i = 1:navg
writematrix(avg_mat{i}, ['avg_val_' num2str(i)]);
R = avg_mat{i}(:,1);
C = avg_mat{i}(:,2);
F = avg_mat{i}(:,5);
for tidx = 1:nt
t = tvals(tidx);
r{tidx} = linspace(min(R), max(R), 1000)*cosd(t);
c{tidx} = linspace(min(C), max(C), 1000)*cosd(t);
[Rg, Cg] = meshgrid(r, c);
Fg{i,tidx} = griddata(R, C, F, Rg, Cg);
result(i,tidx) = trapz(c{tidx}, trapz(r{tidx}, Fg{i,tidx}, 2));
end
end
MS
on 28 Apr 2020
thanks for correction. please look into it.
i got this error.
Error using griddata (line 64)
GRIDDATA no longer requires Qhull-specific options.
Please remove these options when calling GRIDDATA.
Error in average_folder_new (line 49)
Fg{i,tidx} = griddata(R, C, F, Rg, Cg);
Walter Roberson
on 28 Apr 2020
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
navg = numel(avg_mat);
r = cell(navg,1);
t = cell(navg,1);
maxi = 90;
t1 = linspace(0, maxi, navg);
t2 = fliplr(t1(1:end-1));
tvals = [t1, t2];
nt = length(tvals);
results = zeros(navg,nt);
Fg = cell(navg,nt);
for i = 1:navg
writematrix(avg_mat{i}, ['avg_val_' num2str(i)]);
R = avg_mat{i}(:,1);
C = avg_mat{i}(:,2);
F = avg_mat{i}(:,5);
rvec = linspace(min(R), max(R), 1000);
cvec = linspace(min(C), max(C), 1000);
for tidx = 1:nt
t = tvals(tidx);
r{tidx} = rvec*cosd(t);
c{tidx} = cvec*cosd(t);
[Rg, Cg] = meshgrid(r{tidx}, c{tidx});
Fg{i,tidx} = griddata(R, C, F, Rg, Cg);
result(i,tidx) = trapz(c{tidx}, trapz(r{tidx}, Fg{i,tidx}, 2));
end
end
MS
on 28 Apr 2020
Edited: MS
on 28 Apr 2020
Thanks it worked without any error. But, the output resut gives various column and various rows. Please have a look at it.
for eg., the below code gives one column result which is right
results = zeros(numl(avg_mat),1);
for i=1:numel(avg_mat)
writematrix(avg_mat{i}, ['avg_val_' num2str(i)]);
R = avg_mat{i}(:,1);
C = avg_mat{i}(:,2);
F = avg_mat{i}(:,5);
r = linspace(min(R), max(R), 1000);
c = linspace(min(C), max(C), 1000);
[Rg, Cg] = meshgrid(r, c);
Fg = griddata(R, C, F, Rg, Cg);
result(i) = trapz(c, trapz(r, Fg, 2));
end
writematrix(result(:), 'integral_val.txt')
Walter Roberson
on 28 Apr 2020
What is numl and how does it differ from numel ?
MS
on 28 Apr 2020
Edited: MS
on 28 Apr 2020
sorry, it was a typo. it should be numel. can you please Check the error in the results output. It has 23 columns. The first and last column are the results with the same values, remaining columns are NAN values.
Walter Roberson
on 28 Apr 2020
Your original posted code had
result(i,t) = trapz(c(t), trapz(r(t), Fg, 2));
That was apparently intended to produce one result for every element of avg_mat and every t value. So that is what my code does.
More Answers (1)
MS
on 28 Apr 2020
Edited: Walter Roberson
on 28 Apr 2020
I am attaching you the 12 input files for your kind refernce. Every file should give one result(1*1) result. I should get only result as 12 *1 for the 12 files. did i do something wrong? can you please check.
earlier code was
result(i) = trapz(c, trapz(r, Fg, 2));
it gave me 1*1 result for every file.
we are mutiplying the limts[Rg, Cg] by cos (t). it should give the results 1*1 for every file. I might have confused you.
now, it converted to
result(i,t) = trapz(c(t), trapz(r(t), Fg, 2));
can you please check.
32 Comments
Walter Roberson
on 28 Apr 2020
we are mutiplying the limts[Rg, Cg] by cos (t)
Yes, but cos(t) has a different value for each different t value, so you need to calculate for each different t, and that set of calculations has to be done for each different files. For 12 different files, and the t values increasing to a maximum and decreasing again to a minimum according to the number of files (not clear why the step size depends on the number of files), that would be like an 12 x (12*2-1) output.
Note: when you decrease t again from the maximum, you are going to have t values that you have already processed. It seems a waste to re-process them. Why are you doing that? Is it possible that instead you want to process negative values through 0 and then up to positive ??
MS
on 28 Apr 2020
Edited: MS
on 28 Apr 2020
thanks, it is an area integeration(trapz) of tweleve files each file has twelve limts [Rg, Cg] by cos (t)
eg.,
first file trapz(avg_mat(1) has the integeration limits [Rg, Cg]*cos(0)
second file trapz(avg_mat(2))has the integeration limits [Rg, Cg]*cos(90/6)
...........
sixth file trapz(avg_mat(6)) has the integeration limits [Rg, Cg]*cos(90)
...............
last file trapz(avg_mat(12)) has the integeration limts [Rg, Cg]*cos(0)
so every file should give one unique result(1*1). I should get (12*1) results combined for 12 files. But the code is giving (12*23) (navg*tidx) result values fro some reason. Please let me know if you need further clarifications.
Walter Roberson
on 28 Apr 2020
I do not see the rule for the pattern 0, 90/6, ?, ?, ?, 90/6, ?, ?, ?, ?, ?, 0 .
It would make more sense to me if the pattern were like 180/6 * (file number minus 1)
MS
on 28 Apr 2020
Edited: MS
on 28 Apr 2020
thanks, becasue the integeration area is rotating only upto to the maximum 90 deg and comes back to the orginal position(0). We are setting the integeration area which is rotating for every file by cos(t). in a simple terms, each file should be integerated by diffrent limts
with the condition such as
first file , integrated by the first limit
second file, second limit
last file, last limit.
eg.,
first file trapz(avg_mat(1) has the integeration limits [Rg, Cg]*cos(0) , gives one out put(1*1)
second file trapz(avg_mat(2))has the integeration limits [Rg, Cg]*cos(90/6) , gives one out put(1*1)
...........
sixth file trapz(avg_mat(6)) has the integeration limits [Rg, Cg]*cos(90), gives one out put(1*1)
...............
last file trapz(avg_mat(12)) has the integeration limts [Rg, Cg]*cos(0), gives one out put(1*1)
finally, i get results as (12*1). i think, your code is doing 12files*12 limts, instead of first file should be integrated by first limit as well as last file should be integerated by the last limit. Please let me know if you need any calrification to modify the code.
Walter Roberson
on 28 Apr 2020
So the first file is 90/6 * (1-1) = 0, second file is 90/6 * (2-1) = 15, third file is 90/6 * (3-1) = 30, fourth file is 90/6 * (4-1) = 45, fifth file is 90/6 * (5-1) = 60, sixth file is 90/6 * (6-1) = 75, seventh file is 90/6 * (7-1) = 90, eigth file is back to 75, nineth file is back to 60, tenth file is back to 45, 11th file is back to 30, 12th file is back to 15, and a missing 13th file not mentioned is back to 0??
If you want the first file to be 0, and the second to be 15, and the sixth to be 90, then you cannot have an even increment: you would need five steps of 15 to get from 15 to 90, but you are only permitting 4 steps.
We need a formula for the rule, one that will work to give second step of 15 and 6th of 90.
Walter Roberson
on 28 Apr 2020
By the way: wouldn't integration limit of [Rg, Cg]*cos(90 degrees) be 0 to 0 ? Is there any point in doing that calculation if the integration limit is 0 to 0 ?
MS
on 28 Apr 2020
Edited: MS
on 28 Apr 2020
yeah, thanks for very good inputs. i am doing the sinusodial rotaion. I need to consider mutiply by sin(t) instead of cos(t)! i will try to do both (cos, sin) as of now.
since, the integeration area is rotating from (90/6) upto to the maximum 90 deg and comes back to the orginal position(90/6). We are setting the integeration area which is rotating for every file by [rmin, cmin]*cos(t) and [rmax, camx] *sin(t).
I request you please modify the code.
MS
on 28 Apr 2020
Edited: Walter Roberson
on 28 Apr 2020
thanks. I checked the code. I got only one issue. The result should be 12by1 matrix. But, i am getting result as 12by12 matrix. please help me to fix the error.
navg = numel(avg_mat);
r = cell(navg,1);
t = cell(navg,1);
maxi = 90;
t1 = linspace(15, maxi, 6);
t2 = fliplr(t1(1:end));
tvals = [t1, t2];
nt = length(tvals);
results = zeros(navg,nt);
Fg = cell(navg,nt);
%results = zeros(numel(avg_mat),1);
for i = 1:navg
writematrix(avg_mat{i}, ['avg_val_' num2str(i)]);
R = avg_mat{i}(:,1);
C = avg_mat{i}(:,2);
F = avg_mat{i}(:,5);
rvec = linspace(min(R), max(R), 1000);
cvec = linspace(min(C), max(C), 1000);
for tidx = 1:nt
t = tvals(tidx);
r{tidx} = rvec*cosd(t);
c{tidx} = cvec*sind(t);
[Rg, Cg] = meshgrid(r{tidx}, c{tidx});
Fg{i,tidx} = griddata(R, C, F, Rg, Cg);
Fg{i,tidx}(isnan(Fg{i,tidx})) = 0;
result(i,tidx) = trapz(c{tidx}, trapz(r{tidx}, Fg{i,tidx}, 2));%issue only here.
end
end
MS
on 28 Apr 2020
Edited: MS
on 28 Apr 2020
Hi walter,
Thank you very much for all the help. I think, i need to choose only diagonal elements[(1,1),(2,2)...(12,12) ] of the results. Please give your suggestions.
final_result = diag(result); is this correct?
Walter Roberson
on 29 Apr 2020
navg = numel(avg_mat);
r = cell(navg,1);
c = cell(navg,1);
maxi = 90;
%this code is designed to handle odd navg as well as even
t1 = linspace(0, maxi, ceil(navg/2)+1);
t2 = linspace(maxi, 0, navg - length(t1)+2);
tvals = [t1(2:end), t2(2:end)];
results = zeros(navg,1);
Fg = cell(navg,1);
for i = 1:navg
writematrix(avg_mat{i}, ['avg_val_' num2str(i)]);
R = avg_mat{i}(:,1);
C = avg_mat{i}(:,2);
F = avg_mat{i}(:,5);
r{i} = linspace(min(R), max(R), 1000) * cosd(tvals(i));
c{i} = linspace(min(C), max(C), 1000) * sind(tvals(i));
[Rg, Cg] = meshgrid(r{i}, c{i});
tg = griddata(R, C, F, Rg, Cg);
tg(isnan(tg)) = 0;
Fg{i} = tg;
result(i) = trapz(c{i}, trapz(r{i}, Fg{i}, 2));
end
MS
on 29 Apr 2020
Edited: MS
on 29 Apr 2020
I need a help to check the limts. I am making a huge logic mistake here. I want to keep the area (width, height) cosntant and rotate the limts by cos(t). When i try to mutiply by cosd(tvals(i)). It is changing the area values insead the roating the area as shown in the figure1. can you help me to keep the same area and rotate the rectangle limits. We need the figure 2 instead of figure1. current code is giving figure1. ![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/288644/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/288644/image.png)
Walter Roberson
on 29 Apr 2020
Edited: Walter Roberson
on 30 Apr 2020
This is to be expected from the way you are processing. You are creating vector r and vector c and doing meshgrid() based on those. vector r and vector c without the cosd() and sind() would obviously be a rectangular grid. vector r and vector c times cosd() and sind() is just doing a linear scaling on the limits and producing a rectangular grid from what results. In order to be a rotation, you need the transformed r position to relate to c also, and the transformed c position to relate to r also. You need a rotation matrix. https://en.wikipedia.org/wiki/Rotation_matrix or equivalent
navg = numel(avg_mat);
r = cell(navg,1);
c = cell(navg,1);
maxi = 90;
%this code is designed to handle odd navg as well as even
t1 = linspace(0, maxi, ceil(navg/2)+1);
t2 = linspace(maxi, 0, navg - length(t1)+2);
tvals = [t1(2:end), t2(2:end)];
results = zeros(navg,1);
Fg = cell(navg,1);
for i = 1:navg
writematrix(avg_mat{i}, ['avg_val_' num2str(i)]);
R = avg_mat{i}(:,1);
C = avg_mat{i}(:,2);
F = avg_mat{i}(:,5);
r0 = linspace(min(R), max(R), 1000);
c0 = linspace(min(C), max(C), 1000);
[Rg0, Cg0] = meshgrid(r0, c0);
cth = cosd(tvals(i));
sth = sind(tvals(i));
Rg = Rg0 .* cth - Cg0 .* sth;
Cg = Rg0 .* sth + Cg0 .* cth;
tg = griddata(R, C, F, Rg, Cg);
tg(isnan(tg)) = 0;
Fg{i} = tg;
result(i) = trapz(c{i}, trapz(r{i}, Fg{i}, 2));
end
MS
on 30 Apr 2020
Edited: MS
on 30 Apr 2020
thanks, i found the way to rotate the points using your link. I am unable to implelmet the rotation matrix in a loop. Please correct the mistake.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
RO{i} = [cosd(tvals(i)) -sind(tvals(i)); sind(tvals(i)) cosd(tvals(i))];%rotation matrix
point1 = [min(R) min(C)]'
point2= [max(R) max(C)]'
rotpoint1 = RO{i}*point1;
rotpoint2 = RO{i}*point2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Walter Roberson
on 30 Apr 2020
the code I posted already has the rotation coded into it.
MS
on 30 Apr 2020
Edited: MS
on 30 Apr 2020
i got this error while running the code. Can you please check it.
Error using trapz (line 62)
X must be a vector.
Error in average_folder_new2 (line 55)
result(i) = trapz(c{i}, trapz(r{i}, Fg{i}, 2));
i think, c{i} and r{i} is a mistake here. I have no idea. how to implement it. Please help me to resolve it.
Walter Roberson
on 30 Apr 2020
I did forget to change the line
result(i) = trapz(c{i}, trapz(r{i}, Fg{i}, 2));
Logically it should probably be more like
result(i) = trapz(Cg, trapz(Rg, Fg{i}, 2));
Unfortunately, you seem to have deleted the zip of data files, so I cannot test it.
I have to think more about whether it is valid to use trapz on a non-rectangular grid
Walter Roberson
on 30 Apr 2020
Question for you: when you rotate the rectangle, where is the center of rotation ? The simple rotation matrix that I coded above was for the case of rotation around (0,0), and with the simple test I did a moment ago, even with a modest angle, the rotation moved all the sample points outside the valid area, giving nan for everything. I retested with data centered around (0,0) and got more useful results.
You might want to have a look at imtransform()
Otherwise:
R = avg_mat{i}(:,1);
C = avg_mat{i}(:,2);
F = avg_mat{i}(:,5);
r0 = linspace(min(R), max(R), 1000);
c0 = linspace(min(C), max(C), 1000);
[Rg0, Cg0] = meshgrid(r0, c0);
%assuming center of rotation is (Xc, Yc) and angle th is degrees
th = tvals(i);
m = makehgtform('translate',-Xc,-Yc,0, 'zrotate',deg2rad(th), 'translate',Xc,Yc,0);
Coords = [Rg0(:), Cg0(:), zeros(numel(Rg0),1), ones(numel(Rg0),1)];
TrCoords = Coords*m.';
Rg = reshape(TrCoords(:,1),size(Rg0));
Cg = reshape(TrCoords(:,2), size(Cg0));
tg = griddata(R, C, F, Rg, Cg);
tg(isnan(tg)) = 0;
Fg{i} = tg;
%result(i) = trapz(Cg, trapz(Rg, Fg{i}, 2)); %WILL NOT WORK
However, now you encounter the problem that trapz() will not accept grids of data, and your coordinates cannot be reduced to marginals.
Now, you use linear spacing for your r0 and c0, and if you translate, rotate, translate back, with no sheer, then the area of each grid element is the same as the area of each other grid element. Theory says that in that case, you can calculate the area using uniform grid, and then multiply the result by the area of one cell. Rotation without sheer preserves area, so you can calculate area from the pre-rotational marginals:
result(i) = trapz(trapz(Fg{i},2)) * (r0(2)-r0(1)) * (c0(2)-c0(1));
%... I think
MS
on 30 Apr 2020
Edited: MS
on 1 May 2020
Thank you very much for the inputs. It worked like a wonder. you are a great advisor to me. your logic tackled this problem easily. I think, you are right about the area integeral and mutiplying with the same elemental area.
I am attaching the folder for you to test the code. Also, I am attaching the correct result document obtained using other manual tool for your refernce. Please assume that the rectangle is rotating with respect to centroid. Kindly let me know the correct approach to get the correct results.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/288938/image.png)
I have few request, is there a way to plot the translated and rotated points as rectangle as shown in figure2. It will help us to verify the rotated area of the griddata is same.
I do not understand the logic completely. I have commented my understaning on the specific line. Please educate me, rotation and translation logic and the diffrence bewteen Coords and TrCoords. I am getting zero when i subtract them. Please have a look into it.
Xc = 0;%center of rotationx
Yc = 0;%center of rotationy
th = tvals(i);
m = makehgtform('translate',-Xc,-Yc,0, 'zrotate',deg2rad(th), 'translate',Xc,Yc,0);%translation,and rotation
Coords = [Rg0(:), Cg0(:), zeros(numel(Rg0),1), ones(numel(Rg0),1)];
TrCoords = Coords*m.';
result(i) = trapz(trapz(Fg{i},2)) * (r0(2)-r0(1)) * (c0(2)-c0(1)) %double integeration(Fg)*constant elemental area
Walter Roberson
on 1 May 2020
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
navg = numel(avg_mat);
r = cell(navg,1);
t = cell(navg,1);
maxi = 90;
t1 = linspace(0, maxi, navg);
t2 = fliplr(t1(1:end-1));
tvals = [t1, t2];
nt = length(tvals);
results = zeros(navg,nt);
Rg = cell(navg,1);
Cg = cell(navg,1);
Fg = cell(navg,nt);
filenames = cell(navg,1);
for i = 1:navg
filenames{i} = sprintf('avg_val_%d.txt', i);
writematrix(avg_mat{i}, filenames{i});
R = avg_mat{i}(:,1);
C = avg_mat{i}(:,2);
F = avg_mat{i}(:,5);
r0 = linspace(min(R), max(R), 1000);
c0 = linspace(min(C), max(C), 1000);
[Rg0, Cg0] = meshgrid(r0, c0);
Xc = mean(R);
Yc = mean(C);
%assuming center of rotation is (Xc, Yc) and angle th is degrees
th = tvals(i);
m = makehgtform('translate',Xc,Yc,0, 'zrotate',deg2rad(th), 'translate',-Xc,-Yc,0);
Coords = [Rg0(:), Cg0(:), zeros(numel(Rg0),1), ones(numel(Rg0),1)];
TrCoords = Coords*m.';
Rg{i} = reshape(TrCoords(:,1),size(Rg0));
Cg{i}= reshape(TrCoords(:,2), size(Cg0));
tg = griddata(R, C, F, Rg{i}, Cg{i});
tg(isnan(tg)) = 0;
Fg{i} = tg;
results(i) = trapz(trapz(Fg{i},2)) * (r0(2)-r0(1)) * (c0(2)-c0(1));
end
for i = 1 : navg
fig = figure();
ax = axes('Parent', fig);
surf(ax, Rg{i}, Cg{i}, Fg{i}, 'edgecolor','none');
hold(ax, 'on');
h(1) = plot(ax, [Rg0(1,[1 end]),Rg0(end,[end 1]),Rg0(1,1)],[Cg0(1,[1 end]),Cg0(end,[end 1]), Cg0(1,1)],'k');
h(2) = plot(ax, [Rg{i}(1,[1 end]),Rg{i}(end,[end 1]),Rg{i}(1,1)],[Cg{i}(1,[1 end]),Cg{i}(end,[end 1]), Cg{i}(1,1)],'r');
hold(ax, 'off');
title(ax, filenames{i}, 'interpreter', 'none');
end
fig = figure();
ax = axes('Parent', fig);
hold(ax, 'on');
for i = 1 : navg
h(i) = plot(ax, [Rg{i}(1,[1 end]),Rg{i}(end,[end 1]),Rg{i}(1,1)],[Cg{i}(1,[1 end]),Cg{i}(end,[end 1]), Cg{i}(1,1)],'r', 'DisplayName', filenames{i});
end
hold(ax, 'off');
L = legend(h, filenames);
L.Interpreter = 'none';
MS
on 1 May 2020
Edited: MS
on 2 May 2020
Thank you very much. This is an excellent code. I have two doubts. Please clear it.
1,
results(i) = trapz(trapz(Fg{i},2)) * (r0(2)-r0(1)) * (c0(2)-c0(1));
In the above line, we are taking double integeration and mutiplying with the elemental area. In that case, the final result unit may be wrong.
For eg.,
if Fg{i} has the units of 1/seconds, then the area integeration should give the unit as m^2/s.
trapz(Fg{i},2) = Q m/s % integerates along the row
trapz(Q) = XX (m^2/s) % integerates along the column gives the single value
results= XX*Area (m^4/s) . sorry, I may be wrong. Please give your suggestions.
2,
The rectangle area should return to the orginal postion as shown in figure . But, our results shows that the rectangle is not returning to the orginal postion. I will attach you the figure for your refernce.
Baciallly, we need to integerate/pick 12 diffrent (F) values at 12 diffrent rectangle area positions (6 area forward + 6 area reverse).
Since, we have the information of F(R,C) from the data and the 12 diffrent area positions from the code.
Can we mutiply the mean data (dF(i)) at a particular area postion by an area dA(i). will it be an issue? Please share your thoughts and the final code.
Result = (F1 + F2 + F3 + F4 + F5+F6+F7+F8+F9+F10+F11+F12).Area(R,C)
or
Result = ∬(dF1.dA1 + dF2.dA2 +dF3.dA3.......+dF12.dA12)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/289472/image.png)
Walter Roberson
on 2 May 2020
trapz(Fg{i},2) = Q m/s % integerates along the row
No, the dimension for it is "units/s"
trapz(Q) = XX (m^2/s) % integerates along the column gives the single value
and that takes it to units^2/s
results= XX*Area (m^4/s)
The area is m/unit * m/unit, so units^2/s * m^2/units^ = m^2/s
Walter Roberson
on 2 May 2020
Can we mutiply the mean data (dF(i)) at a particular area postion by an area dA(i). will it be an issue?
I do not think I understand the question yet.
If you are talking about taking mean(Fg{i}(:)) times the area of one (rotated) unit times the number of rotated units, then mean Fg{i}(:) is sum(Fg{i}(:))/numel(Fg{i}) and when you multiply that by number of rotated units, which is numel(Fg{i}) then the /numel(Fg{i}) and the *numel(Fg{i})) cancel out, giving you sum(Fg{i}(:)) times area of a single rotated unit, so you might as well not bother taking the means.
The difference between sum(Fg{i}(:)) and trapz(trapz(Fg{i},2)) is in the calculation of the boundary conditions. trapz(trapz{i},2)) works out to already include sum(Fg{i}(2:end-1,2:end-1)) plus (sum(Fg{i}(2:end-1,1)) + sum(Fg{i}(2:end-1,end)) + sum(Fg{i}(1,2:end-1)) + sum(Fg{i}(end,2:end-1))) / 2 + something for corners. I have not worked out yet exactly how the corners fit into the calculations... I think it might be their sum divided by 4.
Mathematically, using sum(Fg{i}(:)) is like the assumption that every grid rectangle is a horizontal flat-topped plate with instantaneous change at the boundaries, whereas using trapz(trapz()) is the assumption that each grid rectangle is a flat-topped plate that tilts to meet the neighbours.
MS
on 2 May 2020
Edited: MS
on 2 May 2020
Thank you very much for the detailed answer and the guidance. I have understood that,the units are correct for the results output and the explanation about the idea of sum(Fg{i}(:)).
We need to integerate the (Fg{i}(:)) of revesing areas from 6th to 12th files. It seems, the rectangle area is not revesing back at the end.
But, our results shows that the rectangle is not returning to the orginal postion.The rectangle area should return to the orginal postion as shown in figure. I request you to have a look at it.
Thank you ![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/289648/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/289648/image.png)
Walter Roberson
on 2 May 2020
Count the number of visible rectangles. You will find that there are 6 visible rectangles. If the positions were not reverting, then the other 6 rectangles would be visible. However, the other 6 rectangles are not visible because they are immediately underneath the last 6 rectangles, which have reverted to the exact positions as the earlier rectangles.
For example look for the yellow rectangle, the one that appears first in the legend. You cannot see it, because it is covered over by one of the reverting rectangles.
If you change the linestyle for the second 6 rectangles to ':' or '-.' then you will be able to see the lines underneath, in the gaps in the lines drawn with ':' or '-.'
MS
on 3 May 2020
Edited: MS
on 3 May 2020
Thank you very much. You are absolutely correct. I am trying to use 'sum(Fg{i}(:))' as you suggested to cross check the answers now. It will be very helpful to make a solid conculsion from the flow field. I request you to let me know if you know the better approach.
results1{i} = sum(Fg{i})* (r0(2)-r0(1)) * (c0(2)-c0(1));
I have used the below line as the inputs degrees of rotation. But, i got the resutls as 12by12 results(except first colums all other colums are zeros) instead of 12by1. I am attaching you the results for you reference. Please let me know if you find any mistake.
navg = numel(avg_mat);
r = cell(navg,1);
t = cell(navg,1);
maxi = 90;
t1 = linspace(15, maxi, 6);
t2 = fliplr(t1(1:end));
tvals = [t1, t2];
nt = length(tvals);
results = zeros(navg,nt);
Rg = cell(navg,1);
Cg = cell(navg,1);
Fg = cell(navg,nt);
filenames = cell(navg,1);
for i = 1:navg
filenames{i} = sprintf('avg_val_%d.txt', i);
writematrix(avg_mat{i}, filenames{i});
R = avg_mat{i}(:,1);
C = avg_mat{i}(:,2);
F = avg_mat{i}(:,5);
r0 = linspace(min(R), 0.1, 1000);
c0 = linspace(min(C), 0.1, 1000);
[Rg0, Cg0] = meshgrid(r0, c0);
Xc = .05;%centroid fo the roation
Yc = .05;
%assuming center of rotation is (Xc, Yc) and angle th is degrees
th = tvals(i);
m = makehgtform('translate',Xc,Yc,0, 'zrotate',deg2rad(th), 'translate',-Xc,-Yc,0);
Coords = [Rg0(:), Cg0(:), zeros(numel(Rg0),1), ones(numel(Rg0),1)];
TrCoords = Coords*m.';
Rg{i} = reshape(TrCoords(:,1),size(Rg0));
Cg{i}= reshape(TrCoords(:,2), size(Cg0));
tg = griddata(R, C, F, Rg{i}, Cg{i});
tg(isnan(tg)) = 0;
Fg{i} = tg;
results(i) = trapz(trapz(Fg{i},2)) * (r0(2)-r0(1)) * (c0(2)-c0(1));
results1{i} = sum(Fg{i})* (r0(2)-r0(1)) * (c0(2)-c0(1));
end
writematrix(results, 'integral_val.txt')
writecell(results1, 'integral_val1.txt')
Thank you very much for your philanthropist efforts to help the boys like me. You have been very kind to me to solve the problem. In this process, i have learned bit of coding from you and the logic to tackle the big problem.
Walter Roberson
on 3 May 2020
sum(Fg{i}(:))
MS
on 3 May 2020
Edited: MS
on 4 May 2020
Thanks. i request you to please look at the result document from the below code. It shows the result zeros of 11by11 matrix. is this a flaw of the choosing the limts wrong? i remeber, you have got the same kind of resutls as zeros while testing. Kindly let me know your inputs.
navg = numel(avg_mat);
r = cell(navg,1);
t = cell(navg,1);
maxi = 90;
t1 = linspace(15, maxi, 6);
t2 = fliplr(t1(1:end));
tvals = [t1, t2];
nt = length(tvals);
results = zeros(navg,nt);
Rg = cell(navg,1);
Cg = cell(navg,1);
Fg = cell(navg,nt);
filenames = cell(navg,1);
for i = 1:navg
filenames{i} = sprintf('avg_val_%d.txt', i);
writematrix(avg_mat{i}, filenames{i});
R = avg_mat{i}(:,1);
C = avg_mat{i}(:,2);
F = avg_mat{i}(:,5);
r0 = linspace(min(R), 0.1, 1000);
c0 = linspace(min(C), 0.1, 1000);
[Rg0, Cg0] = meshgrid(r0, c0);
Xc = .05;%centroid fo the roation
Yc = .05;
%assuming center of rotation is (Xc, Yc) and angle th is degrees
th = tvals(i);
m = makehgtform('translate',Xc,Yc,0, 'zrotate',deg2rad(th), 'translate',-Xc,-Yc,0);
Coords = [Rg0(:), Cg0(:), zeros(numel(Rg0),1), ones(numel(Rg0),1)];
TrCoords = Coords*m.';
Rg{i} = reshape(TrCoords(:,1),size(Rg0));
Cg{i}= reshape(TrCoords(:,2), size(Cg0));
tg = griddata(R, C, F, Rg{i}, Cg{i});
tg(isnan(tg)) = 0;
Fg{i} = tg;
results(i) = trapz(trapz(Fg{i},2)) * (r0(2)-r0(1)) * (c0(2)-c0(1));
end
writematrix(results, 'integral_val.txt')
It seems, the rectangle area is deforming visually instead of keeping the shape constant. Is it due to the axis restricton or just a visual illusion? Please let me know your valuable comments.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/290066/image.png)
Walter Roberson
on 4 May 2020
navg = numel(avg_mat);
r = cell(navg,1);
t = cell(navg,1);
maxi = 90;
t1 = linspace(15, maxi, 6);
t2 = fliplr(t1(1:end));
tvals = [t1, t2];
nt = length(tvals);
results = zeros(navg,1);
Rg = cell(navg,1);
Cg = cell(navg,1);
Fg = cell(navg,nt);
filenames = cell(navg,1);
for i = 1:navg
filenames{i} = sprintf('avg_val_%d.txt', i);
writematrix(avg_mat{i}, filenames{i});
R = avg_mat{i}(:,1);
C = avg_mat{i}(:,2);
F = avg_mat{i}(:,5);
r0 = linspace(min(R), 0.1, 1000);
c0 = linspace(min(C), 0.1, 1000);
[Rg0, Cg0] = meshgrid(r0, c0);
Xc = .05;%centroid fo the roation
Yc = .05;
%assuming center of rotation is (Xc, Yc) and angle th is degrees
th = tvals(i);
m = makehgtform('translate',Xc,Yc,0, 'zrotate',deg2rad(th), 'translate',-Xc,-Yc,0);
Coords = [Rg0(:), Cg0(:), zeros(numel(Rg0),1), ones(numel(Rg0),1)];
TrCoords = Coords*m.';
Rg{i} = reshape(TrCoords(:,1),size(Rg0));
Cg{i}= reshape(TrCoords(:,2), size(Cg0));
tg = griddata(R, C, F, Rg{i}, Cg{i});
tg(isnan(tg)) = 0;
Fg{i} = tg;
results(i) = trapz(trapz(Fg{i},2)) * (r0(2)-r0(1)) * (c0(2)-c0(1));
end
writematrix(results, 'integral_val.txt')
and when you do your plotting,
axis equal
which is not the default.
See Also
Categories
Find more on Data Synthesis in Help Center and File Exchange
Tags
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!An Error Occurred
Unable to complete the action because of changes made to the page. Reload the page to see its updated state.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)