//                                                                      
// Mobius for POVRay 3.5 Mesh2 Object v3.5                              
// Trevor G Quayle, July 18, 2002                                       
//                                                                      
// This is a mobius strip made from mesh2. The strip twists around a    
// neutral-axis defined by a closed ellipse on the x-z plane centred on 
// the origin. The mobius can have either round or square edges.        
// UV mapping can now be used.                                          
//                                                                      
// syntax: Mobius(Ax,d,s,nR,SA,SM,UV,RE)                                
// Ax = <x,z> half-length of x-axis and z axis                          
// d  = <W,T> Width/thickness  of strip (centred on N.A.)               
// s  = <m,n,o> # segments in circumference (m), width (n) &            
//              thickness (o)                                           
// nR = # half-twists of the strip (0=no twist, 1=180deg, etc)          
//     (negative #s will reverse the direction of the twist)            
// SA = angle of twisting at +x axis                                    
// SM = smoothing (yes/no)                                              
// UV = uv_mapping (yes/no)                                             
// RE = round edges (yes/no)                                            
//                                                                      
// ex:                                                                  
// object{Mobius(<20,10>,4,2,<6,50>,1,45,yes,yes,yes) material{MyMat}}  
//                                                                      
// creates a smooth mobius 40x20 units in size, 4 units wide, 2 units   
// thick with one twist (traditional mobius) and round edges. The twist 
// starts at 45° on the +x axis. There are 6 mesh segments in the width 
// and 50 around the circumference of the mobius.                       
// Using uv_mapping: The mesh uses a texture defined from <0,0,0> to    
//                   <2,2,0>. The top and bottom faces use the texture  
//                   from <0,0,0> to <1,2,0> and the side faces use the 
//                   texture from <1,2,0> to <2,2,0>. If there is an    
//                   odd number of twists, there is only one face, and  
//                   the texture will run smoothly around the entire    
//                   face. With an even number of twists, the texture   
//                   will be split between the top and bottom faces and 
//                   likewise for the sides.                            
//                                                                      

#macro Mobius(Ax,d,s,nR,SA,SM,UV,RE)
  #local W=d.u;#local T=d.v;
  #local Lx=Ax.u;#local Lz=Ax.v;
  #local m=int(s.x);#local n=int(s.y);#local o=int(s.z);
  #local UVS=<1/m,1/n,1/o>;
  #if (Lx<=0) #debug "X-axis must be  > 0, assigning 2*W" #local Lx=2*W; #end
  #if (Lz<=0) #debug "Z-axis must be  > 0, assigning 2*W" #local Lz=2*W; #end
  #if (m<1) #debug "m must be an integer > 0, assigning 1\n" #local m=1; #end
  #if (n<4) #debug "n must be an integer > 3, assigning 4\n" #local n=4; #end
  #if (min(Lx,Lz)<W) #debug "A minimum axis smaller than the width may produce unpredictable results!\n" #end
  #if (mod(nR,1)!=0) #debug "nR is not an integer. Strip is not closed!\n" #end
  #macro PE(A,B,C) <B*cos(A),0,C*sin(A)> #end
  #macro NE(P,B,C) vnormalize(<C*P.x,0,pow(B,2)*sqrt(1-pow(P.x/B,2))#if(P.z<0)*-1#end>) #end//ok
  #macro NX(i,j,a,b) vnormalize(vcross(MP[i][b]-MP[i][a],MP[i+1][j]-MP[i-1][j])) #end
  #local MP=array[(m+3)*2][n+o+2]
  #local sgn=0;#while(sgn<=1)
    #local i=0;#while (i<(m+3))
      #local Rd=((i-1)*nR*180/m+SA+180*sgn);
      #local A=radians((i-1)*360/m);
      #local P=PE(A,Lx,Lz);
      #local Nu=NE(P,Lx,Lz);
      #local Nv=vcross(Nu,y);
      #local Pu=(vaxis_rotate(Nu,Nv,Rd));
      #local Pv=(vaxis_rotate( y,Nv,Rd));
      #local j=-n/2;#while (j<=n/2)
        #local MP[i+(m+3)*sgn][j+n/2]=P+Pu*W/n*j+Pv*T/2;
      #local j=j+1;#end
      #local k=-o/2;#while (k<=o/2)
        #if (RE)
          #local MP[i+(m+3)*sgn][(n+1)+k+o/2]=P+vaxis_rotate(Pv*T/2,Nv,k*180/o-90)+Pu*W/2;
        #else
          #local MP[i+(m+3)*sgn][(n+1)+k+o/2]=P+Pv*T/o*k+Pu*W/2;
        #end
      #local k=k+1;#end
    #local i=i+1;#end
  #local sgn=sgn+1;#end
  mesh2{
    vertex_vectors{(m+1)*2*(n+o+2)
      #local k=0;#while(k<=1)
        #local i=1;#while(i<=(m+1))
          #local j=0;#while (j<=n)
            ,MP[i+k*(m+3)][j]
          #local j=j+1;#end
        #local i=i+1;#end
      #local k=k+1;#end
      #local k=0;#while(k<=1)
        #local i=1;#while(i<=(m+1))
          #local j=0;#while (j<=o)
            ,MP[i+k*(m+3)][j+(n+1)]
          #local j=j+1;#end
        #local i=i+1;#end
      #local k=k+1;#end
    }
    #if (SM)
      normal_vectors{(m+1)*2*(n+o+2)
        #local k=0;#while(k<=1)
          #local i=1;#while(i<=(m+1))
            #local j=0;#while (j<=n)
              ,-NX(i+k*(m+3),j,0,n)
            #local j=j+1;#end
          #local i=i+1;#end
        #local k=k+1;#end
        #local k=0;#while(k<=1)
          #local i=1;#while(i<=(m+1))
            #local j=0;#while (j<=o)
            #if (RE)
              #local Nu=NX(i+k*(m+3),j+(n+1),(n+1),o+(n+1));
              #local Nv=vnormalize(MP[i+k*(m+3)-1][j+(n+1)]-MP[i+k*(m+3)+1][j+(n+1)]);
              #local NN=vaxis_rotate(Nu,Nv,(o-j)/o*180-90);
              ,NN
            #else
              ,-NX(i+k*(m+3),j+(n+1),n+1,n+o)
            #end
            #local j=j+1;#end
          #local i=i+1;#end
        #local k=k+1;#end
      }
    #end
    #if (UV)
      uv_vectors{(m+1)*2*(n+o+2)
        #local k=0;#while(k<=1)
          #local i=0;#while(i<=m)
            #local j=0;#while (j<=n)
              ,<k+i*UVS.x,j*UVS.y>
            #local j=j+1;#end
          #local i=i+1;#end
        #local k=k+1;#end
        #local k=0;#while(k<=1)
          #local i=0;#while(i<=m)
            #local j=0;#while (j<=o)
              ,<k/2+i*UVS.x,1+j*UVS.z>
            #local j=j+1;#end
          #local i=i+1;#end
        #local k=k+1;#end
      }
    #end
    face_indices{m*(n+o)*2*2
      #local k=0;#while(k<=1)
        #local i=0;#while(i<m)
          #local j=0;#while (j<n)
            #local T1=<j+(n+1)*i,j+(n+1)*(i+1),(j+1)+(n+1)*i>+k*(n+1)*(m+1);
            ,T1,T1+x*(n+2)
          #local j=j+1;#end
        #local i=i+1;#end
      #local k=k+1;#end
      #local k=0;#while(k<=1)
        #local i=0;#while(i<m)
          #local j=0;#while (j<o)
            #local T1=<j+(o+1)*i,j+(o+1)*(i+1),(j+1)+(o+1)*i>+2*(m+1)*(n+1)+k*(o+1)*(m+1);
            ,T1,T1+x*(o+2)
          #local j=j+1;#end
        #local i=i+1;#end
      #local k=k+1;#end
    }
    #if (UV)  uv_mapping #end
  }
#end


#macro Orient(V,N)
  #local Rz=#if(V.y>0)-#end degrees(acos(vdot(x,vnormalize(V*<1,1,0>))));
  #local Vnz=vrotate(V,z*Rz);
  #local NVz=vrotate(N,z*Rz);
  #local Ry=#if(Vnz.z<0)-#end degrees(acos(vdot(x,vnormalize(Vnz*<1,0,1>))));
  #local NVy=vrotate(NVz,y*Ry);
  #local Rx=#if(NVy.z>0)-#end degrees(acos(vdot(y,vnormalize(NVy*<0,1,1>))));
  rotate <-Rx,-Ry,-Rz>
#end


//                                                              
// MobiusP for POVRay 3.5 v2.1                                  
// Trevor G Quayle, July 16, 2002                               
//                                                              
// This returns the vectors associated with a point on a mobius 
// strip.                                                       
//                                                              
// syntax: MobiusP(Ax,nR,SA,I,P,F,N)                            
// Ax = <x,z> half-length of x-axis and z axis                  
// nR = # half-twists of the strip                              
// SA = angle of twisting at +x axis                            
// I  = <u,v> point to return values for. u=angle in degrees    
//      from the +x axis*, v=deviation from the N.A. (+='right' 
//      - = 'left'                                              
// P  = point on mobius      \  calculated by the macro and     
// F  = 'forward' vector at P |-returned to the scene. They     
// N  = 'up' vector at P     /  must be declared identifiers    
//                                                              
// identifiers: P, F & N may be declared as anything as they    
//              are not used by the macro, they return <x,y,z>  
//              vectors to the scene. See example.              
//                                                              
// *note: the angle represents the parametric angle, NOT the    
//        actual angle from the +x axis (except in the case of  
//        a circle). Each degree increment will not produce a   
//        uniform distance increment, but will be close. (i.e., 
//        45°=approx 1/8 the distance, 90°=exactly 1/4, etc)    
//        As the position nears the major axes of the ellipse,  
//        the change in distance will decrease. This amount     
//        will depend on the eccentricity of the ellipse. (In   
//        the case of a circle, the distances will be equal).   
//        This is the best I could do with getting into complex 
//        formulae. Elliptical arc lengths get quite messy!     
//                                                              
// ex:                                                          
// #declare mP=0;#declare mF=0;#declare mN=0;                   
// MobiusP(<20,10>,1,45,<135,1>,mP,mF,mN)                       
//                                                              
// for a 40x20 mobius with one twist starting at 45°, mP, mF &  
// mN will return the position, 'forward' and 'up' vectors for  
// a point approx 3/8 around the strip (135°), 1 unit 'right'   
// of the N.A.                                                  
//                                                              
// 'forward', 'up' & 'right' are relative to the mobius surface 
//                                                              

#macro MobiusP(Ax,nR,SA,I,P,V,N)
  #local Lx=Ax.u;#local Lz=Ax.v;
  #macro PE(A,B,C) <B*cos(A),0,C*sin(A)> #end
  #macro NE(P,B,C) vnormalize(<C*P.x,0,pow(B,2)*sqrt(1-pow(P.x/B,2))#if(P.z<0)*-1#end>) #end
  #local A=radians(I.u);
  #local T=radians(I.u/360*nR*180+SA);
  #local P=PE(A,Lx,Lz);
  #local N=NE(P,Lx,Lz);
  #local R=<N.x*cos(T),sin(T),N.z*cos(T)>;
  #local P=P+I.v*R;
  #local A0=radians(I.u-1);
  #local T0=radians((I.u-1)/360*nR*180+SA);
  #local P0=PE(A0,Lx,Lz);
  #local N0=NE(P0,Lx,Lz);
  #local P0=P0+I.v*<N0.x*cos(T0),sin(T0),N0.z*cos(T0)>;
  #local A1=radians(I.u+1);
  #local T1=radians((I.u+1)/360*nR*180+SA);
  #local P1=PE(A1,Lx,Lz);
  #local N1=NE(P1,Lx,Lz);
  #local P1=P1+I.v*<N1.x*cos(T1),sin(T1),N1.z*cos(T1)>;
  #local V=vnormalize(P1-P0);
  #local N=vnormalize(vcross(V,R));
#end



//                                                              
// Mobius for POVRay 3.5 from objects v2.1                      
// Trevor G Quayle, July 16, 2002                               
//                                                              
// This is a mobius strip made from a user-defined object. The  
// strip twists around a neutral-axis defined by a closed       
// ellipse on the x-z plane.                                    
//                                                              
// syntax: MobiusO(Ax,W,s,nR,SA,Ob)                             
// Ax = <x,z> half-length of x-axis and z axis                  
// W  = Width of strip (centred on N.A.)                        
// s  = <m,n> # objects in width (m) & circumference (n)        
// nR = # half-twists of the strip (0=no twist, 1=180deg, etc)  
//      (negative #s will reverse the direction of the twist)   
// Ob = User-defined object. In order to have object oriented   
//      correctly, create it at the origin with x='forward',    
//      y='up'.                                                 
//                                                              
// ex:                                                          
// object{Mobius(<20,10>,4,<6,200>,1,45,yes) MyObj}             
//                                                              
// creates a smooth mobius 40x20 units in size, 4 units wide,   
// with one twist (traditional mobius). The twist starts at 45° 
// on the +x axis. There are 6 'Myobj's across and 200 around   
// the circumference of the mobius                              
//                                                              
// This macro uses both the MobiusP and Orient macros included  
// in "Mobius.inc" They must be available to use this macro.    
//                                                              

#macro MobiusO(Ax,W,s,nR,SA,Obj)
  #local s=<int(s.u),int(s.v)>;
  #if (Ax.u<=0) #debug concat("X-axis must be  > 0, assigning 2*W" #local Ax=<2*W,Ax.v>; #end
  #if (Ax.v<=0) #debug concat("Z-axis must be  > 0, assigning 2*W" #local Ax=<Ax.u,2*W>; #end
  #if (s.u<1) #debug "m must be an integer > 0, assigning 1\n" #local s=<1,s.v>; #end
  #if (s.v<4) #debug "n must be an integer > 3, assigning 4\n" #local s=<s.u,4>; #end
  #if (min(Ax.u,Ax.v)<W) #debug "A minimum axis smaller than the width may produce unpredictable results!\n" #end
  #if (mod(nR,1)!=0) #debug "nR is not an integer. Strip is not closed!\n" #end
  #local mP=0;#local mV=0;#local mN=0;
  union{
    #local i=0;#while (i<s.v)
      #local j=0;#while (j<s.u)
          #local D=(0#if(s.u>1)+(j-(s.u-1)/2)*W/(s.u-1)#end);
          MobiusP(Ax,nR,SA,<i*360/s.v,D>,mP,mV,mN)
        object{Obj transform{Orient(mV,mN)} translate mP}
      #local j=j+1;#end
    #local i=i+1;#end
  }
#end