#include "rmannotes.sl"

/*
#define pulse(a,b,fuzz,x) (smoothstep((a)-0.5*(fuzz),(a)+0.5*(fuzz),(x)) - \
                           smoothstep((b)-0.5*(fuzz),(b)+0.5*(fuzz),(x)))
*/


displacement 
trackdisp(
	float     hfreq = 1,
                  vfreq = 1;
	float     rail_width = 0.02,
                  Km = 0.1)
{
  float   surface_mag,
          layer_mag;
  float   fuzz = 0.005;
  float   ss, tt;
  float   row, col;
  point   diff_N;
  float r, theta, angle;
  float d, d0, d1;
  float straight = 0,
        curved = 0;
  float rail, ties;


  ss = repeat(s, hfreq);
  tt = repeat(t, vfreq);
  col = whichtile(s, hfreq);
  row = whichtile(t, vfreq);

  if (snoise2(row + 0.05, col + 0.05) > 0.0)
    straight = 1;
  else if (snoise2(row * 3 + 0.05, col * 7 + 0.05) > 0.0)
    curved = 1;


  /* railroad ties layer */

  if (straight == 1)
    {
      /* straight ties */
      
      d0 = intersection(pulse(0.33, 0.67, fuzz, ss),
			pulse(0.33, 0.67, fuzz, tt));
      if (abs(0.5 - ss) > 0.2)
	{
	  d1 = mod(abs(0.5 - ss) * 100.0 + 5.0, 10.0);
	  d1 = intersection(complement(pulse(2.0, 8.0, fuzz, d1)),
			    pulse(0.33, 0.67, fuzz, tt));
	  d0 = union(d0, d1);
	}
      if (abs(0.5 - tt) > 0.2)
	{
	  d1 = mod(abs(0.5 - tt) * 100.0 + 5.0, 10.0);
	  d1 = intersection(complement(pulse(2.0, 8.0, fuzz, d1)),
			    pulse(0.33, 0.67, fuzz, ss));
	  d0 = union(d0, d1);
	}
      ties = d0;
    }
  else
    {
      if (curved == 1)
	{
	  /* curved ties NW SE */

	  d0 = distance((0, 0, 0), (ss, tt, 0));
	  d1 = distance((1.0 - fuzz, 1.0 - fuzz, 0), (ss, tt, 0));
	  if (d0 <= d1)
	    {
	      d = d0;
	      topolar2d(ss, tt, r, theta);
	    }
	  else
	    {
	      d = d1;
	      topolar2d(1.0 - fuzz - ss, 1.0 - fuzz - tt, r, theta);
	    }
	}
      else
	{
	  /* curved ties NE SW */
	  
	  d0 = distance((0, 1.0 - fuzz, 0), (ss, tt, 0));
	  d1 = distance((1.0 - fuzz, 0, 0), (ss, tt, 0));
	  if (d0 <= d1)
	    {
	      d = d0;
	      topolar2d(ss, 1.0 - fuzz - tt, r, theta);
	    }
	  else
	    {
	      d = d1;
	      topolar2d(1.0 - fuzz - ss, tt, r, theta);
	    }
	}
      
      angle = mod(degrees(theta) + 365.0, 10.0);
      ties = intersection(complement(pulse(2.0, 8.0, fuzz, angle)),
			  pulse(0.33, 0.67, fuzz, d));
    }



  /* rail layer */

  if (straight == 1)
    {
      /* straight rails */

      d0 = union(pulse(0.4 - rail_width, 0.4 + rail_width, fuzz, ss),
		 pulse(0.6 - rail_width, 0.6 + rail_width, fuzz, ss));
      d1 = union(pulse(0.4 - rail_width, 0.4 + rail_width, fuzz, tt),
		 pulse(0.6 - rail_width, 0.6 + rail_width, fuzz, tt));
      rail = union(d0, d1);
    }
  else
    {
      /* curved rails */

      if (curved == 1)
	{
	  d0 = distance((0, 0, 0), (ss, tt, 0));
	  d1 = distance((1.0 - fuzz, 1.0 - fuzz, 0), (ss, tt, 0));
	}
      else
	{
	  d0 = distance((0, 1.0 - fuzz, 0), (ss, tt, 0));
	  d1 = distance((1.0 - fuzz, 0, 0), (ss, tt, 0));
	}
      d = min(d0, d1);
      rail = union(pulse(0.4 - rail_width, 0.4 + rail_width, fuzz, d),
		   pulse(0.6 - rail_width, 0.6 + rail_width, fuzz, d));
    }



  /* surface magnitude calculation */

  surface_mag = max(rail, 0.4 * ties);


  diff_N = normalize(N) - normalize(Ng);
  P += Km * surface_mag * normalize(N);
  N = normalize(calculatenormal(P)) + diff_N;
}