I (reluctantly?) agree with Britton that there's an error. He's going to
fix it, unless somebody speaks now. (Else, holding peace forever is
expected.)
Forwarded conversation
Subject: yt: embarrassing question, gravity
------------------------
From: *david collins* <antpuncher(a)gmail.com>
Date: Sat, Jan 17, 2009 at 5:28 PM
To: Matthew Turk <matthewturk(a)gmail.com>
Hey, Matthew--
I have an embarrassingly dumb question that I'm not going to see until
I ask someone. And I swear it made sense to me once.
In your _IsBound function, you have
pot = 2*G*PointCombine = 2 * G * Sum_i{ Sum_j { m_i * m_j / || r_i - r_j ||
}}
Why the 2 in front? Why isn't there a 1/2 for double counting in the
sum? And isn't the virial thm 2*KE < PE? So Now I'm missing an extra
2, so I'm off by 8.
Thanks,
d.
----------
From: *Matthew Turk* <matthewturk(a)gmail.com>
Date: Sat, Jan 17, 2009 at 5:40 PM
To: Britton Smith <brittonsmith(a)gmail.com>
Begin forwarded message:
*From:* "david collins" <antpuncher(a)gmail.com>
*Date:* January 17, 2009 5:28:31 PM PST
*To:* "Matthew Turk" <matthewturk(a)gmail.com>
*Subject:* *yt: embarrassing question, gravity*
----------
From: *Matthew Turk* <matthewturk(a)gmail.com>
Date: Sat, Jan 17, 2009 at 5:45 PM
To: david collins <antpuncher(a)gmail.com>
Cc: Britton Smith <brittonsmith(a)gmail.com>
Hey,
Britton wrote it so I fwd'd it to him. We don't double count, tho -
iteration is n < m inside PointCombine.c. Hopefully Britton can point out
more?
-Matt
----------
From: *david collins* <antpuncher(a)gmail.com>
Date: Sat, Jan 17, 2009 at 5:49 PM
To: Matthew Turk <matthewturk(a)gmail.com>
> iteration is n < m inside PointCombine.c.]
Todja it was dumb-- I fail both reading and programming on that one... XD
d.
----------
From: *Britton Smith* <brittonsmith(a)gmail.com>
Date: Sat, Jan 17, 2009 at 6:39 PM
To: Matthew Turk <matthewturk(a)gmail.com>
Cc: david collins <antpuncher(a)gmail.com>
After looking at this again, I think the first 2, as in 2 * G * Sum_i{ Sum_j
{ m_i * m_j / || r_i - r_j || }}, may be an error. As Matt said, the sum is
actually i = 0 ... N-1 and j = i+1...N, so there is no need for an
additional 0.5. When I had done this originally, I was somewhat uncertain
about calculating the boundness of N particles as to whether you should do
the full sum or this partial sum mentioned above. I had found the source to
some other code that calculates gravitational potential for particles and it
had this partial sum with a 2 out front, so I went with that. That, I'm
pretty sure now, it not right. The total binding energy is actually the sum
over all N * (N - 1) / 2 pairs of cells or particles, which is exactly the
sum above without that 2.
My impression of the virial theorem is that it applies to the time average
of KE and PE. Any body in a stable orbit will satisfy the virial theorem
when averaged over the entire orbit. A perfectly circular orbit has 2*KE =
PE all the time, but an elliptical orbit will vary in both directions
throughout its orbit. So, gravitational boundness for an object requires
that KE < PE at all times, but when you do a time average of that object you
will find 2 <KE> = <PE>.
In conclusion, there should be no problem with respect to the virial
theorem, but the first 2 in front of the sum should probably be removed. If
you guys agree with this, I will go ahead and remove that 2 from the
calculation. Let me know what you guys think.
Thanks for pointing this out, Dave.
Britton