1 view (last 30 days)
Paramjit Yadav on 18 Jul 2021
Commented: Alan Weiss on 19 Jul 2021
Hi ,
Here when I run the program and 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.
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]
% I have taken few manual values of V1 to V9 to Optmize MF1 (at the bottom) as much I can.
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]
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)
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
Alan Weiss on 19 Jul 2021
It would be much easier to read your question if you would mark up the code portions as code. Use the CODE section of the formatting bar at the top of the editing pane. Code should look like this:
j = 0;
for i = 1:1000
j = j + i;
end
Alan Weiss
MATLAB mathematical toolbox documentation