|
I have a code that allows the user to specify the position of a light
source and an ellipse-shaped obstruction. To keep this post concise,
I've omitted details about how all the equations are solved and merely
put in the results that lead to the IF loop; To summate, the code
specifies the angular extent between the light at XY and the tangent
points of the ellipse, and all points outside this angular extent are
illuminated (This is shown in figure 2 in the code).
To compute the irradiances INSIDE the angular extent (while all points
in and beyond the ellipse are not illuminated, the points before the
ellipse are and this makes the question tricky) the code finds a
circumcentre to the tangent points and XY, and the angular extent of
these points. It then works out the radius of this circumcircle. It
draws a modified arc between the tangent points, with the
circumcentre as the centre and the arc modified to with a cosine
condition to be a straightline (shown in Figure 3 in the code)
These conditions are combined in Figure I, but as you can see, the
most right tangent line is no longer in agreement with the beta
condition..in fact, light "bends" around the ellipse at one point and
I wish to eradicate this error. Can anyone help? I know this is messy
coding but I really need an answer soon. Thanks in advance, David.
---------------------------------------------------------------------------------------
format long
Irad=zeros([1000 1000]);
warning off all
clear
H = 1.75;
P = 100;
gsize=100;
a = .15; %condition for ellipse
b = .05; %condition for ellipse
h = .5; %condition for ellipse
k = .5; %condition for ellipse
A = P/(2*pi*H); %Power rule
X = 0.1; %position of tube X co-ord
Y = 0.1; %position of tube Y-co-ord
%Circumcentre of two Tangent points and tube point
w = 0.45750000000000
v = 0.16250000000000
%Radius of circle
Rcirc = 0.36292216796443;
%Lambda, the relative angle between the circumcentre and tangent
points
cosinelambda = 0.92007972078891;
lambda = 0.40251238126841;
%Gamma, the angle between XY and the bisector of the tangent line
gamma = 0.81660756839543;
%The angle beta, between XY and the tangent points
beta = 0.20125619063420;
%The real angle between the midpoint of the tangent lines and
circumcentre
gammae = 1.46013910562100;
%Loop 1
for n = 1:1000
for m = 1:1000
rp = sqrt( ((n/1000)-(X)).^2 + ((m/1000)-(Y)).^2 ); % distance from
XY to any point
D = sqrt( ((n/1000)-(w)).^2 + ((m/1000)-(v)).^2 ); % distance from wv
to any point circumcentre
eq = (a^-2)*((n/1000)^2) - (a^-2)*(2*h*(n/1000)) + (b^-2)*((m/1000)^2)
- (b^-2)*(2*k*(m/1000)) + ((a^-2)*(h^2) + (b^-2)*(k^2) - 1);
ang = atan2(((m/1000) - Y),((n/1000) - X)); %angle between XY and
point
delta = atan2(((m/1000) - v),((n/1000) - w)); % angle between WV and
point
delta_eff = abs(gammae - atan2(((m/1000) - v),((n/1000) - w))); %
angle between WV and point relative lambda
cosdel = abs(cos(delta_eff));
ratio = cosinelambda / cosdel;
R_e = Rcirc*ratio;
if abs(mod(gamma -ang +pi,2*pi)-pi) < beta && (abs(mod(gammae -delta
+pi,2*pi)-pi) < lambda && D > R_e) || eq <= 0
I = 0;
else
I = A ./ rp;
end
Irad(n,m) = sum(I);
end
end
[XM,YM]=meshgrid(0.1:0.1:gsize , 0.1:0.1:gsize);
figure(1)
pcolor(YM,XM,log(Irad)), title('Log of Intensity versus distance
graph')
xlabel('Centimetres'), ylabel('Centimetres')
shading interp
%Loop 2 Demonstrating the beta confines
for n = 1:1000
for m = 1:1000
rp = sqrt( ((n/1000)-(X)).^2 + ((m/1000)-(Y)).^2 ); % distance from
XY to any point
D = sqrt( ((n/1000)-(w)).^2 + ((m/1000)-(v)).^2 ); % distance from wv
to any point circumcentre
eq = (a^-2)*((n/1000)^2) - (a^-2)*(2*h*(n/1000)) + (b^-2)*((m/1000)^2)
- (b^-2)*(2*k*(m/1000)) + ((a^-2)*(h^2) + (b^-2)*(k^2) - 1);
ang = atan2(((m/1000) - Y),((n/1000) - X)); %angle between XY and
point
delta = atan2(((m/1000) - v),((n/1000) - w)); % angle between WV and
point
cosdel = abs(cos(atan2((m/1000) - v,(n/1000) - w)));
cosdel2 = abs(cos(delta));
R_e = Rcirc*(cosinelambda/cosdel2);
if abs(mod(gamma -ang +pi,2*pi)-pi) < beta
I = 0;
else
I = A ./ rp;
end
Irad2(n,m) = sum(I);
end
end
[XM,YM]=meshgrid(0.1:0.1:gsize , 0.1:0.1:gsize);
figure(2)
pcolor(YM,XM,log(Irad2)), title('Log of Intensity versus distance
graph')
xlabel('Centimetres'), ylabel('Centimetres')
shading interp
%Loop 3 The extent of lambda
for n = 1:1000
for m = 1:1000
rp = sqrt( ((n/1000)-(X)).^2 + ((m/1000)-(Y)).^2 ); % distance from
XY to any point
D = sqrt( ((n/1000)-(w)).^2 + ((m/1000)-(v)).^2 ); % distance from wv
to any point circumcentre
eq = (a^-2)*((n/1000)^2) - (a^-2)*(2*h*(n/1000)) + (b^-2)*((m/1000)^2)
- (b^-2)*(2*k*(m/1000)) + ((a^-2)*(h^2) + (b^-2)*(k^2) - 1);
ang = atan2(((m/1000) - Y),((n/1000) - X)); %angle between XY and
point
delta = atan2(((m/1000) - v),((n/1000) - w)); % angle between WV and
point
delta_eff = abs(gammae - atan2(((m/1000) - v),((n/1000) - w))); %
angle between WV and point relative lambda
cosdel = abs(cos(delta_eff));
ratio = cosinelambda / cosdel;
R_e = Rcirc*ratio;
if abs(mod(gammae -delta +pi,2*pi)-pi) < lambda && D > R_e
I = 0;
else
I = A ./ rp;
end
Irad3(n,m) = sum(I);
end
end
[XM,YM]=meshgrid(0.1:0.1:gsize , 0.1:0.1:gsize);
figure(3)
pcolor(YM,XM,log(Irad3)), title('Log of Intensity versus distance
graph')
xlabel('Centimetres'), ylabel('Centimetres')
shading interp
|
Feed for this Thread