Hey folks; I have written a program that utilises the atan2 function,
only it's giving me some very strange behaviour at certain values; in
essence, I have a grid 1000 x 1000 units, and a circle at the centre
(g,f) = (500,500) with radius 100 units. The program asks the user to
specify an Xvalue between 0 and 1 and a Y  value between 0 and 1,
corresponding to the units. (For example, a user inputting X = 0.3 and
Y = 0.2 would specify the point [300,200]). With this data, the
program uses the atan2 function to work out the angle gamma between
the point (x,y) and the circle centre (g,f). It then calculates the
half angle bisector of the two tangent points, beta. Using this and
some basic trig, it computes the angle of the point to the lower
tangent line (theta) and the upper tangent line (phi). Then it
analyses all 1000 x 1000 points; if they lie between theta and phi,
then they're assigned an I value of 0, otherwise they are assigned an
I value of 1000/r. The code is here;

%Circle tangent code
format long
Irad=zeros([1000 1000]);
warning off all
g = .50; %Circle centre, Xposition metres
f = .50; %Circle centre, Y Position metres
crad = .10; %Circle radius in metres
gsize=100;
A = 1000;
X = input('Xposition (between 0 and 1): ');
Y = input('Yposition (between 0 and 1): ');
tubecent = sqrt((g  X)^2 + (f  Y)^2); %distance from (X,Y) to (g,f)
gamma = atan2((f  Y),(g  X));
beta = asin((crad)/tubecent);
theta = gamma  beta;
phi = gamma + beta;
for j = 1:1000
for k = 1:1000
r = sqrt( ((j/1000)(X)).^2 + ((k/1000)(Y)).^2 );
%r is the distance of the point in question to the point (X,Y)
ang = atan2(((k/1000)  Y),((j/1000)  X));
if ang <= theta  ang >= phi %tests of angle is less than theta
I = A ./ r;
else
I = 0;
end
Ir(j,k) = sum(I);
end
end
[XM,YM]=meshgrid(0.1:0.1:gsize , 0.1:0.1:gsize);
figure(2)
pcolor(YM,XM,log(Ir)), title('Log of Intensity versus distance graph')
xlabel('Centimetres'), ylabel('Centimetres')
shading interp

This works fine, except when you try to run the code at (X,Y) =
(1,0.5). In fact, you get strange results anywhere that satisfies the
following condition.
X > g, and (Y >= f  crad, Y <= f + crad)
Can anyone please help me? I've been trying to weeks to solve this
problem! If you input the same point under symmetry (ie (0,0.5),
(0.5,0), (0.5,1) ) the results are perfect! Any ideas? Thank you so
much in advance.
