#include <maps.h>

/*----------------------------------------------------------------------
 *
 * diaknurl puts an impressed diamond pattern on the surface
 * of a "cylindrical workpiece" with a cylindrical projection.
 *
 * freq - "frequency" of the knurl
 * depth - magnitude of knurl
 * width - width of slots
 * radius - radius of cylinder (so diamonds will be diamond-shaped)
 * zmin,zmax - extent of knurl
 * dampzone - distance from ends to do a smoothstep diminishing of depth
 * maporigin - the bottom point of the center axis of the cylindrical proj.
 * xaxis, yaxis, zaxis - 3 points defining the coordinate system of
 *        the cylinder.  The central axis is along zaxis-cylinder.
 *	  The point on the equator where (uu,vv)==(0,0) is mapped is at xaxis.
 *--------------------------------------------------------------------*/
displacement
diaknurl(float freq = 10, float depth = 0.025, float width = 0.05, float radius = 1/(2*M_PI),
	float zmin = -1e20, float zmax = 1e20, float dampzone = 0,
	point maporigin=point ("shader", 0,0,0),
	point xaxis=point ("shader", 1,0,0), point yaxis=point ("shader", 0,1,0),
	point zaxis=point ("shader", 0,0,1))
{
    point	PP,O,X,Y,Z;
    normal      sN;
    float	uu, vv,
			dist1,dist2,rim,pos1,pos2,S,T,displacement,
			magnitude;

    /* do projection in shader space */
	PP=transform("shader",P);
	O=transform("shader",maporigin);
	X=transform("shader",xaxis);
	Y=transform("shader",yaxis);
	Z=transform("shader",zaxis);
	cylinmap(PP,uu,vv,O,X,Y,Z);

    /* if inside diamond pattern area... */
    if ((vv >= zmin) && (vv <= zmax)) {
    /* calculate new surface point and normal */
	if ((vv >= zmin+dampzone) && (vv <= zmax-dampzone)) {
		magnitude = depth;
	} else {
		magnitude = depth * min(smoothstep(zmin,zmin+dampzone,vv),
			1-smoothstep(zmax-dampzone,zmax,vv));
	}

    /* calculate diamond pattern */
	S = uu * freq * 2 * M_PI * radius;
	T = vv * freq;
	pos1 = mod(S-T,1);
	S = (1-uu) * freq * 2 * M_PI * radius;
	pos2 = mod(S-T,1);

    /* have a smoothstep rim */
	rim = 2*width; /* rim = blur boundary */
	dist1 = abs(pos1 - (width+rim));
	dist2 = abs(pos2 - (width+rim));

	if ( dist2 < dist1 )  /* reverse loops */
		dist1 = dist2;

    /* displace surface */
	displacement = 0.0;
	if (dist1 < width)
		displacement = -1.0;
	else if (dist1 < width+rim)
		displacement = -(1.0-smoothstep(0.0,1.0,(dist1-width)/rim));

	sN = normalize(transform("shader",N));
	P = transform("shader","common",PP + sN*displacement*magnitude);
	N = calculatenormal(P);
    }
}
