simgear/scene/sky/clouds3d/SkyCloud.cpp

   1 //------------------------------------------------------------------------------
   2 // File : SkyCloud.cpp
   3 //------------------------------------------------------------------------------
   4 // SkyWorks : Adapted from skyworks program writen by Mark J. Harris and
   5 //                                              The University of North Carolina at Chapel Hill
   6 //                                      : by J. Wojnaroski Sep 2002
   7 //------------------------------------------------------------------------------
   8 // Permission to use, copy, modify, distribute and sell this software and its
   9 // documentation for any purpose is hereby granted without fee, provided that
  10 // the above copyright notice appear in all copies and that both that copyright
  11 // notice and this permission notice appear in supporting documentation.
  12 // Binaries may be compiled with this software without any royalties or
  13 // restrictions.
  14 //
  15 // The author(s) and The University of North Carolina at Chapel Hill make no
  16 // representations about the suitability of this software for any purpose.
  17 // It is provided "as is" without express or
  18 // implied warranty.
  19 /**
  20  * @file SkyCloud.cpp
  21  *
  22  * Implementation of class SkyCloud.
  23  */
  24
  25 // warning for truncation of template name for browse info
  26 #pragma warning( disable : 4786)
  27
  28 #include <plib/ul.h>
  29
  30 #include "SkyCloud.hpp"
  31 #include "SkyRenderableInstance.hpp"
  32 #include "SkyContext.hpp"
  33 #include "SkyMaterial.hpp"
  34 #include "SkyLight.hpp"
  35 #include "SkyTextureManager.hpp"
  36 #include "SkySceneManager.hpp"
  37 #include <algorithm>
  38
  39 //! The version used for cloud archive files.
  40 #define CLOUD_ARCHIVE_VERSION 0.1f
  41
  42 //------------------------------------------------------------------------------
  43 // Static initialization
  44 //------------------------------------------------------------------------------
  45 SkyMaterial*  SkyCloud::s_pMaterial         = NULL;
  46 SkyMaterial*  SkyCloud::s_pShadeMaterial    = NULL;
  47 unsigned int  SkyCloud::s_iShadeResolution  = 32;
  48 float         SkyCloud::s_rAlbedo           = 0.9f;
  49 float         SkyCloud::s_rExtinction       = 80.0f;
  50 float         SkyCloud::s_rTransparency     = exp(-s_rExtinction);
  51 float         SkyCloud::s_rScatterFactor    = s_rAlbedo * s_rExtinction * SKY_INV_4PI;
  52 float         SkyCloud::s_rSortAngleErrorTolerance = 0.8f;
  53 float         SkyCloud::s_rSortSquareDistanceTolerance = 100;
  54
  55 //------------------------------------------------------------------------------
  56 // Function               : SkyCloud::SkyCloud
  57 // Description      :
  58 //------------------------------------------------------------------------------
  59 /**
  60  * @fn SkyCloud::SkyCloud()
  61  * @brief Constructor.
  62  */
  63 SkyCloud::SkyCloud()
  64 : SkyRenderable(),
  65   _bUsePhaseFunction(true),
  66   _vecLastSortViewDir(Vec3f::ZERO),
  67   _vecLastSortCamPos(Vec3f::ZERO)
  68 {
  69   if (!s_pShadeMaterial)
  70   {
  71     s_pShadeMaterial = new SkyMaterial;
  72     s_pShadeMaterial->SetAmbient(Vec4f(0.1f, 0.1f, 0.1f, 1));
  73     s_pShadeMaterial->EnableDepthTest(false);
  74     s_pShadeMaterial->SetBlendFunc(GL_ONE, GL_ONE_MINUS_SRC_ALPHA);
  75     s_pShadeMaterial->EnableBlending(true);
  76     s_pShadeMaterial->SetAlphaFunc(GL_GREATER);
  77     s_pShadeMaterial->SetAlphaRef(0);
  78     s_pShadeMaterial->EnableAlphaTest(true);
  79     s_pShadeMaterial->SetColorMaterialMode(GL_DIFFUSE);
  80     s_pShadeMaterial->EnableColorMaterial(true);
  81     s_pShadeMaterial->EnableLighting(false);
  82     s_pShadeMaterial->SetTextureApplicationMode(GL_MODULATE);
  83   }
  84   if (!s_pMaterial)
  85   {
  86     s_pMaterial = new SkyMaterial;
  87     s_pMaterial->SetAmbient(Vec4f(0.3f, 0.3f, 0.3f, 1));
  88     s_pMaterial->SetDepthMask(false);
  89     s_pMaterial->SetBlendFunc(GL_ONE, GL_ONE_MINUS_SRC_ALPHA);
  90     s_pMaterial->EnableBlending(true);
  91     s_pMaterial->SetAlphaFunc(GL_GREATER);
  92     s_pMaterial->SetAlphaRef(0);
  93     s_pMaterial->EnableAlphaTest(true);
  94     s_pMaterial->SetColorMaterialMode(GL_DIFFUSE);
  95     s_pMaterial->EnableColorMaterial(true);
  96     s_pMaterial->EnableLighting(false);
  97     s_pMaterial->SetTextureApplicationMode(GL_MODULATE);
  98     _CreateSplatTexture(32); // will assign the texture to both static materials
  99   }
 100 }
 101
 102
 103 //------------------------------------------------------------------------------
 104 // Function               : SkyCloud::~SkyCloud
 105 // Description      :
 106 //------------------------------------------------------------------------------
 107 /**
 108  * @fn SkyCloud::~SkyCloud()
 109  * @brief Destructor.
 110  */
 111 SkyCloud::~SkyCloud()
 112 {
 113 }
 114
 115
 116 //------------------------------------------------------------------------------
 117 // Function               : SkyCloud::Update
 118 // Description      :
 119 //------------------------------------------------------------------------------
 120 /**
 121 * @fn SkyCloud::Update(const Camera &cam, SkyRenderableInstance* pInstance)
 122 * @brief Currently does nothing.
 123 */
 124 SKYRESULT SkyCloud::Update(const Camera &cam, SkyRenderableInstance* pInstance)
 125 {
 126   return SKYRESULT_OK;
 127 }
 128
 129
 130
 131 //------------------------------------------------------------------------------
 132 // Function               : DrawQuad
 133 // Description      :
 134 //------------------------------------------------------------------------------
 135 /**
 136  * @fn DrawQuad(Vec3f pos, Vec3f x, Vec3f y, Vec4f color)
 137  * @brief Draw a quad.
 138  */
 139 inline void DrawQuad(Vec3f pos, Vec3f x, Vec3f y, Vec4f color)
 140 {
 141   glColor4fv(&(color.x));
 142   Vec3f left  = pos;  left   -= y;
 143   Vec3f right = left; right  += x;
 144   left  -= x;
 145   glTexCoord2f(0, 0); glVertex3fv(&(left.x));
 146   glTexCoord2f(1, 0); glVertex3fv(&(right.x));
 147   left  += y;  left  += y;
 148   right += y;  right += y;
 149   glTexCoord2f(1, 1); glVertex3fv(&(right.x));
 150   glTexCoord2f(0, 1); glVertex3fv(&(left.x));
 151 }
 152
 153
 154 //------------------------------------------------------------------------------
 155 // Function               : SkyCloud::Display
 156 // Description      :
 157 //------------------------------------------------------------------------------
 158 /**
 159  * @fn SkyCloud::Display(const Camera &camera, SkyRenderableInstance *pInstance)
 160  * @brief Renders the cloud.
 161  *
 162  * The cloud is rendered by splatting the particles from back to front with respect
 163  * to @a camera.  Since instances of clouds each have their own particles, which
 164  * are pre-transformed into world space, @a pInstance is not used.
 165  *
 166  * An alternative method is to store the particles untransformed, and transform the
 167  * camera and light into cloud space for rendering.  This is more complicated,
 168  * and not as straightforward.  Since I have to store the particles with each instance
 169  * anyway, I decided to pre-transform them instead.
 170  */
 171 SKYRESULT SkyCloud::Display(const Camera &camera, SkyRenderableInstance *pInstance)
 172 {
 173   // copy the current camera
 174   Camera cam(camera);
 175
 176   // This cosine computation, along with the if() below, are an optimization.  The goal
 177   // is to avoid sorting when it will make no visual difference.  This will be true when the
 178   // cloud particles are almost sorted for the current viewpoint.  This is the case most of the
 179   // time, since the viewpoint does not move very far in a single frame.  Each time we sort,
 180   // we cache the current view direction.  Then, each time the cloud is displayed, if the
 181   // current view direction is very close to the current view direction (dot product is nearly 1)
 182   // then we do not resort the particles.
 183   float rCosAngleSinceLastSort =
 184       _vecLastSortViewDir * cam.ViewDir(); // dot product
 185
 186   float rSquareDistanceSinceLastSort =
 187     (cam.Orig - _vecLastSortCamPos).LengthSqr();
 188
 189   if (rCosAngleSinceLastSort < s_rSortAngleErrorTolerance ||
 190       rSquareDistanceSinceLastSort > s_rSortSquareDistanceTolerance)
 191   {
 192     // compute the sort position for particles.
 193     // don't just use the camera position -- if it is too far away from the cloud, then
 194     // precision limitations may cause the STL sort to hang.  Instead, put the sort position
 195     // just outside the bounding sphere of the cloud in the direction of the camera.
 196     _vecSortPos = -cam.ViewDir();
 197     _vecSortPos *= (1.1 * _boundingBox.GetRadius());
 198     _vecSortPos += _boundingBox.GetCenter();
 199
 200     // sort the particles from back to front wrt the camera position.
 201     _SortParticles(cam.ViewDir(), _vecSortPos, SKY_CLOUD_SORT_TOWARD);
 202
 203     //_vecLastSortViewDir = GLVU::GetCurrent()->GetCurrentCam()->ViewDir();
 204     //_vecLastSortCamPos = GLVU::GetCurrent()->GetCurrentCam()->Orig;
 205     _vecLastSortViewDir = cam.ViewDir();
 206     _vecLastSortCamPos = cam.Orig;
 207   }
 208
 209   // set the material state / properties that clouds use for rendering:
 210   // Enables blending, with blend func (ONE, ONE_MINUS_SRC_ALPHA).
 211   // Enables alpha test to discard completely transparent fragments.
 212   // Disables depth test.
 213   // Enables texturing, with modulation, and the texture set to the shared splat texture.
 214   s_pMaterial->Activate();
 215
 216   Vec4f color;
 217   Vec3f eyeDir;
 218
 219   // Draw the particles using immediate mode.
 220   glBegin(GL_QUADS);
 221
 222   int i = 0;
 223   for (ParticleIterator iter = _particles.begin(); iter != _particles.end(); iter++)
 224   {
 225     i++;
 226     SkyCloudParticle *p = *iter;
 227
 228     // Start with ambient light
 229     color = p->GetBaseColor();
 230
 231     if (_bUsePhaseFunction) // use the phase function for anisotropic scattering.
 232     {
 233       eyeDir = cam.Orig;
 234       eyeDir -= p->GetPosition();
 235       eyeDir.Normalize();
 236       float pf;
 237
 238       // add the color contribution to this particle from each light source, modulated by
 239       // the phase function.  See _PhaseFunction() documentation for details.
 240       for (int i = 0; i < p->GetNumLitColors(); i++)
 241       {
 242         pf = _PhaseFunction(_lightDirections[i], eyeDir);
 243         // expand this to avoid temporary vector creation in the inner loop
 244         color.x += p->GetLitColor(i).x * pf;
 245         color.y += p->GetLitColor(i).y * pf;
 246         color.z += p->GetLitColor(i).z * pf;
 247       }
 248     }
 249     else  // just use isotropic scattering instead.
 250     {
 251       for (int i = 0; i < (*iter)->GetNumLitColors(); ++i)
 252       {
 253         color += p->GetLitColor(i);
 254       }
 255     }
 256
 257     // Set the transparency independently of the colors
 258     color.w = 1 - s_rTransparency;
 259
 260     // draw the particle as a textured billboard.
 261     DrawQuad((*iter)->GetPosition(), cam.X * p->GetRadius(), cam.Y * p->GetRadius(), color);
 262   }
 263   glEnd();
 264
 265   return SKYRESULT_OK;
 266 }
 267
 268
 269 //------------------------------------------------------------------------------
 270 // Function               : SkyCloud::DisplaySplit
 271 // Description      :
 272 //------------------------------------------------------------------------------
 273 /**
 274  * @fn SkyCloud::DisplaySplit(const Camera &camera, const Vec3f &vecSplitPoint, bool bBackHalf, SkyRenderableInstance *pInstance)
 275  * @brief The same as Display(), except it displays only the particles in front of or behind the split point.
 276  *
 277  * This is used to render clouds into two impostor images for displaying clouds that contain objects.
 278  *
 279  * @see SkyRenderableInstanceCloud
 280  */
 281 SKYRESULT SkyCloud::DisplaySplit(const Camera           &camera,
 282                                  const Vec3f            &vecSplitPoint,
 283                                  bool                   bBackHalf,
 284                                  SkyRenderableInstance  *pInstance /* = NULL */)
 285 {
 286   // copy the current camera
 287   Camera cam(camera);
 288
 289   Vec3f vecCloudSpaceSplit = vecSplitPoint;
 290
 291   if (bBackHalf) // only sort when rendering the back half.  Reuse sort for front half.
 292   {
 293     // compute the sort position for particles.
 294     // don't just use the camera position -- if it is too far away from the cloud, then
 295     // precision limitations may cause the STL sort to hang.  Instead, put the sort position
 296     // just outside the bounding sphere of the cloud in the direction of the camera.
 297     _vecSortPos = -cam.ViewDir();
 298     _vecSortPos *= (1.1 * _boundingBox.GetRadius());
 299     _vecSortPos += _boundingBox.GetCenter();
 300
 301     // sort the particles from back to front wrt the camera position.
 302     _SortParticles(cam.ViewDir(), _vecSortPos, SKY_CLOUD_SORT_TOWARD);
 303
 304     // we can't use the view direction optimization when the cloud is split, or we get a lot
 305     // of popping of objects in and out of cloud cover.  For consistency, though, we need to update
 306     // the cached sort direction, since we just sorted the particles.
 307     // _vecLastSortViewDir = GLVU::GetCurrent()->GetCurrentCam()->ViewDir();
 308     // _vecLastSortCamPos = GLVU::GetCurrent()->GetCurrentCam()->Orig;
 309
 310     // compute the split distance.
 311     vecCloudSpaceSplit  -= _vecSortPos;
 312     _rSplitDistance = vecCloudSpaceSplit * cam.ViewDir();
 313   }
 314
 315   // set the material state / properties that clouds use for rendering:
 316   // Enables blending, with blend func (ONE, ONE_MINUS_SRC_ALPHA).
 317   // Enables alpha test to discard completely transparent fragments.
 318   // Disables depth test.
 319   // Enables texturing, with modulation, and the texture set to the shared splat texture.
 320   s_pMaterial->Activate();
 321
 322   Vec4f color;
 323   Vec3f eyeDir;
 324
 325   // Draw the particles using immediate mode.
 326   glBegin(GL_QUADS);
 327
 328   // if bBackHalf is false, then we just continue where we left off.  If it is true, we
 329   // reset the iterator to the beginning of the sorted list.
 330   static ParticleIterator iter;
 331   if (bBackHalf)
 332     iter = _particles.begin();
 333
 334   // iterate over the particles and render them.
 335   for (; iter != _particles.end(); ++iter)
 336   {
 337     SkyCloudParticle *p = *iter;
 338
 339     if (bBackHalf && (p->GetSquareSortDistance() < _rSplitDistance))
 340       break;
 341
 342     // Start with ambient light
 343     color = p->GetBaseColor();
 344
 345     if (_bUsePhaseFunction) // use the phase function for anisotropic scattering.
 346     {
 347       eyeDir = cam.Orig;
 348       eyeDir -= p->GetPosition();
 349       eyeDir.Normalize();
 350       float pf;
 351
 352       // add the color contribution to this particle from each light source, modulated by
 353       // the phase function.  See _PhaseFunction() documentation for details.
 354       for (int i = 0; i < p->GetNumLitColors(); i++)
 355       {
 356         pf = _PhaseFunction(_lightDirections[i], eyeDir);
 357         // expand this to avoid temporary vector creation in the inner loop
 358         color.x += p->GetLitColor(i).x * pf;
 359         color.y += p->GetLitColor(i).y * pf;
 360         color.z += p->GetLitColor(i).z * pf;
 361       }
 362     }
 363     else  // just use isotropic scattering instead.
 364     {
 365       for (int i = 0; i < p->GetNumLitColors(); ++i)
 366       {
 367         color += p->GetLitColor(i);
 368       }
 369     }
 370
 371     // set the transparency independently of the colors.
 372     color.w = 1 - s_rTransparency;
 373
 374     // draw the particle as a textured billboard.
 375     DrawQuad((*iter)->GetPosition(), cam.X * p->GetRadius(), cam.Y * p->GetRadius(), color);
 376   }
 377   glEnd();
 378
 379   return SKYRESULT_OK;
 380 }
 381
 382 //------------------------------------------------------------------------------
 383 // Function               : SkyCloud::Illuminate
 384 // Description      :
 385 //------------------------------------------------------------------------------
 386 /**
 387  * @fn SkyCloud::Illuminate(SkyLight *pLight, SkyRenderableInstance* pInstance, bool bReset)
 388  * @brief Compute the illumination of the cloud by the lightsource @a pLight
 389  *
 390  * This method uses graphics hardware to compute multiple forward scattering at each cloud
 391  * in the cloud of light from the directional light source @a pLight.  The algorithm works
 392  * by successively subtracting "light" from an initially white (fully lit) frame buffer by
 393  * using hardware blending and read back.  The method stores the illumination from each light
 394  * source passed to it separately at each particle, unless @a bReset is true, in which case
 395  * the lists of illumination in the particles are reset before the lighting is computed.
 396  *
 397  */
 398 SKYRESULT SkyCloud::Illuminate(SkyLight *pLight, SkyRenderableInstance* pInstance, bool bReset)
 399 {
 400   int iOldVP[4];
 401
 402   glGetIntegerv(GL_VIEWPORT, iOldVP);
 403   glViewport(0, 0, s_iShadeResolution, s_iShadeResolution);
 404
 405   Vec3f vecDir(pLight->GetDirection());
 406
 407   // if this is the first pass through the lights, reset will be true, and the cached light
 408   // directions should be updated.  Light directions are cached in cloud space to accelerate
 409   // computation of the phase function, which depends on light direction and view direction.
 410   if (bReset)
 411     _lightDirections.clear();
 412   _lightDirections.push_back(vecDir); // cache the (unit-length) light direction
 413
 414   // compute the light/sort position for particles from the light direction.
 415   // don't just use the camera position -- if it is too far away from the cloud, then
 416   // precision limitations may cause the STL sort to hang.  Instead, put the sort position
 417   // just outside the bounding sphere of the cloud in the direction of the camera.
 418   Vec3f vecLightPos(vecDir);
 419   vecLightPos *= (1.1*_boundingBox.GetRadius());
 420   vecLightPos += _boundingBox.GetCenter();
 421
 422   // Set up a camera to look at the cloud from the light position.  Since the sun is an infinite
 423   // light source, this camera will use an orthographic projection tightly fit to the bounding
 424   // sphere of the cloud.
 425   Camera cam;
 426
 427   // Avoid degenerate camera bases.
 428   Vec3f vecUp(0, 1, 0);
 429   if (fabs(vecDir * vecUp) - 1 < 1e-6) // check that the view and up directions are not parallel.
 430     vecUp.Set(1, 0, 0);
 431
 432   cam.LookAt(vecLightPos, _boundingBox.GetCenter(), vecUp);
 433
 434   // sort the particles away from the light source.
 435   _SortParticles(cam.ViewDir(), vecLightPos, SKY_CLOUD_SORT_AWAY);
 436
 437   // projected dist to cntr along viewdir
 438   float DistToCntr = (_boundingBox.GetCenter() - vecLightPos) * cam.ViewDir();
 439
 440   // calc tight-fitting near and far distances for the orthographic frustum
 441   float rNearDist = DistToCntr - _boundingBox.GetRadius();
 442   float rFarDist  = DistToCntr + _boundingBox.GetRadius();
 443
 444   // set the modelview matrix from this camera.
 445   glMatrixMode(GL_MODELVIEW);
 446   glPushMatrix();
 447   float M[16];
 448
 449
 450   cam.GetModelviewMatrix(M);
 451   glLoadMatrixf(M);
 452
 453   // switch to parallel projection
 454   glMatrixMode(GL_PROJECTION);
 455   glPushMatrix();
 456   glLoadIdentity();
 457   glOrtho(-_boundingBox.GetRadius(), _boundingBox.GetRadius(),
 458           -_boundingBox.GetRadius(), _boundingBox.GetRadius(),
 459           rNearDist, rFarDist);
 460
 461   // set the material state / properties that clouds use for shading:
 462   // Enables blending, with blend func (ONE, ONE_MINUS_SRC_ALPHA).
 463   // Enables alpha test to discard completely transparent fragments.
 464   // Disables depth test.
 465   // Enables texturing, with modulation, and the texture set to the shared splat texture.
 466   s_pShadeMaterial->Activate();
 467
 468   // these are used for projecting the particle position to determine where to read pixels.
 469   double MM[16], PM[16];
 470   int VP[4] = { 0, 0, s_iShadeResolution, s_iShadeResolution };
 471   glGetDoublev(GL_MODELVIEW_MATRIX, MM);
 472   glGetDoublev(GL_PROJECTION_MATRIX, PM);
 473
 474   // initialize back buffer to all white -- modulation darkens areas where cloud particles
 475   // absorb light, and lightens it where they scatter light in the forward direction.
 476   glClearColor(1, 1, 1, 1);
 477   glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
 478
 479   float rPixelsPerLength = s_iShadeResolution / (2 * _boundingBox.GetRadius());
 480
 481   // the solid angle over which we will sample forward-scattered light.
 482   float rSolidAngle = 0.09;
 483   int i = 0;
 484   int iNumFailed = 0;
 485   for (ParticleIterator iter = _particles.begin(); iter != _particles.end(); ++iter, ++i)
 486   {
 487     Vec3f vecParticlePos = (*iter)->GetPosition();
 488
 489     Vec3f vecOffset(vecLightPos);
 490     vecOffset -= vecParticlePos;
 491
 492     // compute the pixel area to read back in order to integrate the illumination of the particle
 493     // over a constant solid angle.
 494     float rDistance  = fabs(cam.ViewDir() * vecOffset) - rNearDist;
 495
 496     float rArea      = rSolidAngle * rDistance * rDistance;
 497     int   iPixelDim  = sqrt(rArea) * rPixelsPerLength;
 498     int   iNumPixels = iPixelDim * iPixelDim;
 499     if (iNumPixels < 1)
 500     {
 501       iNumPixels = 1;
 502       iPixelDim = 1;
 503     }
 504
 505     // the scale factor to convert the read back pixel colors to an average illumination of the area.
 506     float rColorScaleFactor = rSolidAngle / (iNumPixels * 255.0f);
 507
 508     unsigned char *c = new unsigned char[4 * iNumPixels];
 509
 510     Vec3d vecWinPos;
 511
 512     // find the position in the buffer to which the particle position projects.
 513     if (!gluProject(vecParticlePos.x, vecParticlePos.y, vecParticlePos.z,
 514                     MM, PM, VP,
 515                     &(vecWinPos.x), &(vecWinPos.y), &(vecWinPos.z)))
 516     {
 517       FAIL_RETURN_MSG(SKYRESULT_FAIL,
 518                       "Error: SkyCloud::Illuminate(): failed to project particle position.");
 519     }
 520
 521     // offset the projected window position by half the size of the readback region.
 522     vecWinPos.x -= 0.5 * iPixelDim;
 523     if (vecWinPos.x < 0) vecWinPos.x = 0;
 524     vecWinPos.y -= 0.5 * iPixelDim;
 525     if (vecWinPos.y < 0) vecWinPos.y = 0;
 526
 527     // read back illumination of this particle from the buffer.
 528     glReadBuffer(GL_BACK);
 529     glReadPixels(vecWinPos.x, vecWinPos.y, iPixelDim, iPixelDim, GL_RGBA, GL_UNSIGNED_BYTE, c);
 530
 531     // scattering coefficient vector.
 532     Vec4f vecScatter(s_rScatterFactor, s_rScatterFactor, s_rScatterFactor, 1);
 533
 534     // add up the read back pixels (only need one component -- its grayscale)
 535     int iSum = 0;
 536     for (int k = 0; k < 4 * iNumPixels; k+=4)
 537       iSum += c[k];
 538     delete [] c;
 539
 540     // compute the amount of light scattered to this particle by particles closer to the light.
 541     // this is the illumination over the solid angle that we measured (using glReadPixels) times
 542     // the scattering coefficient (vecScatter);
 543     Vec4f vecScatteredAmount(iSum * rColorScaleFactor,
 544                              iSum * rColorScaleFactor,
 545                              iSum * rColorScaleFactor,
 546                              1 - s_rTransparency);
 547     vecScatteredAmount &= vecScatter;
 548
 549     // the color of th particle (iter) contributed by this light source (pLight) is the
 550     // scattered light from the part of the cloud closer to the light, times the diffuse color
 551     // of the light source.  The alpha is 1 - the uniform transparency of all particles (modulated
 552     // by the splat texture).
 553     Vec4f vecColor = vecScatteredAmount;
 554     vecColor      &= pLight->GetDiffuse();
 555     vecColor.w     = 1 - s_rTransparency;
 556
 557     // add this color to the list of lit colors for the particle.  The contribution from each light
 558     // is kept separate because the phase function we apply at runtime depends on the light vector
 559     // for each light source separately.  This view-dependent effect is impossible without knowing
 560     // the amount of light contributed for each light.  This, of course, assumes the clouds will
 561     // be lit by a reasonably small number of lights (The sun plus some simulation of light reflected
 562     // from the sky and / or ground.) This technique works very well for simulating anisotropic
 563     // illumination by skylight.
 564     if (bReset)
 565     {
 566       (*iter)->SetBaseColor(s_pMaterial->GetAmbient());
 567       (*iter)->ClearLitColors();
 568       (*iter)->AddLitColor(vecColor);
 569     }
 570     else
 571     {
 572       (*iter)->AddLitColor(vecColor);
 573     }
 574
 575     // the following computation (scaling of the scattered amount by the phase function) is done
 576     // after the lit color is stored so we don't add the scattering to this particle twice.
 577     vecScatteredAmount *= 1.5; // rayleigh scattering phase function for angle of zero or 180 = 1.5!
 578
 579     // clamp the color
 580     if (vecScatteredAmount.x > 1) vecScatteredAmount.x = 1;
 581     if (vecScatteredAmount.y > 1) vecScatteredAmount.y = 1;
 582     if (vecScatteredAmount.z > 1) vecScatteredAmount.z = 1;
 583     vecScatteredAmount.w = 1 - s_rTransparency;
 584
 585     vecScatteredAmount.x = 0.50;  vecScatteredAmount.y = 0.60;  vecScatteredAmount.z = 0.70;
 586
 587     // Draw the particle as a texture billboard.  Use the scattered light amount as the color to
 588     // simulate forward scattering of light by this particle.
 589     glBegin(GL_QUADS);
 590     DrawQuad(vecParticlePos, cam.X * (*iter)->GetRadius(), cam.Y * (*iter)->GetRadius(), vecScatteredAmount);
 591     glEnd();
 592
 593     //glutSwapBuffers(); // Uncomment this swap buffers to visualize cloud illumination computation.
 594   }
 595
 596   // Note: here we could optionally store the current back buffer as a shadow image
 597   // to be projected from the light position onto the scene.  This way we can have clouds shadow
 598   // the environment.
 599
 600   // restore matrix stack and viewport.
 601   glMatrixMode(GL_PROJECTION);
 602   glPopMatrix();
 603   glMatrixMode(GL_MODELVIEW);
 604   glPopMatrix();
 605   glViewport(iOldVP[0], iOldVP[1], iOldVP[2], iOldVP[3]);
 606
 607   return SKYRESULT_OK;
 608 }
 609
 610
 611 //------------------------------------------------------------------------------
 612 // Function               : SkyCloud::CopyBoundingVolume
 613 // Description      :
 614 //------------------------------------------------------------------------------
 615 /**
 616  * @fn SkyCloud::CopyBoundingVolume() const
 617  * @brief Returns a new copy of the SkyMinMaxBox for this cloud.
 618  */
 619 SkyMinMaxBox* SkyCloud::CopyBoundingVolume() const
 620 {
 621   SkyMinMaxBox *pBox = new SkyMinMaxBox();
 622   pBox->SetMax(_boundingBox.GetMax());
 623   pBox->SetMin(_boundingBox.GetMin());
 624   return pBox;
 625 }
 626
 627 SKYRESULT SkyCloud::Load(const unsigned char *data, unsigned int size,
 628                          float rScale, /* = 1.0f */
 629                          double latitude, double longitude )
 630 {
 631   unsigned int iNumParticles;
 632   Vec3f vecCenter = Vec3f::ZERO;
 633   //Vec3f vecCenter;
 634   //float rRadius;
 635
 636   //_boundingBox.SetMin(vecCenter - Vec3f(rRadius, rRadius, rRadius));
 637   //_boundingBox.SetMax(vecCenter + Vec3f(rRadius, rRadius, rRadius));
 638
 639   for (unsigned int i = 0; i < size*size*4; i += 4)
 640   {
 641     Vec3f pos;
 642     Vec4f color;
 643     float radius;
 644
 645     color.x = data[i];                  // FIXME: Which unit?
 646     color.y = data[i+1];
 647     color.z = data[i+2];
 648     color.w = data[i+3];
 649
 650     radius = (color.w/255) * rScale;
 651     //radius = (color.x * color.y * color.z * color.w * rScale) /  4096;
 652
 653
 654     pos.x = (i / (size*4)) * 10;        // FIXME: Which unit?
 655     pos.y = (i % (size*4)) * 10;
 656     pos.z = radius / 2;
 657
 658     if (radius > 0)
 659     {
 660       SkyCloudParticle *pParticle = new SkyCloudParticle((pos + vecCenter) * rScale, radius * rScale, color);
 661       _boundingBox.AddPoint(pParticle->GetPosition());
 662
 663       _particles.push_back(pParticle);
 664     }
 665   }
 666
 667   // this is just a bad hack to align cloud field from skyworks with local horizon at KSFO
 668   // this "almost" works not quite the right solution okay to get some up and running
 669   // we need to develop our own scheme for loading and positioning clouds
 670   Mat33f R;
 671   Vec3f  moveit;
 672
 673   R.Set( 0, 1, 0,
 674          1, 0, 0,
 675          0, 0, 1);
 676   // clouds sit in the y-z plane and x-axis is the vertical cloud height
 677   Rotate( R );
 678
 679 // rotate the cloud field about the fgfs z-axis based on initial longitude
 680 float ex = 1.0;
 681 float ey = 1.0;
 682 float ez = 1.0;
 683 float phi = longitude / 57.29578;
 684 float one_min_cos = 1 - cos(phi);
 685
 686 R.Set(
 687 cos(phi) + one_min_cos*ex*ex, one_min_cos*ex*ey - ez*sin(phi), one_min_cos*ex*ez + ey*sin(phi),
 688 one_min_cos*ex*ey + ez*sin(phi), cos(phi) + one_min_cos*ey*ey, one_min_cos*ey*ez - ex*sin(phi),
 689 one_min_cos*ex*ez - ey*sin(phi), one_min_cos*ey*ez + ex*sin(phi), cos(phi) + one_min_cos*ez*ez );
 690
 691 Rotate( R );
 692
 693 // okay now that let's rotate about a vector for latitude where longitude forms the
 694 // components of a unit vector in the x-y plane
 695 ex = sin( longitude  / 57.29578 );
 696 ey = -cos( longitude  / 57.29578 );
 697 ez = 0.0;
 698 phi = latitude / 57.29578;
 699 one_min_cos = 1 - cos(phi);
 700
 701 R.Set(
 702 cos(phi) + one_min_cos*ex*ex, one_min_cos*ex*ey - ez*sin(phi), one_min_cos*ex*ez + ey*sin(phi),
 703 one_min_cos*ex*ey + ez*sin(phi), cos(phi) + one_min_cos*ey*ey, one_min_cos*ey*ez - ex*sin(phi),
 704 one_min_cos*ex*ez - ey*sin(phi), one_min_cos*ey*ez + ex*sin(phi), cos(phi) + one_min_cos*ez*ez );
 705
 706 Rotate( R );
 707 // need to calculate an offset to place the clouds at ~3000 feet MSL  ATM this is an approximation
 708 // to move the clouds to some altitude above sea level. At some locations this could be underground
 709 // will need a better scheme to position clouds per user preferences
 710 float cloud_level_msl = 3000.0f;
 711
 712 float x_offset = ex * cloud_level_msl;
 713 float y_offset = ey * cloud_level_msl;
 714 float z_offset = cloud_level_msl * 0.5;
 715 moveit.Set( x_offset, y_offset, z_offset  );
 716
 717   Translate( moveit );
 718
 719   return SKYRESULT_OK;
 720 }
 721
 722
 723 SKYRESULT SkyCloud::Load(const SkyArchive &archive,
 724                          float rScale, /* = 1.0f */
 725                          double latitude, double longitude )
 726 {
 727   unsigned int iNumParticles;
 728   Vec3f vecCenter = Vec3f::ZERO;
 729   //Vec3f vecCenter;
 730   //float rRadius;
 731   //archive.FindVec3f("CldCenter", &vecCenter);
 732   //archive.FindFloat32("CldRadius", &rRadius);
 733
 734   //_boundingBox.SetMin(vecCenter - Vec3f(rRadius, rRadius, rRadius));
 735   //_boundingBox.SetMax(vecCenter + Vec3f(rRadius, rRadius, rRadius));
 736
 737   archive.FindUInt32("CldNumParticles", &iNumParticles);
 738   iNumParticles = ulEndianLittle32(iNumParticles);
 739
 740   //if (!bLocal)
 741     archive.FindVec3f("CldCenter", &vecCenter);
 742     vecCenter.x = ulEndianLittleFloat(vecCenter.x);
 743     vecCenter.y = ulEndianLittleFloat(vecCenter.y);
 744     vecCenter.z = ulEndianLittleFloat(vecCenter.z);
 745
 746   Vec3f *pParticlePositions = new Vec3f[iNumParticles];
 747   float *pParticleRadii     = new float[iNumParticles];
 748   Vec4f *pParticleColors    = new Vec4f[iNumParticles];
 749
 750   unsigned int iNumBytes;
 751   archive.FindData("CldParticlePositions", ANY_TYPE, (void**const)&pParticlePositions, &iNumBytes);
 752   archive.FindData("CldParticleRadii",     ANY_TYPE, (void**const)&pParticleRadii,     &iNumBytes);
 753   archive.FindData("CldParticleColors",    ANY_TYPE, (void**const)&pParticleColors,    &iNumBytes);
 754
 755   for (unsigned int i = 0; i < iNumParticles; ++i)
 756   {
 757
 758     pParticlePositions[i].x = ulEndianLittleFloat(pParticlePositions[i].x);
 759     pParticlePositions[i].y = ulEndianLittleFloat(pParticlePositions[i].y);
 760     pParticlePositions[i].z = ulEndianLittleFloat(pParticlePositions[i].z);
 761
 762     pParticleRadii[i] = ulEndianLittleFloat(pParticleRadii[i]);
 763
 764     pParticleColors[i].x = ulEndianLittleFloat(pParticleColors[i].x);
 765     pParticleColors[i].y = ulEndianLittleFloat(pParticleColors[i].y);
 766     pParticleColors[i].z = ulEndianLittleFloat(pParticleColors[i].z);
 767     pParticleColors[i].w = ulEndianLittleFloat(pParticleColors[i].w);
 768
 769     SkyCloudParticle *pParticle = new SkyCloudParticle((pParticlePositions[i] + vecCenter) * rScale,
 770                                                        pParticleRadii[i] * rScale,
 771                                                        pParticleColors[i]);
 772     _boundingBox.AddPoint(pParticle->GetPosition());
 773
 774     _particles.push_back(pParticle);
 775   }
 776   // this is just a bad hack to align cloud field from skyworks with local horizon at KSFO
 777   // this "almost" works not quite the right solution okay to get some up and running
 778   // we need to develop our own scheme for loading and positioning clouds
 779   Mat33f R;
 780   Vec3f  moveit;
 781
 782   R.Set( 0, 1, 0,
 783                          1, 0, 0,
 784                          0, 0, 1);
 785   // clouds sit in the y-z plane and x-axis is the vertical cloud height
 786   Rotate( R );
 787
 788 // rotate the cloud field about the fgfs z-axis based on initial longitude
 789 float ex = 0.0;
 790 float ey = 0.0;
 791 float ez = 1.0;
 792 float phi = longitude / 57.29578;
 793 float one_min_cos = 1 - cos(phi);
 794
 795 R.Set(
 796 cos(phi) + one_min_cos*ex*ex, one_min_cos*ex*ey - ez*sin(phi), one_min_cos*ex*ez + ey*sin(phi),
 797 one_min_cos*ex*ey + ez*sin(phi), cos(phi) + one_min_cos*ey*ey, one_min_cos*ey*ez - ex*sin(phi),
 798 one_min_cos*ex*ez - ey*sin(phi), one_min_cos*ey*ez + ex*sin(phi), cos(phi) + one_min_cos*ez*ez );
 799
 800 Rotate( R );
 801
 802 // okay now that let's rotate about a vector for latitude where longitude forms the
 803 // components of a unit vector in the x-y plane
 804 ex = sin( longitude  / 57.29578 );
 805 ey = -cos( longitude  / 57.29578 );
 806 ez = 0.0;
 807 phi = latitude / 57.29578;
 808 one_min_cos = 1 - cos(phi);
 809
 810 R.Set(
 811 cos(phi) + one_min_cos*ex*ex, one_min_cos*ex*ey - ez*sin(phi), one_min_cos*ex*ez + ey*sin(phi),
 812 one_min_cos*ex*ey + ez*sin(phi), cos(phi) + one_min_cos*ey*ey, one_min_cos*ey*ez - ex*sin(phi),
 813 one_min_cos*ex*ez - ey*sin(phi), one_min_cos*ey*ez + ex*sin(phi), cos(phi) + one_min_cos*ez*ez );
 814
 815 Rotate( R );
 816 // need to calculate an offset to place the clouds at ~3000 feet MSL  ATM this is an approximation
 817 // to move the clouds to some altitude above sea level. At some locations this could be underground
 818 // will need a better scheme to position clouds per user preferences
 819 float cloud_level_msl = 3000.0f;
 820
 821 float x_offset = ex * cloud_level_msl;
 822 float y_offset = ey * cloud_level_msl;
 823 float z_offset = cloud_level_msl * 0.5;
 824 moveit.Set( x_offset, y_offset, z_offset  );
 825
 826   Translate( moveit );
 827
 828   return SKYRESULT_OK;
 829 }
 830
 831
 832 //------------------------------------------------------------------------------
 833 // Function               : SkyCloud::Save
 834 // Description      :
 835 //------------------------------------------------------------------------------
 836 /**
 837 * @fn SkyCloud::Save(SkyArchive &archive) const
 838 * @brief Saves the cloud data to @a archive.
 839 *
 840 * @todo <WRITE EXTENDED SkyCloud::Save FUNCTION DOCUMENTATION>
 841 */
 842 SKYRESULT SkyCloud::Save(SkyArchive &archive) const
 843 {
 844   SkyArchive myArchive("Cloud");
 845   //myArchive.AddVec3f("CldCenter", _center);
 846   //myArchive.AddFloat32("CldRadius", _boundingBox.GetRadius());
 847   myArchive.AddUInt32("CldNumParticles", _particles.size());
 848
 849   // make temp arrays
 850   Vec3f *pParticlePositions = new Vec3f[_particles.size()];
 851   float *pParticleRadii     = new float[_particles.size()];
 852   Vec4f *pParticleColors    = new Vec4f[_particles.size()];
 853
 854   unsigned int i = 0;
 855
 856   for (ParticleConstIterator iter = _particles.begin(); iter != _particles.end(); ++iter, ++i)
 857   {
 858     pParticlePositions[i] = (*iter)->GetPosition(); // position around origin
 859     pParticleRadii[i]     = (*iter)->GetRadius();
 860     pParticleColors[i]    = (*iter)->GetBaseColor();
 861   }
 862
 863   myArchive.AddData("CldParticlePositions",
 864                     ANY_TYPE,
 865                     pParticlePositions,
 866                     sizeof(Vec3f),
 867                     _particles.size());
 868
 869   myArchive.AddData("CldParticleRadii",
 870                     ANY_TYPE,
 871                     pParticleRadii,
 872                     sizeof(float),
 873                     _particles.size());
 874
 875   myArchive.AddData("CldParticleColors",
 876                     ANY_TYPE,
 877                     pParticleColors,
 878                     sizeof(Vec3f),
 879                     _particles.size());
 880
 881   archive.AddArchive(myArchive);
 882
 883   return SKYRESULT_OK;
 884 }
 885
 886
 887 //------------------------------------------------------------------------------
 888 // Function               : SkyCloud::Rotate
 889 // Description      :
 890 //------------------------------------------------------------------------------
 891 /**
 892  * @fn SkyCloud::Rotate(const Mat33f& rot)
 893  * @brief @todo <WRITE BRIEF SkyCloud::Rotate DOCUMENTATION>
 894  *
 895  * @todo <WRITE EXTENDED SkyCloud::Rotate FUNCTION DOCUMENTATION>
 896  */
 897 void SkyCloud::Rotate(const Mat33f& rot)
 898 {
 899   _boundingBox.Clear();
 900   for (int i = 0; i < _particles.size(); ++i)
 901   {
 902     _particles[i]->SetPosition(rot * _particles[i]->GetPosition());
 903     _boundingBox.AddPoint(_particles[i]->GetPosition());
 904   }
 905 }
 906
 907
 908 //------------------------------------------------------------------------------
 909 // Function               : SkyCloud::Translate
 910 // Description      :
 911 //------------------------------------------------------------------------------
 912 /**
 913  * @fn SkyCloud::Translate(const Vec3f& trans)
 914  * @brief @todo <WRITE BRIEF SkyCloud::Translate DOCUMENTATION>
 915  *
 916  * @todo <WRITE EXTENDED SkyCloud::Translate FUNCTION DOCUMENTATION>
 917  */
 918 void SkyCloud::Translate(const Vec3f& trans)
 919 {
 920   for (int i = 0; i < _particles.size(); ++i)
 921   {
 922     _particles[i]->SetPosition(_particles[i]->GetPosition() + trans);
 923   }
 924   _boundingBox.SetMax(_boundingBox.GetMax() + trans);
 925   _boundingBox.SetMin(_boundingBox.GetMin() + trans);
 926 }
 927
 928
 929 //------------------------------------------------------------------------------
 930 // Function               : SkyCloud::Scale
 931 // Description      :
 932 //------------------------------------------------------------------------------
 933 /**
 934  * @fn SkyCloud::Scale(const float scale)
 935  * @brief @todo <WRITE BRIEF SkyCloud::Scale DOCUMENTATION>
 936  *
 937  * @todo <WRITE EXTENDED SkyCloud::Scale FUNCTION DOCUMENTATION>
 938  */
 939 void SkyCloud::Scale(const float scale)
 940 {
 941   _boundingBox.Clear();
 942   for (int i = 0; i < _particles.size(); ++i)
 943   {
 944     _particles[i]->SetPosition(_particles[i]->GetPosition() * scale);
 945     _boundingBox.AddPoint(_particles[i]->GetPosition());
 946   }
 947 }
 948
 949
 950 //------------------------------------------------------------------------------
 951 // Function               : SkyCloud::_SortParticles
 952 // Description      :
 953 //------------------------------------------------------------------------------
 954 /**
 955  * @fn SkyCloud::_SortParticles(const Vec3f& vecViewDir, const Vec3f& sortPoint, SortDirection dir)
 956  * @brief Sorts the cloud particles in the direction specified by @a dir.
 957  *
 958  * @vecSortPoint is assumed to already be transformed into the basis space of the cloud.
 959  */
 960 void SkyCloud::_SortParticles(const Vec3f& vecViewDir,
 961                               const Vec3f& vecSortPoint,
 962                               SortDirection dir)
 963 {
 964   Vec3f partPos;
 965         for (int i = 0; i < _particles.size(); ++i)
 966         {
 967                 partPos = _particles[i]->GetPosition();
 968                 partPos -= vecSortPoint;
 969                 _particles[i]->SetSquareSortDistance(partPos * vecViewDir);//partPos.LengthSqr());
 970         }
 971
 972   switch (dir)
 973   {
 974   case SKY_CLOUD_SORT_TOWARD:
 975     std::sort(_particles.begin(), _particles.end(), _towardComparator);
 976     break;
 977   case SKY_CLOUD_SORT_AWAY:
 978     std::sort(_particles.begin(), _particles.end(), _awayComparator);
 979     break;
 980   default:
 981     break;
 982   }
 983 }
 984
 985
 986
 987 //------------------------------------------------------------------------------
 988 // Function               : EvalHermite
 989 // Description      :
 990 //------------------------------------------------------------------------------
 991 /**
 992  * EvalHermite(float pA, float pB, float vA, float vB, float u)
 993  * @brief Evaluates Hermite basis functions for the specified coefficients.
 994  */
 995 inline float EvalHermite(float pA, float pB, float vA, float vB, float u)
 996 {
 997   float u2=(u*u), u3=u2*u;
 998   float B0 = 2*u3 - 3*u2 + 1;
 999   float B1 = -2*u3 + 3*u2;
1000   float B2 = u3 - 2*u2 + u;
1001   float B3 = u3 - u;
1002   return( B0*pA + B1*pB + B2*vA + B3*vB );
1003 }
1004
1005 // NORMALIZED GAUSSIAN INTENSITY MAP (N must be a power of 2)
1006
1007 //------------------------------------------------------------------------------
1008 // Function               : CreateGaussianMap
1009 // Description      :
1010 //------------------------------------------------------------------------------
1011 /**
1012  * CreateGaussianMap(int N)
1013  *
1014  * Creates a 2D gaussian image using a hermite surface.
1015  */
1016 unsigned char* CreateGaussianMap(int N)
1017 {
1018   float *M = new float[2*N*N];
1019   unsigned char *B = new unsigned char[4*N*N];
1020   float X,Y,Y2,Dist;
1021   float Incr = 2.0f/N;
1022   int i=0;
1023   int j = 0;
1024   Y = -1.0f;
1025   for (int y=0; y<N; y++, Y+=Incr)
1026   {
1027     Y2=Y*Y;
1028     X = -1.0f;
1029     for (int x=0; x<N; x++, X+=Incr, i+=2, j+=4)
1030     {
1031       Dist = (float)sqrt(X*X+Y2);
1032       if (Dist>1) Dist=1;
1033       M[i+1] = M[i] = EvalHermite(0.4f,0,0,0,Dist);// * (1 - noise);
1034       B[j+3] = B[j+2] = B[j+1] = B[j] = (unsigned char)(M[i] * 255);
1035     }
1036   }
1037   SAFE_DELETE_ARRAY(M);
1038   return(B);
1039 }
1040
1041
1042 //------------------------------------------------------------------------------
1043 // Function               : SkyCloud::_CreateSplatTexture
1044 // Description      :
1045 //------------------------------------------------------------------------------
1046 /**
1047  * @fn SkyCloud::_CreateSplatTexture(unsigned int iResolution)
1048  * @brief Creates the texture map used for cloud particles.
1049  */
1050 void SkyCloud::_CreateSplatTexture(unsigned int iResolution)
1051 {
1052   unsigned char *splatTexture = CreateGaussianMap(iResolution);
1053   SkyTexture texture;
1054   TextureManager::InstancePtr()->Create2DTextureObject(texture, iResolution, iResolution,
1055                                                        GL_RGBA, splatTexture);
1056
1057   s_pMaterial->SetTexture(0, GL_TEXTURE_2D, texture);
1058   s_pShadeMaterial->SetTexture(0, GL_TEXTURE_2D, texture);
1059   s_pMaterial->SetTextureParameter(0, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
1060   s_pShadeMaterial->SetTextureParameter(0, GL_TEXTURE_MIN_FILTER, GL_LINEAR);
1061   s_pMaterial->SetTextureParameter(0, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
1062   s_pShadeMaterial->SetTextureParameter(0, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
1063   s_pMaterial->SetTextureParameter(0, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
1064   s_pShadeMaterial->SetTextureParameter(0, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE);
1065   s_pMaterial->EnableTexture(0, true);
1066   s_pShadeMaterial->EnableTexture(0, true);
1067
1068   SAFE_DELETE_ARRAY(splatTexture);
1069 }
1070
1071
1072 //------------------------------------------------------------------------------
1073 // Function               : SkyCloud::_PhaseFunction
1074 // Description      :
1075 //------------------------------------------------------------------------------
1076 /**
1077  * @fn SkyCloud::_PhaseFunction(const Vec3f& vecLightDir, const Vec3f& vecViewDir)
1078  * @brief Computes the phase (scattering) function of the given light and view directions.
1079  *
1080  * A phase function is a transfer function that determines, for any angle between incident
1081  * and outgoing directions, how much of the incident light intensity will be
1082  * scattered in the outgoing direction.  For example, scattering by very small
1083  * particles such as those found in clear air, can be approximated using <i>Rayleigh
1084  * scattering</i>.  The phase function for Rayleigh scattering is
1085  * p(q) = 0.75*(1 + cos<sup>2</sup>(q)), where q  is the angle between incident
1086  * and scattered directions.  Scattering by larger particles is more complicated.
1087  * It is described by Mie scattering theory.  Cloud particles are more in the regime
1088  * of Mie scattering than Rayleigh scattering.  However, we obtain good visual
1089  * results by using the simpler Rayleigh scattering phase function as an approximation.
1090  */
1091 float SkyCloud::_PhaseFunction(const Vec3f& vecLightDir, const Vec3f& vecViewDir)
1092 {
1093   float rCosAlpha = vecLightDir * vecViewDir;
1094   return .75f * (1 + rCosAlpha * rCosAlpha); // rayleigh scattering = (3/4) * (1+cos^2(alpha))
1095 }