atan2

[tenochtitlanuk]

Many of us have found atan2, which takes variables y and x and gives an angle, to be very useful. I know I've coded it myself, and seen/used versions by Chris Iverson /Anatoly /B+ /Steelweaver /Rod /etc. Trouble is it is only needed infrequently.

Recently, playing with arcs, I realised there are two common algorithm /definition versions in circulation. Checking on wikipedia- atan2 I find there is little consistency across programming languages. Up to 8 versions can be considered to exist, depending on whether you regard angles as being measured clockwise or anti-clockwise, and from which axis/direction regarded as zero. And you can define atan2 as giving results in the range 0 to 2 pi, or from -pi/2 to +pi/2!!

No wonder I've sometimes got weird results and had to do a bit of fiddling! Just another symptom of my aging brain!

[Rod]

Hmm... Here is a statement I would like to understand better.

You don't have to use atan2 to calculate the angle between two vectors. If you just want the quickest way, you can use dot (v1, v2)=|v1|*|v2|*cos A to get. A = Math.acos ( dot (v1, v2)/ (v1.length ()*v2.length ()) );

In my ball collision code I use atan2 but according to the above statement I dont need to! So how to remove atan2 from this 2D collision code.

  '------------------------------------------
   '2D BILLIARD STYLE COLLISION
   '------------------------------------------
   'converted from a Blitz BASIC demo
   'By Joseph 'Phish' Humfrey
   '
   'converted with the help of
   'Anatoly and Stefan



   nomainwin
   WindowWidth  = 800+themeWidth
   WindowHeight = 600+themeHeight
   UpperLeftX   = (DisplayWidth-WindowWidth)/2
   UpperLeftY   = (DisplayHeight-WindowHeight)/2
   graphicbox #1.g, 0,0,800,600
   open "2D Bumping" for graphics_nf_nsb as #1
   #1 "trapclose [quit]"

   #1.g "down ; fill black ; flush"






   [setup]
   width=800
   height=600
   friction=.99
   maxballs=15
   dim ball(maxballs,6) 'x y delta x delta y radius mass
   b.x=1   'ball x
   b.y=2   'ball y
   b.dx=3  'ball delta x
   b.dy=4  'ball delta y
   b.r=5   'ball radius
   b.m=6   'ball mass

   b.i=0   'ball index



   'set out red balls in a triangle
   ballTriangleSize=5
   for xloop = ballTriangleSize to 1 step -1
       for yloop = 1 to xloop
           b.i=b.i+1
           ball(b.i,b.x) = (5-xloop)*27 + 200
           ball(b.i,b.y) = yloop*31-(xloop*31)/2.0 + 300
           ball(b.i,b.dx)= 0
           ball(b.i,b.dy)= 0
           ball(b.i,b.r) = 15
           ball(b.i,b.m) = 50
       next
   next

   'set out white cue ball
   cue=0
   ball(cue,b.x) = 600
   ball(cue,b.y) = 300 +20
   ball(cue,b.dx)= -20
   ball(cue,b.dy)= 2-int(rnd(0)*4+.5)
   ball(cue,b.r) = 15
   ball(cue,b.m) = 50




   #1.g "when characterInput [setup]"
   #1.g "setfocus"

   timer 17,[updateballs]
   wait



   [updateballs]
   for b.i= 0 to maxballs

       'update positions
       ball(b.i,b.x)=ball(b.i,b.x)+ball(b.i,b.dx)
       ball(b.i,b.y)=ball(b.i,b.y)+ball(b.i,b.dy)

       'gradually slow down
       ball(b.i,b.dx)=ball(b.i,b.dx)*friction
       ball(b.i,b.dy)=ball(b.i,b.dy)*friction

       '------------------------------------------
       'COLLISION CHECKING
       '------------------------------------------
       ' Check each ball against every other

       for b= b.i to maxballs

           ' not against itself
           if b<>b.i then

           collisionDistance = ball(b.i,b.r)+ball(b,b.r)
           actualDistance = sqr((ball(b,b.x)-ball(b.i,b.x))^2 + (ball(b,b.y)-ball(b.i,b.y))^2)

           ' the collision distance is the some of the radius of both balls
           ' if the actual distance appart is less than that then they have collided
           If actualDistance<collisionDistance Then




               ' calculate the angle the colliding ball is travelling at
               collNormalAngle=atan2(ball(b,b.y)-ball(b.i,b.y),ball(b,b.x)-ball(b.i,b.x))

               ' re position the balls as exactly touching, no intersection
               ' this prevents overlapping balls being seen as still colliding
               moveDist1=(collisionDistance-actualDistance)*(ball(b,b.m)/(ball(b.i,b.m)+ball(b,b.m)))
               moveDist2=(collisionDistance-actualDistance)*(ball(b.i,b.m)/(ball(b.i,b.m)+ball(b,b.m)))
               ball(b.i,b.x)=ball(b.i,b.x) + moveDist1*Cos(collNormalAngle+3.14159)
               ball(b.i,b.y)=ball(b.i,b.y) + moveDist1*Sin(collNormalAngle+3.14159)
               ball(b,b.x)=ball(b,b.x) + moveDist2*Cos(collNormalAngle)
               ball(b,b.y)=ball(b,b.y) + moveDist2*Sin(collNormalAngle)


               '------------------------------------------
               'COLLISION RESPONSE
               '------------------------------------------
               'n = vector connecting the centers of the circles.
               'we are finding the components of the normalised vector n
               nX=Cos(collNormalAngle)
               nY=Sin(collNormalAngle)

               'now find the length of the components of each movement vectors
               'along n, by using dot product.
               a1 = ball(b.i,b.dx)*nX  +  ball(b.i,b.dy)*nY
               a2 = ball(b,b.dx)*nX +  ball(b,b.dy)*nY

               'optimisedP = 2(a1 - a2)
               '             ----------
               '              m1 + m2
               optimisedP = (2.0 * (a1-a2)) / (ball(b.i,b.m) + ball(b,b.m))

               'now find out the resultant vectors
               'r1 = c1\v - optimisedP * mass2 * n
               ball(b.i,b.dx) = ball(b.i,b.dx) - (optimisedP*ball(b,b.m)*nX)
               ball(b.i,b.dy) = ball(b.i,b.dy) - (optimisedP*ball(b,b.m)*nY)

               'r2 = c2\v - optimisedP * mass1 * n
               ball(b,b.dx) = ball(b,b.dx) + (optimisedP*ball(b.i,b.m)*nX)
               ball(b,b.dy) = ball(b,b.dy) + (optimisedP*ball(b.i,b.m)*nY)

           end If
           end if
       next

       'Simple Bouncing off walls.
       if ball(b.i,b.x)<ball(b.i,b.r) then
           ball(b.i,b.x)=ball(b.i,b.r)
           ball(b.i,b.dx)=ball(b.i,b.dx)*-0.9
       end if
       if ball(b.i,b.x)>width-ball(b.i,b.r) then
           ball(b.i,b.x)=width-ball(b.i,b.r)
           ball(b.i,b.dx)=ball(b.i,b.dx)*-0.9
       end if
       if ball(b.i,b.y)<ball(b.i,b.r) then
           ball(b.i,b.y)=ball(b.i,b.r)
           ball(b.i,b.dy)=ball(b.i,b.dy)*-0.9
       end if
       if ball(b.i,b.y)>height-ball(b.i,b.r) then
           ball(b.i,b.y)=height-ball(b.i,b.r)
           ball(b.i,b.dy)=ball(b.i,b.dy)*-0.9
       end if

   next


   [drawballs]
   #1.g "discard ; redraw"
   for b.i= 0 to maxballs
       if b.i=0 then
           #1.g "color white ; backcolor white"
       else
           #1.g "color red ;  backcolor red"
       end if
       #1.g "place ";ball(b.i,b.x);" ";ball(b.i,b.y)
       #1.g "circlefilled "; ball(b.i,b.r)
       #1.g "color white ;\";b.i
   next

   wait

   [quit]
   close #1
   end


function atan2(y,x)
   pi = asn(1) * 2
   if x <> 0 then arctan = atn(y/x)

   select case
       case x > 0
       atan2 = arctan

       case y>=0 and x<0
       atan2 = pi + arctan

       case y<0 and x<0
       atan2 = arctan - pi

       case y>0 and x=0
       atan2 = pi / 2

       case y<0 and x=0
       atan2 = pi / -2
   end select
end function





[tsh73]

In my ball collision code I use atan2 but according to the above statement I dont need to!

I fail to see how
A = Math.acos ( dot (v1, v2)/ (v1.length ()*v2.length ()) );
is any better then
collNormalAngle=atan2(ball(b,b.y)-ball(b.i,b.y),ball(b,b.x)-ball(b.i,b.x))
We don't have "dot()" or "length()" functions in LB
so we will end up using more user defined functions
Now, what for?

[Rod]

Well my trig is poor and I think I understand vectors and vector maths a little better. Maybe I should put the time in to better understand atan2.

But the statement did say the vector maths was quicker.

[bplus]

Ha! My trig is far better than my vectors.

I base my x, y coordinates around point cx, cy like this:
x = cx + r * cos(RadianAngle)
y = cy + r * sin(RadianAngle)
r = radial distance between x, y and cx, cy or r = SQR((x-cx)^2 + (y-cy)^2)

When RadianAngle = 0 then Cos(0) = 1, Sin(0) = 0 this graphs due East of point cx, cy
When RadianAngle = Pi/2 = 90 degrees Cos(RadianAngle) = 0, Sin(RadianAngle) = 1 that graphs due South of cx, cy and moving from East to South is Clockwise around cx, cy and the angle increases from 0 to pi/2 (or if you prefer, 0 to 90 Degrees) an onwards around the point in 2*pi Radians or 360 Degrees.


For Atan2, I've used Andy Amaya's code with great results:


Function Atan2(y, x)
   'Atan2 is a function which determines the angle between points
   'x1, y1 and x2, y2. The angle returned is in radians
   'The angle returned is always in the range of
   '-PI to PI radians (-180 to 180 degrees)
   '==============================================================
   'NOTE the position of Y and X arguments
   'This keeps Atan2 function same as other language versions
   '==============================================================
  If x = 0 Then
      If y < 0 Then
           Atan2 = -1.5707963267948967
       Else
           Atan2 = 1.5707963267948967
       End If
  Else
      chk = atn(y/x)
      If x < 0 Then
          If y < 0 Then
               chk = chk - 3.1415926535897932
           Else
               chk = chk + 3.1415926535897932
           End If
      End If
       Atan2 = chk
  End If
  'thanks Andy Amaya
End Function



Sometimes when comparing angles, you add 2*Pi to Radian Angle if < 0 so both angles compared are positive.
Whether you use all pos or positive and negative angles, they graph to same point using the Trig functions I went over above.

And John is right, arc drawing in QB64 is backwards, so I had to create my own arc drawing Sub to draw arcs that match Basic's Trig functions.
And I use this code in JB when translating:


global PI, DEG, RAD  ' trig constants
PI = acs(-1)
DEG = 180 / PI
RAD = PI / 180

'draws an arc with center at xCenter, yCenter, radius from center is arcRadius
sub arc xCenter, yCenter, arcRadius, dAStart, dAMeasure
   'notes:
   'you may want to adjust size and color for line drawing
   'using angle measures in degrees to match Just Basic ways with pie and piefilled
   'this sub assumes drawing in a CW direction if dAMeasure positive

   'for Just Basic angle 0 degrees is due East and angle increases clockwise towards South

   'dAStart is degrees to start Angle, due East is 0 degrees

   'dAMeasure is degrees added (Clockwise) to dAstart for end of arc

   rAngleStart = RAD * dAStart
   rAngleEnd = RAD * dAMeasure + rAngleStart
   Stepper = RAD / (.1 * arcRadius) 'fixed
   lastX = xCenter + arcRadius * cos(rAngleStart)
   lastY = yCenter + arcRadius * sin(rAngleStart)
   #gr "set ";int(lastX);" ";int(lastY)
   for rAngle = rAngleStart+Stepper to rAngleEnd step Stepper
       nextX = xCenter + arcRadius * cos(rAngle)
       nextY = yCenter + arcRadius * sin(rAngle)
       #gr "goto ";int(nextX);" ";int(nextY)    'int speeds things up
   next
end sub

These arcs should start same place the Trig functions would graph them and move clockwise as the angle is increased about the center point, same way as the Trig function graphings go!

[tenochtitlanuk]

Just a nice graphic to confirm I CAN draw arcs reliably! Shades of Sauron's eye, with hgv colouring. This example was created with the arc-ends and radial centre marked. Wihout them and with more arcs drawn you get a smoother but less dramatic pic!

This thread also led me to examine the dot.product method. Thanks, friends..

[tsh73]

Ok I'm still at it

You don't have to use atan2 to calculate the angle between two vectors. If you just want the quickest way, you can use dot (v1, v2)=|v1|*|v2|*cos A to get. A = Math.acos ( dot (v1, v2)/ (v1.length ()*v2.length ()) );

To calculate angle between vectors A, B with ATAN2:

get A angle (-pi..pi) calling ATAN2
get B angle (-pi..pi) calling ATAN2 second time
subtract.
Now, it could end up as (-2pi..2pi)
   so we should normalize say to 0..2pi (or -pi..pi, whatever)

Now, to get angle by this formula:

calculate vectDotProduct(a$, b$)/ (vectLen(a$)*vectLen(b$))
  (that's three user defined functions for me)
then take ACS of it
!! Alas! ACS returns only angles in 0..pi range
  Will it be enough? I dunno, might depend on task at hand.


So I used old mine Vector Lib 2d for vector functions
run angle A from 0 to 2Pi step pi/6 (variable i, 1st column)
run angle B from 0 to 2Pi step pi/6 (variable j, 2nd column)
made random vectors with these angles
3rd column - just subtract (and normalise) A, B angles I have
4th - get angle through ATAN2 (happened to == with 3rd, good sign I think)
5th - get angle through ASC and vector functions.
(and I had to limit ACS argument because it happen to went out of [-1 .. 1] with some 1e-15 reminder, bringing ACS to error - so it does not run without these checks)

Now, all this is highly inefficient - vectors are stored as strings, I use val(word$())) for getting any numbers
Still ACS version works 30% slower. (counted by adding outer loop and removing all else)

'angle between vectors

'using Vector Lib 2d v.01
'vector 2d lib
'vectors as "x y" pairs, to be split by Word$
global pi
pi = acs(-1)   'global for atan2, dePi
''

for i = 0 to 11
   for j = 0 to 11
       a$=vectFromPolar$(1+rnd(0)*5, 2*pi*i/12)
       b$=vectFromPolar$(1+rnd(0)*5, 2*pi*j/12)
       da1=(2*pi*j/12)-(2*pi*i/12) 'b$-a$
       print using ("##",i);" ";using ("##",j);" ";
        print dePiNorm$(da1);" ";
       'now angle by atan2
       a1=vectAngle(a$)
       b1=vectAngle(b$)
       da2=b1-a1
       'print  dePi$(a1);" ";   'so it is from -pi to pi
       'print a$;":";b$;",";
       print  dePiNorm$(da2);" ";
       'print  dePi$(da2);" ";
       'now angle by acos
       'A = Math.acos ( dot (v1, v2)/ (v1.length ()*v2.length ()) );
       dd=vectDotProduct(a$, b$)/ (vectLen(a$)*vectLen(b$))
       'print dd, dd+1;" ";
       'somehow it get over -1
       if dd<-1 then dd=-1
       if dd>1 then dd=1
       da3=acs(dd)
       'print  dePiNorm$(da3)
       print  dePi$(da3)';" ";
       'print dd

   next
next

end
'---------------------------------------------

function vect$(x,y)
   vect$=x;" ";y
end function

function vectX(v$)
   vectX=val(word$(v$,1))
end function

function vectY(v$)
   vectY=val(word$(v$,2))
end function

function vectLen(v$)
   x=val(word$(v$,1))
   y=val(word$(v$,2))
   vectLen=sqr(x*x+y*y)
end function

function vectUnit$(v$)
   x=val(word$(v$,1))
   y=val(word$(v$,2))
   vectLen=sqr(x*x+y*y)
   vectUnit$=x/vectLen;" ";y/vectLen
end function

function vectAdd$(v1$,v2$)
   x1=val(word$(v1$,1))
   y1=val(word$(v1$,2))
   x2=val(word$(v2$,1))
   y2=val(word$(v2$,2))
   vectAdd$=x1+x2;" ";y1+y2
end function

function vectSub$(v1$,v2$)
   x1=val(word$(v1$,1))
   y1=val(word$(v1$,2))
   x2=val(word$(v2$,1))
   y2=val(word$(v2$,2))
   vectSub$=x1-x2;" ";y1-y2
end function

function vectDotProduct(v1$,v2$)
   x1=val(word$(v1$,1))
   y1=val(word$(v1$,2))
   x2=val(word$(v2$,1))
   y2=val(word$(v2$,2))
   vectDotProduct=x1*x2+y1*y2
end function

function vectScale$(a,v$)   'a * vector v$
   x=val(word$(v$,1))
   y=val(word$(v$,2))
   vectScale$=a*x;" ";a*y
end function

function vectTangent$(v$,base$)
   n$=vectUnit$(base$)
   vectTangent$=vectScale$(vectDotProduct(n$,v$),n$)
end function

function vectNorm$(v$,base$)
   vectNorm$=vectSub$(v$,vectTangent$(v$,base$))
end function

function vectAngle(v$)
   x=val(word$(v$,1))
   y=val(word$(v$,2))
   vectAngle=atan2(y,x)
end function

function vectFromPolar$(rho, phi)
   vectFromPolar$=rho*cos(phi);" ";rho*sin(phi)
end function

function vectRotate$(v$,alpha)
   x=val(word$(v$,1))
   y=val(word$(v$,2))
   rho=sqr(x*x+y*y)
   phi=atan2(y,x)+alpha
   vectRotate$=rho*cos(phi);" ";rho*sin(phi)
end function

function dePi$(x)   'pure aestetics
   'dePi$=x/pi;"Pi"
   dePi$=using("##.###",x/pi);"Pi" 'place for sign
end function

function normalizeAngle(x)
   x=x mod (2*pi)
   if x <0 then x=x+(2*pi)
   normalizeAngle = x
end function

function dePiNorm$(x)   'pure aestetics
   'norm z first
   x=x mod (2*pi)
   if x <0 then x=x+(2*pi)
  dePiNorm$=using("##.###",x/pi);"Pi"
end function


'---------------------------
function atan2(y,x)
   'pi = asn(1) * 2    'global
   if x <> 0 then arctan = atn(y/x)

   select case
       case x > 0
       atan2 = arctan

       case y>=0 and x<0
       atan2 = pi + arctan

       case y<0 and x<0
       atan2 = arctan - pi

       case y>0 and x=0
       atan2 = pi / 2

       case y<0 and x=0
       atan2 = pi / -2
   end select
end function