// Macros include file

// 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 (X1,Y1,Z1,X2,Y2,Z2,R)
//    creates box from <X1,Y1,Z1> to <X2,Y2,Z2> and rounds edges by R
// roundcylinder (H,R1,R2)
//    Rounded Cylinder from Y=0 to Y=H with radius R1 and rounded by R2
// roundcylinder2 (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
//    *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



/* 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,x1,y1,z1,c,d,x2,y2,z2)
   #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,x2-x1,1,z2-z1,0,0,1,0,0,0>
      scale <1,y2-y1,1>
      translate <x1,y1,z1>
   }
#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 (X1,Y1,Z1,X2,Y2,Z2,R)
   #local T=0;
   #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 (H,R1,R2)
   merge {
      cylinder { <0,0,0>,<0,H,0>,R1-R2 }
      cylinder { <0,R2,0>,<0,H-R2,0>,R1 }
      torus { R1-R2,R2 translate <0,R2,0> }
      torus { R1-R2,R2 translate <0,H-R2,0> }
      bounded_by { cylinder { 0,H*y,R1 } }
   }
#end

#macro roundcylinder2 (p1,p2,R1,R2)
   #local H=sqrt(pow(p1.x-p2.x,2)+pow(p1.y-p2.y,2)+pow(p1.z-p2.z,2));
   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=sqrt(pow(p1.x-p2.x,2)+pow(p1.y-p2.y,2)+pow(p1.z-p2.z,2));
   #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 (sqrt(pow(vec.x,2)+pow(vec.y,2)+pow(vec.z,2))=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)/sqrt(pow(vec.x,2)+pow(vec.y,2)+pow(vec.z,2)))*180/pi;
   rotate -RotY*y
   rotate RotX*x
   rotate ang*z
   rotate <-RotX,RotY,0>
#end

#macro ascale(vec,scl)
   #if (sqrt(pow(vec.x,2)+pow(vec.y,2)+pow(vec.z,2))=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)/sqrt(pow(vec.x,2)+pow(vec.y,2)+pow(vec.z,2)))*180/pi;
   rotate -RotY*y
   rotate RotX*x
   scale <1,1,scl>
   rotate <-RotX,RotY,0>
#end

#macro reposition(vec1,vec2)
   #if (sqrt(pow(vec1.x,2)+pow(vec1.y,2)+pow(vec1.z,2))=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)/sqrt(pow(vec1.x,2)+pow(vec1.y,2)+pow(vec1.z,2)))*180/pi;
   rotate -RotY*y
   rotate RotX*x
   scale sqrt(pow(vec2.x,2)+pow(vec2.y,2)+pow(vec2.z,2))/sqrt(pow(vec1.x,2)+pow(vec1.y,2)+pow(vec1.z,2))
   #if (sqrt(pow(vec2.x,2)+pow(vec2.y,2)+pow(vec2.z,2))=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)/sqrt(pow(vec2.x,2)+pow(vec2.y,2)+pow(vec2.z,2)))*180/pi;
   rotate <-RotX,RotY,0>
#end