/* rings.sl
 *
 * Started with Steve May's rmannotes ring.sl.  I've attempted to make the shape
*rings dependent on surface properties. 
 *
 */
#define snoise(x)    (2*noise(x) - 1)
surface rings(
   float Kd = .8;
   float Ka = .2;
   float roughness = .01;
   float size = 2;
   float freq = 1;
   float skewscale = 1;
   float time = .5;
   float fuzz = 0.025;
   float radius = .25;           /* radius of ring */
   float half_width = 0.05;       /* 1/2 width of ring */
   color ambientcolor = color(1);
   color specularcolor = color(1);)
{

  color surface_color, layer_color;
  color layer_opac;
  color blue = color (0,0,1);
  point center;
  float d;

  surface_color = Cs;

  layer_color = blue;
  color layer_curve = color (1,0,0);
  float ss = round(freq*s)*size;
  float tt = round(freq*t)*size;
  ss = ss - freq*s*size;
  tt = tt - freq*t*size;

/*Now to skew the circle as the derivative varies.*/
/*Because I'm working with a windmill, I can compute the derivatives.  The
Ashli viewer seems to have numerical difficulties with dPdu.*/


float dPds = -s*s + t*t + 1;
dPds = dPds*dPds + 4*s*s*t*t + 4*s*s;
float dPdt = s*s - t*t + 1;
dPdt = dPdt*dPdt + 4*s*s*t*t + 4*t*t;
dPds = sqrt(dPds);
dPdt = sqrt(dPdt);
/*
float dPds = -ss*ss + tt*tt + 1;
dPds = dPds*dPds + 4*ss*ss*tt*tt + 4*ss*ss;
float dPdt = ss*ss - tt*tt + 1;
dPdt = dPdt*dPdt + 4*ss*ss*tt*tt + 4*tt*tt;
dPds = sqrt(dPds);
dPdt = sqrt(dPdt);

float dPds = -u*u + v*v + 1;
dPds = dPds*dPds + 4*u*u*v*v + 4*u*u;
float dPdt = u*u - v*v + 1;
dPdt = dPdt*dPdt + 4*u*u*v*v + 4*v*v;
dPds = sqrt(dPds);
dPdt = sqrt(dPdt);
*/
 ss = ss;
 tt = tt;

  center = point (.5, .5, 0);  /* position of ring */
  d = sqrt(ss*ss/dPds + tt*tt/dPdt);
  float curvature = 0;

/*  layer_opac = pulse(radius - half_width, radius + half_width, fuzz, d);
  surface_color = blend(surface_color, layer_color, layer_opac);
*/

float a = radius - half_width;
float b = radius + half_width;
  layer_opac = smoothstep(a-fuzz, a, d) - smoothstep(b-fuzz, b, d);
  surface_color = surface_color * (1 - layer_opac) + layer_color * layer_opac;

/*surface_color = d < radius ? layer_color : surface_color;*/

point Nf=normalize(N);
point V=normalize(-I);

  Ci = layer_curve*curvature + surface_color*(Ka * ambientcolor+ Kd * diffuse(Nf)+ Kd * diffuse(-Nf))+
    specularcolor*specular(Nf,V,roughness)+ specularcolor*specular(-Nf,V,roughness);
}