Thousands of intersections

[Rob_Durham]

Hi,

I am trying to get AutoCAD to draw a great many circles of differing radii with the centre point of the circles at a common foci point along the semi-major axis of an ellipse. The radius of all the circles are such that they will intersect the ellipse at some point.

The XYZ coordinates of these intersections are ultimately what I want to extract using AutoCADs data extraction utility.

I have successfully used Lee Macs intersections Lisp routine to draw ‘points’ at these intersections however it doesn’t like it when I have thousands of circles….it just crashes.

;; Intersections in Set  -  Lee Mac
;; Returns a list of all points of intersection between all objects in a supplied selection set.
;; sel - [sel] Selection Set

(defun LM:intersectionsinset ( sel / id1 id2 ob1 ob2 rtn )
    (repeat (setq id1 (sslength sel))
        (setq ob1 (vlax-ename->vla-object (ssname sel (setq id1 (1- id1)))))
        (repeat (setq id2 id1)
            (setq ob2 (vlax-ename->vla-object (ssname sel (setq id2 (1- id2))))
                  rtn (cons (LM:intersections ob1 ob2 acextendnone) rtn)
            )
        )
    )
    (apply 'append (reverse rtn))
)

(defun LM:intersections ( ob1 ob2 mod / lst rtn )
    (if (and (vlax-method-applicable-p ob1 'intersectwith)
             (vlax-method-applicable-p ob2 'intersectwith)
             (setq lst (vlax-invoke ob1 'intersectwith ob2 mod))
        )
        (repeat (/ (length lst) 3)
            (setq rtn (cons (list (car lst) (cadr lst) (caddr lst)) rtn)
                  lst (cdddr lst)
            )
        )
    )
    (reverse rtn)
)

(defun c:interset ( / sel )
    (if (setq sel (ssget))
        (foreach pnt (LM:intersectionsinset sel)
            (entmake (list '(0 . "POINT") (cons 10 pnt)))
        )
    )
    (princ)
)
(vl-load-com) (princ)

Since I don’t ultimately need the circles (only the intersection points they make with the ellipse) I thought I would try using the lisp below to delete the last circle drawn (actually the second to last entity after the ‘point’) which seems to work if I run it manually. This way, the first lisp routine only needs to find one intersection each time.

(defun countback(steps / ms)
(vl-load-com)
(setq ms (vla-get-modelspace (vla-get-activedocument
(vlax-get-acad-object))))
(vlax-vla-object->ename (vla-item ms (- (vla-get-count ms) steps)))
)
 (command ".erase" (countback 2) "")
(princ)

The trouble is I am using a copy paste method from excel to the command line to draw the circles. I thought I would try and also call the above routines from the spreadsheet but I can’t seem to get it to work, it will draw the points at the intersections but then won’t delete the last circle drawn.

Since I’m no expert in Lisp and I found both of these routines on the web, can someone possibly a) help me combine them into 1 lisp adding the erase command to delete the 2nd to last entity (countback 2) or b) suggest a better way of doing this whole process?

Thanks in advance

Rob

[Anonymous]

It should also be noted that command calls slow down programs quite a bit as well.

[Anonymous]

Try this one:

;; Intersections in Set - Lee Mac
;; Returns a list of all points of intersection between all objects in a supplied selection set.
;; sel - [sel] Selection Set

(defun LM:intersections ( ob1 ob2 mod / lst rtn )
(if (and (vlax-method-applicable-p ob1 'intersectwith)
(vlax-method-applicable-p ob2 'intersectwith)
(setq lst (vlax-invoke ob1 'intersectwith ob2 mod))
)
(repeat (/ (length lst) 3)
(setq rtn (cons (list (car lst) (cadr lst) (caddr lst)) rtn)
lst (cdddr lst)
)
)
)
(reverse rtn)
)
(defun c:CreateIntersections ( / sel )
(defun Circle (cen rad)
(entmakex (list (cons 0 "CIRCLE")
(cons 10 cen)
(cons 40 rad))))
;;Updated this portion to select all
(setq CirRad (getreal "\nEnter the radius of the circle: "))
(if (and CirRad
	 (setq sel (ssget "_X" (list(cons 410 (getvar 'ctab)) (cons 0 "LWPOLYLINE")))))
(progn
  (if (setq mycirss (ssget "_X" (list(cons 410 (getvar 'ctab)) (cons 0 "CIRCLE"))))
    (progn
      (setq myVLAcir (vlax-ename->vla-object (setq mycir (cadar (ssnamex mycirss)))))
      (vla-put-center myVLAcir (vlax-3d-point 0 0 0))
      (vla-put-radius myVLAcir CirRad))
    (progn
      (circle '(0 0 0) cirrad)
      (setq myCir (entlast))))
(foreach pnt (LM:intersections (vlax-ename->vla-object myCir) (vlax-ename->vla-object (cadar(ssnamex sel))) acextendnone)
(entmake (list '(0 . "POINT") (cons 10 pnt)))
)
)
)
  (gc) 
(princ)
)
(vl-load-com) (princ)

[Rob_Durham]

Well its running!

I'm now at home using my personal pc and although the process is slower (40 minutes to execute approximately 10,000 points) it seems more stable.

I can maximise and minimise at will and see it working at the command line and it doesn't say "not responding" so fingers crossed!

In case you are interested, these points are positions of the earth on its orbit around the sun every minute from its perihelion position on 4th Jan 2017.

To get a whole year (perihelion to perihelion) I need 394,691 intersections but I'm going to be doing it monthly batches. The point coordinates will be imported into 3d Studio Max via a maxscript so that I can create a scaled 3d model with animation.

Adding the moon will be a whole other set of problems due to its 5 degree inclination to the ecliptic. I hope you don't mind me coming back here again if I need further assistance?

[Anonymous]

I was curious as to what you were doing. I would start a new thread if you come back, but always glad to help.

[john.uhden]

That is truly exciting!

Do we get an honorable mention?  🙂

Is the orbit really a true ellipse?  Since you are using a polyline ellipse, I presume you can tweak it to the true orbit if necessary.

BTW, covering a year in just 40 minutes seems pretty fast to me.

It reminds me of one of my favorite limericks:

"There was a young lady named Bright

Whose speed was much faster than light.

She set out one day

In her usual way

And returned on the previous night."

[Kent1Cooper]

---

@Rob_Durham wrote:

.... these points are positions of the earth on its orbit around the sun every minute .... 

---

If that's what you're doing, here's a better way [he said, not quite understanding exactly how you're doing this with an Ellipse and Circles, that is, what the Circles represent, etc.].  No Circles needed, no intersections to find, not an approximate orbital path with a Polyline, etc.

Orbits are elliptical [well, to be really picky, the orbit of the Earth itself around the Sun wobbles from the true elliptical path because it's the orbit of the center of gravity of the Earth/Moon combination that is truly elliptical, but if you can ignore that....].  And the amount of motion of an orbit in a fixed time period is such that the areas enclosed in the "wedges" swept by the motion over each time increment are equal.

HERE'S THE MAGIC:  In VLA terms, parameter values for an Ellipse work just that way!  Locations at equally-spaced parameter values "sweep out" equal areas, and are closer together along the Ellipse in the vicinity of the major-axis endpoints [any orbiting object moves more slowly out at the extremes], which are the perihelion and aphelion locations, and farther apart along it in the vicinity of the minor-axis endpoints [it moves more quickly when closer in].

A full closed Ellipse's parameter values go from 0 to 2pi [all the way around].  To test this idea in simplified form, I drew an Ellipse, and using a similar approach to the code below, drew Lines from its center to 8 locations at parameter values that were at the multiples of 1/8 of 2pi.  I then used BOUNDARY to put Regions in all the wedges.  The Regions were all of equal area down as far as AutoCAD can display in Units settings [8 decimal places, i.e. to the nearest 100,000,000th of a drawing unit].  To check with a less "clean" situation, I did it again with 11 -- same degree of equality.

So to get the position of the Earth every minute, you need to divide 2pi by how many minutes there are in a year, and put a Point at each multiple of that parameter increment all the way around.  I used 365.25 days for a year, which I think is not quite precisely accurate, but you can fine-tune that if it matters.

Draw an Ellipse [a "real" one, not a Polyline approximation -- parameter values on a Polyline are entirely different] representing the Earth's orbit, then:

(setq ent (entlast)); the Ellipse entity name
(setq inc (/ (* 2 pi) (* 60 24 365.25)))
  ; the parameter-value increment [full orbit divided by (minutes/hour x hours/day x days/year)]
(repeat (setq n (fix (* 60 24 365.25))); same number of Points, converted to integer
  (entmake
    (list
      '(0 . "POINT")
      (cons 10 (vlax-curve-getPointAtParam ent (* (setq n (1- n)) inc)))
    ); list
  ); entmake
); repeat

Your results may vary, but my computer put all 525,960 Points around my sample Ellipse in about 25 seconds.  And since the 0 parameter is at a major-axis endpoint, one of the Points [the last one drawn, actually, because it counts down, though it could be adjusted to start with that one and count up instead] is exactly at the perihelion.  And yes, they're closer together where they should be, and farther apart where they should be [by how much depends on the eccentricity of the Ellipse].

[Rob_Durham]

Thanks for your response. I truly wish I could do it the way you have suggested but I see a couple of problems.

The first issue is that the sun is not at the centre of the ellipse. I would guess your end result is somewhat symmetrical whereas the earth travels faster at the perihelion than at the aphelion making the distance between points very different in January and July.

The second issue is that the perihelion length is different from one year to the next. My ellipse (or polyline) is drawn from the exact minute of the perihelion in 2017 to the exact minute of the aphelion in 2017 using NASA's given eccentricity value of 0.01671022. The perihelion position in January 2018 is a smaller distance to the sun. Therefore I have to draw 2 ellipses and join them together at the aphelion. This means that the ellipse is not a closed entity. It is open ended.

NASA have an ephemeris application on their website that give geocentric (centre-to centre) distances to all the major bodies in the solar system in km to 2 decimal places (I think its 2) for any specified time period of any year in 1 minute increments. If I draw a circle from the sun position at a radius equal to the distance NASA give from the earth to the sun then the point at which the circle intersects the ellipse is the position of the earth at that moment. It seems to be working well so far and is very very accurate but is painfully slow to execute.

I'd love for it to be completed in 25 seconds but I'm going for accuracy. If you can adapt your lisp to account for the above I'd love to know.

The reason for using a polyline was that an ellipse does not import into 3ds max however I'm not doing that now (only importing point data) so I can revert back to using an ellipse in ACAD if necessary.

[pbejse]

Very interesting discussion 

[Rob_Durham]

So the lisp seems to be working however it is taking hours and hours to complete.

It seems to get off to a really fast start, but as time progresses, it slows down to a very slow pace. If I monitor the memory usage in task manager it constantly increases.

I was wondering if there was a command I could add in the spreadsheet at regular intervals to pause/reset and free up some memory. Either that or I break it down manually into smaller batches??

[john.uhden]

Please review post 39.  You could even add in the (gc) after N iterations.

BTW, Kent's idea of using parameters is still valuable.  You must know the full length of the pellipse, so you could set points (circles, whatever) using (vlax-curve-getpointatdist.  Also, I'm pretty sure that points are cheaper to use than circles as points don't have a radius property.  Plus you can set whatever PDMODE and PDSIZE you want, even after the fact.

BTW, I just learned that you can set the size of a circle by area, without doing any math...

(vlax-put CircleObject 'Area 10)

Same for circumference.

[Rob_Durham]

@Kent1Cooper

I really like your idea and I previously posed questions on some maths forums (with no solutions thus far) as I wondered if I could use excel to calculate the XYZ coordinates using the known variables (Eccentricity, Perihelion Length, Aphelion Length and Elapsed Time) thus negating the need for ACAD altogether.

From what I understand, there are 3 angles that must be computed, referred to as the Mean Anomaly, Eccentric Anomaly and the True Anomaly. I believe there is no direct formula that can be derived for the Eccentric Anomaly and it is solved using an iterative method (whatever that means) I'm sure some whizz in excel could come up with something but I have no idea how to do it. 

---

@Kent1Cooper wrote:

---

Locations at equally-spaced parameter values "sweep out" equal areas, and are closer together along the Ellipse in the vicinity of the major-axis endpoints [any orbiting object moves more slowly out at the extremes], which are the perihelion and aphelion locations, and farther apart along it in the vicinity of the minor-axis endpoints [it moves more quickly when closer in]. 

---

Just to clarify,....your statement above would be true if the object being orbited ie. the sun, was at the centre of the ellipse. What they say is that because the sun is not at the centre (it is at a foci point along the major axis) the earth moves faster as it gets closer and gets slower the farther away it gets. This is essentially Kepler's Second Law of Planetary Motion. There is an animated gif on the wiki page that demonstrates it. https://en.wikipedia.org/wiki/Kepler%27s_laws_of_planetary_motion

With this in mind, is it possible to adapt your method and move the origin point off centre? Will it still yield "wedges" of equal area?

@john.uhden

The code provided by zraboin already includes a (gc)

Unless the method Kent has suggest can be modified I think I may be stuck with breaking the circles down into smaller batches to preserve memory.

Thank you so much to everyone who has helped so far and please remember I am not that familiar with lisp so I need to be spoon fed

[Kent1Cooper]

---

@Rob_Durham wrote:

.... 

The first issue is that the sun is not at the centre of the ellipse. I would guess your end result is somewhat symmetrical whereas the earth travels faster at the perihelion than at the aphelion ....

.... I'm going for accuracy. If you can adapt your lisp to account for the above I'd love to know.

The reason for using a polyline was that an ellipse does not import into 3ds max however I'm not doing that now (only importing point data) so I can revert back to using an ellipse in ACAD if necessary.

---

Well, duh on me for not thinking it through enough to realize that the Earth doesn't move faster at both ends and slower at the sides, but fastest at one end and slowest at the other.  I don't know of a way to map that with parameters.

But I have to say, the accuracy issue and the use of a Polyline [also from Post 32: "The ellipse was drawn with the pellipse system variable set to 1 so its actually a polyline...."] are mutually contradictory.  You don't have to Zoom in all that far to see where a true Ellipse and a PELLIPSE=1 Polyline approximation don't coincide -- quite close in a lot of areas, but off in others by what I imagine are unacceptable amounts.  Given that the Circles are going to hit whatever the orbital path object is at very shallow angles of incidence, a small degree of difference in path position makes a difference in intersection position several times the amount of the path variation.  In a pair of such Ellipses defined by exactly the same axis endpoints, with slight eccentricity a little wider than they are tall, and with a Circle centered slightly off-center from the Ellipses, you get this kind of situation in the vicinity of 24 degrees up from the right major-axis endpoint:

Dividing that orbit into the 525,960 minutes'-worth increments, that difference amounts to almost 12 hours off  between the true and the approximated Ellipses!

[Anonymous]

I can cut down the time by 1/3:

Different setup though. You type the command once and then enter all of the rads in order and then the command runs. When I tested it on 2000 points it took 52 seconds the first way and less than 18 with the new lisp.

(defun c:CreateIntersections ( / sel )
(defun Circle (cen rad)
(entmakex (list (cons 0 "CIRCLE")
(cons 10 cen)
(cons 40 rad))))
;;Updated this portion to select all
(while (setq rad1 (getreal "\Enter a new radius: "))
  (setq allrads (cons rad1 allrads)))
(if (and (> (length allrads) 0)
	 (setq sel (ssget "_X" (list(cons 410 (getvar 'ctab)) (cons 0 "LWPOLYLINE")))))
(progn
  (circle '(0 0 0) (car allrads))
  (setq mycir (entlast))
    
    (while (> (length allrads) 0)
      (setq myVLAcir (vlax-ename->vla-object mycir))
      (vla-put-radius myVLAcir CirRad)
      (entmake (list '(0 . "POINT") (cons 10 (vlax-invoke (vlax-ename->vla-object myCir) 'intersectwith (vlax-ename->vla-object (cadar(ssnamex sel))) acextendnone))))
      (setq allrads (cdr allrads))
      (gc)(gc))

)
)
 
(princ)
)
(vl-load-com) (princ)

[Anonymous]

Since my original command wasn't placing thousands of circles, it was unlikely that it was that bogging it down. It was the fact that we were running the command so many times. By cutting out the command call it should improve it drastically. 

[john.uhden]

If the constant creation and erasure of the circles is slowing things down,
I have a few functions for calculating intersection points of lines and
arcs without having to create any new entities.

[Anonymous]

The program isn't creating and erasing circles. (Read above) If anything, the creation of thousands of points is bogging it down over time.

[Rob_Durham]

@Kent1Cooper

It is my understanding (rightly or wrongly) that, geometrically, the ellipse and the Pellipse are the same however they are rendered on screen differently.

The Pellipse seems to be constructed on screen from lots of straight lines at slightly different angles. In practice, when the intersections are made and you zoom in, they do not fall exactly on the Pellipse, moreover they appear in the position where they would have been if it were a true ellipse.

I just had a look at the Autodesk help page and it doesn't really clarify either way.

Regardless, I'll change it to a true ellipse as we know for sure that this is correct.

I assume the code will need to change to reflect this?

[Rob_Durham]

Thanks again zraboin, but I'm not quite getting what you mean by "enter all the rads in order"?

Where / how do I enter these?

And I assume I can just substitute 'LWPOLYLINE' for 'ELLIPSE' if I want to use a true ellipse?

[Anonymous]

Yes, you can replace it with ellipse.

Try the command. Easiest way to figure it out.

Enter the command once, then copy and paste only the rad values into autocad. 

[Anonymous]

How have you been entering in the command and values in the past?