Radioactive Fallout Model (Radfo) - Part II

Stretched to the Limit: 06/24/10

1. Expanding Surface Cylinder Application:

Due to the request of NOAA's first Radioactive Fallout model developer, Jerome Heffter, could we keep the shape of a circle instead of the equivalent area rectangle. We stated, it wasn't worth the effort and the added CPU time. A new look at this central core of the main 1+ hr cranking program RADFORUG.c , we now extend the theory of overlapping circles for this effort. Previously, we saw how to rotate local wind box corners into the absolute grid and generate xg,yg's by the equations of the box lines to be re-rotated back again to the x,y local wind box grid.

Here we have an improved diagram and explanation in the combined all grid coordinate systems Grids.dia . We have included a blowup of the local wind box in Rotate.dia . Earlier, we assumed rotation around the upwind front center of the box but now we are more realistic and rotate around the center of the bottom of the surface cylinder, inward by cap/stem radius like a flagpole. We also put an extension upwind of the box so it is symmetrical for rotation about this landing base cylinder xgo,ygo center so the 4 corners can be separated upwind-downwind. Then the upwind extension corners can be lopped off.

       A. Xo,Yo Quadratic Roots Centers for Downwind Cylinder:

This previous explanation was fair but now with an expanded diagram and definite formulas for overlapping, more efficient code was developed based upon Quad.dia . A printout sample and application of these roots is Roots.txt .

       B. Stacked Cylinders - Partial & Max:

Previously, we had used the term overlap loosely. What was meant was that this referred to max overlap where all the infinite cylinders are stacked vertically where we have a disgr=1 for that area of max overlap. With any wind, we have a downwind dwn term that separates the top from bottom of the surface cylinder centers. All this is linked together with more diagrams and explanation details in NonOvrl.dia . We noticed that this max overlapped distance is csdia-dwn and everything else is still overlapped but not max, so disgr < 1 but still has a plateau.

            (1) Graphics Support:

Even with the above diagrams, for the purpose of my review and to comprehend this in 3-D better, we generated some diagrams without labels for easy viewing. They say that a picture is worth 1000 words - I hope so here.

                (a) Return of the Bracelets:

Even though we had excellent displays of bracelets such as Brace #2 , lets imagine a multitude, not an infinite number, of more cylinders settling between the top and bottom of the surface cylinder downwind as seen in Multi Brace #3 . The center of the bottom cylinder (brown o) and the top last cylinder center (purple o) are separated by downwind dwn. From the left side to the left arc of the purple o, we have the linear separation of circles along the axis of all the centers -oooooooo-. The vertical stacking is therefore linear also. Off from this axis, an xg,yg point follows along it's arc to -oooooooo-. It's root center is calculated by Xo,Yo Quadratic Roots explained above. The max overlap (disgr=1) is small indicated by the purple ( and brown ). Notice the tight gradient from top or bottom inward to -oooooooo- as the Quadratic Equation will calculate that region also.

                (b) Max Overlaps:

Now that we have seen the x-y display, we now flip it to show the x-z display. However, we have a somewhat 3-D effect here but don't display the y axis but tilted the sub cylinders a bit as seen in Dwn=R . Here we see the tilt caused by the wind shear from the top layer of the settling cylinder. As it lands, there would be more separation due to the surface wind displacement. So let's assume surface wind nil or calm to show the stacking principle. Again, we have linear stacking between the left aqua vertical lines. This top cylinder's left circumference had it's center shift dwn=r downwind. Again, we have our max overlap (disgr=1) between the bottom aqua ) to the top green ( or between their respective color vertical lines.

                (c) Minimum Max Overlaps:

We shift the dwn another r units so now dwn=r+r = d diameter units as seen in Dwn=D . Again, we have the same linear region but now the bottom aqua ) and top yellow ( are at the same downwind x value shown by the vertical yellow line so there is max overlap (disgr=1) only at this x=r value. This also shows that if dwn < d, we have max overlap (disgr=1) but if dwn > d, no max overlap but still have an overlap plateau region with disgr < 1 as depicted in NonOvrl.dia and below.

                (d) No Max Overlap:

Now we see this non max overlap in Dwn > D . Here we see a longer plateau region but all disgr < 1. There is a linear increase between the left and the right aqua line but plateaus after there to the top yellow (. Since that is to the right, not left of the bottom aqua ), disgr < 1. The means no change in the number of vertically stacked cylinders because the number of overlaps is constant. This is really seen there as we go to the right of x=r, as we lose the lower cylinders on the left but gain the upper ones on the right so net gain is -1+1=0. However, to the right of the first yellow line, we just lose cylinders so 0-1=-1, the inverse left of the left aqua line where we just gain 0+1=1. This matches the 2nd & 3rd diagrams in NonOvrl.dia where we denote ( ) as bottom, { } as top. We have the following regions:
linear ( disgr< 1 {
max { disgr=1 )
linear ) disgr< 1 }
With the max overlap stacking, we lose none and gain none so 0+0=0.

       C. Surface Cylinder Zones

We saw where we explained the linear up & downs and the plateaus but now how do we calculate the disgr for Dwn > D? By using NonOvrl.dia again, the middle 1/3 diagrams progress to Dwn > D. Below this, we look at the x-y plane again from above and break these areas into 3 zones for volume considerations. With this elongation, disgr < 1 but we consider we see the top of the Zone 3 rectangle as vertical height (rat)=1 to be scaled down later.

To consider Zone 1, look at Dwn=R again. We want to sum the volume of the vertical cylinders between the 2 aqua lines or from -r -> 0, or rather use 0 -> r, just flipping the area to the opposite direction but same interval for integration. We also set r=1.

We use VolCyl.c
for the integration of volumes. Reviewing, x^2 + y^2 = r^2, y^2 = x^2 - r^2, y=(r^2 - x^2)^1/2 since we can have -,+ square roots. We want y dx or Int [(r^2 - x^2)^1/2] dx = r/2 [x*sqrt(r^2 - x^2) + r^2*asin(x/r)] with lim x as x1-> x2. Note: The limits x1 & x2 are the distance from center o. Now we have the partial area for each circle but the cylinder height dz must be chosen. We let dz=dx but to see the convergence, we try dx=.1,.01 and .001. So from 0 to 1.000 height, that's 10, 100, 1000 cylinders respectively. As x -> 0 to r, x1 remains 0 but x2=r-x.

One little problem is the sqrt(r^2 - x^2) term. We know that 0^1/2=0 but the computer dislikes it and gives a domain range error. What we do is set quad=(r*r - x*x). If (quad < 0) quad=0 and if(quad > 0) quad=sqrt(quad). This way, quad is always >= 0 without having to take sqrt(0). For each vertical cylinder, the center moves +dx so the -dx brings us back to the 2nd vertical aqua line again, thus keeping Zone 1 in the same area.

With the integration formula above, we are integrating only one quadrant of a circle. Later, it is doubled on both ends to complete the volumes in Zones 1. The results are shown in VolCyl.txt . We chose 100 cylinders to compare the volume for each dx. The y for each x is also calculated but not needed. The limits x1,x2 are also given then the cumulated vertical height of the cylinders, then it volume and it's cumulated volume. Notice the first area is a=0.7854. That is pi*r^2/4 so the integration code checks out with r=1. Notice that the cum vol (cvol1)=0.4560 and dwn=r=1.00.

The 2nd table below in VolCyl.txt is quite interesting for vol2. View Dwn=D again. This time, we double the cylinder's tilt slope from Dwn=R to Dwn=D. Notice the x limits, x1 & x2 are both changing by -.01. We start with x1=0.000 & x2=1.000 and end with x1=-1.000 x2=0.000 but the distance interval (x2-x1) between the two is always 1=radius r. Since the limits x1 & x2 are the distance from center o, each vertical dz and horizontal dx cylinder, the center also shifts by dx, so the -dx brings us back between the aqua and yellow lines again for Zone 2.

This time we notice that vol2 is increasing instead of decreasing up to a max (0.00957) at n=50, 1/2 the height or the light gray cylinder, then decreases. Finally, we note that cvol2=0.9120 is double cvol1=0.4560 because Zone 2 was always integrating more area of the circles above the base circle than Zone 1.

       D. Accumulated Volumes vs. Dwn

Now that we have the methods for the accumulated volumes within Zones 1 & 2, lets increase the dwn to 1.25r for Zones 1 as indicated by the 3rd table. Note that now cvol1=0.3648. Now increasing dwn to 2r, cvol1=0.2280 as we see in the 4th table. Finally, the 5th table has cvol1=0.1824 @ dwn=2.50.

Since we have an idea of the table values, let's see if the total cvol1's converge with increase # of circles. This is all given below the tables in VolCyl.txt . First we backtrack to n=10 cvol1=0.1962, then we had cvol1=0.1824 from the n=100 table. Then with n=1000, we had cvol1=0.1810 or near a .1% decrease after the near 1% decrease from n=10 -> n=100. We don't need more accuracy than this so we will now use n=1000 for all dwn's. So we backtracked with n=1000 to cvol1=0.2011 for dwn=2.25, cvol1=0.1616 @ dwn=2.80 and cvol1=0.1508 @ dwn=3.00.

Now, we run the unity case (dwn=1.00 r) again @ n=1000 and we have cvol1=0.4525, where we first had cvol1=0.4560 with n=100. This is still < 1% decrease. Obviously notice the decrease with increase dwn. Are they related? Yes, they have to be because the downwind ratio dwr means that 1/dwr as many cylinders are in Zone 1 because the cylinder tilt is dwr times as much as it is for dwr=1r radius. So, instead of using the integration formula 1000 times, since we did it once, all other accumulated cylinder volumes are inversely proportional to dwr so that cvol1=0.4525/dwr for Zone 1 & Zone 2 (2*Zone 1).

       D. New Local Box Heights for DWR > 2:

Why do we need these values? If we know the relative volumes of the local wind box by summing the area zones, then it's actual volume has to match the original volume of any surface cylinder before deposition. With no wind, it would be area O = pi r2. See the middle one-thirds diagram NonOvrl.dia again for these zones and the derivation of rat = pi r^2/[12 x .4525/dwr + 2r^2(dwr-2)] underneath. For a range of dwn's, see the last table in VolCyl.txt again. Notice the smooth transition from 1 -> >1 for dwn > 2. Note, rat > 1 at dwr=2. This is due to the lower linear values from 0 to 1 up the circular sides of the beginning & end cylinders so the rectangle has to make up for those heights < 1 for the total vol to be 1. When dwr=2, we have end to end touching circles with some rectangular areas, but not much > 1 (rat=1.1571).

2. Application in RADFORUG.c :

Referring to the last 2 diagrams in NonOvrl.dia , again we would analyze the complete end circles by the Quadratic Roots method explained in Quad.dia linearly to the Zone 3 rectangle but now disgr is scaled back from 1 by rat so disgr=disgr*rat. When x,y is inside the rectangle, then disgr=1/[2r^2(dwr-2)] or disgr=1/[csdia(dwn-csdia) and disgr=disgr*rat again.

       A. DISGR by CPU or Files?: 08/24/10

As we note that much cpu time is required to run RADFORUG.c , then would it be faster to calculate DISGR values outside the program by another program and store the results in binary form and access these by file position in RADFORUG.c ? We went ahead and took that chance and since this update is 2 months later, we have a lot to report on.

Review - DISGR refers to DIStance Gradient Ratio, the shifting of cylinder centers downwind determined by the roots explained above and seen in Multi Brace #3 . Since we had ratios near 1 but not exactly 1 for overlaps and correction factors (rat) just mentioned above, wouldn't it be better to apply these correction ratios with each DISGR for RADFORUG.c because it is also pre calculated? Actually, to apply these corrections while running RADFORUG.c, all the DISGR would have to be in another grid and the correction factor calculated after all DISGR were done, using more CPU time and a mess of more programming. If we had equal dwr's at times, why spend more cpu time when it could be done only once outside the main program. The program responsible for this is OVERLAP.c . Actually, it includes dwr < 2 as we had before and now dwr >= 2.

Skip below and go to latest update Stretched to the Limit: 06/24/10 & 03/02/14?

            (1) Range & Precision:

We know we need a dwr range to < 2 and the new formulas applied to dwr >= 2 say to dwr <= 50. With a stem radius of 3.62 mi, this extend to 181 miles. If the top & bottom of the surface cylinders took that long to settle out, then they were small particles and their landing centers are even further out and probably out of grid range. The cap radius top would be 543 miles, out of grid range but maybe partially in the grid. If dwr > 50, we set dwr = 50.

The spacing between the grid points is .1 mi to .1414 mi. What's that .1414 mi? We will explain that later below. So we go for 1/20 or .05 mi for 0 < dwr < 2. That is we round off to nearest .05 mi so .05/2 = .025 mi as max possible error. Since the grid point assumed .1, then rounded .1/2 = .05 mi, greater precision than the grid. For dwr >= 2, we cut back on the precision and use every 1/4 or .25 mi. Rounded by +,- .25/2=.125 mi almost matches .1 but less than .1414. So worse error would be near 2.0 or .125/2 = 6 1/4%. Actually, it's the difference in the disgr values that would be the error and not the distance itself.

            (2) DISGR files Part 1:

Using all the formulas described above, OVERLAP.c generates 2 sets of files. The first set matches the dwr with the needed correction ratio for stem and cap for dwr < 2 and dwr >= 2. First, we have Stem.l2r and Cap.l2r . Since each line differs by dwr 1/20 or 1/4, then we can round off any dwr and since each line is forced to have the same character length, the rounded off dwr finds it's associated file line position rapidly by it's appropriate formula. Notice the correction ratio becomes less accurate with dwr=2 to about 3%. This factor is applied to each DISGR by division disgr/rat. The 'filp' shown in the table is the file size of the actual disgr file and the ix2 refers to the downwind row length of each dwr file. Then the # of rows yi is filp/ix2 where xi goes from 1 -> ix2, and yi from 1 -> # of rows. These conversions are seen later.

Notice the larger file sizes filp for the larger cap over the stem radius. For dwr >= 2, we have Stem.g2r and Cap.g2r . This time, we have even larger files but the increasing rectangle area portion of the whole area keeps the ratio correction factor to less than .1% from 100% for most of the dwr > 2.

            (3) DISGR files Part 2:

For the actual binary files, we use the dwr name with 's' for stem and 'c' for cap and use the dwr value with it's decimal to denote the suffix part of the file name. The dwr's < 2 & dwr >= 2 are combined because their dwr values are different to distinguish each file. No sense listing the actual binary characters here but view the directory list of Stem.dir and Cap.dir . Notice each has 232 files of 14,581,093 bytes and 130,682,062 bytes respectively. This 145 mb of data of 464 files is outside the program but larger than the 40 mb absolute grid file. A little sacrifice of space for CPU time would be worth it.

            (4) Graphics Plots:

In order to verify the DISGR values, we create a smaller program to read and plot the data to verify both. This program is called POVERLAP.c . We read the character c, idisgr = (int) c then disgr=idisgr/255 then idisgr=disgr*16 for a 16 color plot. We chose a cap range of dwr's as seen by following:

Cap .50
Cap 2.00
Cap 2.50
Cap 3.25
Cap 4.50
Cap 5.25
The final radfo plot is 1000's of these single plots superimposed or scattered and rotated throughout the absolute grid. We now isolated a one case Cap 3.25 to test a new rotate order which becomes one new hairy topic.

Update Stretched to the Limit: 06/24/10 & 03/02/14: