[from Pawel] Our tasks, part 1

Hi guys, I propose to do those things:
1. precede your post title with [from Author] bit, so we know who's talking.
2. if a post is a continuation of something previously discussed by you, give it a number like
[from Anthony] Particle motion, part 3. that will really help to find our way around here.
3. If you update your earlier post, which you shold do if you wish so, note that like (updated 01-01-2011)
at the end of the post title.

* * *

now, to the first point I want to talk about: simulations. I would like all of you to have minimal experience in particle dynamics before we move on, and I have some evidence that Jeffrey (experienced in this) and Anthony and perhaps Josh (though I'd like to see his resuklts posted here and shown to us in real-time openGL one day) are already close to being able to move to the next level. I know Robert was busy with hardware/system software as was Josh, so I can understand. Still, this is not so complicated and I expect that you'll post. James, are you going to post your 1-Body simulation here soon?
If any of you want to team up with others, please do so, dividing tasks between yourself so you know exactly everythign about what the other(s) have done, and report together.

on level 2, we will improve and build out your simulations to include physical interactions such as radiation pressure with shadow casting. Any ideas on how to do this? Obviously, we need to synchronize the calculations for all the particles so the forces they exert are computed at the same times. It is possible to
multi-time the calculations and have let's say most particles for which we update forces at intevals dt,
and the rest which don't need it that often, and hence are updated in force only every (2*dt).
It's not necessary and it would complicate our life in this instance.

You have to device a scheme for optical thickness comptation at each particle's location.
I think you need to youse a fixed Eulerian polar grid of points, enough of then so you have at least as much cells as particles. Each particle will contribute a certain contribution dtau to tau, and that will have to vary as 1/r^2. Let's see why. dtau by definition is dH/H, or relative increase/decrease of flux of radiation coming from the star. Radiation transfer equation is describing the change of flux in a parallel beam of radiation reads
dH = -H d(tau)
where -dH is the absorbed/scattered flux.
However, our beam of radiation is spreading since we have a spherical symmetry, so it's better to
multiply both sides of this equation by 4pi r^2, and call L the total flux or total luminosity of the radiation *going out radially* (not scattered or absorbed yet). Then
dL = -L d(tau)
and now we have the comfortable result that if 0% of radiation is scattered/absorbed in a layer (optical thickness dtau=0), then the total luminosity at all radii is constant and equal L=L_* (stellar luminosity).

from this, dtau due to one particle is simply the fraction of particle cross section pi*s^2 to the
sphere of radius r (r=distance from the star, s = radius of the dust particle),
dtau = pi s^2/ (4 pi r^2) = (1/4) (s/r)^2

a good way to handle the calculation is to have a procedure for evaluation of a map of dtau asa function of (r, theta) coordinates. I would take 360 cells in theta and 1000 in radius between r=0.1 (or so) to r=3 or so.
I assume that inital ring of particles will be chosen between r=1 and r=2. If there are a billion instead of million particles, I'd take 3600 cells in angle and 10000 in radius, for a total of 0.36e8 cells. So ok, there will be on average >1 particle per cell, that's fine, as long it's not >>1.

After that, you integrate along radial paths
tau(r,theta) = integral _0 ^r dtau(r,theta),
which can be though of as adding consecutive dtau values to the running sum called tau, going from
inner to outer radii, separately for every theta. integral^r means integral with upper limit r.

On the first usage of this routine, you set A=1 and dtau = A/r^2, neglecting all the precise constants that would be needed. That is done, because no matter what s^2 we take and how accurate with constants we are,
we simply don't know and don;t want to deal with the number of particles of dust.
We prefer to specify global properties of the disk: as we travere the disk in its midplane, the total tau
or tau(infinity), should accumulate to a given number, which is our input parameter for the simulation.
So in the first pass we divide the **average** (over theta) tau(r=r_outer) by the wanted final tau_infinity,
and store that number as the average optical thickness parameter per particle, A. Then I'd recalculate
the tau map with the updated A.

I'd recommend to try the following choices of radial optical depth if a disk:
tau_infinity = 0, 2, 4, 12
* * *
when the table of dtau is ready and the integrated tau(r,theta) is ready, you can then compte the table of
beta = F_rad/F_grav. that parameter, as you know, does depend on the position only via tau:
beta = beta0 * exp(-tau)

since, however, exp(-tau) happens to be the solution of differential equation
d(tau) = - tau * dn
where dn=1 is by how much the index of the cell jumps from one cell to another, you might as well
save the computer some time on evaluation of math function exp() like so:
you set the initial value of function f, which is your numerical approximation to exp(-tau), to 1 at the inner radius r0. Moving out from one cell to another, you multiply the previous value of f by (1-dtau) in the current cell, store the result on move on to the outer radius.
Why? well, exp(-tau) = exp (-integral^r dtau) = product over all cells of factors exp(-dtau), because
an exp of a sum is a product of the exp's.
and that is very nearly equal to the product of 1-dtau factors, because exp(-dtau) = 1 - dtau/1! + dtau^2/2! - dtau^3/3! +...., and our dtau <<1 (optical thickness of one cell <<1).