// Macros include file

// Sign(value)
//    Returns the sign of a number
// SuperCone (Xr1,Yr1,Xr2,Yr2)
//    Elliptical cone from z=0 to z=1
// supercone (R1,R2,X1,Y1,Z1,R3,R4,X2,Y2,Z2)
//    Elliptical cone anywhere you want it
// PointTo (p)
//    reorients object so that former +Z axis now points along vector p
//    does not scale objects
// ellipsetorus (Xr,Zr,R,Sa,Ea,N)
//    Elliptical torus (o so slow)
//    Xr=major radius along x
//    Zr=Major radius along z
//    R=minor radius
//    Sa=starting angle (clockwise from +z)
//    Ea=ending angle
//    N=blob resolution
// roundbox (P1,P2,R)
//    creates box from P1 to P2 and rounds edges by R
// roundcylinder (p1,p2,r1,r2)
//    cylinder { p1,p2,r1 } rounded by r2
// roundcone (p1,r1,p2,r2,r3)
//    rounded Cone equivalent to cone { p1,r1 p2,r2 } but rounded by r3
// roundprism (h1,h2,array[points],r)
//    rounded prism
//    equivalent to prism { h1,h2,dimension_size[points] points[0]...points[n] }
//    but rounded by r
//    points must go clockwise (from +y looking towards -y) around prism
//    *do not copy first point over to last
//    *"holes" are not allowed
// arotate (vec,ang)
//    rotate an object by ang degrees about arbitrary axis along vector vec
//    rotation direction follows left-handed coordinate system
// ascale (vec,scl)
//    scales an object along arbitrary axis vec by ammount scl
// reposition(vec1,vec2)
//    reorients an object so it points along vec2 instead of vec1
//    scales object according to vector length ratio
// reposition2(vec1,vec2)
//    reorients an object so it points along vec2 instead of vec1
//    does not scale objects
// arotatev (Ivec,vec,ang)
//    same as arotate but does transformations on input vector
// ascalev (Ivec,vec,scl)
//    same as ascale but does transformations on input vector
// repositionv(Ivec,vec1,vec2)
//    same as reposition but does transformations on input vector
// reposition2v(Ivec,vec1,vec2)
//    same as reposition2 but does transformations on input vector
// SphereConnect(p1,r1,p2,r2)
//    makes two spheres and connects them with a tangent cone
// TorusConnect(p1,maj1,min1,p2,maj2,min2)
//    makes two torii and connects them with a tangent cone with an inner hole
// RGB(hls)
//    converts an HLS color to RGB, where 0<=H<=360 and 0<=(l,s)<=1
// HLS(rgb)
//    converts an RGB vector to HLS
// vangle(vec1,vec2)
//    finds the angle (in radians) between two vectors



#macro Sign(Num)
   #if (Num=0) 0 #else Num/abs(Num) #end
#end

/* found this on pov-ray advanced users newsgroup :)
Larry Fontaine <lfontean@isd.net> wrote:
: I'm
: thinking a quadric will do the trick, but I need help making one.

  You need a quartic for that. Here you are:*/

//------------------------------------------------------------------------
#macro SuperCone(a,b,c,d)
  intersection
  { quartic
    { <0, 0,  0,  0, 0,   0,   0,  b*b-2*b*d+d*d, 2*(b*d-b*b), b*b,
       0,  0,   0,  0,   0,  0, 0,  0,  0, 0,
       0, 0,  0, a*a-2*a*c+c*c, 2*(a*c-a*a), a*a, 0,  0,  0, 0,
       -(a*a-2*a*c+c*c)*(b*b-2*b*d+d*d),
       -(2*((b*d-b*b)*(a*a-2*a*c+c*c)+(a*c-a*a)*(b*b-2*b*d+d*d))),
       -(b*b*(a*a-2*a*c+c*c)+4*(a*c-a*a)*(b*d-b*b)+a*a*(b*b-2*b*d+d*d)),
       -(2*(b*b*(a*c-a*a)+a*a*(b*d-b*b))), -a*a*b*b>
      sturm
    }
    cylinder { 0, z, max(max(a,b),max(c,d)) }
    bounded_by { cone { 0, max(a,b), z, max(c,d) } }
  }
#end
//------------------------------------------------------------------------


/*  It creates a cone from <0,0,0> to <0,0,1> with the ends being ellipses,
one with radiuses 'a' and 'b' and the other with radiuses 'c' and 'd'.
  Example:

camera { location -z*10 look_at 0 angle 35 }
light_source { -z*1000,1 }
light_source { y*1000,1 }

object
{ SuperCone(2, .5, .5, 2)
  pigment { rgb x } finish { specular .5 }
  translate -z*.5 scale <1,1,2>
  rotate z*90 rotate -x*45
}

  If someone is interested in the mathematics behind those quartic parameters,
I can explain (although I think nobody will ask... :) ). */

#macro supercone(a,b,p1,c,d,p2)
   #local s=min(a,b);
   #local s=min(s,c);
   #local s=min(s,d);
   #local s=s/10;
   object {
      SuperCone(a/s,b/s,c/s,d/s)
      rotate -90*x
      scale <s,1,s>
      matrix <1,0,0,p2.x-p1.x,1,p2.z-p1.z,0,0,1,0,0,0>
      scale <1,p2.y-p1.y,1>
      translate p1
   }
#end

#macro PointTo(p)
   #if (0+p.x=0 & 0+p.y=0 & 0+p.z=0)
      #local RotX=0;
   #else
      #local RotX=-atan2(p.y,sqrt(pow(p.x,2)+pow(p.z,2)))*180/pi;
   #end
   #if (0+p.x=0 & 0+p.z=0)
      #local RotY=0;
   #else
      #local RotY=atan2(p.x,p.z)*180/pi;
   #end
   rotate <RotX,RotY,0>
#end

#macro ellipsetorus(Xr,Zr,R,N)
   #local C = 0;
   blob {
   #while (C < N)
   sphere { <Xr*sin((C/N)*2*pi),0,Zr*cos((C/N)*2*pi)>,R,1 }
   #declare C=C+1;
   #end
   threshold 0.01
   }
#end

#macro roundbox(P1,P2,R)
   #local T=0;
   #local X1=P1.x;#local Y1=P1.y;#local Z1=P1.z;
   #local X2=P2.x;#local Y2=P2.y;#local Z2=P2.z;
   #if (X1>X2) #declare T=X2; #declare X2=X1; #declare X1=T; #end
   #if (Y1>Y2) #declare T=Y2; #declare Y2=Y1; #declare Y1=T; #end
   #if (Z1>Z2) #declare T=Z2; #declare Z2=Z1; #declare Z1=T; #end
   merge {
      cylinder { <X1+R,Y1+R,Z1+R>,<X1+R,Y1+R,Z2-R>,R }
      cylinder { <X1+R,Y1+R,Z1+R>,<X1+R,Y2-R,Z1+R>,R }
      cylinder { <X1+R,Y1+R,Z1+R>,<X2-R,Y1+R,Z1+R>,R }
      cylinder { <X2-R,Y1+R,Z1+R>,<X2-R,Y1+R,Z2-R>,R }
      cylinder { <X1+R,Y1+R,Z2-R>,<X2-R,Y1+R,Z2-R>,R }
      cylinder { <X2-R,Y1+R,Z1+R>,<X2-R,Y2-R,Z1+R>,R }
      cylinder { <X1+R,Y1+R,Z2-R>,<X1+R,Y2-R,Z2-R>,R }
      cylinder { <X2-R,Y1+R,Z2-R>,<X2-R,Y2-R,Z2-R>,R }
      cylinder { <X1+R,Y2-R,Z1+R>,<X1+R,Y2-R,Z2-R>,R }
      cylinder { <X1+R,Y2-R,Z1+R>,<X2-R,Y2-R,Z1+R>,R }
      cylinder { <X2-R,Y2-R,Z1+R>,<X2-R,Y2-R,Z2-R>,R }
      cylinder { <X1+R,Y2-R,Z2-R>,<X2-R,Y2-R,Z2-R>,R }
      sphere { <X1+R,Y1+R,Z1+R>,R }
      sphere { <X1+R,Y1+R,Z2-R>,R }
      sphere { <X1+R,Y2-R,Z1+R>,R }
      sphere { <X1+R,Y2-R,Z2-R>,R }
      sphere { <X2-R,Y1+R,Z1+R>,R }
      sphere { <X2-R,Y1+R,Z2-R>,R }
      sphere { <X2-R,Y2-R,Z1+R>,R }
      sphere { <X2-R,Y2-R,Z2-R>,R }
      box { <X1,Y1+R,Z1+R>,<X2,Y2-R,Z2-R> }
      box { <X1+R,Y1,Z1+R>,<X2-R,Y2,Z2-R> }
      box { <X1+R,Y1+R,Z1>,<X2-R,Y2-R,Z2> }
      bounded_by { box { <X1,Y1,Z1>,<X2,Y2,Z2> } }
   }
#end

#macro roundcylinder(p1,p2,R1,R2)
   #local H=vlength(p2-p1);
   merge {
      cylinder { <0,0,0>,<0,0,H>,R1-R2 }
      cylinder { <0,0,R2>,<0,0,H-R2>,R1 }
      torus { R1-R2,R2 rotate 90*x translate <0,0,R2> }
      torus { R1-R2,R2 rotate 90*x translate <0,0,H-R2> }
      PointTo(p2-p1)
      translate p1
      bounded_by { cylinder { p1,p2,R1 } }
   }
#end

#macro roundcone(p1,r1,p2,r2,r3)
   #local H=vlength(p2-p1);
   #local a=atan2(H,r1-r2);
   merge {
      torus { r1-tan(.5*pi-.5*a)*r3,r3 translate r3*y }
      torus { r2-tan(.5*a)*r3,r3 translate (H-r3)*y }
      cone { <0,0,0>,r1-tan(.5*pi-.5*a)*r3 <0,H,0>,r2-tan(.5*a)*r3 }
      cone { <0,r3+sin(.5*pi-a)*r3,0>,r1-tan(.5*pi-.5*a)*r3+cos(.5*pi-a)*r3
             <0,H-r3+sin(.5*pi-a)*r3,0>,r2-tan(.5*a)*r3+cos(.5*pi-a)*r3 }
      rotate 90*x
      PointTo(p2-p1)
      translate p1
   }
#end

#macro roundprism(h1,h2,p,R)
   #local ctr=0;
   #local n=dimension_size(p,1);
   #local ap=array[n]
   #local ams=array[n]
   #local ang=array[n]
   #while (ctr<n)
      #local ps=p[mod(ctr-1+n,n)];
      #local pm=p[ctr];
      #local pe=p[mod(ctr+1,n)];
      #local a1=mod(atan2(pm.x-ps.x,pm.y-ps.y)*180/pi+360,360);
      #local a2=mod(atan2(pe.x-pm.x,pe.y-pm.y)*180/pi+360,360);
      #local am=mod(a2-a1+360,360);
      #if (am<180) #local am=180-am; #end
      #local ab=(mod(a1+180,360)+a2)/2;
      #if (floor((max(mod(a1+180,360),a2)-min(mod(a1+180,360),a2))/180) != floor(am/180))
         #local ab=mod(ab+180,360); #end
      #local bv=<sin(ab*pi/180),cos(ab*pi/180)>;
      #local sam=am; #if (sam>180) #local sam=sam-180; #end
      #local l=1/sin((sam/2)*pi/180)*R;
      #local ap[ctr]=pm+bv*l;
      #local ams[ctr]=am;
      #local ang[ctr]=a1;
      #local ctr=ctr+1;
   #end
   #local ctr=0;
   prism { h1,h2,n+1            
   #while (ctr<n)
      ap[ctr]
      #local ctr=ctr+1;
   #end
   ap[0] }
   #local ctr=0;
   #while (ctr<n)
      sphere { <ap[ctr].x,h1+R,ap[ctr].y>,R }
      sphere { <ap[ctr].x,h2-R,ap[ctr].y>,R }
      cylinder { <ap[ctr].x,h1+R,ap[ctr].y>,<ap[ctr].x,h2-R,ap[ctr].y>,R }
      cylinder { <ap[ctr].x,h1+R,ap[ctr].y>,<ap[mod(ctr+1,n)].x,h1+R,ap[mod(ctr+1,n)].y>,R }
      cylinder { <ap[ctr].x,h2-R,ap[ctr].y>,<ap[mod(ctr+1,n)].x,h2-R,ap[mod(ctr+1,n)].y>,R }
      box { <-R,h1+R,0>,<0,h2-R,sqrt(pow(ap[mod(ctr+1,n)].x-ap[ctr].x,2)
                                    +pow(ap[mod(ctr+1,n)].y-ap[ctr].y,2))>
            rotate atan2(ap[mod(ctr+1,n)].x-ap[ctr].x,ap[mod(ctr+1,n)].y-ap[ctr].y)*180/pi*y
            translate <ap[ctr].x,0,ap[ctr].y> }
      #if (ams[ctr]>180)
         difference {
            intersection {
               cylinder { <p[ctr].x,h1-R,p[ctr].y>*2-<ap[ctr].x,h1-R,ap[ctr].y>,
                          <p[ctr].x,h2+R,p[ctr].y>*2-<ap[ctr].x,h2+R,ap[ctr].y>,R inverse }
               prism { h1,h2,4 ap[ctr],ap[mod(ctr-1+n,n)],ap[mod(ctr+1,n)],ap[ctr] }
               plane { -z,0 rotate ang[ctr]*y translate <p[ctr].x,0,p[ctr].y>*2-<ap[ctr].x,0,ap[ctr].y> }
               plane { -z,0 rotate (180+ang[mod(ctr+1,n)])*y translate <p[ctr].x,0,p[ctr].y>*2-<ap[ctr].x,0,ap[ctr].y> }
            }
            intersection {
               torus { R*2,R translate <p[ctr].x,h1+R,p[ctr].y>*2-<ap[ctr].x,h1+R,ap[ctr].y> inverse }
               cylinder { <p[ctr].x,h1+R,p[ctr].y>*2-<ap[ctr].x,h1+R,ap[ctr].y>,
                          <p[ctr].x,h1-R*2,p[ctr].y>*2-<ap[ctr].x,h1-R*2,ap[ctr].y>,R*2 }
            }
            intersection {
               torus { R*2,R translate <p[ctr].x,h2-R,p[ctr].y>*2-<ap[ctr].x,h2-R,ap[ctr].y> inverse }
               cylinder { <p[ctr].x,h2-R,p[ctr].y>*2-<ap[ctr].x,h2-R,ap[ctr].y>,
                          <p[ctr].x,h2+R*2,p[ctr].y>*2-<ap[ctr].x,h2+R*2,ap[ctr].y>,R*2 }
            }
         }
      #end
      #local ctr=ctr+1;
   #end
#end

#macro arotate(vec,ang)
   #if (vlength(vec)=0)
      #render "Warning: arotate vector length is zero.\n"
   #end
   #if (vec.x=0 & vec.z=0)
      #local RotY=0;
   #else
      #local RotY=atan2(vec.x,vec.z)*180/pi;
   #end
   #local RotX=asin((vec.y)/vlength(vec))*180/pi;
   rotate -RotY*y
   rotate RotX*x
   rotate ang*z
   rotate <-RotX,RotY,0>
#end

#macro ascale(vec,scl)
   #if (vlength(vec)=0)
      #render "Warning: ascale vector length is zero.\n"
   #end
   #if (vec.x=0 & vec.z=0)
      #local RotY=0;
   #else
      #local RotY=atan2(vec.x,vec.z)*180/pi;
   #end
   #local RotX=asin((vec.y)/vlength(vec))*180/pi;
   rotate -RotY*y
   rotate RotX*x
   scale <1,1,scl>
   rotate <-RotX,RotY,0>
#end

#macro reposition(vec1,vec2)
   #if (vlength(vec1)=0)
      #render "Warning: reposition initial vector length is zero.\n"
   #end
   #if (vec1.x=0 & vec1.z=0)
      #local RotY=0;
   #else
      #local RotY=atan2(vec1.x,vec1.z)*180/pi;
   #end
   #local RotX=asin((vec1.y)/vlength(vec1))*180/pi;
   rotate -RotY*y
   rotate RotX*x
   scale vlength(vec2)/vlength(vec1)
   #if (vlength(vec2)=0)
      #render "Warning: reposition final vector length is zero.\n"
   #end
   #if (vec2.x=0 & vec2.z=0)
      #local RotY=0;
   #else
      #local RotY=atan2(vec2.x,vec2.z)*180/pi;
   #end
   #local RotX=asin((vec2.y)/vlength(vec2))*180/pi;
   rotate <-RotX,RotY,0>
#end

#macro reposition2(vec1,vec2)
   #if (vlength(vec1)=0)
      #render "Warning: reposition2 initial vector length is zero.\n"
   #end
   #if (vec1.x=0 & vec1.z=0)
      #local RotY=0;
   #else
      #local RotY=atan2(vec1.x,vec1.z)*180/pi;
   #end
   #local RotX=asin((vec1.y)/vlength(vec1))*180/pi;
   rotate -RotY*y
   rotate RotX*x
   #if (vlength(vec2)=0)
      #render "Warning: reposition2 final vector length is zero.\n"
   #end
   #if (vec2.x=0 & vec2.z=0)
      #local RotY=0;
   #else
      #local RotY=atan2(vec2.x,vec2.z)*180/pi;
   #end
   #local RotX=asin((vec2.y)/vlength(vec2))*180/pi;
   rotate <-RotX,RotY,0>
#end

#macro arotatev(InpV,vec,ang)
   #if (vlength(vec)=0)
      #render "Warning: arotatev vector length is zero.\n"
   #end
   #if (vec.x=0 & vec.z=0)
      #local RotY=0;
   #else
      #local RotY=atan2(vec.x,vec.z)*180/pi;
   #end
   #local RotX=asin((vec.y)/vlength(vec))*180/pi;
   #local InpV=vrotate (InpV,-RotY*y);
   #local InpV=vrotate (InpV,RotX*x);
   #local InpV=vrotate (InpV,ang*z);
   #local InpV=vrotate (InpV,<-RotX,RotY,0>);
   InpV
#end

#macro ascalev(InpV,vec,scl)
   #if (vlength(vec)=0)
      #render "Warning: ascalev vector length is zero.\n"
   #end
   #if (vec.x=0 & vec.z=0)
      #local RotY=0;
   #else
      #local RotY=atan2(vec.x,vec.z)*180/pi;
   #end
   #local RotX=asin((vec.y)/vlength(vec))*180/pi;
   #local InpV=vrotate (InpV,-RotY*y);
   #local InpV=vrotate (InpV,RotX*x);      
   #local InpV=<InpV.x,InpV.y,InpV.z*scl>;
   #local InpV=vrotate (InpV,<-RotX,RotY,0>);
   InpV
#end

#macro repositionv(InpV,vec1,vec2)
   #if (vlength(vec1)=0)
      #render "Warning: repositionv initial vector length is zero.\n"
   #end
   #if (vec1.x=0 & vec1.z=0)
      #local RotY=0;
   #else
      #local RotY=atan2(vec1.x,vec1.z)*180/pi;
   #end
   #local RotX=asin((vec1.y)/vlength(vec1))*180/pi;
   #local InpV=vrotate (InpV,-RotY*y);
   #local InpV=vrotate (InpV,RotX*x);
   #local InpV=InpV*vlength(vec2)/vlength(vec1);
   #if (vlength(vec2)=0)
      #render "Warning: repositionv final vector length is zero.\n"
   #end
   #if (vec2.x=0 & vec2.z=0)
      #local RotY=0;
   #else
      #local RotY=atan2(vec2.x,vec2.z)*180/pi;
   #end
   #local RotX=asin((vec2.y)/vlength(vec2))*180/pi;
   #local InpV=vrotate (InpV,<-RotX,RotY,0>);
   InpV
#end

#macro reposition2v(InpV,vec1,vec2)
   #if (vlength(vec1)=0)
      #render "Warning: reposition2v initial vector length is zero.\n"
   #end
   #if (vec1.x=0 & vec1.z=0)
      #local RotY=0;
   #else
      #local RotY=atan2(vec1.x,vec1.z)*180/pi;
   #end
   #local RotX=asin((vec1.y)/vlength(vec1))*180/pi;
   #local InpV=vrotate (InpV,-RotY*y);
   #local InpV=vrotate (InpV,RotX*x);
   #if (vlength(vec2)=0)
      #render "Warning: reposition2v final vector length is zero.\n"
   #end
   #if (vec2.x=0 & vec2.z=0)
      #local RotY=0;
   #else
      #local RotY=atan2(vec2.x,vec2.z)*180/pi;
   #end
   #local RotX=asin((vec2.y)/vlength(vec2))*180/pi;
   #local InpV=vrotate (InpV,<-RotX,RotY,0>);
   InpV
#end

#macro SphereConnect(p1,r1,p2,r2)
   merge {
      sphere { p1,r1 }
      sphere { p2,r2 }
      #local Dist = vlength(p2-p1);
      #local Ang = asin((r1-r2)/Dist);
      cone { p1+vnormalize(p2-p1)*r1*sin(Ang),r1*cos(Ang)
             p2+vnormalize(p2-p1)*r2*sin(Ang),r2*cos(Ang) }
   }
#end

#macro TorusConnect(p1,maj1,min1,p2,maj2,min2)
   #local Dist=vlength(p2-p1);
   #local Dir=vnormalize(p2-p1);
   #local Ang1=atan2((maj1-maj2)/Dist,1);
   #local Ang2=asin((min1-min2)/sqrt(pow(Dist,2)+pow(maj1-maj2,2)));
   #local OWid1 = maj1+min1*cos(Ang1+Ang2);
   #local OWid2 = maj2+min2*cos(Ang1+Ang2);
   #local ODst1 = min1*sin(Ang1+Ang2);
   #local ODst2 = min2*sin(Ang1+Ang2);
   #local IWid1 = maj1-min1*cos(Ang1-Ang2);
   #local IWid2 = maj2-min2*cos(Ang1-Ang2);
   #local IDst1 = min1*sin(Ang1-Ang2);
   #local IDst2 = min2*sin(Ang1-Ang2);
   #local ConP1 = p1+Dir*(ODst1-(ODst1+IDst1)*(OWid1/(OWid1-IWid1)));
   #local ConP2 = p2+Dir*(ODst2-(ODst2+IDst2)*(OWid2/(OWid2-IWid2)));
   #local ConP3 = p1+Dir*ODst1;
   #local ConP4 = p2+Dir*ODst2;
   #if (vdot(vnormalize(ConP1-ConP3),Dir)>0) #local Dir1=1; #else #local Dir1=-1; #end
   #if (vdot(vnormalize(ConP2-ConP4),Dir)>0) #local Dir2=1; #else #local Dir2=-1; #end
   merge {
      difference {
         cone { p1-Dir*min1,OWid1+(OWid1-OWid2)*((min1+ODst1)/(Dist-ODst1+ODst2))
                p2+Dir*min2,OWid2-(OWid1-OWid2)*((min2-ODst2)/(Dist-ODst1+ODst2)) }
         cone { p1-Dir*(min1+.0000001),IWid1+(IWid1-IWid2)*((min1-IDst1)/(Dist+IDst1-IDst2))
                p2+Dir*(min2+.0000001),IWid2-(IWid1-IWid2)*((min2+IDst2)/(Dist+IDst1-IDst2)) }
         #if (Dir1=-1)
            difference {
               plane { Dir,-.0000001 translate ConP3 }
               cone { ConP3,OWid1 ConP1,0 }
            }
         #end
         #if (Dir1=1)
            cone { ConP3,OWid1 ConP1,0 }
            plane { Dir,.0000001 translate ConP3 }
         #end
         #if (Dir2=-1)
            cone { ConP4,OWid2 ConP2,0 }
            plane { -Dir,.0000001 translate ConP4 }
         #end
         #if (Dir2=1)
            difference {
                plane { -Dir,-.0000001 translate ConP4)
                cone { ConP4,OWid2 ConP2,0 }
            }
         #end
      }
      torus { maj1,min1 reposition(y,Dir) translate p1 }
      torus { maj2,min2 reposition(y,Dir) translate p2 }
   }
#end

#macro RGB(hls)
   #local hue=hls.x/60;
   #local lig=hls.y;
   #local sat=hls.z;
   #local MixColor=<1,0,6-hue>;
   #if (hue<5) #local MixColor=<hue-4,0,1>; #end
   #if (hue<4) #local MixColor=<0,4-hue,1>; #end
   #if (hue<3) #local MixColor=<0,1,hue-2>; #end
   #if (hue<2) #local MixColor=<2-hue,1,0>; #end
   #if (hue<1) #local MixColor=<1,hue,0>; #end
   #if (lig<.5)
      #local MixColor=MixColor*(lig*2)+<0,0,0>*(1-lig*2);
   #end #if (lig>.5)
      #local MixColor=MixColor*(2-lig*2)+<1,1,1>*(lig*2-1);
   #end
   #local MixColor=MixColor*sat+lig*<1,1,1>*(1-sat);
   MixColor
#end

#macro HLS(RGBv)
   #local Max=max(max(RGBv.x,RGBv.y),RGBv.z);
   #local Min=min(min(RGBv.x,RGBv.y),RGBv.z);
   #local Lig=(Min+Max)/2;
   #local Sat=(Min-Max)/(2*abs(Lig-.5)-1);
   #if (RGBv.x=Max)
      #if (RGBv.y=Min)
         #local Hue=5;
         #local Med=RGBv.z;
      #end #if (RGBv.z=Min)
         #local Hue=0;
         #local Med=RGBv.y;
      #end
   #end #if (RGBv.y=Max)
      #if (RGBv.x=Min)
         #local Hue=2;
         #local Med=RGBv.z;
      #end #if (RGBv.z=Min)
         #local Hue=1;
         #local Med=RGBv.x;
      #end
   #end #if (RGBv.z=Max)
      #if (RGBv.x=Min)
         #local Hue=3;
         #local Med=RGBv.y;
      #end #if (RGBv.y=Min)
         #local Hue=4;
         #local Med=RGBv.x;
      #end
   #end
   #local Hue=Hue+mod(Hue,2)-2*(mod(Hue,2)-0.5)*(Med-Min)/(Max-Min);
   <Hue*60,Lig,Sat>
#end

#macro vangle(vec1,vec2)
   #local len1 = vlength(vec1);
   #local len2 = vlength(vec2);
   #local len3 = vlength(vec2-vec1);
   acos((pow(len3,2)-pow(len1,2)-pow(len2,2))/(-2*len1*len2))
#end