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

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

[Walter Roberson]

V = randn(1,9) * 10

V = 1×9

-14.8698 0.8622 -1.0558 -3.8018 -0.1254 6.4482 -9.0944 -3.8677 11.6803

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)));
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) 

MAX_Firstsurface = 1×18

3030 5694 5445 4431 5337 3573 2745 4005 3132 2841 2757 3144 3333 4374 5250 6039 4488 2905

MIN_First_surface= min(W.Values)

MIN_First_surface = 1×18

0 0 0 0 0 0 0 9 0 27 9 0 0 0 0 0 0 0

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) 

MAX_Second_surface = 5022

MIN_Second_surface= min(h.Values,[],'all')

MIN_Second_surface = 0

Ratio = MIN_Second_surface/MAX_Second_surface

Ratio = 0

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

Okay, so we have that middle plot... but what does it mean to have a ratio of 1 there? Ratio of 1 implies that two quantities should be equal, but which two quantities should be equal for your purposes?

[Paramjit Yadav]

Hi Walter Roberson,

Can you please help wit optimization.

Here when I run the program I get MF1 as 36. Ideally I must have 1 but near to it is also fine.

Inoder to achieve this we have to set the value of V1 to V9 variables in a way that we MF1 = 1 or near to 1.

Also I have bit modified the programme. Please help.

tic
clear

Nsample = 600; % Number of sample points in x and in y direction
Windowsample =60; % Size of samling window [µm], sampling points to be distributed evenly across 
FFbins =60; % Number of bins in histogram of FF

x = linspace( 1/2*Windowsample/Nsample, Windowsample-1/2*Windowsample/Nsample, Nsample); 
y = linspace( 1/2*Windowsample/Nsample, Windowsample-1/2*Windowsample/Nsample, Nsample); 

[X,Y] = meshgrid(x,y);

z= zeros(Nsample,Nsample);

R2x =  zeros(Nsample,Nsample);
R2y =  zeros(Nsample,Nsample);
R2z =  zeros(Nsample,Nsample);

R1x = zeros(Nsample,Nsample);
R1y = zeros(Nsample,Nsample);
R1z = zeros(Nsample,Nsample);

Theta2x = zeros(Nsample,Nsample);
Theta2y = zeros(Nsample,Nsample);
%Angle2z = zeros(Nsample,Nsample); Not needed

Theta3x = zeros(Nsample,Nsample);
Theta3y = zeros(Nsample,Nsample);
%Angle4z = zeros(Nsample,Nsample); Not needed

% 0.030208


Target_Value = ones(FFbins,FFbins)*(Nsample/FFbins)^2;
% z = V1*sin(2*pi/20*1*X)  + V2* sin(2*pi/20*1*Y)...
%      + -V3*sin(2*pi/20*1*X).*sin(2*pi/20*1*Y) + V4*sin(2*pi/20*3*X)  +  V5*sin(2*pi/20*3*Y)...
%      + V6*sin(2*pi/20*5*X) + V7*sin(2*pi/20*5*Y) + -V8*sin(2*pi/20*3*X).*sin(2*pi/20*1*Y)...
%      + -V9*sin(2*pi/20*1*X).* sin(2*pi/20*3*Y); % z in [µm]        %MARKED

% I have taken few manual values of V1 to V9 to Optmize MF1 (at the bottom) as much I can.  %MARKED

z = 3.7019*sin(2*pi/20*1*X)  + 3.700191* sin(2*pi/20*1*Y)...
     + -0.015*sin(2*pi/20*1*X).*sin(2*pi/20*1*Y) + 0.2000*sin(2*pi/20*3*X)  +  0.200*sin(2*pi/20*3*Y)...
     + 0.042016*sin(2*pi/20*5*X) + 0.0553*sin(2*pi/20*5*Y) + -0.10001*sin(2*pi/20*3*X).*sin(2*pi/20*1*Y)...
     + -0.200005*sin(2*pi/20*1*X).* sin(2*pi/20*3*Y); % z in [µm]    %MARKED

nx1 = zeros(Nsample,Nsample); % Normal vector of first surface
ny1 = zeros(Nsample,Nsample);
nz1 = zeros(Nsample,Nsample);

[nx1, ny1, nz1] = surfnorm(X,Y,z);

S1 = [0 0 1]; % Incident light being collimated
 
NX2=  0; % Normal vector of second surface
NY2=  0;
NZ2=  -1;
 
for c = 1:600 
    for r = 1:600
    
        %%%%%%%%%%%%% 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));
        %Angle4z(r,c) = Theta4(3); Not needed             
    end
end
toc
f = figure(1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% First  surface Plot  %%%%%%%%%%%%%%%%%%%%%%%%%%        
        
subplot(1,3,1)       
Xedges = [-30:1:30];
Yedges = [-30:1:30];
W1 = histogram2(Theta2x,Theta2y,Xedges,Yedges,'DisplayStyle','tile','ShowEmptyBins','on');%,'DisplayStyle','tile','ShowEmptyBins','on'

counts_First_surface = W1.Values ;%%%%%%%%%%%%%% Data in each Bin  %%%%%%%%%%%%%%%%%%
Bincounts_First_surface = W1.BinCounts;

xlabel('x Angle with yz plane [°]');
ylabel('y Angle with xz plane [°]');
colormap(jet(256));
colorbar;
axis equal    
xlim([-35 35]);
ylim([-35 35]);    
title('FF after first surface')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Second  surface Plot  %%%%%%%%%%%%%%%%%%%%%%

subplot(1,3,2) 
Xedges = [-30:1:30];
Yedges = [-30:1:30];

W2 = histogram2(Theta3x,Theta3y,Xedges ,Yedges,'DisplayStyle','tile','ShowEmptyBins','on'); %,'DisplayStyle','tile','ShowEmptyBins','on'

counts_Second_surface = W2.Values; %%%%%%%%%%%%%% Data in each Bin  %%%%%%%%%%%%%%%%%%
Bincounts_Second_surface = W2.BinCounts;
MAX_Second_surface =max(W2.Values,[],'all');%max(h.Values) 
MIN_Second_surface= min(W2.Values,[],'all');

Ratio = MIN_Second_surface/MAX_Second_surface; %%%%%%%%%%%%%%%%%%%% Ratio of Max and Min  %%%%%%%%%%%%%%%%%%

Real_Value = (W2.Values);                       %%%%%%%%%%%%%% Merit Function rms of differences  %%%%%%%%%%%%%%%%%%
SamplesinFF = sum(Real_Value,'all')                   % Test if all rays are getting into FF
Merit_Function = Real_Value - Target_Value;
MF1 =  sqrt( sum(Merit_Function.*Merit_Function,'all')/FFbins/FFbins)   %MARKED

xlabel('x Angle with yz plane [°]');
ylabel('y Angle with xz plane [°]');
colormap(jet(256));
colorbar;
axis equal    
xlim([-35 35]);
ylim([-35 35]);          
title('FF after second surface')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Contour surface Plot  %%%%%%%%%%%%%%%%%%%% 
subplot(1,3,3) 
contourf(X,Y,z)
xlabel('x [µm]');
ylabel('y [µm]');
title('Contour of first surface')
colorbar;
axis equal
xlim([-5 65]);
ylim([-5 65]);   
toc

[Walter Roberson]

You are calculating with histograms, and because you can get the same histogram for conditions that are not (exactly) the same, then your calculation is not differentiable -- and so you cannot use fmincon() or any optimizer that relies upon differentiation or upon gradient, or Newton-Raphson type techniques.

I have it going calculating based upon ga(). It is improving relative to where it started, but it is taking a while. It is not especially close to your designed 1, though.

Current best is near

-0.067418510588264 -0.0929171153276874 1.29078839763246 0.368323624975687 1.068752729548 0.458787283054125 -0.029959986835674 1.5907867966087 -0.172926961924405

which has a MF1 in the range of 38

[Walter Roberson]

Current best is near

-0.067052299650764 -0.0929071613520099 1.29054425700746 0.368615210897928 1.068752729548 0.45903336083008 -0.0298500758707082 1.59079219082905 -0.172926961924405

but it is not a lot better, only about 0.1 better MF1. It looks likely that it is going to stall out.

[Walter Roberson]

It stalled at

-0.067052299650764 -0.0930292316645099 1.29054425700746 0.368615210897928 1.068752729548 0.45903336083008 -0.0298500758707082 1.59079219082905 -0.172926961924405

I will attach the code. You will want to modify the P matrix to have values close to the above.

[Paramjit Yadav]

Dear Walter Roberson,

Thank you so much for such a nice help. It is really helpful for me to understand how the optimization will work in our case.

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.

I am thankful to you once again.

[Paramjit Yadav]

Dear Walter Roberson,

Thank you so much for providing on how the Genetic algortihm can be applied.

I have tried my best to find out the solution near 1 or 0 for MF1. I have iterated it over the night. But I am really sorry I am not able to get it. Can you please help to look in to the matter and help me to achieve something near to 0 or 1 for MF1. I have tried the whole week.

I have also modified the function as below to get better answers for MF1 near to 1 or 0.

z = V1*cos(2*pi/20*1*X)  + V2* cos(2*pi/20*1*Y)...
     + V3*cos(2*pi/20*1*X).*cos(2*pi/20*1*Y) + V4*cos(2*pi/20*3*X)  +  V5*cos(2*pi/20*3*Y)...
     + V6*cos(2*pi/20*5*X) + V7*cos(2*pi/20*5*Y) + V8*cos(2*pi/20*3*X).*cos(2*pi/20*1*Y)...
     + V9*cos(2*pi/20*1*X).* cos(2*pi/20*3*Y); % z in [µm]

I am very thankfull to you, Please help me to get something near to 1 or 0 for MF1. Please.

Kind Regards

Paramjit

[Paramjit Yadav]

Also I have achieved uptil 34 MF1 on the variable(V) as below.

3.7024 3.8881 -0.3572 -0.1537 -0.2530 0.0414 0.0571 -0.0760 0.2001

Thanking you,

Kind Regards

Paramjit

[Paramjit Yadav]

Dear Walter Roberson,

I have also noticed one thing that the upper and lower boundary are set to -10 to 10 but it does not takes the value from it. It only revolves near to first intial supplied matirx. I am not able to understand why it is doing so.

Please help to clear the doubts.

Thank you so much Dear Walter Roberson, for your help.

Kind Regards

Paramjit

[Walter Roberson]

If you want the full range from -10 to +10 to be examined, then do not supply an InitialPopulation, or make the population size larger than the size of the InitialPopulation .

If you do not supply an initial population, then you should expect a minimum of 90 iterations of it just looking around in areas that are not at all close to the solution, and another 150 after while it moves towards better solutions.

[Walter Roberson]

As an outsider, I have no reason to expect that MF even close to 1 is possible.

[Paramjit Yadav]

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

Answer 1

[Paramjit Yadav]

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.

Answer 2

[Walter Roberson]

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]

Nsample = 600; % Number of sample points in x and in y direction

FFbins =60; % Number of bins in histogram of FF

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

W2 = histogram2(Theta3x,Theta3y,Xedges ,Yedges,'DisplayStyle','tile','ShowEmptyBins','on'); %,'DisplayStyle','tile','ShowEmptyBins','on'

The default Normalization for histogram2 is 'count'; see https://www.mathworks.com/help/matlab/ref/matlab.graphics.chart.primitive.histogram2-properties.html#buvhqe4-Normalization

Real_Value = (W2.Values);                       %%%%%%%%%%%%%% Merit Function rms of differences  %%%%%%%%%%%%%%%%%%

When the Normalization property is 'count', as it is by default, then the Values property is the histogram counts; see https://www.mathworks.com/help/matlab/ref/matlab.graphics.chart.primitive.histogram2-properties.html#buvhqe4-Values

Merit_Function = Real_Value - Target_Value;

As discussed, above, Target_value is a matrix whose values are all about 7.7e-10 . And as discussed above, Real_Value is actual histogram counts, which are integers.

MF1 =  sqrt( sum(Merit_Function.*Merit_Function,'all')/FFbins/FFbins)   %MARKED

Does this really make sense to you? To subtract 7.7e-10 from counts, and calculate the sum, and divide that all by 3600, and hope that the result is 0??

Consider what it would mean for there to be a perfect fit, MF1 exactly 0: all of the counts would have to be 7.7e-10 . That isn't going to happen; it can't happen.

The lowest hypothetically achievable MF1 would be if the counts were all 0, in which case the merit function entries would be -7.7e-10, and assuming that you have FFbins by FFbins of those, the MF1 would end up as 7.7e-10 . Sounds promising in terms of target, right? But for the counts to be 0, that would mean that the X or Y values would have to be out of range or complex valued. Are you sure that is what you want to optimize for, maximum avoidance of the area?

Are you aware that your code calculates and displays the first histogram, but otherwise never uses the results of the calculation? And that therefore to improve performance you should comment out that histogram creation? And that once you have done that, then to improve performance you should also comment out calculation of Theta2x and Theta2y ?

You can also vectorize the calculation of Theta3x and Theta3y, which are

Theta3x = -(180.*asin(((3.*NX1.*(1 - (4.*NY1.^2)./9 - (4.*NX1.^2)./9).^(1./2))./2 + NX1.*NZ1)./(abs((9.*(NX1.*conj((1 - (4.*NY1.^2)./9 - (4.*NX1.^2)./9).^(1./2)) + (2.*NX1.*NZ1)./3).*(NX1.*(1 - (4.*NY1.^2)./9 - (4.*NX1.^2)./9).^(1./2) + (2.*NX1.*NZ1)./3))./4 + (9.*(NY1.*conj((1 - (4.*NY1.^2)./9 - (4.*NX1.^2)./9).^(1./2)) + (2.*NY1.*NZ1)./3).*(NY1.*(1 - (4.*NY1.^2)./9 - (4.*NX1.^2)./9).^(1./2) + (2.*NY1.*NZ1)./3))./4 - 1) + abs((3.*NX1.*(1 - (4.*NY1.^2)./9 - (4.*NX1.^2)./9).^(1./2))./2 + NX1.*NZ1).^2 + abs((3.*NY1.*(1 - (4.*NY1.^2)./9 - (4.*NX1.^2)./9).^(1./2))./2 + NY1.*NZ1).^2).^(1./2)))./pi
Theta3y = -(180.*asin(((3.*NY1.*(1 - (4.*NY1.^2)./9 - (4.*NX1.^2)./9).^(1./2))./2 + NY1.*NZ1)./(abs((9.*(NX1.*conj((1 - (4.*NY1.^2)./9 - (4.*NX1.^2)./9).^(1./2)) + (2.*NX1.*NZ1)./3).*(NX1.*(1 - (4.*NY1.^2)./9 - (4.*NX1.^2)./9).^(1./2) + (2.*NX1.*NZ1)./3))./4 + (9.*(NY1.*conj((1 - (4.*NY1.^2)./9 - (4.*NX1.^2)./9).^(1./2)) + (2.*NY1.*NZ1)./3).*(NY1.*(1 - (4.*NY1.^2)./9 - (4.*NX1.^2)./9).^(1./2) + (2.*NY1.*NZ1)./3))./4 - 1) + abs((3.*NX1.*(1 - (4.*NY1.^2)./9 - (4.*NX1.^2)./9).^(1./2))./2 + NX1.*NZ1).^2 + abs((3.*NY1.*(1 - (4.*NY1.^2)./9 - (4.*NX1.^2)./9).^(1./2))./2 + NY1.*NZ1).^2).^(1./2)))./pi

Those have a number of common sub-expressions so the calculation can be made even more efficient. The idea is that you would be able to take the nx1, ny1, nz1 and directly calculate the arrays of Theta3x, Theta3y values, instead of looping doing that. A function for these could be generated ahead of time, saving calculation time compared to the current looping.

[Paramjit Yadav]

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]

Dear Walter Roberson,

I am not able to understand completely-

Theta3x = -(180.*asin(((3.*NX1.*(1 - (4.*NY1.^2)./9 - (4.*NX1.^2)./9).^(1./2))./2 + NX1.*NZ1)./(abs((9.*(NX1.*conj((1 - (4.*NY1.^2)./9 - (4.*NX1.^2)./9).^(1./2)) + (2.*NX1.*NZ1)./3).*(NX1.*(1 - (4.*NY1.^2)./9 - (4.*NX1.^2)./9).^(1./2) + (2.*NX1.*NZ1)./3))./4 + (9.*(NY1.*conj((1 - (4.*NY1.^2)./9 - (4.*NX1.^2)./9).^(1./2)) + (2.*NY1.*NZ1)./3).*(NY1.*(1 - (4.*NY1.^2)./9 - (4.*NX1.^2)./9).^(1./2) + (2.*NY1.*NZ1)./3))./4 - 1) + abs((3.*NX1.*(1 - (4.*NY1.^2)./9 - (4.*NX1.^2)./9).^(1./2))./2 + NX1.*NZ1).^2 + abs((3.*NY1.*(1 - (4.*NY1.^2)./9 - (4.*NX1.^2)./9).^(1./2))./2 + NY1.*NZ1).^2).^(1./2)))./pi

Theta3y = -(180.*asin(((3.*NY1.*(1 - (4.*NY1.^2)./9 - (4.*NX1.^2)./9).^(1./2))./2 + NY1.*NZ1)./(abs((9.*(NX1.*conj((1 - (4.*NY1.^2)./9 - (4.*NX1.^2)./9).^(1./2)) + (2.*NX1.*NZ1)./3).*(NX1.*(1 - (4.*NY1.^2)./9 - (4.*NX1.^2)./9).^(1./2) + (2.*NX1.*NZ1)./3))./4 + (9.*(NY1.*conj((1 - (4.*NY1.^2)./9 - (4.*NX1.^2)./9).^(1./2)) + (2.*NY1.*NZ1)./3).*(NY1.*(1 - (4.*NY1.^2)./9 - (4.*NX1.^2)./9).^(1./2) + (2.*NY1.*NZ1)./3))./4 - 1) + abs((3.*NX1.*(1 - (4.*NY1.^2)./9 - (4.*NX1.^2)./9).^(1./2))./2 + NX1.*NZ1).^2 + abs((3.*NY1.*(1 - (4.*NY1.^2)./9 - (4.*NX1.^2)./9).^(1./2))./2 + NY1.*NZ1).^2).^(1./2)))./pi

Can you please give some refence from where I can find this.

Thanking you,

Kind Regards

Paramjit

[Walter Roberson]

Ah, you are right; I am not sure how I multiplied instead of dividing.

[Walter Roberson]

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'