Radioactive Fallout Model (Radfo) - Part II

by Buddy McCloskey     02/02/2009     Contact

5. Use of Radiosonde Soundings:

We have a sample sounding from Pittsburgh, PA - 12z 3 June 2002. These soundings are taken every 12 hours nationally at 12z and 00z daily. We chose a rather poor example of a typical weather day to perform a rugged test of the program. We will be using more realistic data later. However, we ran this data to show how the older RADFO versions plots appeared and to exhibit the errors with this older method used.

       A. Wind Shear on a Dust Layer Cylinder:

            (1) Garbage In   =   Garbage Out:

In the earlier models, recent attempts by some NATO countries developing a Radfo model, EPA Air Pollution Models, and my work with the City of Phila, it was discovered that modeling dust layers was a funny experience. The concept is to examine each wind layer with a particle distribution. As shown above, we have a distribution of particle sizes in each layer. With each particle size, we have a different fall velocity in considering a wind layer cylinder, each particle size cylinder falls at different rates. We have the same cylinder layer of air but each group of particles with its different travel times land independently of it different neighbor sized particles. Thus, it's the same physical cylinder landing over and over again in longer downwind locations due to longer travel time of the higher layered particles and/or smaller particles.

This is why we needed many particles so that the difference fall velocities are less apparent. Another parameter is the thickness of the wind layer. The less the thickness, the earlier it passes through the layer. With fewer but thicker wind layers, a specific dust cylinder would land further down wind. What's then is the problem? It's the smoothing of this real time data. We need an overlap of settling dust cylinders. For instance, in the real world, consider a snowstorm. In a given region of generally flat terrain, we have many different sized flakes and unlimited wind layers. The results show that over the same region, the final snow depth is generally the same. This is the natural way in the nature of things, we would have very smooth contours of snow depth if a dense grid reported their correct measured depths.

            (2) Early Examples:

If we examine the settling plot of each dust layer, we could then determine what a smoothing contour program has to work with. We assigned each different particle size a different color. For our first example, we present Radfo #0. It's obvious the wind shear and different particle sizes make for a difficult analysis without a smoothing technique. Our next example Radfo #1 shows a similar pattern. In a better plot, see Radfo #2. A different plot is seen in the final plot Radfo #3.

            (3) A New Method:

From the examples above, smoothing must be applied. However, do we do this after the dust has settled or before? How do we do this before? Let us examine the 1st layer of wind next to the ground. In all the previous models, we would take the mid level height to apply the fall velocity for a representative particle center. At that time, the entire cylinder is shown on the surface as a circle indicating that the cylinder collapsed or falls to the surface.

In the real world, the bottom of the cylinder is dragged across the ground and deposits its lower levels of dust. The upper level of the cylinder continues to fall from greater heights, greater travel time to land further downwind. If we represent the top layer as dark blue and the bottom layer as gray, then a wind from 240 degrees is shown by 240 Shear . Notice the overlap area shown in lighter blue.

Here we call this the top and bottom layer. To show our method, let's divide the cylinder into 15 layers that can be represented by 15 colors later (cylinder.dia) . Instead of using the midpoint + to represent the average height of the particle cylinder, we use 15 layers instead. When the top layer falls to the surface, we have the blue outer circle and the lighter blue overlap seen in 240 Shear . Within the intersecting circles, all 15 layers have deposited all its particles size. This would mean that here in the lighter blue area, we have our highest accumulation of dust stacked vertically.

                (a) Dimensions Needed?:

In order to develop our final method, we needed to find the width of these colored circles. Therefore, we first found the distance between the 1st and last layered circle. In our 1st approach, it thought that a distance d was needed. For a wind direction of 300 degrees, we calculate some other inner dimensions as seen and explained in WVCP300 . This was an early attempt to calculate areas. The next paragraph partly explains the formulas we thought we might need.

Now we use a wind direction of 300 degrees to show the rest of this explanation. Let's temporarily subdivide this cylinder into 9 layers and show the plot of all the layers in a 9 layer plot . Here we see the colored centers corresponding to its colored ring or circle and the mid of the overlaps by (+) as the cyan (baby blue) center. Notice the gray lines intersect into the overlap center there. If we have a grid point G anywhere on any gray line inside the circles and extend it to the overlap center (+), all the distances between the circles appear to be in proportion. We can use proportions to see how far G is away from the brown circle towards the white circle.

However, this was not the correct method. How much error, about 13% to high. As we have higher winds and the circles separate more with less overlap area, the error increase to over 30%. After connecting the line from G to origin o instead of the overlap, there was a slight improvement. These two methods was a complex solution substituting the formula for a line in the circle formula and then finding its roots boundaries on the circles. Before abandoning this, we used it to plot the gray circle overlay on the blue circle as we searched for the beginning and ending pixels by their roots across the gray-blue area and coloring the blue and light blue area by comparing pixel colors. Thus we have 240 Shear . At least we used this to count these lighter blue pixels (100% overlap) which came to 48% of the original circle area.

                (b) Determining the Grid:

In the previous dust models, we can determine the concentration of dust related to its deposition height here. We let it fall where it may then run the whole grid to determine where the layers of dust land. Rather here, we choose a sub boundary section of the grid that is only involved. So in our program, we need to find the minimum and maximum range of the 2 circles. This is indicated by the red border in the diagram 240 Shear again. This is simply the +,- radius applied from the centers of the 2 circles. The center of the blue circle has shifted to +dx, -dy as seen in WVCP300 .

                (c) Another Attempt:

Starting from the simple to the complex, let's assume the wind is west or consider that we have a grid based on downwind and crosswind. This would be a mess to calculate and then a new grid for another wind direction and cylinder and then convert each grid point to the absolute (S-N,E-W) grid field. We worry about that later. We refer to a new set of diagrams but explain in detail using Bracelets #1 . That's exactly what these would represent and actually borrowed 9 different colored ones from my granddaughter's Christmas present to visually see this preliminary concept's end results before proceeding.

With the use of the full screen by an expanded scale, we see the lines A through F intersecting to the original center (most left center light brown). Definitely, all circle arcs intersecting line A are in proportion as each height interval, has its corresponding travel time interval the same, all headed downwind on A. Even along B, we notice no spatial differences. At C, there is a difference near the periphery and then along D, the difference is noticed and is obvious along E. So even to the original centers, the spacing is not linear away from the downwind x-axis. We must examine another method (# 3).

Consider a grid point in the upper left (green circle) of Bracelets #1 . Now we only have the grid coordinate and have no idea that it is closest to the green contour yet. How do we determine which circle it's nearest to? The equation of a circle is (x-Xo)^2 + (y-Yo)^2 = r^2 where Xo,Yo is the coordinate's center of the circle. Here x,y is really a grid point xg,yg and the center of the shifted circle Xo,Yo represents x,y from the original center (0,0). Instead of finding a point on the circle solving for x, we have a grid point on one of infinite circles so we already know the x,y coordinate. This circle and all the circles have their center on the x axis or y=0. That is, the circle has its center 0 units from Yo or Yo= 0. We also know the radius r (the cap or stem of the nuclear cloud). We sub Yo=0 to simplify immediately before squaring that term so now the equation of that unknown circle is:

(x-Xo)^2 + (y)^2 = r^2, squaring we have
x^2 - 2Xo*x + Xo^2 + y^2 = r^2 putting in quadratic form for Xo
Xo^2 -2x*Xo + x^2 + y^2 - r^2 = 0 coefficients for the quadratic
a = 1, b = -2x, c = x^2 + y^2 - r^2 satisfying the form
ax^2 + bx + c= 0 for x = [-b +- sqrt(b^2-4ac)]/2a.

We note the roots, x1,x2 where x1 is the left side on the circle, x2 the right. The left side is applied when the grid point is to the left of the overlap center (cyan o with +) or xfc/2. This root on the x-axis (line A) is shown as the small white o. Here it is seen as just to the left of the green circle. As a grid point, it was further up along that circle. This value at the grid point is calculated by the ratio of the root x1 to the length between the 2 centers, which is the same as the ratio between its 2 circumferences on the x axis. Just to the left of the brown circle, the ratio rat >= 0 up to 1.0 for the last (white) circle. So for x1, rat=x1/xfc and for x2, rat=(xfc-x2)/xfc or computed from the furthest circle to the right side crossing the x-axis towards the white circle. We cannot have a ratio > 1 in the white area.

Other examples of this are seen in Bracelets #2 , Bracelets #3 and Bracelets #4 . Notice that it appears that the root o, not matching its nearest color circle this time so a color contrast is visible, appears to be further away from its nearest circle especially in the last example where grid and root are both white. Is this an error? The answer is no, because it is placed at the ratio distance between the 2 circles and in the tight complex region, it is still halfway between the red and purple circles as it is on the corresponding intersect on the x-axis showing ratio, not distance consistency. Also, there is a value there as the original circle advances downwind to the shifted vector circle. At the same relative position in 240 Shear , it would appear as out of the region. It is when only looking at the initial grey and final blue circles but Bracelets #4 is more accurate as it shows the transition between the two which is in the region of intermittent circles.

                (d) Infinite Layers:

How many points are there in a line, or how many arcs of circles can pass through a line segment? Even though we showed 9 here, each layer would take 1/9 of a distance and G could be rounded off to the nearest 1/9 or nearest circle. If we had 10 layers, then to the nearest .1. But why round off? We can have an infinite number of layers by just using the actual ratio. If we had a ratio of 4.026, then that would be the same as considering that we are between the 4th and 5th circle out of 10 or 4th and 5th thousandth's circle out of 10,000 or at circle # 4026. Therefore, we smooth the data at its source by considering an infinite number of layers producing a infinite number of circles-cylinders landing on the surface as they fall.

                (e) Grid System Abandoned?:

Now that we have a value for a point, which is true for the wind direction (local) coordinate system. Let's assume the real world, our wind direction of 300 degrees. Well temporarily, let's keep Xo and the final Xo (xfc) downwind along the x axis, ignoring direction. All we have to do is bring our absolute grid system temporarily into our local wind coordinate system by rotation of the angle between them. So we consider that the 300 degree wind blows across the x axis (270 degree) by 270-300 = -30 degree. To track a minimal of absolute grid points rotated 30 degrees counter clockwise, we show Rotate 30 . The green circle is the absolute grid (xg,yg) and it is rotated along the gray arc until it stops at the red circle or the temporary wind coordinate system x,y. Each arc shown is 30 degrees cc rotation. Now we apply our root calculations from above to this temporary rotated grid but we never show its plot but save its (x,y) ratio value to apply to the absolute grid xg,yg where it was rotated from. The formula for the rotation is:
x=xg*cos(deg) + yg*sin(deg) and y=yg*cos(deg) - xg*sin(deg) where deg = -30. Our final result is seen in Contour 300 . Notice the overlap center in red, the original circle center in white and the wind vector shifted circle center as blue mini-circles o.

                (f) Original Volume   =   Final Volume :

Now let's consider the volume of this cylinder. We may think of this as a high stack of pancakes. We push from the top and all the pancakes are moved as shown by example in cylinder.dia. Then we push down vertically until they are all squashed and there is no air space between any. Even though they are squashed, the total weight is the same. We may have made it denser due to the squashing, but ignoring that, the mass is the same. Assume we didn't make it more dense, the volume is then also the same. A sketch of these is seen in the 3 diagrams of sideview.dia. A simple concept but now to show its application. We have calculated a ratio value rat for the grid point G shown in the bottom fig. as over 1/3 the way up the steps.

                (g) Application of the Ratio:

Instead of pancakes, we imagine these as dust layers whose thickness cumulated one upon another. Of course, this would be really loaded with dust but shows the concept. Or we can imagine that a real pile of dust is just scaled up or down in the vertical. Or, consider in the real world, the snow example, but is now dust instead. We now represent each layer by a color or 15 colored layers and view this from the top as just seen in Contour 300 . Notice we have a grid overlaid and we just calculated the height or value of each grid point. For our new method, the correction factor of the weighted values is now 0.998 and at worst is .981 for the Kissing Circles . In any event, we will use this method with < 2% maximum error.

            (4) Multiple Plots:

                (a) Common Initial Cylinder Centers:

In order to extend this method for multiple dust cylinders, we put each cylinder's dust deposit into an absolute grid. We will have different cylinder position landings based on wind direction, wind speed, particle size, cylinder height and cylinder thickness between two wind levels. We also use the nuclear cloud cap or stem for the cylinder diameter. For a possible cylinder, let its base height at the surface be a thickness of 2000' with a particle fall velocity of 2.7'/sec. We used a stem radius for 1 mega as stemr=3.62 mi. Here, we combined a set of 4 wind combinations as dir=300 deg @ vel=15 mph, dir=315 @ vel=30, dir=240 @ vel=25 and dir=210 @ vel=20.

For our test case, we use a small absolute grid min and max denoted by agminx, agmaxx, agminy and agmaxy in mi. If we define the dx and dy increment as dxy, we have idxg=(agmaxx-agminx)/dxy and idyg=(agmaxy-agminy)/dxy # of xg and yg grid points along each axis. We have a subdivision of this absolute grid by the red rectangle around the plots above defined by xmin,ymin,xmax and ymax where the number of grid points in this grid sub system was idx=(xmax-xmin)/dxy and idy=(ymax-ymin)/dxy, the same method as in the absolute grid system. For all our test cases here and above, we let dxy=.05 - .08 mi. After initializing the grid file with 0's, we summed the 4 grid values and put the grid sub system values into the absolute grid by ixg=(xmin-agminx)/dxy+ix and iyg=(ymin-agminy)/dxy+iy where ix,iy represent the column and row grid # of the sub section grid and ixg,iyg the absolute.

Now we convert ixg,iyg into a one-dimensional file position by filp=iyg*idxg+ixg. We read what was previously stored there, add its new value and rewrite back to the same position in the file. With the sub section grid, we don't have to bother with any other section of the much larger absolute grid but go directly to its sub section. In our test case, that's roughly 80 mi^2 compared to (100mi)^2 = 10,000 mi^2. In our first multi cylinder plot, we have all 4 original cylinders centers starting from 0,0. Thus, it's still has a kind of max center but depositing in different directions as seen in Multi Plot # 3 .

                (b) Shifted Initial Cylinder Centers:

Now we need the real word, cylinders moving by wind vectors and depositing. In this next case, Multi Vector Plot # 1 , each of the 4 cylinders has a new center by assuming that each cylinder has moved by the upper level wind vectors as it travelled down through the lower levels before it landed within the surface (0-2000') layer. This shows a slightly different pattern than having 4 common centers with no translation. This translation is just a grid shift (dx,dy) by the mean wind vector between the reported wind layers. These are added to the coordinates for the ixg,iyg term above.

            (5) Cylinders Aloft:

Above, we concentrated on cylinders near the surface we called the 1st layer above the ground. The bottom part of the cylinder dragging the ground has no wind direction and velocity. Thus the top of the cylinder moved at the mean wind speed within this layer relative to the wind speed of 0 at the bottom layer of the cylinder in contact with the earth's surface. Now extending this principle to upper levels, the tops also move relative to its lowest wind level. That is, it's the wind shear between the two levels that reshapes this cylinder layer just as what happened at the surface. However, it's not as sudden as at the surface but the cylinder may become even more distorted with the higher wind velocities and direction shifts aloft.

Another difference is, as the cylinder is deposited in the 1st layer at the surface, it was deposited and we have our flattened pancakes. Aloft, the cylinder keeps its thickness until it lands on the surface but becomes twisted laterally. In our examples, we assumed the cylinder was vertical with no tilt. In the real world, these are twisted or tilted cylinders. This means that with the tilt, we have multiple centers if we divide this surface cylinder into layers. This is what we shall do and would appear as seen in Surface Cylinder .

                (a) Attempted Displays:

It would be nice to show the trajectory and twisted cylinders as they pass through the different wind layers. However, much time would be spent in using the appropriate 3D graphics program although it could be done later. For now, we could present the top and bottom layers of the cylinders again by our top view looking down to the ground as the shear reacts on each cylindrical layer. Since we want to keep them as pairs, we might as well indicate the translation of both by the mean wind vector also.

                (b) Wind Vector Formulas: TRAJ.c

For now, we will use some dummy test wind data and develop the formulas needed. Regardless of a particle's fall velocity, the particle passes through increasing or decreasing wind velocities changing its horizontal downwind speed during its descent through the layer. The particle starts at the top with a vector speed of V1 and ends up with a speed V2 at the bottom level. So the average speed through the layer it transits is (V1+V2)/2. An example of this is explained best and depicted in Fall Vectors .

In our first example, we have Vector Plot #3 . The white line indicates the track of the cylinder as it passes through each wind layer using the mean velocity derivation above. We stop at the last 2 circles, which we call the Surface Cylinder, indicated by the top (white) and the bottom (blue) as used before. Now we can apply our new depositing method seen above. However, before we do that, we used the Pittsburgh, PA strange wind data to show its vector path. We set each wind layer at 1000' intervals just for demonstration purposes. Since there are more levels, higher heights with more velocity, we just were able to squeeze its display on the screen as seen in Vector Plot #2 with smaller circles.

We now extend Vector Plot #3 to Vector Plot #4 . Now this looks appealing like this is working but we applied the new method to the Surface Cylinder and it fit into the right location after juggling coordinate systems. However, that was a demo to match our original system. The difference is that previously, we assumed the Surface Cylinder was vertical and the bottom layers deposited while the top layers moved downwind before it deposited. Earlier we showed the blue (top) as it settled but positioned from the same center as the bottom of the cylinder because it was vertical.

Now we back up a bit to look at Vector Plot #3 again and noticed we just mentioned that the Surface Cylinder was indicated by the white top and the blue bottom, the last 2 plotted circles. So now we see that the last cylinder is the Surface Cylinder and is not vertical. Since we can only apply our new deposition method to vertical cylinders, we make Surface Sub Cylinders and assume they are vertical. The more, the greater accuracy of applying this theory. Review of (cylinder.dia) again, a straight line through the midpoint of the steps along the tilt of the many vertical cylinders approximate the main tilted cylinder. Here we tried 5 Sub Cylinders and as indicated by the last blue rings in Vector Plot #5 .

Now we have newer results in Vector Plot #7 , similar to Vector Plot #4 but this greater accuracy is evident as the upwind cylinder layers are depositing more dust as they travel downwind over the center position of the tilted Surface Cylinder and hence, a wider area of higher contours upwind, less downwind. Also, the cumulated weighted grid values were corrected for 1/5 for each of the 5 layers and the value came up from .998 to .999 out of a 1.0. We might as well ignore this correction factor now.

                (c) How Many Surface Sub Cylinders?:

This resulting plot shows the advantage of using these sub cylinders while keeping the high correlation. For our final model, we want to predict radiation from the dust fall probably out to 300 miles. Our guess is that we want a grid interval dxy about .1 mi. We used < .1 mi above to show more details for confirmation that the method works. In order for the sub cylinders to be of use, we need the grid points to fall on the fringe area, between the colored contours and the white area seen in Vector Plot #7 again or in the blue rings of Vector Plot #5 .

In the plots, we had the parameters given above in (4)(a) Common Initial Cylinder Centers: We have the travel time of 2000'/2.7 '/sec = 740.7 sec = 0.206 hrs. Now x 15 mph = 3.086 mi. A grid spacing of .1 mi would have 3.086/.1 = 30.86 grid points/sub cylinder or at least a grid point within each deposited sub cylinder. For a diagonal spacing of the grid, we have 3.086/.1414 = 21.82 grid points/sub cylinder. So when we chose 5 sub layers, there would be at least 1 grid for each of the 5 sub cylinders. Actually it would be 3.086/5/.1 = .617/.1 = 6.17 grid points/sub cylinder at its maximum downwind deposit wind axis, less off the downwind axis. What we want is at least 1 grid point/sub cylinder so we won't miss the sub cylinders. That being reworded means that a fringe area of a sub cylinder doesn't lie between two grid points. That is if dh/fps*vel < .1.

The 1st reportable wind layer above the surface by the Pittsburgh, PA sounding was 449-373 m = 76 m= 249'. If we use 250', a minimum wind speed of 2 kt, use mph, then 250'/fps/3600sec/hr x 2 mph = .1mi, 2 x .069/fps = .1, fps=2 x .069/.1=1.39'/sec. Using a minimum cylinder thickness (height @ surface) and the lowest reportable wind speed, a particle with fall velocity of > 1.39'/sec will just touch the grid to apply sub cylinders. Now twice this fall velocity will half the distance it has to travel so this larger particle size has a 50/50 chance of landing on a grid point. So 2 x 1.39'/sec = 2.78'/sec corresponds to a particle size of near > 130 u. However, the main sub cylinder circles will not miss a grid point. Even if the fringe is missed, the bulk of the particles are within the white zone anyhow with these larger particle falling through a lower wind speed surface layer, reverting back to a near vertical cylinder because of such little tilt. To save cpu time and grid size, we may keep the 5 sub cylinders and the grid spacing of .1 mi.

                (d) Only One Sub Cylinder?: 05/07/2009

The increase in the number of sub cylinders multiplies the CPU time proportionally. For a near surface dust source like WTC 9/11, we would apply this. In the real world weather, using the 1st two lowest wind levels from Pittsburgh, PA the surface wind shear on appears as shown in Real Grid Plot #2 and with smaller particles or higher winds as seen in Real Grid Plot #5 . With our diameter of the cylinder over 10 mi, this shows little effect here. Another reason is that we have examined real wind shears that overwhelm this surface effect as explained below.

                (e) Don't Forget the Shear Dear:

As the wind layers below a wind level adds to the translation, the difference in the wind vectors adds to the distortion or tilt of the vertical dust cylinder by the shearing effect just as the surface does. Many have seen this in the side view of thunderstorms. This means a cylinder will change it's shape by it's decent through all the wind levels. However, the wind levels are not fixed, except in the PIBAL (PPBB) section of the original rawinsonde code. The significant levels are determined by a 1 degree or a 10% relative humidity departure from a linear line between the 2 levels below. Too bad we don't have significant levels also determined by wind shear vectors instead. The PIBAL section winds are determined from all these significant levels anyhow so we will use them directly.

Above, we explained the fall velocity translation vectors by Fall Vectors . However, to explain the shear effect and how different wind layer thicknesses affect the particles moving through it, we use part of the same diagram for review but expand to add this new feature as seen in Vector Shear 1 . So our conclusion now is that we have the shear components, des=(velb*cos(degb)-velt*cos(degt))/2 and dns=(velb*sin(degb)-velt*sin(degt))/2 where t=top, b=bottom layers. That is, where translation is half the vector sums, vector shear is half the vector differences. Also, regardless of the thickness of a lower wind layer, the factor that adjusts for this and links it all together is the travel time tt through the lower layer.

Update 12/08/09
A better explanation and simpler method has been developed. We now refer to a new diagram that explains the shear - tilt of a descending dust cylinder as seen in Vector Shear 2 . If identical particles were dropped from the same point and time, would have parallel and simultaneous paths, which is obvious. Here we have 7 layers, the top of a vertical dust cylinder and decreasing winds through each layer. The top level (t) of the cylinder is always moving further downwind than the bottom portion (b) through all the layers, thus the wind shear causes the tilt by the relative positions of t & b. From levl 6 down, they both have parallel paths, so the only difference is the initial starting position of t when it arrives at the top of layer 6. That is, b7 --> t6, (holding b7 constant until t7 arrives at t6). That is tt[V7+V6]/2 = (14+12)/2 = 26/2= 13 assuming tt=1 unit of time.

Our other previous approach was to accumulate the tilt as the top cylinder layer passes through each wind layer below. Since the tilt is 1/2 the shear caused by relative travel vector positions of top vs. bottom, we summed up all the expressions and using numbers, we arrived at the same results as seen at the bottom of Vector Shear 2 . The advantage is now we don't need to keep track of all these shear terms in pttc.cor, just the difference between the top and bottom of the upper wind layer where the top layer of a particle dust cylinder resides before descent. End Update 12/08/09

                (f) Stretched to the Limit:

If the wind velocity is high near the surface and the particles are small or the surface reportable wind layer is thick, the travel time could be large so the top of the surface cylinder may drift well past it's bottom as seen at the bottom of Vector Shear 1 again. Also, if there is much vector wind shear from the aloft layers where these particles fall, the cylinder could be further elongated. With vector wind shear, direction or velocity changes causes a shear. As the distance between the centers increase, the correlation of 1.000 decreases. Therefore, we set a limit where to apply the overlapping circle method as downwind diameters dwd=2 (two circles end to end = rectangle length) as indicated by the Kissing Circles again. After they kiss, the top circle departs downwind as the area swept out by the top approaches a rectangle in shape. Actually, in the fig, the gap, below minimal values detected, has almost separated them already whose value would be as a minimal concentration.

However, we really lose the shape of the circle to a rectangle by diffusion. So if the area swept out is the diameter D > 10 miles, then dwd = 2D = 20 miles travel should diffuse to within the rectangle. However, with dwd = 1, area O = pi r^2, rectangle R area would be >= 2r x 2r = 4r^2. The ratio is O/R = pi r^2/4r^2 = pi/4 = 78.5%. That is for each grid point that is in this larger area rectangle, the concentration would be reduced by a factor of 78.5% from being within a circle.

If a circle is translated down wind 5 diameters (10 radii) as such, we would have (----------), our rectangle equivalent would be [---------]. The true area would be a 4 diameter rectangle plus 1 diameter circle. R = d x 4d = 2r x 4 x 2r = 16r^2. Now plus area O = pi r^2 so we have (16+pi)r^2 = (16+3.14)r^2 = (19.14)r^2. To assume a true rectangle, we would have 2r x 5 x 2r = 20r^2. So O/R = (19.14)r^2/20r^2 = 19.14/20 = 95.7%. So with greater downwind distance, the error becomes less. Actually, we have no error in our total, we just distorted the elongated circle paths into a rectangle which is extending the perimeter to the 4 corners. A preview of our next topic shows this in 2553 O . Here, dwd = 1.6, used for testing purposes, and the 4 blank corners become filled in 2553 R .

If the rectangle becomes longer, what are the values for all the blue grid points within the rectangle? Using unity as the original value or z height of the cylinder, its volume is Rect 1 area x z-value = Rect 2 area x z-value = original Circle area x z-value. In the expanding rectangle area, the z-value is inversely proportional to it's base area. This is seen and explained best in Non Overlapping Cylinder Method.

                (f) Stretched to the Limit: Update 06/24/10

                (g) Tight Rectangle Formulas:

Since the calculation of overlapping circles was seen to be complex, we had used an estimated oversized box to make sure the overlapping circles would be detected by the downwind root equations as explained earlier. This caused much extra CPU time. Now we have a better method, just fill in the rectangle with only the required grid points as they all will have the same value as just explained previously. Since we determine the width of the box by the cap/stem diameter and have the downwind by the cumulated wind shears aloft and the surface layer vector, then we will have the exact dimensions of the box to use it for the rectangle calculations and the overlapping circles that the rectangle surrounds. See this tight fit in 2553 O again.

                (h) Rectangle Rotation:

In 5 (3) (e) Grid System Abandoned?: above, we showed the rotation method in Local Wind.dia and an example in Rotate 30 . We did not test it for all directions and rotate a box of grids. Couldn't we just use the simpler local downwind coordinates (x,y) and rotate them to the local grid portion of the absolute grid points xg,yg? Yes, but the math and rounding off errors, we may not have a unique grid point xg,yg for each x,y. Also the possibility exists that some xg,yg's may be missed while others may be counted twice. This makes for a ragged receptor grid while we are trying to smooth the process. See the update on (g) Tight Rectangle Formulas: a few paragraphs down.

However, this idea brings up a possible solution. We must have an incremented xg,yg only which ones? Lets take this rejected idea about each x,y rotation to xg,yg by just considering the downwind box boundary corners and just rotate this box and calculate all the xg,yg's within the box and rotate them back to the local downwind coordinates (x,y). This is best indicated in Rotate Rect .

The original downwind box is in blue with the x,y mins & maxes at the corners. The box is rotated 'a' degrees clockwise around it's center 0,0 to the red box. We labeled the corners after the rotation with t=top, b=bottom, left=l and right=r. Note that we said after the rotation. The x,y mins & maxes identify the original corners and the x,y l&r t&b after the rotation may not be the same corners because the box may be rotated up to nearly 360 degrees. Since we will have new x,y mins & maxes, lets identify the original corners as 0,1,2,3 clockwise with 0 being the original lower left corner or xmin,ymin colored white. The diameter is D and the cross wind cwy = sqrt(D^2 + cwx^2) where cwx/D = tan(a). However, due to rounding off errors after all this, it is simply better to work with the limits of all 4 red lines (top and sides) of the rotated box.

We need the equations (y=bx+a) for all the red lines. The slope (b) of the sides are equal because they are parallel, as well as the top and bottom. That is br=bl=(ytl-ybl)/(xtl-xbl) for the side slopes and bt=(ytr-ytl)/(xtr-xtl) or bt=-1/bl for the top and bottom slope. Now intercept (a) left is al=ytl-bl*xtl or al=ybl-bl*xbl and right side intercept (a) is ar=ytr-br*xtr or ar=ybr-br*xbr. The top line intercept (a) is at=ytl-bt*xtl or at=ytr-bt*xtr and bottom line intercept (a) is ab=ybl-bb*xbl or ab=ybr-bb*xbr. The better estimate of cwy is the distance between the 2 lines, so at x=0, it's the difference of the intercepts cwy=at-ab and remains the same all along the top line.

We now have new max and min grid xg's & yg's. Depending on the side slope, we see a + slope has the xgmin at the bottom left, and xgmax at the top right. With a - slope, xgmin is at the top left and xgmax at the bottom right. We don't really need the ygmin & ygmax as we have the line equations for them as xg goes from xgmin -> xgmax. Now the # of ix dxy (.1) to scan is rni=(xgmax-xgmin)/dxy and the initial ix is ix0=xgmin/dxy. So our loop is ix = ix0 -> ni and xg=ix*dxy. With xg given, we chose to use the upper line equation to find yg2 by subbing for xg in yg2=bt*xg+at. For the bottom line yg1, we also sub xg into it's equation yg1=bb*xg+ab or could have used yg1=yg2-cwy also. Now we have the yg limits and the corresponding iy limits are from iy1=yg1/dxy to iy2=yg2/dxy.

We have swept the inner rectangle for all the grid points. Now we must include the corners. Looking at Rotate Rect again, we see that left of the cwy line (xg < xtl), we use the left side instead of the top line equation, namely yg2=bl*xg+al, for a +bl slope. At the right corner, we start yg1 on the right side instead of the lower line by if(xg > xbr) yg1=br*xg+ar. For corresponding -bl side slopes, if(xg < xbl) yg1=bl*xg+al and yg2=br*xg+ar. For small particles, the total travel time may be too long to land within the action box's absolute grid limits. Hence, that particle cylinder is ignored so we loop to the next particle size id. However, part of the cylinder may land within the box even if the upper parts of the cylinder lands past the absolute grid box limits. Thus, the total concentration (100%) from that cylinder will be less than 1.000 within the box, but no problem since we have to cut the downwind limit somewhere.

                (i) Comparing Plots:

The summated multiple shear effect from the winds aloft overwhelms the surface shear effect by the ground so the number of surface cylinders is meaningless and time consuming. Pictured is a 5-sub cylinder surface cylinder with no multiple wind shear effect in Flash Non Shear 5 Cyl. Notice the details but the ragged effect. It probably would smooth with some math, but with the cumulated wind shear aloft effect and especially the expanding rectangles, we reproduce the same pic with this smoothed effect in Flash Shear 5 Cyl . As for the number of cylinders, now we use only one as seen in the next pic Flash Shear 1 Cyl . Can the difference really be seen and if so, is it worth it? The most difference I saw is in the yellow solid contour along the downwind axis.

                (g) Tight Rectangle Formulas: Update 08/24/10

       B. Extracting the Radiosonde Data:

In order to simplify and not overload the main RAM memory in the RADFO program and not have variable name conflicts, we will run a series of batched programs that create files for subsequent programs. This is also done in the dust inventory process mention in the very beginning of RADFO I. In the original Pittsburgh, PA sounding, we strip off the excess header and trailer info that the storm chasers use. The first program is called RAOBDHM.c = RAOB Delta Height' from Mega. That is, while extracting data and converting input units, let's use the data to calculate some input parameters also. We only need pressure, height, temperature, wind direction and speed. However, we generate some other mean layer parameters also. Actually, this Pittsburgh, PA file has also been extracted from the original radiosonde code by the Storm Machine Web Site.

Update 12/08/09
The Storm Machine Web Site gives us the original soundings at 00Z & 12Z. However, it is made available a few hours later. By that time, it is better to consider predicted soundings and not just for 1 station, but a grid of stations. This is a trajectory path explained in more detail in our next update in section (2) below. End Update 12/08/09

            (1) Pre Processor Programs #1: RAOBDHM.c -> RUCDHM.c

                (a) Nuclear Cloud Heights:

Using the formulas from 3. Nuclear Cloud Growth: (subroutine megap) and the 2. Bomb Yield Parameters:, 3. A. through D. above (dhb), we scale and determine the cloud cap top and cap bottom (stem top) in feet above ground. We convert heights from feet to meters msl to find the wind velocities then convert the plume rises back to heights above ground for the cap/stem heights. The diagram showing this msl correction procedure is seen in Cloud Msl.dia .

                (b) Cloud Layer Parameters:

Our new sounding input file is now 03060312.pit. We see the 1st line give some info about the bomb size and the 2nd line the plume rise info CapTop/CapBot which is used for the dust inventory programs mentioned earlier. The 3rd line is the header. The pressure is still in mb, heights remained in m above msl, temperature units unchanged as deg C as well as the wind direction and speed still in knots. However, new parameters are layer average temperature, dew points, lapse rate in deg C/100m. Also layer average density and viscosity are used in the fall velocity calculations.

            (2) Pre Processor Programs #2: PTTC.c

                (a) Landing Particle Coordinates:

The name of this program is PTTC.c = Particle Travel Times Coordinates. This is unique for each new sounding due to different profiles of wind direction, speed and temperature. We simply run the fall velocity from the cloud's cap height to the surface for 1 to 449 microns. In order to keep the RAM memory free, we put this to a external file called pttc.cor using the so called binary type file, but is actually a base 256 number of characters instead of a base 10. Reading 03060312.pit into an array of levels, we find the transport and shear components explained above from each level cumulated to the surface. For each level, we must calculate the fall velocity for the 449 particles. The coordinates are stored as the number of .1 mi components east and north.

Since we can only store characters as positive numbers, we set a min east as -100 mi and a min north of -250 miles and we add a correction factor of 500 miles. We then store the cumulative delta east (cde) + 500 east as an example, but broken up into 2-256 characters. We do the same for the north. For the cumulative shear, we add the constant of 500 miles also. When translating these character back in the RADFO program, we must subtract off the 500 miles, then multiply by dxy (.1 mi). We also use 2 characters to identify the total cumulated travel time (tt=dz/fps) for each layer in regard for the half life decay factor described below. The fps is '/sec fall velocity for each particle size id for each wind layer.

Update 10/30/09
The +,- 500 miles was too small. If the winds aloft shift from the north component to the south component, for just a few layers, the smaller particles could return back to mid grid very soon which is incorrect tracking. However, we don't want to extend the grid greater than 400 mi x 500 mi or 4000 x 5000 grid points of .1 mi, as it fairly large enough now. How do we handle this problem? Considering just the east component, we cumulate it by variable cde regardless of the grid limit but when it comes to output, we use another variable cdeo that we truncate to a new larger upper limit that will still be a 2 byte value. This is +32,767 but we use 32,700 (3,270 mi) instead.

The grid S-N is 500 miles or +,- 250 miles, but here we can now accumulate components to +,- 3,270 mi instead in case the wind direction shifts with height. So for say a final component of -3000 miles (south), to keep the output file as + 2 byte value, we now add 3270 = 270. If we had +3000 miles (north) and we add 3270, we have 6270 temporary dummy miles, or 62700 .1 mi. This is still under the limit of 65535 (256^2-1) of a 2 byte value. When reading these 2 byte values in the subsequent programs below, we must subtract off the 32,700, then divide by 10 for the transport component miles of the dust cylinder center. End update 10/30/09

STAPTTCT.c<--->PTTCT.c = trSTApt.c <-- Click 3/07/10

                (b) Surface Cylinder:

We use the mean velocity integration formula, vel.dia, as explained above for the surface layer instead of a linear approach. First, we solve for pex in velz=v0*(z2/zp)^pex where v0 = surface wind @ zp=30' where z2=h[1]-h[0]. Here velz= vel[1], the 1st level above the surface. Taking logs, ln velz = ln [v0*(z2/zp)^pex] = pex ln [v0*(z2/zp)], pex=ln velz/ln [v0*(z2/zp)]. Now use pex in the long velm expression 3.E. Integration of Wind Profile Power Law: in RADFO 1 or p in vel.dia. We then use vel[0]=velm, instead of (V1+V2)/2, times the surface direction dir[0] for the surface translation wind components.

            (3) Pre Processor Programs #3: RADFOPT.c

This program, RADFOPT.c = RADFO Particle Tons, combines cylinder height levels into common wind layers. Again, to cut down on CPU time, we shorten up the number of 449 particle distributions by combining these non-flash multiple levels into a new ton distribution for a given height of a wind layer. This doesn't change much accuracy but gives us less levels to process. Above, we had the file nfhdpd1.NYC explained in 3. P. (2) Non Flash Mode % Dust:. This is the result from the remaining inflow into the nuclear cloud stem after the flash cools down. Notice from inflow.dia again, we have the dry particles start near 52,000' in the cap, distributed to the cap top near 70,000', then down the cap into the stem to near the surface. This is the order they entered the cloud from the updraft.

Each of these levels, along with's it's 449 particle ton distribution, belong within a wind level but is expanded to the whole reportable wind layer. There is no major error here since we still have the same amount of tons, but just expanding it's vertical distribution a little to fill a wind layer. We need to scan up then down to see what wind layers they can be combined into. This seems simple but this was a pesky little program. First, we read the height hc, find what 2 wind levels it lies between, then round it up or down into the appropriate wind layer (l).

We store the count in the winds tons level wtl[L][I]=nfl, where L is the wind level # and I the count # within this layer and nfl is the non flash level (record) # from this tons of dust file being read. We have to read the tons per particle (tonp) for all 449 particles to advance to the next nfl record. Could have used file positions but the number of tons vary and have decimal points so simplier to just read them all.

Update 12/12/09
After all levels are read and stored, we examine the wtl[L][I] array for hc levels or nfl stored there by index I. Now we 1st loop through all the L wind levels again, initialize an array of 449 tons (ton[n]=0) for the combo particle sizes then in the 2nd inner loop, the I's in wtl[L][I] again until the I is reached where we have no more count (wtl[L][I]=0). We reread the nfl file again until the nfl record number matches wtl[L][I]. If so, then the nfl is written to the output combo file with a -, showing the combined last level nfl. For the matched nfl, we add the tonp @ particle id into ton[id]=ton[id]+tonp for id=1,449 for the wind level L. We repeat this again for the next 'I' in wtl[L][I] until we find the last nfl, when we write out all of the 449 ton[id] after it was accumulated for each id.

Now with all the combos done, we write out the non combos as they are transferred from the nfhdpd1.NYC file. We reread this file again, count the nfl records again, continue reading until wtl[L][I]=0 (no nfl's stored) with only the first nfl stored @ I=0. We write out the nfl level, then we read & output consecutive tonp's for each of the 449 particles (id) as explained above.

When to skip reading the hc levels, reading them again and exiting from reading them is controlled by a few logical if statements. As mentioned above, the case that wtl[L][I]=0 means no hc at count nfl has any records. In the first case, when I=0, wtl[L][0]=0, we skip and search for the next wind layer. If I=0 that wtl[L][0]=nfl, where nfl>0, then we have at least one hc level at wind level L. This may or may not be a stand alone. Simple test then is with I=1, that wtl[L][1]=nfl, we have multiple hc's at wind level L, but if wtl[L][1]=0, we have a stand alone that is written out to the output file as +nfl for the first variable in the record, then hc and finally ton[id], id=1,449 or simply just the tonp. Combos are -nfl. Otherwise, we read 449 tonp's to summate ton[id], id=1,449, increment I & nfl then read the next hc record.

We repeat this until wtl[L][I]=0, then output the last -nfl, hc then the summated ton[id], id=1,449. For a given level L, nfl's are in in order of sequence because of nfl increments so the hc file doesn't need to be rewound. However, when wtl[L][I]=0 with I > 0, then we reset nfl=i=0 again and the hc file rewinds to start from the beginning again but L remains constant. Even though nfl increments, the corresponding hc values can be up or down the cap or stem so can appear anywhere.

The finished output file is nf1.NYC (Non Flash 1 Mega for NYC). The original file, (nfhdpd1.NYC) , had 41 original levels. We now have 20 levels, 9 with combos and 11 single stand alones. The 9 combo wind layers are placed within 30 hc layers. Considering that we have tons for each of the 449 particles at each level, this major non flash section cpu time was reduced by 50%. With the new RUC, we have 61 layers and should have nearly the same results depending on the Nuclear Cloud height. As an example, the last section of the wtl[L][I] is seen in radfopt.sam . The wind levels L increase down by rows and the # of hc's [I] to the right by columns with the nfl stored value in wtl[L][I]. Here, we have combined 28 hc levels into 8 leaving 21 stand alones. So we reduced 49 levels to 29 or a (49-29)/49 = 41% reduction. End Update 12/12/09

6. Hierarchy of Grid Systems: 03/02/2009

It may not appear obvious, but we transit back and forth between grid systems before we have a computer graphics format. The first coordinate system in the math standard for sines, cosines, tangets, etc. we will call:

       A. Mathematical Coordinates:

All the trigonometric angle functions and signs is according to Math.dia . The cosines on the right side are + x's, left side -x's while the sines above the x axis are + y's, below -y's. All math operations must be converted to this standard. The mathematicians and navigators never got together because, the math 0 degrees points east. The other point is that the degrees must be converted to radians (57.29578 degrees). The pi factors in the diagram refer to the angle in pi radians, note 2 pi radians = 360 degrees, before using the math sin,cos in the program language.

       B. Map - Wind Direction Coordinates:

At least, we now have wind directions matching N,E,S,W as most maps showing 0 = 360 = North at the top. We consider all Easterly components as +, Westerly as -. Northerly components as +, Southerly as - seen in Map.dia . To convert from map to math, we apply deg=(270-dir)/rad where we use the deg in math, dir is degrees from a map.

       C. Local Wind Direction Coordinates:

Here we plot the downwind direction on the right side of the X axis (easterly) for winds with a westerly component, and to the left (westerly) for winds with a easterly component. We would have crosswinds to the top -, to the bottom + because we want to match the local downwind axis to rotate the centers of the cylinders plots along the downwind x-axis. This is explained above but seen also by WindC.dia . The Counter Clockwise Rotation is -, Clockwise Rotation +. Also given above was the rotation formula x=xg*cos(deg) + yg*sin(deg) and y=yg*cos(deg) - xg*sin(deg) where deg = -45 in our display. The reason for this was to have the center of the dust cylinders (yfc) = 0 or on the x axis so we can calculate the xfc for any grid point Xg,Yg on the circles.

       D. Local and Absolute Grid Coordinates:

Grid Local & Abs.dia shows that we have + xg for East, + yg for North. We could use 0,0 as for a local cylinder plot (top) and an absolute bomb grid of 100 mi to the west of 0,0 as the real 0 for a map system and 300 mi to the east for the eastern end of the absolute grid. We have 250 miles both north and south to complete this absolute grid as seen in the (bottom) diagram. We can work easier with the smaller local cylinder grid when needed rather than working with the larger absolute grid then superimpose these local grids to the absolute by applying a DX,DY rounded off local grid to absolute grid values.

       E. UTM Absolute Grid Coordinates:

Thus, UTM stands for Universal Transverse Mercator Projection coordinate system. It is a rectangular system but because the earth is curved, there are different zones (every 6 degrees Longitude) that we set up another rectangular system again because of the converging of the Longitudes towards the poles. A brief sketch is shown in UTM.dia More info and zones can be seen in UTM wikipedia.org . The transition from our absolute bomb grid to the local UTM zone isn't too difficult. However, not knowing our final based map to display the contour map layer, we will hold off on this.

       F. Polar Grid Coordinates:

This is simply the Latitude - Longitude grid system of the earth showing its curvature. This will be converted from the UTM system above. A crude representation can be seen in Polar.dia . This system can be developed by using a common earth's radius from an imaginary N-S pole through the earth. However, because the earth is elliptical in shape, there is a much more complicated procedure, actually a computer program, that will convert to polar from UTM. Again, not knowing our final based map to display the contour map layer, we will hold off on this. I have a web site that combines Topo Maps to make T-Shirts. To choose which Topos to combine, view samples of this polar plot of 128 topos in this nice regional map which may be used as a general overall base map showing the nuclear effects from multiple city bomb bursts. These base maps are seen in 40N 075W and 39N 077W .

       G. Graphics Plot Coordinates:

This will include the fallout pattern layer and/or any base maps. Simply put, the graphics plot of any compressed bit map is the pixels of a video card - crt monitor. A slight modification of a map, pixels increase from top (yp=0) to bottom (yp=py), # of vertical screen pixels. This is seen in Graphics.dia . Sometimes it's easier to work from the center, xpc,ypc = px/2,py/2.

To keep equal area graphics perspectives, by working with distances, we want dp/dx and dp/dy to be equal where dp = # of pixels/unit x and y and then multiplied by your monitor's aspect ratio given as ar=(swi/pax)/(shi/py) or the factor to change x-pixels distance to y-pixels distance. Here pax,py is total x-pixels, total y-pixels of the monitor screen and swi,shi is screen width, screen height in inches that cover all the pixels displayed.

7. Radfo Grid Plots: 05/10/2009

With all the pre processor programs mentioned above, we test the accumulation of multiple cylinders to make sure the grid is picking them up. Each dust particle cylinder is weighted equally with no decay factor, % probability distribution or particle ton distribution, which will be incorporated in the final RADFO program. The basic program flow is given below with a few added features.

       A. Setting Grid Limits:

The absolute grid's maximum east distance is 300 mi, but starting from a minimum of -100 mi. This choice is because the upper winds are usually westerlies up to say 100,000'. A low pressure circulation may cause some westerly components but unless it's a major closed low system, we hope the westerlies cumulated east components would prevail. As for north and south, we go 250 miles north and south of the east axis because there is more of a chance of this than the easterlies aloft. With the scale resolution at .1 mile, then the grid dimensions are idx=maxde/dxy times idy=maxdn/dxy = 400/.1 x 500/.1 = 4000 x 5000 = 20^6 grid points or a 20 mega-pixel camera. We initialize this grid output file radfo.grd with the 2 null characters '\0', which has a 0 value, per grid point. The result is that the whole grid has 0 values and is fixed in size to a 2 x 20^6 = 40 mb external file. After the program is run, another program will plot panned sections of this grid file. Finally, each grid point value will become a colored pixel making a 20 megabyte uncompressed 256 color picture file.

       B. Radfo Application of Decoded RAOB: RADFOGF.c -> RADFORUG.c

Update 12/08/09
The 03060312.pit file was extracted from RAOBDHM.c but now we use the RUC format RUCDHM.c in our updates mentioned above and explained more in the next topic RadfoGeo . We change the heights (h[n]) to feet and also store in arrays dir[n] & vel[n]. Next, we perform a special surface cylinder wind profile routine. In PTTC.c above, we found the mean integration velocity mean within the surface layer (velm). Here, we do it again, thinking that we would use the # of cylinders (nc) for the surface cylinder up to 5, while we only needed 1 cylinder, top & bottom circles there and decided we only wanted one for here too. Now we use this method again, as before, but apply an algorithm to adjust the surface layer direction by the shear associated with velm. First, we apply surface shear dd=dir[1]-dir[0]. Now if one direction is just west of north, while the other is just east of north, dd > 180, or dd < -180, we correct dd by the -360 or +360 respectively.

Now the proportion (rat) velm is between vel[0] to vel[1], is applied to the wind shear as dirhl=dir[0]+dd*rat. This will be the surface shear components de & dn. Since our first flash particle size is 202u, we find de202 & dn202 by multiplying de & dn by the travel time (tth) in hours. This is obtained from a subroutine called cortt(), where we access the particle coordinate and travel time file pttc.cor created in PTTC.c , is now pttct.mb[m] created in PTTCT.c . We previously found the file position by filp=(id-1)*ln*10+lev*10 where id is particle size, ln = total # of wind layers, and lev = actual wind layer starting at surface= 0. The 10 is the 5 pairs of characters, de,dn trans cor ,tt, ded,dnd shear cor. So with lev=0 for surface, id=220, we find the de,dn translation coordinate for the surface layer. Knowing the surface thickness, z2=h[1]-h[0], we find all the 449 id particle fall velocities pfps[id]=z2/tt for this surface layer in '/sec. However, with our update, we don't need the last 2 pairs (summated vector differences ded,dnd) anymore due to our new method to find cumulated shear-tilt shown above in Vector Shear 2 so we use filp=(id-1)*ln*6+lev*6 instead.

We then call the subroutine sfcpbox() and add the surface shear components de202,dn202 to the cumulated wind shear components cded,cdnd. All subsequent particle (flash: 201->1u & non-flash: 1->449u) and calls to sfcpbox() are multiplied by the fps ratio (fpsrat=fps202/pfps[id]) to expand (pfps[id] < fps202) or contract (pfps[id] > fps202) the 202u trans & shear coordinates to new de,dn because the travel time tt is inversely proportional to fall velocity pfps[id]. The rest of sfcpbox() prepares the surface cylinder box explained above in (h) Rectangle Rotation: and seen in Rotate Rect . Now we can define the max & mins. The downwind xmin in the fig is xmin=-csrad, where csrad is the cap or stem radius set. Now xmax=csrad+dwn, where the downwind dwn=sqrt(de^2+dn^2) distance from center (0,0) of bottom surface cylinder base circle. Now ymin=-csrad and ymax=csrad. In other words, the radius out along the axes from the circle center (ground zero) 0,0 determines the surface rectangle limits except for the downwind dwn which is added to the x direction. The rest of the subroutine sfcpbox() deals with the box rotation explained above in 5. A. (5) (h) Rectangle Rotation:.

            (1) Main Crunching Subroutine:

This grid's xg,yg's are rotated into the downwind box local coordinates and the calculated values for them are in the subroutine chigrid(). The values are determined by the overlapping circles or the expanded rectangle explained above. Now since we have many small probability values and small total tons that may accumulate, we set the minimum value low so not to miss the accumulated effect. We set 3-256 character values for each grid point based upon the value. This gives us a good range with much precision, maybe too much, as 1 -> 256^3 = 16,777,216. This is the value used on your monitor color settings when it states "Millions of Colors or High Color" settings (24 bit). For now, we said above that each original cylinder value is unity before being stretched across the surface while landing instead of using it's actual probability and tons. With our update, we now only need 2-256 chars or 65,536 which is the "Medium Color" settings (16 bit) on your monitor

Above, we explained that ix,iy generated xg,yg's (local grid) to be rotated into the local wind grid (x,y). We have to link dx,dy increments (integer ix,iy) to the absolute grid output file radfo.grd to store the values. We use file positions again to read in, accumulate values and write back this new 2 character chi value. We have a fixed (40 mb) xg,yg grid so file positions is the way. We have an absolute ix,iy as ixa=ixo+ix, iya=iyo-iy. The -iy is because during file reading & plotting, the y-pixel numbers increase from top (0) to bottom pix where we start the iy vertical local grid scan from the minimum y values. Now ixo,iyo is the rounded off coordinate to the nearest .1 value (dxy) of the absolute grid as seen by ixo=(xgo-mine)/dxy and iyo=(maxn-ygo)/dxy. Now xgo,ygo or the # of absolute grid's xg,yg (ix,iy) increments is also at the cylinder center +. So the absolute grid's file position for reading and writing is filp=iya x 2idx+2ixa. As usual, this is best explained by a diagram that combines all the grid coordinate systems in Grids.dia .

            (2) Pan Grid Plots: PANGRID.c

Finally, we see the results from all this by another program PANGRID.c = PAN GRID. It reads the 40 mb radfo.grd file and displays a simple 640 x 480 pixels section of the 4000 x 5000 absolute grid. This resolution could be expanded, with overlapping 1/2 screen heights and widths. It runs searching for values then displays a screen, hit space bar advances to next pic or c to capture the screen display. The first pic is the 202 particle in Flash Only . The next four pics are all the flash & non-flash combined, but with no half life or probabilities and cylinders = 1.0 only. The difference in pics, is the number of cycles of 16 contrasting colors. Here we have, Too Many Color Cycles , then Many Color Cycles also Acceptable Color Cycles ? and Classic Color Cycles . Although not labeled, the dimensions here are the screen width of 640 x dxy=.1 = 64 miles horizontal (west-east) and 480 x dxy=.1 = 48 miles vertical (north-south). It may be best to apply a log scale of color values later.

8. Emission Rates: 05/14/2009

We have to assign a unit of emission from the nuclear blast. A review in air pollution, we have an emission rate in lbs/hr or say g/sec. In a stack, we have a flow rate of m^3/sec. So the concentration in the stack is emission rate/flow rate or (g/sec)/(m^3/sec) or g/m^3. In the diluting atmosphere, we have the same units of flow rate in m^3/sec, but a much larger value. So if all 3 terms stay constant, we have a steady state ambient air concentration in g/m^3 and also in the stack.

For the nuclear, we will use the suggested normalized factor R1 value of 2900 (rads/hr)/1 mi^2 for 1 kiloton. Back to the stack analogy, we didn't have (g/sec)/m^2 but just g/sec. We call this a point source but it really is a small area source. Now we could say per m^2 also, a large stack. That would be a strange unit if we had a tiny stack with just a few inches or cm, which we will do. So with the same emission rate (g/sec)/cm^3, the concentration is the same even though the volume in 1 sec is 1/1000 the 1m^3.

If we doubled the emission rate and kept everything else constant, the concentration would double. The point is, it's a point as we don't mention area when it's small enough, we just call it in g/sec being emitted. This is what they probably meant with this 1 mi^2 as seen by the Cap Radius sizes in Mega. Here, if we double the bomb yield, and the same with the dust area it creates, we would have 5800 (rads/hr)/2 mi^2 = 2900 (rads/hr)/1 mi^2 again. However, keeping the unit area the same, we would have 5800 (rads/hr)/1 mi^2, like the stack instead of the expanding cloud if we keep the flow rates equal. However, since the cloud doubled in size, then we would still have 2900 (rads/hr) over twice the area or twice the total amount of radiation when considering the whole cloud size doubled.

It wouldn't matter which you used if a dust cylinder stretched out to 4 times the area when it landed. So for the original 1 mi^2, we would have 1/4 x 5800 = 1450 but for the 2 mi^2, we have 1/2 x 2900 the same 1450. We will stay with the smaller area and would use 2 x 2900 (rads/hr)/1 mi^2. For our program, we have the expanded rectangle approach for rval in Non Overlapping Cylinder Method where we have d/ds as 1/4 there. For the overlapping circles, we use the calculated root values explained above as seen in Grid 300l.

Since we are spewing out radioactive dust, let's just assume unit weight like tons or g instead of rads/hr. Pretend that there is no rads/hr until it lands and it's just a dust particle. Almost true as the dust takes on the radiation during the flash time, not emits it. However, in the air during post flash, it emits but we are mainly concerned of it's radiation from a cumulated ground dust layer. So the link between g and rads/hr is just some constant of proportionality or rads/hr = k g, a portion of g. So if we double g, the amount of dust from the crater, then k 2g = 2(k g) = 2 rads/hr. Keeping the unit area as 1 mi^2, we have an emission rate Q, as in air pollution, as 2900 (rads/hr)/kiloton = 2900 (rads/hr)/kt x 1000 kt/mt = 29e5 (rads/hr)/mega.

Reviewing this procedure, we have an r1 based on the total rads/hr, name it as temporary, because it is proportional to tons, lbs or grams of dust per area. This is concentration, not in ambient air, but as a layer of dust per mi^2. It is the same concentration or depth over the entire settled cloud source area but the total tons = tons/mi^2 x mi^2 (cloud area). It is based upon the kiloton yield with a proportional radiation rate R1. If we double kiloton, the concentration (tons/mi^2) is the same assuming the cloud area doubled also, so tons/mi^2 x 2mi^2 = 2tons. We have the dust from the crater produced by the 1kt, which we won't calculate here but could, emitting 2900 (rads/hr), by a portion of it's weight or proportionality. In that same proportion, 2 kt spreading out over twice the area still yields the same dust concentration.

However, using our cap radius formula, it isn't actually double. According to the cap radius formula, capr = 10.86 x mega^.4197. Using double the megatons, doubles the overall dust. The ratio of radii for these would be 10.86 x 2^.4197/10.86 x 1^.4197 = 2^.4197. The ratio of areas would be pi x (2^.4197)^2/pi x 1 = 2^(2 x.4197) = 2^.8304 = 1.789 < 2.0. Why is this? The explosion has a vertical component too because of no lid. This would be (2-1.789)/2= 10.5%. Now we must assume that the winds aloft have not acted on the cloud or assume completely calm winds aloft because we want the 2 mega dust to settle back down into a circular area of 1.789 times as much as the 1 megaton cloud area. This means we have a greater (1.12) times as much concentration per mi^2 than the 1 megaton because it didn't spread out to twice the area but only 89.5% of it.

The best way to consider this then is total tons / original cloud cap area. Original means no winds yet to tilt dust cylinders. This is why we don't need an extra term like per nm^2 but treat it as a pseudo stack in 2 stages. Putting this to the stack analogy, at first stage, keep the inner area the same, like a stack, then double the dust by doubling the blast. In the model, we adjust by dividing by the cloud base area later. The 2nd stage is the depositing of the dust buy the spreading out of overlapping particle cylinders as they fall, a type of 2nd dilution technique. For a given bomb yield of 1 megaton, we have 6.1977e+6 tons of crater dust. Proportionately, we have megaton's x 6.1977e+6 tons of dust. Then 1st stage dilution = 1/cap area then the 2nd stage dilution is circle overlap or rectangle stretching by winds, the rval factor. Now we can apply the radiation rate R1 (rads/hr) emitted from the settled dust pile.

       A. Q Application for Dust Particles:

We read a total ton inventory file 'tton1.nyc' produced from another inventory-met pre processor program FNFPARH.c described above. This is given as CraterTons= 6.1977e+006 Building:FlashTons= 7.1443e+006 Non-FlashTons= 1.2843e+008. Now the R1 value of 2900 (rads/hr)/1 mi^2 for 1 kiloton was derived from Nevada test sites, etc. where we only had ground crater dust, not the buildings. So for 1 megaton, we have R1=29e5 (rads/hr) here for the 6.1977e+6 crater dust only. Recall that the building dust flash tons 7.1443e+6 + non flash tons 1.2843e+8 = 1.3557e+8. So divided by the crater dust is 1.3557e+8/6.1977e+6=21.875 times as much building dust than crater dust. This means because of a large city instead of an open plane, it has boosted the 1 megaton bomb to a 20+ megaton when considering fallout in the same areas. Of course, the concrete and asphalt will shield and reflect the crater dust effect a little, but being a few feet deep at the most, it will be broken up and vaporized too.

       B. Final Steps for CHI Values:

            (1) Flash Tons Q Value:

For the FlashTons, we read flhdpd1.nyc also from FNFPARH.c . Here we read the % of activity by decimal form (pdr) for the flash particles 202 -> 1 u. So we apply that times the total flash which is CraterTons 6.1977e+6 + Building FlashTons 7.1443e+6 = 13.342e6. Now divided by the crater tons x R1 x % particle distribution = pdr x 13.342e6/6.1977e+6 x 29e5 (rads/hr) = Q. Finally, we have a Q value for each particle size - now we dilute this over the area of the appropriate cloud cap or stem radius (csrad) then times it's expanded rectangle factor or overlapped circle ratio which are both indicated by the value val (0-1.000). This value is the grid point's (rads/hr) which is now integrated by the Decay Factor Formula explained below.

            (2) Non-Flash Tons Q Value:

For the Non-Flash Tons, we do not have a % particle distribution, by choice, but the tons directly (tons) for each of the 449 particles. So it's emission rate is simply tons divided by the crater tons x R1 = tons/6.1977e+6 x 20e5 (rads/hr) = Q.

9. Half Life Theory: 02/25/2009

As first applied, it's used to define the radioactive decay time for any quantity of a substance to reduce its radioactivity 50% or 1/2 its original activity. The time is usually expressed in hours. Used in carbon dating, and assuming a constant rate of decay for a certain element, we can reverse the formula and solve for the time T.

Its extended use can be applied to air pollution. In modeling SO2 for the City of Phila, EPA suggested that SO2 has a 3 hr half life in urban areas. Actually adjusting the half life, it was determined that it was less. SO2 combining with O2 made SO3 (sulphates) that deposited instead of remaining airborne. We noticed humidity also had an effect and helped reduce the half life to 2 hours.

       A. Half Life Formula:

The simplest expression is Df=.5^nhl where Df is the decimal or % remaining. We define nhl as # of half lives determined by T/hlh, where hlh = half life in hours. Using SO2 as an example, let's try the hlh as 2 hrs, 3 hrs, 5 hrs etc. Now if hlh=3, and we have 3 hrs of time T (travel time), nhl = T/hlh = 3/3=1, Df=.5^nhl = .5^1 =.5. Now if the SO2 travels another 3 hrs aloft, nhl = T/hlh = 6/3=2, Df=.5^2 = .25 or 50% again of the remaining 50%. Also note, if T = 0 hours, nhl = T/hlh = 0/3=0, Df=.5^0 = 1 or 100% of its original.

To extend this for any time T, where nhl (the exponent) is not a whole number, as well as to display the results for hlh=2 & 5 hrs also, see Decay # 1. We have a time T1 for T,
nhl2= # of half lives @ hlh=2 hrs for df2,
nhl3= # of half lives @ hlh=3 hrs for df3,
nhl5= # of half lives @ hlh=5 hrs for df5,
The first to check is the corresponding 2, 3 & 5 hrs where all nhl's = 1.0 and all corresponding df's = .5. For T1=6 hrs, nhl2=3 df2=0.125 & nhl3=2 df3=0.25. So we would have 12.5% or 25% SO2 left for 2 or 3 hr half life respectively. Also note, at T1=10.0 hrs, nhl2=5 df2=0.0313 & nhl5=2 df5=0.25 or a factor of 1/8 using the 2 hr half life over a 5 hr half life.

       B. Other Fractional Life Formulas:

For radioactive decay, if we have Df=.5^nhl, we can define the constant .5 as the dfr or the decay factor rate. Here, we changed the term half-life term dfr to .1 or 10%, since a combination of 1/2 lives can't be represented by a simple Df=.5^nhl, but we have a reduction of 10% each integer power of 7. Wait for loading the PDF file and scroll down to pg. 391 for explanation and pg 392 for the graph in Matches SRFM - Fig 13 Pg 3-13 - US Navy - March 1981.

We now examine the values of the Df=T^-1.2 plot in Decay # 2. Here we note at T1 = 7, 7^2 = 49 and 7^3 = 343 hrs, we have a df2 & df3 of 1.00, .01 and .001 as seen in the bottom table T1 being made multiples of 7 hrs. However, what is this df1 value? It is the plot of Df=T^-1.2 which appears slightly off. So for accuracy, we have a choice. Do we want to keep Df=T^-1.2 or modify it to fit the .1 life every 7^ntl where ntl = # of .1 life. Let's keep the 7^ntl rule and see what the new exponent exp will be. That is, 1/7^1.2 should be 1/10 but 7^1.2 = 10.33 > 10. So df1 = 1/10.33 = .0968. Now in 7^exp=10, ln 7^exp = ln 10, exp ln 7 = ln 10, exp = ln10/ln7 = 1.1833. So our new formula for dfr=1/10 is Df=1/T^1.1833=T^-1.1833 shown in the table by df2. Also notice, we have a column of ntl's where 7^ntl = T. So we calculate ntl by ln 7^ntl = ln T, ntl ln7 = ln T, ntl=lnT/ln7 or the same method used to find the new exp above. Now we can use this as an alternative df3=dfr^ntl=.1^ntl which matches the df2=T^-1.1833 column.

       C. Accumulated Dose:

Suppose we want the accumulated radiation at a grid location Xg,Yg between 2 time intervals T1,T2. Then we would have to integrate Df with respect to time. If we stayed with the radioactive cylinder or it had landed, the decay would decrease with time. Say between 1-3 hrs, it is seen by df2 or df3 in Decay # 2, we decay from df2=1.00 -> df2=0.2725, or an average of .63625/hr. If we sum each df2 for every half hour, then the average is .53398 a better estimate. Then the total dose from the deposited cylinder during 2 hours would be 2 x .63625 = 1.2725 for our 1st estimate but 2 x .53398 = 1.068 for our better 2nd estimate. So if we average the Df over a time interval then take this average times the time interval, we have the accumulated dose. For units, Df and average Df are radiation dose rate or rads/hr. So average Df x hrs = rads/hr x hr = rads - the accumulated dose.

For the mathematical true average of 1 -> 3 hrs, we use the Calculus Integration Formula. This will in effect change the .5 hr interval to .1, .01, .001, .0001 hr or almost the smallest limit (0) hour increments. So we denote this by INT Df. Then the average Df is still INT Df/2 hrs but here it is the exact average for 2 hrs. So then INT Df/2 hrs x 2 hrs = INT Df alone.

Now, if the authors isolated each fission products, which they could have done, but because fractionalization with different particle size distribution may now possibly be done. With the original half life formulas, Df=.5^nhl, or Df=.5^T/hlh. This is the Integral form b^ax dx = b^ax/(a ln b). Here we have b=.5, a=1/hlh, x=T & dx=dT so subbing [(.5^T/hlh)/(1/hlh ln.5)] = [(hlh .5^T/hlh)/ln.5)] lim t1 -> t2 but t1 >=1 hr.

Since the authors combined all complex mixtures of fission products into a decay factor with time T, we will be forced to use it's modified form Df=T^-1.1833. Actually, this is a simpler integral in the basic form x^b, where INT x^b dx = x^(b+1)/(b+1). Here b= -1.1833, x=T, so we have INT Df = T^-1.1833 = -T^-.1833/.1833 lim t1 -> t2 but t1 >=1 hr. So for t1=1 t2=3 hrs, our final Accumulated Dose = INT Df = 0.99506 rads compared to 1.2725 & 1.068 above.

To prove somewhat that the INT Df is correct, we chose smaller increments of t1 -> t2 as .01 hr instead of 1/2 hr. Seen here below:

df=0.875 @ 1.12 hrs rad=0.0087 @ 1.12 -> 1.13 hrs rad/hr=0.870
df=0.865 @ 1.13 hrs rad=0.0086 @ 1.13 -> 1.14 hrs rad/hr=0.861
The average df for the two times 1.12 - 1.13 hrs is df=.870. Notice, this is also the average INT Df/hr shown on the right as rad/hr in INT Df/(t2-t1) = 0.870 for t1=1.12 -> t2=1.13. Therefore, the integrated formula is correct. Again, if we multiply the average INT Df by the time interval, we have INT Df/(t2-t1) x (t2-t1) where the time interval cancels and we are left with the real INT Df (ri), the accumulated dose (rads).

       D. Application to Dust Cylinders:

Each particle's cylinder has to be weighted by its own respective travel time T. A trajectory will give us the path of a particle cylinder from its level aloft all the way to its ground deposit. During this time, the travel time through it original layer and all the subsequent layers to the surface will be accumulated as the total travel time tt or T based on its particle fall velocity. In other models, we may be just interested in the present rads/hr at a given time rather than the Accumulated Dose INT Df. For our absolute grid, explained below, we will have an Xg,Yg coordinate with its deposited dust mass.

It's possible that even though we have the final deposited dust, it may be that for our time of interest, say t1, travel time, tt or T > t1 meaning the dust particle cylinder hasn't settled yet and is still airborne for a particular region of Xg's, Yg's. Therefore, it can't be included in the radiation calculation. Also, for Accumulated Doses, if T > t1, t1 is set to T since at t1, the radiation t1 has occurred upwind somewhere else. Examples for cylinder travel time tt or T is seen below:

Ex 1: tt = T = 2 hrs, Df @ t1=3 hrs so Df for cylinder t1=3 still used.
Ex 2: tt = T = 4 hrs, Df @ t1=3 hrs so Df for cylinder t1=3 can't be used.
Ex 3: tt = T = 2 hrs, INT t1=3 hrs -> t2=6 hrs, so t1=3, t2=6 still used.
Ex 4: tt = T = 4 hrs, INT t1=3 hrs -> t2=6 hrs, so t1=4, t2=6 now used.

Finally, we use Real Integral Dose Rate (ri) = -T^-.1833/.1833 lim t1 -> t2 where t1 = tt, t2 = time after burst as 12 hours? This would be ri=[t2^(exp+1)-t1^(exp+1)]/(exp+1) where exp=-1.1833. So the final cumulated rads at a grid point is chi=Q x ri for 12 hrs, relating this chi to radiation health effects.

10. Final Radfo Plots: GRIDPCX.c - GRIDPCX.c 05/18/2009

Using PANGRID.c , we have hopefully included all the parameters above and present another 64 x 48 mile section of the complete fallout in rads. This is seen first in Radfo #4. We see that the few cycles of 16 colors were used and near the center line. The increase is so rapid that we have no detail. For Radfo #5, we switched to a 2 color cycle log scale which is better but shows the details of the rugged edges. Notice, even though we chose a wild wind profile Pittsburgh, PA, it shows in detail what may have been missed with more uniform wind vectors. For example, the transition from circle overlap to extended rectangle deposits is clearly seen in the wind shifts.

We will now show the complete plot, without coordinates and color value labels for now using GRIDPCX.c Plots . See Radfo #6. Click the [Download] tab above the pic, put in your chosen directory - it is safe. View in your favorite graphics software that uses zoom & pan for details. If not, use Start All Programs ==> Accessories > Paint. Move the sidebar's until you can find the radioactive plume. After changing the color sequence and adding the half life decay factor, we have a newer display as seen in Radfo #7.

       A. Framing by Labels on GRIDPCX.c Plots : 06/19/2009

As you can see by the dates, it took a month to add the labels to the graph. This would never take that long with a graphic artists using a graphics art program. However, we developed programs for each of these labels then programs that combine these labels. This way, the graphics are added automatically within seconds when running these programs in a batch file which will be added to the other batch file programs. Now, no other graphics art programs are required. The 2 sides and the bottom labels are fixed and don't change after combining with axes but the Title.txt file will be updated with each run. The date in the title is given by mo-day-year and times in EST (Eastern Standard Time) rather than Zulu time, which the public in the US would have to convert, so EST is better for the US. The flow diagram of these programs is seen in Labels.fld. The spacing for the 50 mile markers and legend values are depicted in xy50mm.dia.

All the real time fallout plots are put in their respective folders in the MSN Sky Drive. Here we put the final sample in the NYC folder NYC06030307. Suggestion is not to print this but only view it on the screen. For printing, use it's color reverse color palette NYC06030307R. The reason being, your printer will run out of black ink or toner soon so we switch to a white paper background plot. One note here, the S96 in the site's file name represents that both original pics were reduced to 96% of the original size in order for MSN to post it.

Note the cross hairs (+) at ground zero that match the X & Y labels - used for a coordinate check of the plot but decided to keep. Also, the plot concentrations have been doubled to show the full range of the legend values, showing the egg yolk maximum. However, if we doubled the 1 megaton, the plot would be similar but displaced a little upwind. Another simulation for doubling would be halfing the wind speed instead, causing an upwind shift in double concentration. Notice the maximum fallout concentration is just over 100 miles east and 40 miles north of NYC. Instead of referring back to here during real time operation, we will present the health effects table from accumulated RADS in each city folder.

       B. Single Contour Plots:

We could superimpose major highways on top of these plots or we can overlay single contour templates on top of road, state or topo maps. Actually it was too difficult to write a program that would show a continuous 1 pixel contour line thickness without going to a graphics art program. However, we modified the solid contour program somewhat for less colors but variable thickness contours that can be seen in Radfo Singlular Contours that can be improved later.

        C. Time Plots: PPTTC.c 07/31/2009

A request by the pioneers of the early models suggest a 1st approximation of the plume before the detailed model run. Since it appears that the main model run RADFOGF.c will take a minimum of 75 minutes, the input to the model can also be used to time radioactive fallout plots. We will be plotting travel time rather than radioactive dust concentration. However, since this avoids the probability, half life decay integration and overlap-stretching of cylinder areas probabilities, this computer run is immediate and uses the output file pttc.cor from the program PTTC.c explained above. The purpose is for the traffic to stay ahead of or turn away from the earliest fallout time lines to escape any radiation.

            (1) Output Description:

Here we assume each of 449 particles are distributed in each of the multiple wind layers to the height and diameter of the Nuclear Cloud Cap only. The stem will have particles fall in closer but we want the most outward extent of the radiation. Now since the binary file pttc.cor means Particle Travel Time & Coordinates, it has all this info. We plot the minimum time it takes each particle to fall on the grid. This minimum means that if we have 2 particles land at the same grid point, we show the colored time plot of the particle that lands first. As an example, say we have a smaller particle at 40000 feet that takes 3 1/4 hours to land at a grid point. We have another much larger particle at 60000 feet. This particle falls faster and even though it falls another 20000 ft, it landed at the grid point in 2 3/4 hours or 1/2 hour earlier. The associated color plot will be with the larger particle and plotted at the grid location with the 2-3 hr color range.

            (2) Sample Plots:

                (a) 15-Minute Segments:

The first plot is shown in timev4an.gif. The legend is given in the upper left corner. Here we have contrasting color intervals by hours with the alpha-numeric characters representing 15 min sub divisions of the same color 1-hr range. We use 0-9, 15 minute segments up to 2 1/2 hours, then we use the alphabet (capitals) until all the letters are used. Then at 9 hours, we use the lower case alphabet up to 15 hours. However, in the real concentration plot, we have the half life for only 12-hr dose radiation. The gray circle represents the cloud cap diameter and the + is ground zero. The particles are taken from this outer ring as it has already advanced to the cloud cap diameter before the wind transports it. This plot has a strange wind profile and shows a gap which really may not exist. For a directional wind shear with height or a wind shift at the surface, it must pass between it's starting and ending point but won't show here because we use the discrete mean wind vector layers. Therefore, care should be taken if connecting like color values and should show no gap when joining hourly plots in those two distinct regions.

                (b) 1-hr Line Rays:

The second plot is shown in timev4px.gif. The same legend is given in the upper left corner. We kept the colored 4 alpha-numeric in case the rays would obliterate the hour's color. Also, it's a reference if you want to revert back to the 15 min alpha-numeric. Notice in this plot, we use the pixels rather than the 8 x 16 pixels that the alpha-numeric's require. We don't have the 15 min subdivisions but eyeball estimates along a color ray should suffice and may even be a better estimate than rounding off the 8 x 16 pixel block.

If you compare (a) and (b), you will notice the rays are longer. That shouldn't be for a given time value. Surprisingly, this is a flaw by MSN's in it's Visual C++ software that came back to haunt me again after 10 years. For the maximum pixel resolution, they failed to assign the correct number of columns. I'm sure they corrected this with their new Windows controlled Dynamic Link Library programs but I requested a patch to correct for this DOS executable program but no word. Therefore, we will probably use this pixel pic instead, which is referred to as the colored broom sweep. These plots have no axis labels or scales. The reason being is that the pic will hopefully be overlaid on a Google Earth map and viewed similar to the Weather Channel's Radar seen on the internet. Therefore, their scale is built into their maps but more important, the towns and roads are also visible. However, the scale here is the same as the main 12 hr dose radiation plot but condensed. The left border is -100 miles west from ground zero, the right border 300 miles east. The top is 250 miles north, bottom 250 miles south.

This program could be modified to show hourly time lapse advances of the leading edge fallout plume just as the advancing various color radar echoes are displayed on the internet and TV weather. In this later case, the legend would need not be shown as a time clock projection box could be seen as in animated radar.

11. RadfoGeo Radfo-Google Earth Overlays: