simple nested loops error

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

Answer 1

[Walter Roberson]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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.

[MS]

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]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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]

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]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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]

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]

What is numl and how does it differ from numel ?

[MS]

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]

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.

Answer 2

[MS]

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.

[Walter Roberson]

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]

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]

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]

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]

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.

[MS]

thanks, you are correct. The first file should be (90/6) deg instead of 0 and it comes back to 90/6 in twelve files.

[Walter Roberson]

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]

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]

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]

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]

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]

Thank you very much. The code worked really well.

[MS]

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.

[Walter Roberson]

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]

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]

the code I posted already has the rotation coded into it.

[MS]

Awesome.thanks. Sorry, I did not see it.

[MS]

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]

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]

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]

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.

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]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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]

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)

[Walter Roberson]

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]

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]

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

[Walter Roberson]

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]

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]

sum(Fg{i}(:))

[MS]

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.

[Walter Roberson]

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.

[MS]

Thank you very much.