6

6/4/2007

Polytropes Part I [working toward a polytrope sim step by step]

[links]

Astrophysics with a PC [contains useful and interesting treatment on polytropes, but not collisions. The simple polytrope model in this book is explained below]
Dr. Joshua E. Barnes' (University of Hawaii) n-particle polytrope sims [you can see colliding stars here. One of the most inspiring sites on the web, as far as I'm concerned].

Summary of forces

A polytrope is a time independent model of a star which computes dependent quantities density ρ(r), mass M(r), and pressure P(r) as functions of r, the distance from the center of the star. As per usual with infinitesimals, the cylinder below represents a small 3 dimensional volume within the star itself. The dimensions of the cylinder are arbitrary but very small, and in fact, since the only obtainable solutions for most interesting polytropes are purely numeric, they will be in fact small rather than infinitesimal as they would be for analytic solutions.

The cylinder is oriented and defined as follows.

It rides a vector passing through the center of the (theoretically) spherical star. The bottom (the part of the cylinder nearest the center) is at a distance r from the center.
The height of the cylinder is dr (an arbitrarily small value)
The surface area of the circular each cap is dA

The relationship between force and pressure: F = Pressure * Area or F = P(r) dA

Upward force on the cylinder (from the center outward, pushing the cylinder away from the center of the star) = P(r) dA

Downward force on the cylinder (pushing the cylinder toward the center of the star) = -P(r+dr) dA.

Downward force of Gravity on the cylinder (pulling the cylinder toward the center of the star) =

= -G * thin-shell mass carved inside the star * mass of cylinder divided by r²

The factors in the last equation are understood by the relations that

1) The force of gravity is calculated as

where G = gravitational constant, m₁ = mass1, m₂ = mass2, and r = distance between the masses

In the case of the polytrope, the mass acting on the cylinder pulling it inward is computed as a function of r = M(r) = m₁ in 1) and includes the mass enclosed by a shell just inside where the cylinder is (the shell enclosed by radius = r).

The mass of the cylinder itself is computed by using the relation mass = density * volume and thus: The volume of the cylinder is dr * dA and the density is a continues function ρ(r) of r; Therefore the mass of the cylinder = ρ(r) * dr * dA = m₂ in 1)

The cylinder is essentially stationary, so that the sum of forces acting on it = 0.

r >> dr, meaning in general, r is much larger than dr.

P', G and P are all acting on the cylinder. They have been pulled out of the cylinder because it's too confusing with them overlaying the drawing.

The final equation is P + P' + G = 0 where P' and G are negative (opposite of P in any case), so that if you're computing magnitudes (positive values) of the forces then the equation is

P-P'-G = 0

After some manipulation with the above, we get

and finally, Euler's Equation for hydrostatic equilibrium (Hellings p. 113)

From here, the book takes an interesting mathematical turn. At this point, we have three functions of r ρ(r), P(r), M(r) but only two equations relating them, and therefore some additional presumptive relations are required.

Kurt binghamgham

--KB 6/5/2007

_______________________________________________________________________

6/7/2007 Polytropes Part II

We have the following relationships

1) Mass differential

2) Pressure differential (Euler's Equation)

And the Polytrope relation, an additional equation in order to make the system solvable, and basically it suggests a relationship between the density and the pressure based upon (I guess) whatever the star is made of, among other factors, but it's specific to the star and is constant throughout the star (Hellings, p.113). This mystery constant is K.

P(r) = K ρ(r) ^(n+1)/n

3) Polytrope Relation, basically what's above but without the mystery variable K and instead, a ratio of density and pressure

Notice that n suddenly appears: n is a constant between 0 and 5 and is discussed below, requiring we look at the boundary conditions of the functions. The equations are only sensible on certain [closed] intervals

These lead to the following relationships and resultant substitutions

4a)

4b)

4c)

4d)

In 4d), when differentiated wrt r we get

With 2) we get:

Using 2), 4c) and 4d) and more substitutions (differentiating 4d), we obtain

Now, differentiating both sides (will show all steps, although Hellings does not), using the chain rule and substitutions above, we get

Create a variant of r to be computed at the end (when the equations are solved)

Create a new variable x to be referenced in the equations (the actual value of which will be worked out later, and will give us a real value r)

And finally, the Lane-Emden Equation

Next will be analysis of the boundary conditions and generating a numerical procedure. He actually has the quick basic pasted in the book. Naturally, I will use my own functions, but analyzing his solution will be an interesting exercise. In my case, I will be using a class Polytrope which inherits my standard 'body' class (which computes its own point set and handles motion etc), and use this to generate a single polytrope after which, I will create two and have them interact with one another.

--KB 6/7/2007

_______________________________________________________________________

Polytropes continued...

6/18/2007

Taking the original Lane-Emden equation and solving it is fairly straightforward, but as with all realistic problems, it offers a couple of problems, like a singular point at x = 0, the center of our star, which is circumvented via power series expansions of F and F'.

I used Euler for the first class rather than the midpoint method as he used in the book. I want to do everything from the ground up since it's more fun that way. Even though the Euler method isn't as accurate, it's still pretty close to his results.

Once again, the Lane-Emden Equation

Converted into all too familiar DE notation:

Use these Taylor expansions to get the initial values for y and y' of epsilon. Our seed value (e) is very small will produce our initial conditions for y and y'.

(Taylor series to get your seed values)

A and B usually denote matrices. Here, they are simply scalar known quantities computed from the above expansions. We know A and B when we start our recursion.

In order to setup an iteration, we need to manipulate the equations a bit, and figure out what our i+1 terms look like as well as our first 1 terms. This is done below.

What we're solving for is the point at which F or y becomes zero, where the pressure/density crosses the x-axis (thereby becoming zero), in which case you've reached the surface of the star.

[A very brief overview of solving a pair of 1ODEs via Euler and Heun is given here. This online book is concise and easy to follow. In practice, more computationally intensive solvers will be needed, and the online document gives and example of Heun at the end of that chapter.]

First, the original equation in familiar y, y' notation:

Setting up two FODEs based on the above we get

And given our initial conditions, we get the following computed values for a region very near the center of the star (e is very small)

Now substituting that back into the above equations, we get the following iterations in i, each followed by the first iteration.

The program below is the result of the above steps (which is by no means simple or obvious). It took some fiddling around to assemble these steps on paper, and then programming them was straightforward. The online book does go over this very well as I said above.

I created a basic 2nd order differential equation solver (Second_O_DiffyQ_Eqn) and then incorporated that into a polytrope class. This solves by either Euler or Heun methods. It also has some sample code that can be uncommented to test the solver itself (this way you can add different Runge-Kutta methods and test them). The function is the same as that used in the online DiffyQ book mentioned above.

Polytrope Differential Eqn solver

Public Class hC

Public Const nf = ".00000000"

End Class

Public Structure xyzS

Public x As Double

Public y As Double

Public z As Double

Public Sub New(ByVal ix As Double, ByVal iy As Double, ByVal iz As Double)

x = ix

y = iy

z = iz

End Sub

End Structure

Public Class Second_O_DQ_Solver

'solves 2nd order diffy q using euler by default but can be overriden for F

'to change the function, and compute to change the method (to heun or other runge-kutta method)

Public XnYn() As xyzS

Public Sub Print()

For k As Integer = 0 To UBound(XnYn)

Console.WriteLine(k & " : x=" & Format(XnYn(k).x, hC.nf) & ",y=" & Format(XnYn(k).y, hC.nf) & ",z=" & Format(XnYn(k).z, hC.nf) & ": (nevermind this)" & Format(-XnYn(k).x ^ 2 * XnYn(k).z, hC.nf))

End Sub

Public Sub PrintLast()

Dim k As Integer = UBound(XnYn)

Console.WriteLine(k & " : " & Format(XnYn(k).x, hC.nf) & ", " & Format(XnYn(k).y, hC.nf) & ", " & Format(XnYn(k).z, hC.nf) & ": " & Format(-XnYn(k).x ^ 2 * XnYn(k).z, hC.nf))

End Sub

Public ReadOnly Property F()

Get

Return XnYn(UBound(XnYn)).z

End Get

End Property

Public ReadOnly Property FX()

Get

Return XnYn(UBound(XnYn)).x

End Get

End Property

Private Function f2(ByVal xi As Double, ByVal yi As Double, ByVal zi As Double, ByVal n As Double, ByVal h As Double) As Double

Dim zk As Double

zk = (-yi ^ n - 2 * zi / xi)

'zk = (Math.E ^ (-xi) - 2 * zi - yi)

Return zk

End Function

Protected Overridable Function Compute_Euler(ByVal p0 As xyzS, ByVal h As Double, ByVal n As Double) As xyzS

Dim p1 As xyzS = New xyzS

'yi+1 = yi + h*zi

'yi+1 = yi + h*f1

p1.y = p0.y + h * p0.z

'zi+1 = zi + h*f2(xi,yi,zi)

p1.z = p0.z + h * f2(p0.x, p0.y, p0.z, n, h)

p1.x = p0.x + h

Return p1

End Function

Protected Overridable Function Compute_Heun(ByVal p0 As xyzS, ByVal h As Double, ByVal n As Double) As xyzS

Dim p1 As xyzS

Dim k1Y As Double = p0.z

Dim k1Z As Double = f2(p0.x, p0.y, p0.z, n, h)

Dim k2Y As Double = p0.z + h * k1Z

Dim k2z As Double = f2(p0.x + h, p0.y + h * k1Y, p0.z + h * k1Z, n, h)

p1.y = p0.y + (k1Y + k2Y) * h / 2

p1.x = p0.x + h

p1.z = p0.z + (k1Z + k2z) * h / 2

Return p1

End Function

Public Sub New(ByVal x0 As Double, ByVal y0 As Double, ByVal z0 As Double)

ReDim XnYn(0)

XnYn(0).x = x0

XnYn(0).y = y0

XnYn(0).z = z0

End Sub

Public Sub Compute(ByVal h As Double, ByVal max_itns As Integer, ByVal n As Double, ByVal method As String)

For k As Integer = 0 To max_itns - 1

ReDim Preserve XnYn(UBound(XnYn) + 1)

Select Case method

Case "euler"

XnYn(UBound(XnYn)) = Compute_Euler(XnYn(k), h, n)

Case ("heun")

XnYn(UBound(XnYn)) = Compute_Heun(XnYn(k), h, n)

End Select

If XnYn(UBound(XnYn)).y <= 0 Then

ReDim Preserve XnYn(UBound(XnYn) - 1)

Exit Sub

End If

End Sub

End Class

Public Class Poly_funcs

Public Shared Function Ye0(ByVal e As Double, ByVal n As Double) As Double

Return 1 - (1 / 6) * e ^ 2 + (n / 120) * e ^ 4

End Function

Public Shared Function YeP0(ByVal e As Double, ByVal n As Double) As Double

Return -(1 / 3) * e + (n / 30) * e ^ 3

End Function

End Class

Public Class PolyTrope

Private n As Double

Private M As Double

Private R As Double

Private Equn As Second_O_DQ_Solver

Private Pc As Double

Private RoM As Double

Private RoC As Double

Private L As Double

Private Rn As Double

Private Xn As Double

Private Hn As Double

Private Const Sp = 904800000000000.0

Public Sub New(ByVal ni As Double, ByVal Mi As Double, ByVal Ri As Double)

n = ni

M = Mi

R = Ri

Dim x0 As Double = 0.0001

Equn = New Second_O_DQ_Solver(x0, Poly_funcs.Ye0(x0, n), Poly_funcs.YeP0(x0, n))

End Sub

Public Sub Compute(ByVal h As Double, ByVal max_its As Integer, ByVal method As String)

Equn.Compute(h, max_its, n, method)

Xn = Equn.FX

Hn = Equn.F

Pc = Sp * M ^ 2 / (n + 1) / Hn ^ 2 / R ^ 4

RoM = 1.42 * M / R ^ 3

RoC = -RoM * Xn / 3 / Hn

L = -Xn ^ 2 * Hn / M

Rn = R / Xn

End Sub

Public Sub Print()

Console.WriteLine("Pc=" & Format(Pc, hC.nf))

Console.WriteLine("RoM=" & Format(RoM, hC.nf))

Console.WriteLine("RoC=" & Format(RoC, hC.nf))

Console.WriteLine("L = " & Format(L, hC.nf))

Console.WriteLine("Rn =" & Format(Rn, hC.nf))

Console.WriteLine("-----------")

End Sub

End Class

Module Main_Module

Sub Main()

'this is a test equation you can use for both euler and heun based on the

'linked document, for their equation; in the f2 function, simply comment

'out the polytrope equn and use the other one, then you can create

'new runge-kutta solvers and test the results

'Dim p As Second_O_DQ_Solver

'p = New Second_O_DQ_Solver(0, 1, 2)

'p.Compute(0.25, 4, 0, "heun")

'p.Print()

'************ end test section for DQ_Solver class

Dim p As New PolyTrope(1.5, 2, 3)

p.Compute(0.05, 300, "heun")

p.Print()

Console.ReadLine()

End Sub

End Module

kurt bingha

KB 6/19/2007

______________________________________________________________________

Verifying the Polytrope I program with a GUI

6/26/2007

In the past, I have raced toward the graphical portion of a program too quickly and ended up with my TI calculator in hand, sifting through console outputs to verify results because some object didn't look right or because various numbers seemed off. I spent a lot of time writing down or printing data and then converting it manually to a myriad of different units in order to compare it with known values. This was especially difficult in programs where I had, in order to simplify the frustum, modified constants like G to be in terms of earth AU^3/Earth_mass, and made standard distances and velocities in AU or even in some cases, arbitrary large numbers I chose for different reasons.

Expressing members of a solar system in Earth masses is succinct and simplifies things for the human, especially someone who wants to tinker with things like putting a larger version of Earth where Mercury is and setting the system in motion to see what happens. However, if there's a problem and the programmer can't verify the result, then he needs to go back and put in standard data (simulate the Earth/Moon/Sun system for example), and he will need to convert the known data from its standard units (km/m/kg) to his custom units (AU,earth_radii, etc).

So, before I start putting polytropes on the screen, I made some design commitments

Input data should be in a form easy for humans to create (like 3 Earth Masses and 2 Moon Radii for a new object)
Data coming from NASA or astronomy texts expressed in km/m and other units can be quickly entered directly without needing to do any conversion.
The programmer will want to see the data in easily understood formats, maybe even MPH or in miles or km/h, and that should be available.
Internally the program will use its own units, and the programmer will need a standard not only for defining constants (like G), but also just for his own understanding of things. I choose a constant and stick with it. gr/cm^3 really isn't very useful to me as a measure of density. However, it's what Hellings is using and it allows me to directly compare our data before I get too deep into a process (to catch mistakes before they propagate). What I need then is to make converting to and from different units simple, and to do that, the data types holding these values must be classes.

Units I found useful to define as classes

Force, Mass, Pressure, Density, and Distance

The first step was to simply create my enums (enumerated data). Their values are not important (they are just default integer values). But what is important is that I see a drop down of distance_units enums when a method has been defined to accept this argument type. I didn't create every possibility for all of these though, these are just convenient units I intend to use in the near term.

Classes then accept these enums as arguments along with one or more double [double = .net numeric datatype] values. This tells the class to create, say, a new Mass with the incoming value 200 measured in kg.

Outputting the mass in different units is done by referencing a different property. For example, given a particular mass object M, M.kg() will return the mass in kg whereas M.earths() will return it in earth masses. Likewise, it is easy to create a new mass by calling the constructor and passing the appropriate arguments: a double value for the mass, and an enum specifying the units:

_Mr = New Mass(_LMr / _LFinal, Mass_Units.suns)

_Lmr and _Lfinal are double variables coming from the differential equation solution, and the units are in Sun masses (M₀). Thus this line creates a mass object with a Mass in Suns whose value is _Lmr / _LFinal. Once the_Mr mass object is instantiated, the value cannot be referenced except through an appropriate property: kg, grams, earths, suns, etc. Each of these ReadOnly properties converts the mass to the units you specify.

Using enum to make look up lists for sub and method arguments

Creating a Mass object with its constructor and readonly methods

Public Enum Distance_Units As Integer

miles

Sun_Radii

Earth_Radii

End Enum

Public Enum Volume_Units

cm3

End Enum

Public Enum Mass_Units

earths

suns

grams

End Enum

Public Enum Force_Units

ftlb

dyn

SunGs

End Enum

Public Enum Pressure_Units

pascal

bar

atm

psi

End Enum

Public Class Mass

'internally all mass will be stored in grams

Private M As Double

Public Shared Operator *(ByVal f As Double, ByVal i_F As Mass) As Mass

Return New Mass(i_F.grams * f, Mass_Units.grams)

End Operator

Public Sub New(ByVal iM As Double, ByVal unit As Mass_Units)

Select Case unit

Case Mass_Units.kg

M = iM * 1000

Case Mass_Units.earths

M = iM * EarthMass.grams

Case Mass_Units.suns

M = iM * SunMass.grams

Case Mass_Units.grams

M = iM

Case Else

Err.Description = "Unknown incoming mass unit"

Err.Raise(0, Err)

End Select

End Sub

Public ReadOnly Property kg()

Get

Return M / 1000

End Get

End Property

Public ReadOnly Property earths()

Get

Return M / EarthMass.grams

End Get

End Property

Public ReadOnly Property suns()

Get

Return M / SunMass.grams

End Get

End Property

Public ReadOnly Property grams()

Get

Return M

End Get

End Property

End Class

With intellisense on, just create the enum defs, add the unit as a method argument, and when it comes time to choose which unit you want, the dropdown list will appear.

I also overloaded * in a couple of the classes, so that a line like this in my program (to run a quick test) can be executed without having to create a new mass object

MessageBox.Show(Force.G_Force(2 * Constants.SunMass, New Mass(1, Mass_Units.kg), 3 * Constants.SunRadius).N)

SunMass is a Mass object defined in Constants based on the Sun's mass in kg, which is easy to look up. Ordinarily of course, you can't just multiply two arbitrary objects because the * operator may not even be defined (and isn't for a non-inherited custom class). With VS 2005, VB allows operator overloading.

Public Shared Operator *(ByVal f As Double, ByVal i_F As Mass) As Mass

Return New Mass(i_F.grams * f, Mass_Units.grams)

End Operator

In the end, the GUI version looks like the screen shot below. The idea is to create a very large collection of points in space inside the polytrope, whereby each point has various properties: a force in G's, N, dyn, or whatever; a pressure, a density etc. From this continuum of data points a DirectX mesh or vertexbuffer of points or lines can be created that will represent the object in 3space.

When putting in Helling's test data of M = 2 Suns and R = 3 Sun Radii, it's easy to verify my results. However, given that I want to compute different values, like the G's (G = force of earth's gravity on the surface, not the gravitational constant above) at a particular point or the force in Newtons, I need to be able to verify all of the different quantities against known values before I carry on with the next stage (graphics). This will become esp important when looking at strange mass/radius combos. The program has to be as correct as it can be, given the model. This simple polytrope program is limited in its ability to model a real star. Later I will take on the more complex polytrope type in the following chapter, which includes temp, luminosity, radiative pressure etc, and will attempt to get through it completely with a find DirectX representation and explicit GUI modeling interface. I will follow that with more from his book because there's plenty in there as time/interest permits.

In the test below, I have created a polytrope with Earth's Mass and Earth's Radius. The final G value on the surface should be 1 G, and as you can see ( blue row, bottom grid), it basically is.

[Earth as a polytrope, where we expect the surface gravity to be = 1 G, and where we expect the radius (measured in earth_radii) to be = 1 and for the mass contained inside the shell at the point r to be = Earth's total mass, we verify this is true in the 2nd datagridview, blue row. This data checks out fairly well with the expected quantities.]

Download the source (VB.net 2005)

Code Files | Constants classes | Polytrope classes | DiffyQ Solver | Form

______________________________________________________________________________________

Testing the interaction with the polytrope and another mass (or collection of masses) [Again, prior to beginning the graphical portion and the effects of masses on the polytrope]

' simple 2-d prog to test out the motion of particles

' in a gravitational field generated by many objects

' this program sums the forces acting on the particle (3 objects, but you

'can add as many as you like) and then

' moves the particles, using euler integration for the first 2 points, and verlet 'thereafter. I prefer to use something simple like this in 2d before trying for 'the whole tamale in the main prog (which is already plenty complicated

' Kurt Bingham 7/2007

Public Class Form1

Dim p(2) As pnt

Dim vp As VertPnt

Public Const G = 0.00000000006673 'kg-m-s

Public Const dt = 0.000001

' a simple vector2 class with x,y, a norm func, and a normalize function

' and overloads + and * as will be needed later on

Public Structure vec2

Public x As Single

Public y As Single

Public Sub New(ByVal ix As Single, ByVal iy As Single)

x = ix

y = iy

End Sub

Public Function norm() As Single

Return Math.Sqrt(x ^ 2 + y ^ 2)

End Function

Public Function normalize() As vec2

Return New vec2(x / norm(), y / norm())

End Function

Public Shared Operator *(ByVal v As vec2, ByVal k As Single) As vec2

Return New vec2(v.x * k, v.y * k)

End Operator

Public Shared Operator +(ByVal v As vec2, ByVal k As Single) As vec2

Return New vec2(v.x + k, v.y + k)

End Operator

Public Shared Operator +(ByVal v As vec2, ByVal v2 As vec2) As vec2

Return New vec2(v.x + v2.x, v.y + v2.y)

End Operator

End Structure

'*************************

'the pnt class really represents objects with mass and therefore

'objects contributing to the gravitational field.

Public Class pnt

Public p As vec2

Public m As Single

Public Sub New(ByVal iv As vec2, ByVal im As Single)

p = iv

m = im

End Sub

End Class

'**************************

'the vertpnt class is a mass negligible object that

'models the sum of forces in the region and moves accoringly using verlet integration

Public Class VertPnt

Public p() As vec2

Public ReadOnly Property C_p() As vec2

Get

Return p(UBound(p))

End Get

End Property

Public Sub New(ByVal p0 As vec2, ByVal p1 As vec2)

ReDim p(1)

p(0) = p0

p(1) = p1

End Sub

Public Sub ComputeXY(ByVal dt As Single, ByVal F As vec2)

'verlet integration based on force F and dt

ReDim Preserve p(UBound(p) + 1)

Dim n As Integer = UBound(p)

p(n).x = p(n - 1).x + p(n - 1).x - p(n - 2).x + dt ^ 2 * F.x

p(n).y = p(n - 1).y + p(n - 1).y - p(n - 2).y + dt ^ 2 * F.y

End Sub

End Class

Public Function Gf(ByVal m1 As Single, ByVal m2 As Single, ByVal r As Single) As Single

'compute gravitational force

Return -G * m1 * m2 / r ^ 2

End Function

Public Function GVec(ByVal pos As vec2, ByVal obj2 As pnt) As vec2

'this computes the gravitational vector at this point

Dim v As vec2 = New vec2(pos.x - obj2.p.x, pos.y - obj2.p.y)

Dim unitv As vec2 = v.normalize

Return unitv * Gf(1, obj2.m, v.norm)

End Function

Private Function SumForces(ByVal p0 As vec2) As vec2

'this sums all of the p gravity objects in play and generates a single

'sigma vector representing the sum of these forces

Dim g As vec2 = New vec2(0, 0)

For k As Integer = 0 To UBound(p)

g += GVec(p0, p(k))

Return g

End Function

Private Function ComputeEuler(ByVal dt As Single, ByVal p0 As vec2) As vec2

Dim g As vec2 = New vec2(0, 0)

g = SumForces(p0)

Dim v As vec2 = g * dt

Dim p1 As vec2 = p0 + v

Return p1

End Function

Private Sub Form1_Load(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles MyBase.Load

p(0) = New pnt(New vec2(10, 100), 10000000000.0)

p(1) = New pnt(New vec2(350, 100), 1.0E+20)

p(2) = New pnt(New vec2(200, 400), 1.0E+21)

Dim vp1 As vec2 = ComputeEuler(dt, New vec2(100, 50))

Dim vp2 As vec2 = ComputeEuler(dt, vp1)

vp = New VertPnt(vp1, vp2)

End Sub

Private Sub form1_paint(ByVal sender As Object, ByVal e As PaintEventArgs) Handles MyBase.Paint

vp.ComputeXY(dt, SumForces(vp.C_p))

Dim g As Graphics = e.Graphics

Dim pen As New Pen(Color.Blue)

For k As Integer = 0 To UBound(p)

g.DrawEllipse(pen, p(k).p.x, p(k).p.y, 10, 10)

For k As Integer = UBound(vp.p) - 1 To UBound(vp.p)

g.DrawEllipse(pen, vp.p(k).x, vp.p(k).y, 10, 10)

'Console.WriteLine("drawing")

Me.Invalidate()

End Sub

Private Sub form1_unload(ByVal sender As Object, ByVal e As FormClosedEventArgs) Handles MyBase.FormClosed

RemoveHandler Me.Paint, AddressOf form1_paint

End Sub

End Class

While I am still tinkering with the display of the polytrope itself, I spent a few minutes creating a form application that allows me to test in simple 2-d, the effectiveness of having particles react in various fields and accelerate. So as not to excessively complicate thing with a DirectX windowed application, I just used the graphics object of the main form and threw some circles onto it. The main point is, using Verlet integration and several forces, does the object move as expected.

In this simple test scenario, we have several masses in the region and one particle. We sum the gravity force vectors of all of the masses at the test point (mass = 1 unit for the test point). Doing this amounts basically to computing the G force at the point, creating a vector between the mass and the given point, normalizing that vector, multiplying it by the G force magnitude and repeating this process for all of the masses in the region. The final vector sum is then applied using Verlet to the point and it will move accordingly.

This application can be pasted into a new VB.net class file and run directly. It can be modified as needed for more bodies or more test particles (they're arrays). It doesn't look pretty, but you're really just looking to see that the particles move as expected (toward the most massive objects but with effects from all masses affecting its path).

7/13/2007
_______________________________________________________

Colored Polytrope 1

Program colors points based on the pressure from blue to red with blue being the highest. Each rendered pixel is a set of points distance r from the center of the polytrope (from the GUI program). The program is also capable outputting the values as an equatorial slice. I still tinkering with the colors, and of course, the program currently doesn't handle collisions or other forces aside from those that went into creating the polytrope (again, GUI program). Ignore the FPS--the computer the caps were taken from isn't the computer that will be running the sims. This program reads output files (.dat) files output by the GUI. It can also save Polytrope files (which contain all of the polytrope data, whereas the .dat files contain enough information to create a polytrope).

Color and BW graphical representations of a polytrope.

With this done, the points can be put in motion by using extensions of the above verlet program functions/classes to calculate motion of the particles in response to forces.

KB 7/16/2007

_______________________________________________________________

7/21/2007 The first moving simulation

This is the first moving simulation. It features one polytrope, the large mass to the right. The secondary object is simply a massive white dwarf (Sun mass or more) drawing the polytrope toward it. The simulation itself needs more work. I am not summing all the force vectors just yet. I just want to test the verlet integration and the DirectX simulation.

Currently, the only force applied to the polybit elements of the polytrope is that from the small companion object. The polytrope wouldn't normally be ripped asunder as shown here. The points on the screen are negligible mass objects which ride the gravity field created by the companion. This is the approach I plan to use.

Kurt Bingham

12/2007

-------------------

Update 08

Click here for a simple program capable of rendering the .dat files output by the GUI program.

Which generated the capture below

Home