Vectors in the Plane


Dim A() As Double, B() As Double
Dim C() As Double, D() As Double
Dim R() As Double

A = pt(4, 1, 0)
B = pt(2, -3, 0)
C = pt(-2, 2, 0)
D = pt(-3, -2, 0)

vec_draw1 A, "A"
vec_draw1 B, "B"
vec_draw1 C, "C"
vec_draw1 D, "D"

R = v_add(A, B)
vec_draw1 R, "A+B"

R = v_add(A, C)
vec_draw1 R, "A+C"

R = v_add(A, D)
vec_draw1 R, "A+D"

R = v_add3(A, B, C)
vec_draw1 R, "A+B+C"

R = vec_add_2(1, A, 2, B)
vec_draw1 R, "A+2B"

R = vec_add_2(1, A, -3, C)
vec_draw1 R, "A-3C"

R = v_add3(A, C, v_m(1 / 3, D))
vec_draw1 R, "A+C+1/3D"


Sub vec_draw1(pt1() As Double, Optional str_text As Variant)
'line from origin to pt1
Dim lineobj As AcadLine
Set lineobj = acadDoc.ModelSpace.AddLine(pt(0, 0, 0), pt1)

    If Not IsMissing(str_text) Then
         txt1 CStr(str_text), pt1, 0.25
    End If
End Sub

'add 2 vectors, return vector
Function v_add(pt1() As Double, pt2() As Double) As Double()
v_add = pt(pt1(0) + pt2(0), pt1(1) + pt2(1), pt1(2) + pt2(2))
End Function

'multiply scalar by vector, return vector
Function v_m(m As Double, pt1() As Double) As Double()
v_m = pt(m * pt1(0), m * pt1(1), m * pt1(2))
End Function

'add 3 vectors, return vector
Function v_add3(pt1() As Double, pt2() As Double, pt3() As Double) As Double()
v_add3 = pt(pt1(0) + pt2(0) + pt3(0), pt1(1) + pt2(1) + pt3(1), pt1(2) + pt2(2) + pt3(2))
End Function

'add 2 vectors each multiplied by a scalar, A-B would be vec_add_2(1, A, -1, B)
Function vec_add_2(m As Double, pt1() As Double, n As Double, pt2() As Double) As Double()
Dim temp1() As Double,  temp2() As Double
temp1 = v_m(m, pt1)
temp2 = v_m(n, pt2)

vec_add_2 = v_add(temp1, temp2)
End Function

it would be possible to do this on one line but it would be hard to read.

Dim temp1() As Double, temp2() As Double

temp1 = vec_add_2(2, A, -3, B)
temp2 = vec_add_2(-3, C, 4, D)
R = v_add(temp1, temp2)
vec_draw1 R, "2A-3B-3C+4D"

Length and Direction –

Function vec_len(pt1() As Double) As Double
Dim x As Double, y As Double, z As Double
x = pt1(0): y = pt1(1): z = pt1(2)
vec_len = (x ^ 2 + y ^ 2 + z ^ 2) ^ (1 / 2)
End Function

Dividing a vector by the length (multiplying by the inverse of the length) produces a vector with length one in the same direction –

The angle of a 2D vector is a directed line measured from the positive X axis. It takes a value between 0 and 360 degrees not including 360 or 0 and 2 pi not including 2 pi. The inverse tangent is used with y/x as argument. Tangent returns an angle for an undirected line between -90 to +90 or -pi/2 to +pi/2. The code is a bit tedious but its pretty straightforward. We have to trap out zero values of x, to avoid divide by zero, then check for zeros of y and interpret, then retrieve the value for the inverse tangent and interpret according to which quadrant the head of the vector is in.

Function vec_ang(pt1() As Double) As Double
'returns angle in radians
'check for zero length vector return 0 for angle
Dim x As Double, y As Double

x = pt1(0)
y = pt1(1)

If x = 0 And y = 0 Then
    vec_ang = 0
    Exit Function
End If

'get axis directions
If x = 0 Or y = 0 Then
    If y = 0 And x > 0 Then vec_ang = 0
    If x = 0 And y > 0 Then vec_ang = Pi / 2
    If y = 0 And x < 0 Then vec_ang = Pi
    If x = 0 And y < 0 Then vec_ang = 3 * Pi / 2
Exit Function
End If

'calculate m tangent
Dim m As Double
m = y / x

If x > 0 And y > 0 Then  'First Q
'Debug.Print rad2deg(Atn(m))
vec_ang = Atn(m)
End If

If x < 0 And y > 0 Then  'Second Q
'Debug.Print rad2deg(Atn(m))
vec_ang = Pi + Atn(m)
End If

If x < 0 And y < 0 Then  'Third Q
'Debug.Print rad2deg(Atn(m))
vec_ang = Pi + Atn(m)
End If

If x > 0 And y < 0 Then  'Fourth Q
'Debug.Print rad2deg(Atn(m))
vec_ang = 2 * Pi + Atn(m)
End If

End Function


Sub vec_draw3(startpt1() As Double, vectorpt2() As Double, Optional str_text As Variant)
'line from startpoint at vector distance and angle
Dim lineobj As AcadLine
Dim pt3() As Double
pt3 = v_add(startpt1, vectorpt2)
Set lineobj = acadDoc.ModelSpace.AddLine(startpt1, pt3)

If Not IsMissing(str_text) Then
         txt1 CStr(str_text), pt1, 0.25
    End If
End Sub

triangles with random vertexes and medians – (I moved them after the program drew them)


Sub vec_draw2(pt1() As Double, pt2() As Double, Optional str_text As Variant)
'simple line from pt1 to pt2
Dim lineobj As AcadLine
Set lineobj = acadDoc.ModelSpace.AddLine(pt1, pt2)

If Not IsMissing(str_text) Then
         txt1 CStr(str_text), midpt2(pt1, pt2), 0.25
    End If
End Sub
Advertisements

3D Line – part 1

A straight line is completely determined by the coordinates of its endpoints. A straight line having a definite length and direction but no definite location in space is a vector. In its most basic form its a single point. A vector and a point are both a 3-array of doubles. Make them a dynamic array.

The functions to calculate direction, midpt, cosines, addition and scalar multiplication of vectors all return an array of 3 doubles. The line subroutines work with arrays of 3 doubles. Where I have one point argument, assuming 0,0,0 for the other, I number the routine_1. Where I have two point arguments, I number the routine_2. Passing and assigning points cleans up the code to almost logo like clarity.

dim pt1() as double, pt2() as double
pt1 = pt(2, 4, 6)
pt2 = pt(1, 9, 10)
line1 pt1
line1 pt2
line2 pt1, pt2
line2 midpt1(pt1), midpt1(pt2)

Function pt(x As Double, y As Double, z As Double) As Double()
Dim pnt(0 To 2) As Double
pnt(0) = x: pnt(1) = y: pnt(2) = z
pt = pnt
End Function

Sub line1(pt1() As Double)
'line from origin to pt1
Dim lineobj As AcadLine
Set lineobj = acadDoc.ModelSpace.AddLine(pt(0, 0, 0), pt1)
End Sub

Sub line2(pt1() As Double, pt2() As Double)
'line from pt1 to pt2
Dim lineobj As AcadLine
Set lineobj = acadDoc.ModelSpace.AddLine(pt1, pt2)
End Sub

Sub line3(startpt1() As Double, vectorpt2() As Double)
'line from startpoint at vector distance and angle
Dim lineobj As AcadLine
Dim pt3() As Double
pt3 = add_vectors(startpt1, vectorpt2)
Set lineobj = acadDoc.ModelSpace.AddLine(startpt1, pt3)
End Sub
Function dist1(pt1() As Double) As Double
Dim x As Double, y As Double, z As Double
x = pt1(0): y = pt1(1): z = pt1(2)
dist1 = (x ^ 2 + y ^ 2 + z ^ 2) ^ (1 / 2)
End Function

Function dist2(pt1() As Double, pt2() As Double) As Double
Dim x As Double, y As Double, z As Double
x = pt2(0) - pt1(0)
y = pt2(1) - pt1(1)
z = pt2(2) - pt1(2)
dist2 = (x ^ 2 + y ^ 2 + z ^ 2) ^ (1 / 2)
End Function

Function midpt1(pt1() As Double) As Double()
Dim x As Double, y As Double, z As Double
x = pt1(0): y = pt1(1): z = pt1(2)
midpt1 = pt(x / 2, y / 2, z / 2)
End Function

Function midpt2(pt1() As Double, pt2() As Double) As Double()
Dim x1 As Double, y1 As Double, z1 As Double
x = (pt1(0) + pt2(0)) / 2
y = (pt1(1) + pt2(1)) / 2
z = (pt1(2) + pt2(2)) / 2
midpt2 = pt(x, y, z)
End Function


Function dir_cosines1(pt1() As Double) As Double()
Dim cos_alpha As Double, cos_beta As Double, cos_gamma As Double
Dim d As Double

d = dist1(pt1)

If d = 0 Then
dir_cosines1 = pt(0, 0, 0)
Exit Function
End If

cos_alpha = pt1(0) / d
cos_beta = pt1(1) / d
cos_gamma = pt1(2) / d

dir_cosines1 = pt(cos_alpha, cos_beta, cos_gamma)
End Function


Function dir_cosines2(pt1() As Double, pt2() As Double) As Double()
Dim cos_alpha As Double, cos_beta As Double, cos_gamma As Double
Dim d As Double

d = dist1(pt1)

If d = 0 Then
dir_cosines2 = pt(0, 0, 0)
Exit Function
End If

cos_alpha = (pt2(0) - pt1(0)) / d
cos_beta = (pt2(1) - pt1(1)) / d
cos_gamma = (pt2(2) - pt1(2)) / d

dir_cosines2 = pt(cos_alpha, cos_beta, cos_gamma)
End Function

Function dir_angle(dir_cosine As Double) As Double
dir_angle = WorksheetFunction.Acos(dir_cosine)
dir_angle = rad2deg(dir_angle)
End Function

Function add_vectors(pt1() As Double, pt2() As Double) As Double()
add_vectors = pt(pt1(0) + pt2(0), pt1(1) + pt2(1), pt1(2) + pt2(2))
End Function

Function mult_vector(n As Double, pt1() As Double) As Double()
mult_vector = pt(n * pt1(0), n * pt1(1), n * pt1(2))
End Function

Sub test5()
Call Connect_Acad

pt1 = pt(2, 4, 6)
pt2 = pt(1, 9, 10)
line1 pt1
line1 pt2
line2 midpt1(pt1), midpt1(pt2)

Dim cosines() As Double
cosines = dir_cosines1(pt1)
pt3 = mult_vector(10, cosines)

'makes a line 10 in length along same line as pt1
line1 pt3

'alpha angle
Debug.Print dir_angle(cosines(0))
End Sub

Sub test6()
Call Connect_Acad

pt1 = pt(5, 1, 2)
pt2 = pt(1, 2, 3)
pt3 = add_vectors(pt1, pt2)

line1 pt1
line3 pt1, pt2
line1 pt3

acadApp.Update
End Sub

The Dot Product

Autocad does not have a 3D arc command and its kind of wild how hard it is to draw an arc on a sphere. I think we can do it, but there is a lot of theory involved, which is good because its fun and what else were you going to do? Googling the problem, the best solution I saw (which i won’t use here) was “ProjectGeometry” where a line drawn from point to point on the sphere is projected onto the solid sphere object and converted to an arc. I did not spend a lot of time looking for it, but I am not sure if that method is available in activeX. Its not an object type. Its a command, a written program in autocad, so its a manual technique. Although possibly you could use sendcommand or something similar.

To draw an arc on a sphere – two endpoints, center and radius will be the givens. Two vectors from 0,0,0 to the endpoints are drawn. The dot product, also called the scalar product of two vectors, has all the theory we need to construct a new ucs in the plane of the vectors and the arc we want to draw. Its a math topic in a beginning vectors class. The dot product is a calculation that projects a perpendicular from one vector to the other. The point that it touches becomes the origin of the new ucs. The new x-axis is defined by the endpoint of the vector, and the new y-axis is defined by the projecting vector endpoint. The 3 points must form a right angle, which they do. The dot product also gives us the included angle between the vectors, just from the endpoint coordinates. very nice.

There is no way to draw an arc that is not parallel to the current ucs, so you have to change the ucs. However all the other data that goes into creating an arc stays tied to the World ucs, as far as drawing in activeX is concerned. If you are drawing manually, autocad accomodates your desire to draw in the new ucs. 0,0,0 typed in at the command line is interpreted to be the new ucs origin. Activex generally ignores the ucs, except as it defines the drawing plane.

Draw an arc with any method then look at the arc in the properties window to see what data is read/write and what is grayed out. You can edit the center xyz coordinates, the radius, the start and end angle. Thats all. Those are the exact same prompts the activex addarc requires. Moving the center moves it parallel to its plane of creation. If the other data are changed, the shape of the curve is changed and new endpoints are calculated by autocad, but it stays parallel to the plane it was made in.

Startangle is a hard problem. With the dot product, we can calculate the included angle from the endpoints. We know the center and the radius. The startangle for the arc is not the angle from the current ucs. It is not the angle projected onto the world ucs. With 2D vectors, there is only one reference angle, the angle from the vector to the positive X axis. In the 2D plane autocad measures angles from the positive x-axis. The angle of a 3D vector is defined by the angle from the vector directly to each positive axis. The angle to the X-axis is called the Alpha angle. In the tilted ucs I did not see what the logic of startangle was. I drew an arc with startangle zero, and the arc touched down on the world x-axis. The startangle is the vector direction angle Alpha from the new x-axis directly to the world x-axis. But it is not that simple either. My vectors were 5,5,5 and -5,5,5. They were symmetrical across the z-axis, which made an alpha projection hit the x-axis. If the ucs has another twist, the alpha angle will extend to where it breaks the xy plane apparently.

I need to find the intersection of the ucs with the world xy plane. The only way i know how to do it right now is to draw an arc with startangle 0 and see where it starts. Pretty sure its a solvable problem.

(actually have ny and la backwards labeled. my latitude longtitude converter is not handling east west correctly, the arclength came out correct)

Sub A_scalar_dot_B(x1 As Double, y1 As Double, z1 As Double, x2 As Double, y2 As Double, z2 As Double)

Dim ptA() As Double
Dim ptB() As Double

Dim A_dot_B As Double
Dim len_A As Double
Dim len_B As Double
Dim cos_theta As Double
Dim theta As Double
Dim theta_deg As Double

Dim A_unit_vector() As Double
Dim A_unit_x As Double
Dim A_unit_y As Double
Dim A_unit_z As Double

Dim B_unit_vector() As Double
Dim B_unit_x As Double
Dim B_unit_y As Double
Dim B_unit_z As Double

Dim scalar_B_on_A As Double
Dim scalar_A_on_B As Double

Dim pt_scalar_B_on_A() As Double

ptA = pt(x1, y1, z1)
ptB = pt(x2, y2, z2)

A_dot_B = x1 * x2 + y1 * y2 + z1 * z2
len_A = (x1 ^ 2 + y1 ^ 2 + z1 ^ 2) ^ (1 / 2)
len_B = (x2 ^ 2 + y2 ^ 2 + z2 ^ 2) ^ (1 / 2)
cos_theta = A_dot_B / (len_A * len_B)
theta = WorksheetFunction.Acos(cos_theta)
theta_deg = rad2deg(theta)

A_unit_x = x1 / len_A
A_unit_y = y1 / len_A
A_unit_z = z1 / len_A

B_unit_x = x2 / len_B
B_unit_y = y2 / len_B
B_unit_z = z2 / len_B

scalar_B_on_A = len_B * cos_theta
scalar_A_on_B = len_A * cos_theta

pt_scalar_B_on_A = pt(A_unit_x * scalar_B_on_A, A_unit_y * scalar_B_on_A, A_unit_z * scalar_B_on_A)

ucs pt_scalar_B_on_A, ptA, ptB, "wow"

Debug.Print theta_deg

public_theta_deg = theta_deg

End Sub

Sub testconv()
Call Connect_Acad
  pt0 = pt(0, 0, 0)
 
  set_ucs 0, 0, 0, 1, 0, 0, 0, 1, 0, "world"
 
 pt1 = conv(40.713, 74.006, 3963)  'NEW YORK
 pt2 = conv(34.052, 118.244, 3963)  'LOS ANGELES
  'CONV IS PRELIMINARY
  
 line1 pt0, pt1
 line1 pt0, pt2
 
  Call A_scalar_dot_B(pt1(0), pt1(1), pt1(2), pt2(0), pt2(1), pt2(2))
  
  Dim r As Double
  r = 3963  'RADIUS IN MILES
     
     Dim alpha As Double
     alpha = 85.727
    ' alpha = 0
    'try 0, measure and input value
     
  arc1 pt0, r, alpha, alpha + public_theta_deg

Debug.Print pt1(0)
Debug.Print pt1(1)
Debug.Print pt1(2)

Debug.Print pt2(0)
Debug.Print pt2(1)
Debug.Print pt2(2)

acadApp.Update

End Sub

Spherical Coordinates

My goal is to draw on the sphere converting latitude longitude coordinates to xyz values. Once coordinates are calculated, every line on a sphere is an arc. To draw an arc with known endpoints, known radius and known center, you would think would be simple, but consider what if the points are inaccurate and an arc cannot be drawn. That is why the AddArc method requires a Center, a Radius, a Start Angle and an End Angle. It calculates the points. But notice, there is no opportunity to tilt the arc up into 3D space. It will lie in the XY plane. To draw an arc on a sphere from a center, a radius and two known points, a user coordinate system must be created containing the 3 points, and the angle between points has to be calculated. A UCS point set requires three points that form a right angle. the origin, one of the endpoints, and a normal point. There might be another problem – methods from VBA (ActiveX) generally ignore user coordinate systems and work only in World coordinate system, Autocad’s base of operation. However the help for addarc mentions the center is in WCS, but not the angles. I assume the angles are calculated from UCS. (update a few days later – they are not). Then the angle between points would be calculated from the X axis of the user coordinate system. I have not tried it yet.

Generally if an ActiveX method requires a point, its a WCS point. AddLine is a fully functional 3D command. It takes two 3D points as input. It does not matter what UCS is current. It ignores it. If VBA has the point 0,0,0 as the first point of a line, it starts the line at WCS 0,0,0 no matter where the origin is for the current UCS. That is not so at the command line, 0,0,0 typed at the command line to the line prompt is interpreted to be the origin of the current UCS.

Circle is a little different. you can draw circles at different Z elevations without changing the UCS from the command line, because the center can be any 3D point. You can do the same from VBA, but subsequent prompts and input do not allow you to tip the circle so that it is not parallel to the current UCS.

we need to survey our basic tools and see how they act and what they require.

we need some ucs tools and some viewpoint tools. viewpoint is somewhat more involved than ucs. it has more variables and finally whatever view is created has to be applied to the current viewport. for now i am going to rely on some simple macros and my 3Dconnexion mouse, which works great, the simple one.

The old viewpoint command, before the dialog box, just took a direction vector. such as “1,-1,1” for a southeast viepoint. We can throw that into a program quick if we need it. mostly i will probably rely on the mouse. but here is the macro sub using sendcommand. if you type VP at the command line and get the dialog, type -VP. if that doesnt work look at your pgp which is re-directing. this code works if -vp works at your command line.

Sub vp(str_vp As String)
Dim str As String
'these are commented out for reference 
'acadDoc.SendCommand "-vp" & vbCr & "-1,-1,1" & vbCr  'SW
'acadDoc.SendCommand "-vp" & vbCr & "1,-1,1" & vbCr  'SE
'acadDoc.SendCommand "-vp" & vbCr & "1,1,1" & vbCr  'NE
'acadDoc.SendCommand "-vp" & vbCr & "-1,1,1" & vbCr  'NW
'acadDoc.SendCommand "-vp" & vbCr & "0,0,1" & vbCr  'TOP
'acadDoc.SendCommand "-vp" & vbCr & "0,0,-1" & vbCr  'BOTTOM
'acadDoc.SendCommand "-vp" & vbCr & "-1,0,0" & vbCr  'LEFT
'acadDoc.SendCommand "-vp" & vbCr & "1,0,0" & vbCr  'RIGHT
'acadDoc.SendCommand "-vp" & vbCr & "0,-1,0" & vbCr  'FRONT
'acadDoc.SendCommand "-vp" & vbCr & "0,1,0" & vbCr  'BACK

 Select Case UCase(str_vp)
    Case "TOP"
    str = "0,0,1"
    
    Case "FRONT"
    str = "0,-1,0"
    
    Case "RIGHT"
    str = "1,0,0"
    
    Case "LEFT"
    str = "-1,0,0"
    
    Case "BACK"
    str = "0,1,0"
    
    Case "TOP_FRONT"
    str = "0,-1,1"
     
    Case "SW"
    str = "-1,-1,1"
    
    Case "SE"
    str = "1,-1,1"
    
    Case "NE"
    str = "1,1,1"
     
    Case "NW"
    str = "-1,1,1"
     
    Case "BOTTOM"
    str = "0,0,-1"
      
    Case Else
    Exit Sub
  End Select
 
acadDoc.SendCommand "-vp" & vbCr & str & vbCr
End Sub

next we need to survey our basic tools, line, circle, arc, polyline, and ucs with regard to how they work in 3D space. we will need some text labels, we will want to dimension in 3D. we can download the formulas to convert from latitude longitude to cartesian but we want to show how they are developed. The diameter of the earth is not constant. People who map for a living know exactly what the bulge is. we are going to assume a perfect sphere, but we will use a variable for the Radius so results will be general for a bowling ball or the earth.

once we have the autocad tools, we should be able to download some coordinates into excel and draw arcs between them from vba. assuming i can draw 3D arcs.

Automating Autocad

Years ago an electrical engineer stopped by my desk, leaned over the partition and said, I want you to automate my ladder diagrams. I want to push a button and the drawing comes out. He had seen what we had done, what others had done, with an automated program that spit out a family of parts that were going down our only area of the plant with a production line. I answered him that I could not write that program, which he interpreted I would not write the program. You don’t know how much time that button would save me, he said. In my only perfect response in this lifetime I said, you don’t know how long it would take me to write that button.

Our successful automation program had all the necessary ingredients. The product line had existed a long time. It had sold a lot. The variations in the family of parts were considerable but considering the volume and the number of times they had been reproduced, a system had developed to name and describe them. The product line had several smart people good with computers working in that area for a long time. High volume, good people doing it manual for a long time, discrete number of variations. That is what you must have to automate – a manual system.

The first step to automation, besides an actual product, is a five minute refresher, or introduction, to databases, specifically table relationships, and specifically what is a key field, and what is a foreign key. This is how data is stored. This is how the product is described. MS Access used to be a useful tool, but version changes, file version problems, 32-64 bit problems, getting it funded, being replaced by the big do-all database, all these caused it to drop in my work environment. But the principles of product structure, a main table and parts tables connected by relations, are applied in excel.

I will come back to this. This is a large topic not really suited to the blog format. I need to think up a sample project. It will be math related, not with obvious commercial applications.

Polygon Circle

Regular polygons are called convex and star polygons that cross themselves are non-convex. For both types the vertexes (vertices) are on a circle that has a center and radius.

The turning angle the turtle uses from line to line is more generally called the exterior angle of a polygon. Its the angle from the line extended to the next line.

The sum of the exterior angles of a regular polygon = 360.

nA=360

To find the radius of a regular polygon, divide it into triangles. each wedge has a central angle of 360/n.

Bi-Secting one of these triangles creates a right triangle with a known angle and a known side length.

This formula works for regular polygons and for star polygons. If a line is drawn from the center to each vertex, and a line drawn perpendicular to the edge, right triangles with known sides and angles are created.

now that we have the radius, the center point can be found. Think of the turtle on a vertex, heading turned to the next line. It is on a triangle we already solved to find the radius. we use the same angles. The angle that the turtle needs to turn to look at the center is 90-A/2. I use acad utility polarpoint, not the turtle, to draw the circle. I wrote a new turtle class function to return the current turtle position as a point using a dynamic array and used it directly in PolarPoint.

Sub poly_1(angle As Double, n As Integer, len_side As Double)
Dim inc As Integer
Dim rad As Double
Dim A As Double
Dim ctr() As Double
Dim ang2ctr As Double
Dim acadcirc As acadcircle

A = ang2rad(angle)
rad = len_side / (2 * (Sin(A / 2)))

ang2ctr = ang2rad(turtle1.heading - (angle / 2) + 90)
ctr = acadDoc.Utility.PolarPoint(turtle1.pt1, ang2ctr, rad)

For inc = 1 To n
turtle1.fd len_side
turtle1.left angle
Next inc

Set acadcirc = acadDoc.ModelSpace.AddCircle(ctr, rad)

txt_h "A = " & angle, turtle1.x1, turtle1.y1, 0.125
txt_h "n = " & n, turtle1.x1, turtle1.y1 - 0.25, 0.125
txt_h "R = " & (angle * n / 360), turtle1.x1, turtle1.y1 - 0.5, 0.125
txt_h "Radius = " & Round(rad, 4), turtle1.x1, turtle1.y1 - 0.75, 0.125

End Sub

'this is part of the turtle class module to return current position as a point array
Public Function pt1() As Double()
Dim pnt(0 To 2) As Double
pnt(0) = Me.x1: pnt(1) = Me.y1: pnt(2) = 0
pt1 = pnt
End Function

I am not going to try to draw these triangles, but all the graphics check out, it seems like my formulas derived from simpler polygons gives me the correct center and radius.


Sub turtle_demo_18()
 init_turtle
 Dim n As Integer, R As Integer
 Dim A As Long

 A = 164
 n = LCM(A, 360) / A
 R = LCM(A, 360) / 360
 
 Call poly_1(CDbl(A), n, 1)
End Sub

polygons

The inputs for this example problem are the angle 174 and 360. When a multiple of 174’s add up to a multiple of 360, the program stops. In advance we do not know how many 174s will add up to how many 360s, but we want the first one. This is the Least Common Multiple of Angle A (174) and 360.
LCM (A, 360)

174 has prime factors of 2 * 3 * 29
360 has prime factors of 2 * 2 * 2 * 3 * 3 * 5
the least common multiple is 360 * 29 = 10440.

we can substitute that into either side of the equation
nA=360R
n*174 = 10440
10440 = 360*R
to find that n = 60 and R = 29

The LCM can be figured different ways. There is a famous method of Euclid, and Euler, (no Eugene). There is an excel function to do it. Sample VBA code can be downloaded. The simplest method conceptually is the same way the turtle does it, by adding angles one by one and testing the result.

Here is brute force way one with no advance knowledge of when the program will stop. The program stops when the turtle heading returns to its original position. I also have a second emergency stop at 361 which never comes into play usually. (*footnote – the examples in the group photo that look like folded ribbon were polygons that did not close until 360 lines were drawn but were forced into an early exit with ” If inc > 130 Then Exit Do”) (**footnote – if the turn angle A is an integer, the maximum number of lines to bring the heading back to start is 360.)

Sub poly_360(len_side As Double, angle As Double)

Dim inc As Integer
Dim heading As Double
heading = turtle1.heading

Do
turtle1.fd len_side
turtle1.left angle

inc = inc + 1
If inc > 361 Then Exit Do
Loop While turtle1.heading <> heading

End Sub

Here is a primitive LCM function based on the same method, not intended to be the final version. The angles are added one by one and the result divided by 360 looking for a remainder, breaking out when the remainder is zero. The VBA mod operator works accurately only with integers. I had some overflows on the multiplication. The function seems to work better when all is type Long.

Function LCM(A As Long, B As Long)
 Dim n As Long
 Dim result As Long
 Dim remainder As Long
 n = 0
 
 Do
    n = n + 1
    result = n * A
    remainder = result Mod B
 Loop While remainder <> 0

 LCM = result
End Function

Now the sub to draw the polygon can be taken back to its roots. The loop calculations can be removed, because we will know in advance how many lines will be drawn.

Text information labels are added after the drawing is complete.

Sub poly_1(angle As Double, n As Integer, len_side As Double)
Dim inc As Integer

For inc = 1 To n
turtle1.fd len_side
turtle1.left angle
Next inc

txt_h "A = " & angle, turtle1.x1, turtle1.y1, 0.125
txt_h "n = " & n, turtle1.x1, turtle1.y1 - 0.25, 0.125
txt_h "R = " & (angle * n / 360), turtle1.x1, turtle1.y1 - 0.5, 0.125
End Sub

The sub to call the poly can be fancy or plain. It can draw families of polygons. It can loop and draw a range of turning angles.
This particular one will draw all the polygons with total turns (R) = 29 of angles between 1 and 180.

Sub turtle_demo_16()
init_turtle

Dim inc As Integer
Dim n As Integer, R As Integer
Dim A As Long, B As Long

B = 360
 
For inc = 1 To 180

 A = inc
 n = LCM(A, B) / A
 R = LCM(A, 360) / 360

If R = 29 Then

Debug.Print "LCM of " & A; " and " & B; " = " & LCM(A, B)
Debug.Print "A = " & A
Debug.Print "n = " & n
Debug.Print "R = " & R
Debug.Print " "

Call poly_1(CDbl(A), n, 1)

turtle1.x1 = turtle1.x1 + 5

End If
Next inc

End Sub

Sub turtle_demo_16()
init_turtle

Dim inc As Integer
Dim n As Integer, R As Integer
Dim A As Long, B As Long
B = 360
 
For inc = 160 To 181
 A = inc
 n = LCM(A, B) / A
 R = LCM(A, 360) / 360
Call poly_1(CDbl(A), n, 1)
turtle1.x1 = turtle1.x1 + 1.1
Next inc

End Sub

R<20