#include "rmannotes.sl"

light
prism(
      float     intensity = 10;
      point     from = point "shader" (0, 0, 0),
                to = point "shader" (0, 0, 1),
                up = point "world" (0, 1, 0);
      float     coneangle = radians(45),
                conedeltaangle = radians(3),
                beamdistribution = 2;
      float     petals = 5;
      float     ftime = 0.0;
      )
{
  /* set up a left handed coordinate system from the 
     perspective of the light source (the light source
     is at the origin):
     lightZ = positive direction from "from" to "to"
     lightY = lightZ rotated 90 degrees towards the world "up"
     lightX = positive direction to the right of Y and Z */

  uniform point   lightZ = (to - from) / length(to - from),
                  lightX  = normalize(lightZ ^ up),
                  lightY = normalize(lightX ^ lightZ);

  uniform float   cosoutside = cos(coneangle),
                  cosinside = cos(coneangle - conedeltaangle),
                  numpetals;

  float   atten,
          cosangle,
          angle,
          hue,
          cosX, cosY,
          xbeta, ybeta,
          xtheta, ytheta,
          r, theta,
          rneg, thetaneg,
          x, y;
  point   Lnorm, Lxzplane, Lyzplane, Lplanenorm;
  float   rr, rr1, th,
          curve;


  if (even(petals))
    numpetals = 0.5 * petals;
  else
    numpetals = petals;
           

  illuminate(from, lightZ, coneangle)
    {
      Lnorm = normalize(faceforward(L, -lightZ));

      /* xtheta = the angle between the lightYZ plane and Lnorm, in radians */
      /* ytheta = the angle between the lightXZ plane and Lnorm, in radians */

      cosX = Lnorm . lightX;
      cosY = Lnorm . lightY;
      xtheta = 0.5 * PI - acos(cosX);
      ytheta = 0.5 * PI - acos(cosY);

      /* Lxzplane = projection of Lnorm onto lightXZ plane */
      /* Lyzplane = projection of Lnorm onto lightYZ plane */

      Lxzplane = Lnorm - sin(ytheta) * lightY;
      Lyzplane = Lnorm - sin(xtheta) * lightX;
      
      /* xbeta =  the angle between lightZ and Lxzplane, in radians */
      /* ybeta =  the angle between lightZ and Lyzplane, in radians */

      xbeta = acos(lightZ . normalize(Lxzplane));
      if (cosX < 0)
	xbeta = -xbeta;

      ybeta = acos(lightZ . normalize(Lyzplane));
      if (cosY < 0)
	ybeta = -ybeta;

      /* x = lightX component of Lnorm; y = lightY component of Lnorm */

      x = sin(xbeta);
      y = sin(ybeta);
      r = sqrt(x*x + y*y);
      theta = atan(y, x);

      rneg = -r;
      thetaneg = theta + PI;


      

      /* positive theta's */
      
      curve = 0;
      th = numpetals * theta + 2 * PI * ftime;
      rr = 0.35 * sin(th);
      if (abs(r - rr) < 0.05)
	curve = 1;

      th += 2 * PI;
      rr1 = 0.35 * sin(th);
      while (abs(rr - rr1) > 0.001)
	{
	  if (abs(r - rr1) < 0.05)
	    curve = 1;
	  th += 2 * PI;
	  rr1 = 0.35 * sin(th);
	}


      /* negative theta's */
      
      th = numpetals * thetaneg + 2 * PI * ftime;
      rr = 0.35 * sin(th);
      if (abs(rneg - rr) < 0.05)
	curve = 1;

      th -= 2 * PI;
      rr1 = 0.35 * sin(th);
      while (abs(rr - rr1) > 0.001)
	{
	  if (abs(rneg - rr1) < 0.05)
	    curve = 1;
	  th -= 2 * PI;
	  rr1 = 0.35 * sin(th);
	}


      hue = (theta - 2 * PI * ftime) / (2 * PI);


      cosangle = L.lightZ / length(L);
      atten = pow(cosangle, beamdistribution) / (L.L);
      atten *= smoothstep(cosoutside, cosinside, cosangle);


      if (curve == 1)
	Cl = 0.0;
      else
/*	Cl = intensity * color "rgb" (x/0.3*0.5+0.5, y/0.3*0.5+0.5, 0); */
 	Cl = intensity * color "hsv" (hue, r/0.424, 1);
    }
}