Final Project: Solar Corona
Daniel Chang [dhchang at stanford dot edu], Ming Jiang [mingj at stanford dot edu]
We wanted to model the surface of the sun exhibiting coronal loops and coronal mass ejection. Coronal loops are results of a magnetic effect where loops of magnetic flux rise above the surface and fill with plasma. Coronal mass ejections are massive bursts on the surface of the sun which eject plasma particles from the solar corona.
Figure 1: reference image, 821 x 691 (credit: NASA)
Our solution is divided into rendering the sun's surface and coronal loops.
Sun's Surface (Ming)
In reality, the sun's outer part consists of the photosphere, chromosphere and corona. All these three parts are filled with emissive gas and plasma. The photosphere is the most bright and dense part. So it appears as a solid ball in the center. The chromosphere and corona are much thinner than the photosphere, so they act like the atmosphere of the sun. The chromosphere is relatively shallow, and separates the corona from the photosphere. Therefore, the concept of the sun's surface equals photosphere plus chromosphere in our project.
We model the photosphere as a uniformly emissive perfect sphere. This is simply done by declaring the following area light source:
AreaLightSource "area" "color L" [100 27 0] "integer nsamples"  Translate -55 -355 0 Shape "sphere" "float radius" 
The chromosphere should be very shallow and much less emissive than the photosphere. In the reference image, we can see that it is almost like dark clouds covering the lighting ball. So we extend PBRT to provide a new volume region class to model the chromosphere.
CloudRegion, the new volume region class, is inherited from DensityRegion. The main extension is to overwrite Density() function to return a procedurally generated density value. Also, its region is defined as the space between two homocentric spheres surrounding an inner spherical object which is our photosphere.
There are totally 11 new parameters added for CloudRegion class:
Shape-related parameters: innerradius, outerradius
- Both of them are of floating-point type. The chromosphere fills the space between two homocentric spheres defined by innerradius and outerradius.
Turbulence-related parameters: roughness, octaves, scale, shift
The basic method in our procedural generation is Perlin Noise. PBRT has a built-in function called Turbulence() which implements the basic perlin noise with parameters "roughness" and "octaves" to control the quality of the noise. We added another two parameters: [float]scale, [Vector]shift to control the size of the perlin noise granules, because the built-in Turbulence() always return the value based on the current position's coordinates. So in the actual fact "scale" and "shift" only affect the geometry transformation and then indirectly the granule size.
Here are the comparison of three pictures with different sizes of the granule.
Figure 2: median, large and small sized granule (from left to right)
Cutoff-related parameters: coverage, sharpness, amplitude
It is far from cloud-like if we only use perlin noise. According to an article introducing cloud cover, we applied an exponential cutoff transformation to make it more cloud-like. The transformation function is basically as follows:
F(x) = amplitude * (1 - exp(-(x - coverage) * sharpness))
By "basically" we mean that F(x) will be forced to equal zero if the calculated result is negative.
The following two picture illustrates the difference between no cutoff and after cutoff transformation.
Figure 3: without and with cutoff transformation
Texture-related parameters: biasmap, peak
After applying turbulence and cutoff transformation, we get the following picture:
Figure 4: rendering without texture map
However, the biggest problem with this rendering is that it is unable to simulate the reference picture and thus match the positions of all its coronal loops. To solve this problem, we added support for an external texture map to specify the surface's energy distribution, that is, to add offset to the already procedurally generated noise. Here the "biasmap" parameter specifies such a texture, and "peak" indicates the scale and thus numerically equals the largest offset amount. To add such support, we must modify some of the PBRT's core files as well, because the original VolumeRegion only accepts ParamSet and is unable to deal with textures. So we modified api.cpp, dynload.h and dynload.cpp to enable accepting TextureParam for this new class while keeping compatible with the normal volume region classes.
Here is the texture map used in our final rendering.
Therefore, simply put, our final density return value is: texture offset + procedurally generated noise.
Coronal Loops (Dan)
Coronal loops are resultant from the twisted solar magnetic flux within the solar body , and are a complex physical formation due to the magnetic polarity and differential rotation of the sun. Because modeling the magnetic field acting behind the loop formations was beyond my knowledge of physics, I chose a phenomenological approach to modeling the shapes.
The basic idea is to procedurally generate many such curves so that the general ordered shape of the entire loop body was preserved, but providing enough randomness to make the loops look natural. I chose Bezier curves for the shape of each individual strand of the coronal loops because it gave enough control over the shape of the loop (with 4 knots / control points). The parametric equation for Bezier curves is given as follows:
B(t) = (1 - t)3 P0 + 3(1 - t)2(t) P1 + 3(1-t)(t2) P2 + (t3)P3 for t E (0, 1)
P0 and P3 form the base of the curve, and protrude from the sun's surface. P1 and P2 control the shape of the curve. Because the curve is in 3-space, P1 and P2 can also be moved opposite to each other in Z to provide extra complexity to the curves.
I implemented my volume region generator outside of PBRT, and made the output of the densities compatible with PBRT's volumegrid format. I chose to do so because it was crucial to visualize the results of the generated densities, and using OpenGL in a standalone program would make this a lot easier. For every curve, I traversed every voxel in the volumegrid and found its shortest distance to the curve, then determined if the shortest distance was within the curve's density falloff to receive density contributions from the curve. However, solving the shortest distance analytically for a quartic equation proved unfeasible. For those interested, refer to  for the analytical solution to the quartic formula.
Instead, I estimate the shortest distance from a voxel to the curve by stepping t from 0 to 1, and for each step, finding its distance to the voxel. Undersampling resulted in blocky curves with noticeable aliasing, so t needed to be stepped at around .0025 intervals (400 steps from 0 to 1). With 300x300x300 voxels, 400 curves, and 400 steps, this would require 4.3 * 1012 comparisons, each one solving for a distance between two 3-space points, which has square roots and is thus slow. Thus optimization was needed in finding the shortest distance.
The first optimization was to create bounding boxes for each curve. Voxels outside the bounding box were not compared, which could reduce 300x300x300 voxel comparisons to, say, 290x250x50 voxels, which is a significant improvement.
The second optimization (illustrated above) was to iteratively split the curve into segments, find the segment which had the shortest distance to the point, then split that segment into segments and repeat. Doing this 3 times, the number of comparisons could be reduced from 400 to about 5*cuberoot(400). For instance, we split the curve initially into 14 segments. At each split point, we find the distance to the point. The two segments around the split point with shortest distance are then split into 14 segments, etc.
Pictured (above left) are the first results of the coronal loops modeled from Bezier curves. The left image has the most basic bezier curves in the volume region. I realized after forming the first curves that to model the thin "threads" of plasma that comprise the coronal loop, a lot more had to be done. First of all, each "thread" would need to be thinner and lighter. Second of all, a striking property of plasma that makes it so visually distinct is its "filamentation" property (above right). Plasma is highly conductive and contains many free electrons, and as charged particles rapidly move in plasma, "a ring of magnetic fields forms around the current that can pinch it into filamentary current strands" . From the original reference picture, we can clearly see that the coronal loops have a magnificent amount of filamentation. To simulate this phenomenon, I chose to use many thin strands of bezier curves, and give them initial parameters that will cause random bunching near the top of the coronal loops. The second image above illustrates how this creates the effect of filamentation (from ).
Because the loops are thin layers of plasma, the driving factor behind their appearance is emission. Plasma contains many free and highly excited electrons, which can fall back into equilibrium state, releasing photons as a result of dropping in energy level. Another mechanism is stimulated emission where speeding electrons collide with ions, exciting its electrons which then release more photons when they fall back in state. We suspected the reference image was taken in a spectrum other than the visible light spectrum, with special filters, otherwise the brightness of the surface of the sun should saturate the film. As such, inter-reflection between coronal loops and the sun's surface was ignored due to 1) coronal loops are mostly emissive ionized gas, 2) the effect of interreflection on the emissive surface of the sun is negligible in comparison to the photosphere's own emissivity.
However, the nature of the plasma is not uniformly distributed over the loop. As evinced from the reference picture, the top or bottom of some macro-loops appear much brighter than the rest. Reference  explains that coronal loops have varying degrees of brightness along the loop based on a random sequence of macroflare heating:
This causes the uneven brightness - and brilliant highlights - on the loop bodies, often near the top or bottom of the loop. To model this phenomenon, I mapped density to emissivity, and used falloff functions based on vertical height on the loop to procedurally generate the densities. Example falloffs are linear grow per height, linear fade, logarithmic, and inverse squared. Pictured below are the effects of linear falloff in the openGL visualizer, one with growing emissivity and the other fading:
To generate the translucent, high-contrast brilliance of the loops, some post-processing on the loop densities were performed. An important operation was blurring the densities to create the "glow" effect around the loops, likely due to a thinner periphery of plasma than the bright center of the loop. Operations such as blur and sharpen were implemented in 3-space with the following masks (each mask is 5x5x5, with cross sections of the rows displayed. Rows 0, 4 and 1, 2 are symmetrical and shown once):
Combining the brilliance falloff with some post processing, major improvements were made to the first draft of the loops. You can see that the picture at right appears more translucent and looks more similar to the coronal loops than the first implementation:
Loop bodies were broken into macro-loops, each of which can hold hundreds of threads. The shapes were procedurally generated with parameters such as defining the loop center and the tilt from center, dtilt/ddistance from center, curve height, width, sharpness, density falloff function, height falloff function, base density, etc. This gave intuitive control over each macro loop, while procedurally generating the macrostructure based off of the starting parameters with random noise in its shape (noise coefficients such as random displacement from root, random height, random tilt also parameters) to make it look natural.
The sun has its own atmosphere of plasma, known as the corona, which extends millions of kilometers into space. To generate the inhomogeneous volume for the corona, I loosely based the implementation off of Ken Perlin's Hypertexture paper and Image Synthesizer papers [7, 8], but with a lesser complexity. Perlin used random offsets to starting voxels on a sphere to generate coronas, tribbles, and other inhomogeneous volumes. However, the spline gradients at each voxel which he interpolates need to be predetermined for the volume to look correct. Instead of using spline knots at each voxel, I gave my growing volume initial gradients and random offsets to slope, position, and density, and approximated each as lines. From the reference picture, we can see that the atmosphere seems to extend linearly (or at least with low change in gradient) into space in bundles, likely also due to filamentation on a grander scale (much like the aurora phenomenon). Additionally, from a reference picture of the coronal atmosphere, it is appropriate to approximate the bundles of plasma as lines at this camera position (as the bundles extend millions of kilometers into space with relatively low curvature) (image from ):
A more diffuse density representation of the corona was combined with sharper "bundles" of filamentation to create an effect similar to "god rays" seen under water. Except instead of volumetric caustics, these are caused by plasma filamentation extending from the sun way out into space. The high density of coronal plasma near the photosphere, seen through a grazing angle, was also procedurally generated and added to the horizon of the sun, seen as the fiery streaks off in the horizon. Emissive stars were modeled as sphere densities and randomly placed in the background. The final result looks as follows:
Division of Work
Dan worked on the coronal loops and coronal atmosphere, Ming made the sun's surface.
Descriptions: CoronalDensities.rar contains the pbrt render file, coronal.pbrt, and the coronal loop densities. Very large when unpacked (~600 MB). VolumeGenerationCode contains the .cpp files that create the coronal loops, atmospheric lines, stars, etc.
The following package contains all the pbrt files that are modified as well as the extended class file and its Makefile:
 Perlin Noise http://freespace.virgin.net/hugo.elias/models/m_perlin.htm
 Cloud Cover http://freespace.virgin.net/hugo.elias/models/m_clouds.htm
 F. Reale, Physics of Coronal Loops. Memorie della Societa Astronomia Italiana, Vol. 69, p.669. Viewed 6/8/2009. http://adsabs.harvard.edu/abs/1998MmSAI..69..669R
 K. Perlin, Hypertexture. ACM SIGGRAPH Computer Graphics, Volume 23, Issue 3 (July 1989). http://portal.acm.org/citation.cfm?id=74359&coll=GUIDE&dl=GUIDE&CFID=36533152&CFTOKEN=45777804&ret=1#Fulltext
 K. Perlin, An Image Synthesizer. ACM SIGGRAPH Computer Graphics, Volume 19, Issue 3 (July 1985). http://portal.acm.org/citation.cfm?id=325247