import three;
import graph3;

pen defaultbackpen=linetype("4 4",4,scale=false);

// A solid geometry package.

// Try to find a bounding tangent line between two paths.
real[] tangent(path p, path q, bool side) 
{
  static real fuzz=1.0e-5;

  if((cyclic(p) && inside(p,point(q,0)) || 
      cyclic(q) && inside(q,point(p,0))) &&
     intersect(p,q,fuzz).length == 0) return new real[];

  for(int i=0; i < 100; ++i) {
    real ta=side ? mintimes(p)[1] : maxtimes(p)[1];
    real tb=side ? mintimes(q)[1] : maxtimes(q)[1];
    pair a=point(p,ta);
    pair b=point(q,tb);
    real angle=angle(b-a,warn=false);
    if(abs(angle) <= sqrtEpsilon || abs(abs(0.5*angle)-pi) <= sqrtEpsilon)
      return new real[] {ta,tb};
    transform t=rotate(-degrees(angle));
    p=t*p;
    q=t*q;
  }
  return new real[];
}

path line(path p, path q, real[] t) 
{
  return point(p,t[0])--point(q,t[1]);
}

// Return the projection of a generalized cylinder of height h constructed
// from area base in the XY plane and aligned with axis.
path[] cylinder(path3 base, real h, triple axis=Z,
                projection P=currentprojection) 
{
  base=rotate(-colatitude(axis),cross(axis,Z))*base;
  path3 top=shift(h*axis)*base;
  path Base=project(base,P);
  path Top=project(top,P);
  real[] t1=tangent(Base,Top,true);
  real[] t2=tangent(Base,Top,false);
  path p=subpath(Base,t1[0]/P.ninterpolate,t2[0]/P.ninterpolate);
  path q=subpath(Base,t2[0]/P.ninterpolate,t1[0]/P.ninterpolate);
  return Base^^Top^^line(Base,Top,t1)^^line(Base,Top,t2);
}

// The three-dimensional "wireframe" used to visualize a volume of revolution
struct skeleton {
  struct curve {
    path3[] front;
    path3[] back;
  }
  // transverse skeleton (perpendicular to axis of revolution)
  curve transverse;
  // longitudinal skeleton (parallel to axis of revolution)
  curve longitudinal;
}

// A surface of revolution generated by rotating a planar path3 g
// from angle1 to angle2 about c--c+axis.
struct revolution {
  triple c;
  path3 g;
  triple axis;
  real angle1,angle2;
  triple M;
  triple m;

  static real epsilon=10*sqrtEpsilon;
  
  void operator init(triple c=O, path3 g, triple axis=Z, real angle1=0,
                     real angle2=360) {
    this.c=c;
    this.g=g;
    this.axis=unit(axis);
    this.angle1=angle1;
    this.angle2=angle2;
    M=max(g);
    m=min(g);
  }
  
  // Return the surface of rotation obtain by rotating the path3 (x,0,f(x))
  // sampled n times between x=a and x=b about an axis lying in the XZ plane.
  void operator init(triple c=O, real f(real x), real a, real b, int n=ngraph,
                     interpolate3 join=operator --, triple axis=Z,
                     real angle1=0, real angle2=360) {
    operator init(c,graph(new triple(real x) {return (x,0,f(x));},a,b,n,
                          join),axis,angle1,angle2);
  }

  revolution copy() {
    return revolution(c,g,axis,angle1,angle2);
  }
  
  triple vertex(int i, real j) {
    triple v=point(g,i);
    triple center=c+dot(v-c,axis)*axis;
    triple perp=v-center;
    triple normal=cross(axis,perp);
    return center+Cos(j)*perp+Sin(j)*normal;
  }

  // Construct the surface of rotation generated by rotating g
  // from angle1 to angle2 sampled n times about the line c--c+axis.
  // An optional surface pen color(int i, real j) may be specified
  // to override the color at vertex(i,j).
  surface surface(int n=nslice, pen color(int i, real j)=null) {
    return surface(c,g,axis,n,angle1,angle2,color);
  }

  path3 slice(real position, int n=nCircle) {
    triple v=point(g,position);
    triple center=c+dot(v-c,axis)*axis;
    triple perp=v-center;
    if(abs(perp) <= epsilon*max(abs(m),abs(M))) return center;
    triple v1=center+rotate(angle1,axis)*perp;
    triple v2=center+rotate(angle2,axis)*perp;
    path3 p=Arc(center,v1,v2,axis,n);
    return (angle2-angle1) % 360 == 0 ? p&cycle : p;
  }
  
  triple camera(projection P) {
    triple camera=P.camera;
    if(P.infinity) {
      real s=abs(M-m)+abs(m-P.target);
      camera=P.target+camerafactor*s*unit(P.vector());
    }
    return camera;
  }

  // add transverse slice to skeleton s;
  // must be recomputed if camera is adjusted
  void transverse(skeleton s, real t, int n=nslice,
		  projection P=currentprojection) {
    skeleton.curve s=s.transverse;
    path3 S=slice(t,n);
    triple camera=camera(P);
    int L=length(g);
    real midtime=0.5*L;
    real sign=sgn(dot(axis,camera-P.target))*sgn(dot(axis,dir(g,midtime)));
    if(dot(M-m,axis) == 0 || (t <= epsilon && sign < 0) ||
       (t >= L-epsilon && sign > 0))
      s.front.push(S);
    else {
      path3 Sp=slice(t+epsilon,n);
      path3 Sm=slice(t-epsilon,n);
      path sp=project(Sp,P);
      path sm=project(Sm,P);
      real[] t1=tangent(sp,sm,true);
      real[] t2=tangent(sp,sm,false);
      if(t1.length > 1 && t2.length > 1) {
        real t1=t1[0]/P.ninterpolate;
        real t2=t2[0]/P.ninterpolate;
        int len=length(S);
        if(t2 < t1) {
          real temp=t1;
          t1=t2;
          t2=temp;
        }
        path3 p1=subpath(S,t1,t2);
        path3 p2=subpath(S,t2,len);
        path3 P2=subpath(S,0,t1);
        if(abs(midpoint(p1)-camera) <= abs(midpoint(p2)-camera)) {
          s.front.push(p1);
          if(cyclic(S))
            s.back.push(p2 & P2);
          else {
            s.back.push(p2);
            s.back.push(P2);
          }
        } else {
          if(cyclic(S))
            s.front.push(p2 & P2);
          else {
            s.front.push(p2);
            s.front.push(P2);
          }
          s.back.push(p1);
        }
      } else {
        if((t <= midtime && sign < 0) || (t >= midtime && sign > 0))
          s.front.push(S);
        else
          s.back.push(S);
      }
    }
  }

  // add m evenly spaced transverse slices to skeleton s
  void transverse(skeleton s, int m=0, int n=nslice,
		  projection P=currentprojection) {
    if(m == 0) {
      int N=size(g);
      for(int i=0; i < N; ++i)
	transverse(s,(real) i,n,P);
    } else if(m == 1)
      transverse(s,reltime(g,0.5),n,P);
    else {
      real factor=1/(m-1);
      for(int i=0; i < m; ++i)
	transverse(s,reltime(g,i*factor),n,P);
    }
  }

  // return approximate silhouette based on m evenly spaced transverse slices;
  // must be recomputed if camera is adjusted
  path3[] silhouette(int m=64, projection P=currentprojection) {
    if(is3D())
      warning("2Dsilhouette",
              "silhouette routine is intended only for 2d projections");
    path3 G,H;
    int N=size(g);
    int M=(m == 0) ? N : m;
    real factor=m == 1 ? 0 : 1/(m-1);
    int n=nslice;
    
    real tfirst=-1;
    real tlast;
    for(int i=0; i < M; ++i) {
      real t=(m == 0) ? i : reltime(g,i*factor);
      path3 S=slice(t,n);
      triple camera=camera(P);
      path3 Sp=slice(t+epsilon,n);
      path3 Sm=slice(t-epsilon,n);
      path sp=project(Sp,P);
      path sm=project(Sm,P);
      real[] t1=tangent(sp,sm,true);
      real[] t2=tangent(sp,sm,false);
      if(t1.length > 1 && t2.length > 1) {
	real t1=t1[0]/P.ninterpolate;
	real t2=t2[0]/P.ninterpolate;
	if(t1 != t2) {
	  G=G..point(S,t1);
	  H=point(S,t2)..H;
	  if(tfirst < 0) tfirst=t;
	  tlast=t;
	}
      }
    }
    int L=length(g);
    real midtime=0.5*L;
    triple camera=camera(P);
    real sign=sgn(dot(axis,camera-P.target))*sgn(dot(axis,dir(g,midtime)));

    skeleton sfirst;
    transverse(sfirst,tfirst,n,P);
    triple delta=this.M-this.m;
    path3 cap;
    if(dot(delta,axis) == 0 || (tfirst <= epsilon && sign < 0)) {
      cap=sfirst.transverse.front[0];
    } else {
      if(sign > 0) {
	if(sfirst.transverse.front.length > 0)
	  G=reverse(sfirst.transverse.front[0])..G;
      } else {
      if(sfirst.transverse.back.length > 0)
	G=sfirst.transverse.back[0]..G;
      }
    }
    
    skeleton slast;
    transverse(slast,tlast,n,P);
    if(dot(delta,axis) == 0 || (tlast >= L-epsilon && sign > 0)) {
      cap=slast.transverse.front[0];
    } else {
      if(sign > 0) {
	if(slast.transverse.back.length > 0)
	  H=reverse(slast.transverse.back[0])..H;
      } else {
	if(slast.transverse.front.length > 0)
	  H=slast.transverse.front[0]..H;
      }
    }

    return size(cap) == 0 ? G^^H : G^^H^^cap;
  }

  // add longitudinal curves to skeleton;
  // must be recomputed if camera is adjusted
  void longitudinal(skeleton s, int n=nslice, projection P=currentprojection) {
    real t, d=0;
    // Find a point on g of maximal distance from the axis.
    int N=size(g);
    for(int i=0; i < N; ++i) {
      triple v=point(g,i);
      triple center=c+dot(v-c,axis)*axis;
      real r=abs(v-center);
      if(r > d) {
        t=i;
        d=r;
      }
    }
    path3 S=slice(t,n);
    path3 Sm=slice(t+epsilon,n);
    path3 Sp=slice(t-epsilon,n);
    path sp=project(Sp,P);
    path sm=project(Sm,P);
    real[] t1=tangent(sp,sm,true);
    real[] t2=tangent(sp,sm,false);
    transform3 T=transpose(align(axis));
    real Longitude(triple v) {return longitude(T*(v-c),warn=false);}
    real ref=Longitude(point(g,t));
    real angle(real t) {return Longitude(point(S,t/P.ninterpolate))-ref;}
    triple camera=camera(P);
    void push(real[] T) {
      if(T.length > 1) {
	path3 p=rotate(angle(T[0]),c,c+axis)*g;
	path3 p1=subpath(p,0,t);
	path3 p2=subpath(p,t,length(p));
	if(length(p1) > 0 &&
	   abs(midpoint(p1)-camera) <= abs(midpoint(p2)-camera)) {
	  s.longitudinal.front.push(p1);
          s.longitudinal.back.push(p2);
	} else {
	  s.longitudinal.back.push(p1);
	  s.longitudinal.front.push(p2);
	}
      }
    }
    push(t1);
    push(t2);
  }
  
  skeleton skeleton(int m=0, int n=nslice, projection P=currentprojection) {
    skeleton s;
    transverse(s,m,n,P);
    longitudinal(s,n,P);
    return s;
  }
}

surface surface(revolution r, int n=nslice, pen color(int i, real j)=null)
{
  return r.surface(n,color);
}

// Draw on picture pic the skeleton of the surface of revolution r.
// Draw the front portion of each of the m transverse slices with pen p and
// the back portion with pen backpen. Rotational arcs are based on
// n-point approximations to the unit circle.
void draw(picture pic=currentpicture, revolution r, int m=0, int n=nslice,
	  pen frontpen=currentpen, pen backpen=frontpen,
	  pen longitudinalpen=frontpen, pen longitudinalbackpen=backpen,
	  light light=currentlight, projection P=currentprojection)
{
  pen thin=is3D() ? thin() : defaultpen;
  skeleton s=r.skeleton(m,n,P);
  begingroup3(pic);
  if(frontpen != nullpen) {
    draw(pic,s.transverse.back,thin+defaultbackpen+backpen,light);
    draw(pic,s.transverse.front,thin+frontpen,light);
  }
  if(longitudinalpen != nullpen) {
    draw(pic,s.longitudinal.back,thin+defaultbackpen+longitudinalbackpen,light);
    draw(pic,s.longitudinal.front,thin+longitudinalpen,light);
  }
  endgroup3(pic);
}

revolution operator * (transform3 t, revolution r)
{
  triple trc=t*r.c;
  return revolution(trc,t*r.g,t*(r.c+r.axis)-trc,r.angle1,r.angle2);
}

// Return a right circular cylinder of height h in the direction of axis
// based on a circle centered at c with radius r.
revolution cylinder(triple c=O, real r, real h, triple axis=Z)
{
  triple C=c+r*perp(axis);
  axis=h*unit(axis);
  return revolution(c,C--C+axis,axis);
}

// Return a right circular cone of height h in the direction of axis
// based on a circle centered at c with radius r. The parameter n
// controls the accuracy near the degenerate point at the apex.
revolution cone(triple c=O, real r, real h, triple axis=Z, int n=nslice)
{
  axis=unit(axis);
  return revolution(c,approach(c+r*perp(axis)--c+h*axis,n),axis);
}

// Return an approximate sphere of radius r centered at c obtained by rotating
// an (n+1)-point approximation to a half circle about the Z axis.
// Note: unitsphere provides a smoother and more efficient surface.
revolution sphere(triple c=O, real r, int n=nslice)
{
  return revolution(c,Arc(c,r,180,0,0,0,Y,n),Z);
}