// Persistence of Vision Ray Tracer Include File
// File: Catenary.inc
// Desc: Proper catenary chain. Version 1.0
// Date: 2005.07.30    
// Auth: PM 2Ring
//
// Catenary parameter calculations thanks to Zdislav V. Kovarik. See below 
//
// For instructions, please see "Catenary.txt"
//
//-------------------------------------------------------------------------

#ifndef(Catenary_Inc_Temp)
#declare Catenary_Inc_Temp=version;
#version 3.6;

#ifdef(View_POV_Include_Stack)
  #debug "including Catenary.inc\n"
#end
               
//-------------------------------------------------------------------------
//                                                                             
// From: kovarik@mcmail.cis.McMaster.CA (Zdislav V. Kovarik)                   
// Subject: Re: Catenary                                                       
// Date: 5 Nov 1999 14:33:22 -0500                                             
// Newsgroups: sci.math                                                        
// Keywords: fitting a catenary to match a suspended cable                     
//                                                                             
// Upright catenary:  y - y_0 = a * cosh((x - x_0)/a) ,  a > 0.                
// The vertex is (x_0, y_0+a).                                                 
// The parameter a turns out to be the radius of curvature at the vertex.      
// Remark: The radius of curvature at a point above x is                       
//   a * (cosh((x-x_0)/a)^2.                                                   
//                                                                             
// Problem: Parameters a, x_0, y_0 to be found so that the catenary arc passes 
// through (x_1, y_1), (x_2, y_2) and has length L between these points        
// (provided L > sqrt((x_2 - x_1)^2 + (y_2 - y_1)^2) )                        
//                                                                             
// The equations to be solved are (after some symmetrizing manipulation        
// of the arclength integral)                                                  
//                                                                             
//   y_1 - y_0 = a * cosh ((x_1 - x_0)/a)                                      
//   y_2 - y_0 = a * cosh ((x_2 - x_0)/a)                                      
//   2*a * cosh((x_1 + x_2 - 2*x_0)/(2*a)) * sinh((x_2 - x_1)/(2*a)) = L       
//                                                                             
// and it can be reduced to solving for one unknown at a time.                 
//                                                                             
// First, we introduce an auxiliary unknown z which is to satisfy              
//                                                                             
//   sinh(z) / z = sqrt(L^2 - (y_2 - y_1)^2) / abs(x_2 - x_1)  ,   z > 0       
//                                                                             
// (the only transcendental non-elementary equation)                           
// and then the unknowns pop out:                                              
//                                                                             
//     a = abs(x_2 - x_1) / (2*z)                                              
//   x_0 = (1/2)*(x_1 + x_2 - a * ln ((L + (y_2 - y_1)) / (L - (y_2 - y_1)))   
//   y_0 =  (y_1 + y_2)/2 - (L/2) * coth(z).                                   
//                                                                             
//-------------------------------------------------------------------------

// Find B such that A=sinh(B)/B. Named in parallel to sinc()
#macro asinch(A)
  #local B = sqrt(6*max(1e-4,A-1));   //1st approx, from sinh(B) = B + B^3/3! + ...             
  #local B = asinh(A*B);              //2nd approx, from sinh(B) = A*B                  
                                                                                      
  #local MI = 8;                                                                       
  #local I = 0;                                                                       
  #while(I<MI)                         //Newton's method                                 
    #local S = sinh(B);                                                            
    #local C = B * cosh(B);                                                        
    #local M = (A*B + C - 2*S) / (C - S);                                                                     
    #local B = B * M;                                                              
    #local I = I + 1;                                                                
    #if(abs(M-1)<1e-12)               //bailout                                      
      #local I = MI;                                                                  
    #end                                                                           
  #end                                                                             
  B                                   //return value                                 
#end 

//Find width of an object's bounding box
#macro BBWidth(A) (max_extent(A) - min_extent(A)).x #end

// Make a catenary. Parameters: Link object, Start point, End Point, 
// Slackness of the chain >1, Link overlap, Extra twist on whole cable,
// Alternate 90 degree twist on X-axis, Starting phase of Alternation, 
// Links to Skip at start, Links to Skip at end, Make number of links odd
#macro Chain(Link, StartA, EndA, Slack, Overlap, Alternate, AltPhase, Twist, StartSkip, EndSkip, OddLinks)
  //Temporarily translate Start point to origin & work in XY plane 
  #local End = EndA - StartA;
  #local TH = degrees(atan2(End.z, End.x));
  #local End = vrotate(End, y*TH);              //Rotate to positive x-axis                                 
  
  //Find required chain length and number of links
  #local Len = vlength(End) * Slack;            //Basic chain length
  #local LL = BBWidth(Link) / Overlap;          //Link Length adjusted for link overlap  
  
  #local Steps = ceil(Len/LL);                  //Number of links minus one, to account for half-links at each end 
  #if(OddLinks)                                 //Round up to next even integer
    #local Steps = 2*ceil(.5*Steps); 
  #end   
  #local Len = Steps * LL;                      //Adjusted chain length
  
  //Find vertex of catenary that connects <0,0,0> & End, with length Len.
  #local P = sqrt(Len*Len - End.y*End.y) / End.x;
  #local Q = asinch(P);
  #local A = End.x / (2*Q);                     //Catenary curvature parameter. 
  
  #local X = A * ln((Len + End.y) / (Len - End.y));
  #local Y = Len/tanh(Q);  
  #local V = (End - <X, Y, 0>)*.5;              //Vertex - a*y
  #local S1 = A*sinh(V.x/A);                    //Arclength at vertex

  //Step evenly along catenary parametrized by arclength. 
  union {
    #local I = StartSkip;
    #while (I <= Steps - EndSkip)
      #local S = LL * I - S1;                   //Arclength from Vertex
      #local M = S / A;                         //Slope of tangent
      #local X = asinh(M);    
      #local Y = sqrt(1 + M*M);
      
      object{
        Link       
        rotate (Alternate*90*mod(I+AltPhase,2) + Twist*360*I/(Steps-1))*x
        rotate z*degrees(atan(M))               //Rotate parallel to tangent
        translate V + A*<X, Y, 0>
      }
      #local I = I + 1;
    #end 
  
    //Transform back  
    rotate -y*TH
    translate StartA
  }
#end 

#version Catenary_Inc_Temp;
#end

//-------------------------------------------------------------------------