# I am trying to obtain multiple variable for equation out of the desired result.

1 view (last 30 days)
Paramjit Yadav on 11 Jul 2021
Commented: Walter Roberson on 27 Jul 2021
Hello Sir/Ma'am,
I need help.
I have one equation which gives answer in matrix of 60x 60.
I want to find the 9 variables existing in the equation. But this varaibles must be set based on result. That is ration must be one.
So we have to ask matlab to select variable in such a way that Ratio is 1.
Ratio value is found in subplot(1 3 2)
Below is my equation.
x = 0:0.1:60;
y = 0:0.1:60;
[X,Y] = meshgrid(x,y)
z = (V(1).*sin(2.*pi./20.*1.*X)) + (V(2).* sin(2.*pi./20.*1.*Y))
+ (V(3).*(sin(2.*pi./20.*1.*X).* sin(2.*pi./20.*1.*Y)))
+ (V(4).*(sin(2.*pi./20.*2.*X))) + (V(5).*sin(2.*pi./20.*2.*Y))
+ (V(6).*sin(2.*pi./20.*3.*X)) + (V(7).*sin(2.*pi./20.*3.*Y))
+ (V(8).*(sin(2.*pi./20.*2.*X).*sin(2.*pi./20.*1.*Y)))
+ (V(9).*(sin(2.*pi./20.*1.*X).* sin(2.*pi./20.*2.*Y)));
Here V(1)......V(9) = are variables
nx = zeros(601,601);
ny = zeros(601,601);
nz = zeros(601,601);
[nx, ny, nz] = surfnorm(X,Y,z);
S1 = [0 0 1];
NX2= 0;
NY2= 0;
NZ2= -1;
for c = 1:601
for r = 1:601
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% First surface %%%%%%%%%%%%%%%%%%%%%%%%%%
NX= -nx(r,c);
NY= -ny(r,c);
NZ= -nz(r,c);
Nsurf = [NX NY NZ]; %dot([0.67 0.67 0.67],
Fa = (0.67*(cross(Nsurf,cross(-Nsurf,S1))));
Fb = (Nsurf*sqrt(1-(dot(0.44,
(dot(cross(Nsurf,S1),cross(Nsurf,S1)))))));
S2 = (Fa-Fb);
R1x(r,c) = S2(1);
R1y(r,c) = S2(2);
R1z(r,c) = S2(3);
Theta2 = asind(S2/norm(S2));
Angle2x(r,c) = Theta2(1);
Angle2y(r,c) = Theta2(2);
Angle2z(r,c) = Theta2(3);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Second surface %%%%%%%%%%%%%%%%%%%%%%%%%%
Nsurf2 = [NX2 NY2 NZ2];
Fa2 = (1.5*(cross(Nsurf2,cross(-Nsurf2,S2))));
Fb2 = (Nsurf2*sqrt(1-(dot(2.25,(dot(cross(Nsurf2,S2),cross(Nsurf2,S2)))))));
S3 = (Fa2-Fb2);
%S3(S3~=real(S3)) = NaN;
R2x(r,c) = S3(1);
R2y(r,c) = S3(2);
R2z(r,c) = S3(3);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Theta4 = asind(S3/norm(S3));
Angle4x(r,c) = Theta4(1);
Angle4y(r,c) = Theta4(2);
Angle4z(r,c) = Theta4(3);
end
end
f = figure();
subplot(1,3,1)
W = histogram2(Angle2x,Angle2y,'DisplayStyle','tile','ShowEmptyBins','on');%,'DisplayStyle','tile','ShowEmptyBins','on'%[60 60]
counts_First_surface = W.Values;
DATA_First_surface = W.Data;
Bincounts_First_surface = W.BinCounts;
UserData_First_surface = W.UserData;
MAX_Firstsurface= max(W.Values)
MIN_First_surface= min(W.Values)
xlabel('Angle with x');
ylabel('Angle with y');
colormap(jet(256));
colorbar;
axis equal
% xlim([-20 20]);
% ylim([-20 20]);
title('First Surface')
subplot(1,3,2)
% h = histogram(x,30,'BinLimits',[-0.3 0.3])
h = histogram2(Angle4x,Angle4y,'DisplayStyle','tile','ShowEmptyBins','on'); %,'DisplayStyle','tile','ShowEmptyBins','on'
counts_Second_surface = h.Values;
DATA_Second_surface = h.Data;
Bincounts_Second_surface = h.BinCounts;
UserData_Second_surface = h.UserData;
MAX_Second_surface =max(h.Values,[],'all')%max(h.Values)
MIN_Second_surface= min(h.Values,[],'all')
Ratio = MIN_Second_surface/MAX_Second_surface
xlabel('Angle with x');
ylabel('Angle with y');
colormap(jet(256));
colorbar;
axis equal
% xlim([-40 40]);
% ylim([-40 40]);
title('Second Surface')
subplot(1,3,3)
contourf(X,Y,z)
title('Countour Plot for surface')
axis equal
end
Paramjit Yadav on 27 Jul 2021
Dear Walter Roberson,
Thank you so much. The vector calculation of Theta3x and Theta3y is extremely good.
I have also noted one point in -
Target_Value = ones(FFbins,FFbins)*(Nsample/FFbins)^2;
(Nsample*FFbins)^2 is (60*600)^2 = 1296000000 . So Target_Value is going to be an array that is all 1/1296000000, roughly 7.7e-10
I am sorry but here I think the (Nsample*FFbins)^2 is (60*600)^2 = 1296000000 is mistaken.
Because it is (Nsample/FFbins)^2 = (600/60)^2 = 100.
Thank you so much Dear Walter Roberson, for such a helping hand.
Kind Regards
Paramjit

Paramjit Yadav on 12 Jul 2021
Hi Walter,
Thank you so much for answering the question.
Taking it further.
Here the Ratio is (Ratio = MIN_Second_surface/MAX_Second_surface)
Now the MIN_Second_surface and MAX_Second_surface is the count of values with in a single bin of histogram2.
Hence, if the MAX_Second_surface = 100 and MIN_Second_surface = 100 i.e every bin has 100 values in it.Then the Ratio of them will be 1.
Hence, this is the condition to take variables(V) such that every bin must have equal values in it.
value in histogram2 was found by h.values and W.values where h and W is the place were histogram2 was created.
h = histogram2(Angle4x,Angle4y,'DisplayStyle','tile','ShowEmptyBins','on');
W = histogram2(Angle2x,Angle2y,'DisplayStyle','tile','ShowEmptyBins','on')
Also in our case the idle value in each bin should be 100 for max and min surface.
Also the bins size can be made equal using
Xedges = [-30:1:30];
Yedges = [-30:1:30];
h = histogram2(x,y,Xedges,Yedges)
Hence I have request to please setting the variable in a way that Ration is 1 or around it.

Walter Roberson on 18 Jul 2021
I have one question. Is their a way by which we can set a condition such that everytime it iterates it should show only the lowest value of MF1 as compared to previous one.
Modified version attached.
I incorporated the best point I found before, so it is likely that it will stall out after the first generation. To reduce that, you could remove P0 (by itself) from the P matrix. If you want to see it even less guided (but possibly taking over 100 iterations, remove the P0+randn()/100 line from the P matrix.
Walter Roberson on 27 Jul 2021
To arrive at those vectors, make small changes to your code:
syms nx1 ny1 nz1 real %NEW
S1 = [0 0 1]; % Incident light being collimated
NX2= 0; % Normal vector of second surface
NY2= 0;
NZ2= -1;
for c = 1:1 %changed
for r = 1:1 %changed
%%%%%%%%%%%%% Normal of First surface surface %%%%%%%%%%%%%%%%%%
NX1= -nx1(r,c);
NY1= -ny1(r,c);
NZ1= -nz1(r,c);
Nsurf1 = [NX1 NY1 NZ1];
%%%%%%%%%%%%%%%%%%% Light vector S2 after first surface %%%%%%%%%%%%%%%%%%%%%
Fa = (1/1.5*(cross(Nsurf1,cross(-Nsurf1,S1))));
Fb = (Nsurf1*sqrt(1-(dot(1/2.25,(dot(cross(Nsurf1,S1),cross(Nsurf1,S1)))))));
S2 = (Fa-Fb);
R1x(r,c) = S2(1);
R1y(r,c) = S2(2);
R1z(r,c) = S2(3);
%%%%%%%%%%%%%%%%%%% Angle of S2 after first surface %%%%%%%%%%%%%%%%%%%%
% Theta2 = asind(S2/norm(S2)); Not needed
Theta2x(r,c) = asind(S2(1)/norm(S2));
Theta2y(r,c) = asind(S2(2)/norm(S2));
%Angle2z(r,c) = Theta2(3); Not needed
%%%%%%%%%%%%% Normal of Second surface surface %%%%%%%%%%%%%%%%%%
Nsurf2 = [NX2 NY2 NZ2];
%%%%%%%%%%%%%%%%%%% Light S3 after 2nd surface %%%%%%%%%%%%%%%%%%%
Fa2 = (1.5*(cross(Nsurf2,cross(-Nsurf2,S2))));
Fb2 = (Nsurf2*sqrt(1-(dot(2.25,(dot(cross(Nsurf2,S2),cross(Nsurf2,S2)))))));
S3 = (Fa2-Fb2);
% S3(S3~=real(S3)) = 0; not needed
R2x(r,c) = S3(1);
R2y(r,c) = S3(2);
R2z(r,c) = S3(3);
%%%%%%%%%%%%%%%%%%% Angle of S3 after 2nd surface %%%%%%%%%%%%%%%%
%Theta4 = asind(S3/norm(S3)); not needed
Theta3x(r,c) = asind(S3(1)/norm(S3));
Theta3y(r,c) = asind(S3(2)/norm(S3));
end
end
then
simplify(Theta3x)
ans =
simplify(Theta3y)
ans =
So by letting nx1, ny1, nz1 be symbolic templates, you can get to expressions for Theta3x and Theta3y . You can then
vectorize(simplify(Theta3x))
ans = '-(180.*asin((nx1.*(2.*nz1 - (9 - 4.*ny1.^2 - 4.*nx1.^2).^(1./2)))./(4.*abs(- (3.*nx1.^2.*(2.*nz1 - (9 - 4.*ny1.^2 - 4.*nx1.^2).^(1./2)).*(3.*abs((4.*nx1.^2)./9 + (4.*ny1.^2)./9 - 1) - (2.*nz1.*(9 - 4.*ny1.^2 - 4.*nx1.^2).^(1./2))./3))./(4.*(9 - 4.*ny1.^2 - 4.*nx1.^2).^(1./2)) - (3.*ny1.^2.*(2.*nz1 - (9 - 4.*ny1.^2 - 4.*nx1.^2).^(1./2)).*(3.*abs((4.*nx1.^2)./9 + (4.*ny1.^2)./9 - 1) - (2.*nz1.*(9 - 4.*ny1.^2 - 4.*nx1.^2).^(1./2))./3))./(4.*(9 - 4.*ny1.^2 - 4.*nx1.^2).^(1./2)) - 1) + nx1.^2.*abs((9 - 4.*ny1.^2 - 4.*nx1.^2).^(1./2) - 2.*nz1).^2 + ny1.^2.*abs((9 - 4.*ny1.^2 - 4.*nx1.^2).^(1./2) - 2.*nz1).^2).^(1./2)))./pi'
vectorize(simplify(Theta3y))
ans = '-(180.*asin((ny1.*(2.*nz1 - (9 - 4.*ny1.^2 - 4.*nx1.^2).^(1./2)))./(4.*abs(- (3.*nx1.^2.*(2.*nz1 - (9 - 4.*ny1.^2 - 4.*nx1.^2).^(1./2)).*(3.*abs((4.*nx1.^2)./9 + (4.*ny1.^2)./9 - 1) - (2.*nz1.*(9 - 4.*ny1.^2 - 4.*nx1.^2).^(1./2))./3))./(4.*(9 - 4.*ny1.^2 - 4.*nx1.^2).^(1./2)) - (3.*ny1.^2.*(2.*nz1 - (9 - 4.*ny1.^2 - 4.*nx1.^2).^(1./2)).*(3.*abs((4.*nx1.^2)./9 + (4.*ny1.^2)./9 - 1) - (2.*nz1.*(9 - 4.*ny1.^2 - 4.*nx1.^2).^(1./2))./3))./(4.*(9 - 4.*ny1.^2 - 4.*nx1.^2).^(1./2)) - 1) + nx1.^2.*abs((9 - 4.*ny1.^2 - 4.*nx1.^2).^(1./2) - 2.*nz1).^2 + ny1.^2.*abs((9 - 4.*ny1.^2 - 4.*nx1.^2).^(1./2) - 2.*nz1).^2).^(1./2)))./pi'