// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7 ======= 8 ======= 9 ======= 10

// TubularCurves.pov
// For drawing a simple tube like triangle mesh around parametric functions in the xy-plane
// By Bill Walker and Tor Olav Kristensen, December 2022

#version 3.7;

global_settings { assumed_gamma 1.0 }

#include "colors.inc"
#include "shapes3.inc" // for Segment_of_Torus()

#declare TAU = 2*pi;
#declare EPS = 1e-6;

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7 ======= 8 ======= 9 ======= 10

#declare pCameraPosition = -5*z; // front view
#declare pCameraLookAt  = 0*z;

// zoom for the orthographic view: 4 zooms in 4x, 8 zooms in 8x, etc.
#declare Zoom = 35;

camera {
    orthographic
    right x*image_width/Zoom
    up y*image_height/Zoom
    location pCameraPosition
    look_at pCameraLookAt
    // rotate y*60
}

light_source {
    100*< -5, 5, -25>
    color White*2
    shadowless
}

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7 ======= 8 ======= 9 ======= 10

#macro DrawCentralFunction(FnX, FnY, MinU, MaxU, StepU, Radius)

    #local SpanU = MaxU - MinU;
    #local NoOfSteps = int(SpanU/StepU);
    #local DU = SpanU/NoOfSteps;

    union {
        #local U = MinU;
        #local pC_Prev = <FnX(U), FnY(U), 0>;
        sphere { pC_Prev, Radius }
        #for (I, 1, NoOfSteps)
            #local U = U + DU;
            #local pC = <FnX(U), FnY(U), 0>;
            #if (vlength(pC - pC_Prev) > EPS)
                cylinder { pC_Prev, pC, Radius }
            #end // if
            #local pC_Prev = pC;
        #end // for
    }

#end // macro DrawCentralFunction


#macro CalculatePhi(D_FnX, D_FnY, U)

    #local Phi = atan2(D_FnY(U), D_FnX(U));
    // #debug concat("U: ", str(U, 0, 12), "    Phi: ", str(Phi, 0, 12), "\n")

    Phi

#end // macro CalculatePhi


#macro DrawTube(FnX, FnY, D_FnX, D_FnY, MinU, MaxU, StepU, TubeRadius, NoOfSides)

    // create mesh vertices
    #local NoOfRings = int((MaxU - MinU)/StepU);
    #local Mesh = array[NoOfRings+1][NoOfSides+1];

    #local SpanU = MaxU - MinU;
    #local NoOfSteps = int(SpanU/StepU);
    #for (I, 0, NoOfSteps)
        #local U = MinU + I/NoOfSteps*SpanU;
        #local pC = <FnX(U), FnY(U), 0>;
        #local Phi = CalculatePhi(D_FnX, D_FnY, U);
        // #debug concat("Phi: ", str(degrees(Phi), 0, 5), "\n")
        #for (J, 0, NoOfSides)
            // #debug concat("I (Ring): ", str(I, 0, 0), ", J: ", str(J, 0, 0), "\n")
            #local Theta = J/NoOfSides*TAU;
            // These are the coordinates of the default circle
            #local QY = TubeRadius*cos(Theta);
            #local QZ = TubeRadius*sin(Theta);
            // the circle is in the zy-plane
            // apply a rotation around the z-axis
            #local MX = -QY*sin(Phi);
            #local MY = +QY*cos(Phi);
            #local MZ = QZ;
            #local pM0 = pC + <MX, MY, MZ>;
            #local Mesh[I][J] = pM0;
        #end // for
    #end // for

    union {
        // draw tube mesh
        #for (I, 1, NoOfRings)
            #for (J, 1, NoOfSides)
                // #debug concat("I: ", str(I, 0, 0), ", J: ", str(J, 0, 0), "\n")
                triangle { Mesh[I  ][J-1], Mesh[I-1][J  ], Mesh[I-1][J-1] }
                triangle { Mesh[I  ][J  ], Mesh[I-1][J  ], Mesh[I  ][J-1] }
            #end // for
        #end // for
     }

#end // macro DrawTube


#macro DrawTubeCircles(FnX, FnY, D_FnX, D_FnY, U, NoOfSides, TubeRadius, Radius)

    #local pC = <FnX(U), FnY(U), 0>;
    #local vB = <0, 0, 1>;
    #local vT = vnormalize(<D_FnX(U), D_FnY(U), 0>);
    #local Phi = CalculatePhi(D_FnX, D_FnY, U);
    // #debug concat("Phi: ", str(Phi, 0, 5), "\n")
    #local vN = <sin(Phi), -cos(Phi), 0>;

    union {
        // draw line parallel to y-axis, from x-axis to point on curve
        union {
            sphere { <FnX(U), 0, 0>, Radius }
            #if (vlength(pC - <FnX(U), 0, 0>) > 0)
                cylinder { <FnX(U), 0, 0>, pC, Radius }
            #end // if
            #if (pC.y <= 0)
                cylinder { pC, pC - 6*TubeRadius*y, Radius }
            #end // if
            #if (pC.y >= 0)
                cylinder { pC, pC + 6*TubeRadius*y, Radius }
            #end // if
            pigment { color Green }
        }
        // draw line parallel to x-axis, from y-axis to point on curve
        union {
            sphere { <0, FnY(U), 0>, Radius }
            #if (vlength(pC - <0, FnY(U), 0>) > 0)
                cylinder { <0, FnY(U), 0>, pC, Radius }
            #end // if
            #if (pC.x <= 0)
                cylinder { pC, pC - 6*TubeRadius*x, Radius }
            #end // if
            #if (pC.x >= 0)
                cylinder { pC, pC + 6*TubeRadius*x, Radius }
            #end // if
            pigment { color Red }
        }
        union {
            // point on curve
            sphere {
                <0, 0, 0>, Radius
                pigment { color Red }
            }
            // - - - - - - - - DRAW CIRCLES  - - - - - - 
            // draw default circle
             union {
                #for (I, 0, NoOfSides)
                    #local Theta = I/NoOfSides*TAU;
                    #local QY = TubeRadius*cos(Theta);
                    #local QZ = TubeRadius*sin(Theta);
                    #local pQ = <0, QY, QZ>;
                    sphere {
                        pQ, Radius
                        pigment { color Orange }
                    }
                #end // for
            }
            // draw rotated circle
            union {
                #for (I, 0, NoOfSides)
                    #local Theta = I/NoOfSides*TAU;
                    // the coordinates of the default circle
                    #local QY = TubeRadius*cos(Theta);
                    #local QZ = TubeRadius*sin(Theta);
                    // the circle is in the zy plane
                    // apply rotation around the z-axis
                    #local MX = -QY*sin(Phi);
                    #local MY = +QY*cos(Phi);
                    #local MZ = QZ;
                    #local pM = <MX, MY, MZ>;
                    sphere {
                        pM, Radius
                        pigment { color White }
                    }
                #end // for
            }
            // - - - - - - - - DRAW LINES  - - - - - - - 
            // draw tangent vector of curve
            cylinder {
                <0, 0, 0>, 4*TubeRadius*vT, Radius
                pigment { color Cyan }
            }
            // draw normal vector of curve
            cylinder {
                <0, 0, 0>, 4*TubeRadius*vN, Radius
                pigment { color Yellow }
            }
            // - - - - - - - - DRAW ARCS - - - - - - - - 
            // make a torus segment for angle Phi
            object {
                Segment_of_Torus(
                    2*TubeRadius, // radius major
                    Radius, // radius minor
                    degrees(Phi) // segment angle
                )
                pigment { color Cyan }
                rotate 90*x
            }
            // show circle rotation angle by Phi
            object {
                Segment_of_Torus(
                    2*TubeRadius, // radius major
                    Radius, // radius minor
                    degrees(Phi) // segment angle
                )
                pigment { color Yellow }
                rotate 90*x
                rotate -90*z
             }
            // draw a perpedicular symbol between the normal and the tangent vectors
            union {
                sphere {
                    TubeRadius*vN + TubeRadius*vT, Radius
                }
                 cylinder {
                    TubeRadius*vN, TubeRadius*vN + TubeRadius*vT, Radius
                }
                cylinder {
                    TubeRadius*vT, TubeRadius*vN + TubeRadius*vT, Radius
                }
                pigment { color Gray40 }
            }
            translate pC
        }
    }

#end // macro DrawTubeCircles


#macro GraphPhi(FnX, FnY, D_FnX, D_FnY, MinU, MaxU, StepU, Radius)

    #local Font = "timrom.ttf"
    #local SpanU = MaxU - MinU;
    #local NoOfSteps = int(SpanU/StepU);
    #local DU = SpanU/NoOfSteps;

    // graph function
    union {
        #local U = MinU;
        #local Phi = CalculatePhi(D_FnX, D_FnY, U);
        // scale Phi for graphing
        #local ScaledPhi = Phi/(TAU/2);
        #local p0_Prev = <U, ScaledPhi, 0>;
        sphere { p0_Prev, Radius }
        #for (I, 1, NoOfSteps)
            #local U = U + DU;
            #local Phi = CalculatePhi(D_FnX, D_FnY, U);
            // scale Phi for graphing
            #local ScaledPhi = Phi/(TAU/2);
            #local p0 = <U, ScaledPhi, 0>;
            sphere { p0, Radius }
            #if (vlength(p0 - p0_Prev) > EPS)
                cylinder { p0_Prev, p0, Radius }
            #end // if
            #local p0_Prev = p0;
        #end // for
        cylinder { <MinU,  1, 0>, <MaxU,  1, 0>, Radius }
        cylinder { <MinU,  0, 0>, <MaxU,  0, 0>, Radius }
        cylinder { <MinU, -1, 0>, <MaxU, -1, 0>, Radius }
        text {
            ttf Font, " Plot of curve tangent angle Phi vs u", 0.02, 0.0
            scale <1, 1, 1>/4
            translate <MinU, +1.125, 0>
        }
        text {
            ttf Font, " +pi", 0.02, 0.0
            scale <1, 1, 1>/4
            translate <MinU - 0.35, +1, 0>
        }
        text {
            ttf Font, "  0 ", 0.02, 0.0
            scale <1, 1, 1>/4
            translate <MinU - 0.35,  0, 0>
        }
        text {
            ttf Font, " -pi", 0.02, 0.0
            scale <1, 1, 1>/4
            translate <MinU - 0.35, -1, 0>
        }
        texture {
            pigment { color rgb <2, 1, 0>/3 }
            finish { emission 1 }
        }
     }

#end // macro GraphPhi

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7 ======= 8 ======= 9 ======= 10

#declare NoOfFunctions = 7;

#declare FunctionsX = array[NoOfFunctions];
#declare FunctionsY = array[NoOfFunctions];

#declare DerFunctionsX = array[NoOfFunctions];
#declare DerFunctionsY = array[NoOfFunctions];


#declare I = 0;
#declare FunctionsX[I] = function(u) { u };
#declare FunctionsY[I] = function(u) { pow(u, 2) + 1 };
#declare DerFunctionsX[I] = function(u) { 1 };
#declare DerFunctionsY[I] = function(u) { 2*u };

#declare I = 1;
#declare FunctionsX[I] = function(u) { u };
#declare FunctionsY[I] = function(u) { 0.4*cos(u)*(-pow(u, 2) + 3) - 1 };
#declare DerFunctionsX[I] = function(u) { 1 };
#declare DerFunctionsY[I] = function(u) { 0.4*(pow(u, 2) - 3)*sin(u) - 0.8*u*cos(u) };

#declare I = 2;
#declare FunctionsX[I] = function(u) { u };
#declare FunctionsY[I] = function(u) { pow(u, 3) };
#declare DerFunctionsX[I] = function(u) { 1 };
#declare DerFunctionsY[I] = function(u) { 3*pow(u, 2) };

#declare I = 3;
#declare FunctionsX[I] = function(u) { u };
#declare FunctionsY[I] = function(u) { 4*tanh(u) + 4 };
#declare DerFunctionsX[I] = function(u) { 1 };
#declare DerFunctionsY[I] = function(u) { 16*pow(cosh(u), 2)/pow(cosh(2*u) + 1, 2) };

#declare I = 4;
#declare FunctionsX[I] = function(u) { u };
#declare FunctionsY[I] = function(u) { 8*u/(abs(u) + 1) };
#declare DerFunctionsX[I] = function(u) { 1 };
#declare DerFunctionsY[I] = function(u) { 8/pow(abs(u) + 1, 2) };

/*
#declare I = I + 1;
#declare FunctionsX[I] = function(u) { u };
#declare FunctionsY[I] = function(u) { u*sin(TAU/u) + 1 }; // sind
#declare DerFunctionsX[I] = function(u) { 1 };
#declare DerFunctionsY[I] = function(u) { sin(TAU/u) - (TAU*cos(TAU/u))/u };

#declare I = I + 1;
#declare FunctionsX[I] = function(u) { u };
#declare FunctionsY[I] = function(u) { floor(5* sin(u))/5 + 1 };
#declare DerFunctionsX[I] = function(u) { 1 };
#declare DerFunctionsY[I] = function(u) { sin(5*pi*sin(u)) };
*/

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7 ======= 8 ======= 9 ======= 10

// choose function
#declare I = 1;

#declare NoOfSides = 60;
#declare TubeRadius = 0.1;
#declare LineRadius = 0.02;
#declare MinU = -10.0;
#declare MaxU = +10.0;
#declare StepU = 1/100;

#if (true)
    object {
        DrawCentralFunction(
            FunctionsX[I],
            FunctionsY[I],
            MinU, MaxU, StepU,
            LineRadius*2
        )
        pigment { color White }
    }
#end // if

#if (true)
    object {
        DrawTube(
            FunctionsX[I],
            FunctionsY[I],
            DerFunctionsX[I],
            DerFunctionsY[I],
            MinU, MaxU, StepU,
            TubeRadius,
            NoOfSides
        )
        texture {
            pigment { color rgbf <0.2, 0.3, 1.2, 0.5> }
            // pigment { color White }
            // finish { specular 0.4 }
        }
    }
#end // if

#if (true)
    union {
        #for (U,  MinU, MaxU, 0.50)
            DrawTubeCircles(
                FunctionsX[I],
                FunctionsY[I],
                DerFunctionsX[I],
                DerFunctionsY[I],
                U,
                NoOfSides,
                TubeRadius,
                LineRadius
            )
        #end // for
    }
#end // if

#if (true)
    object {
        GraphPhi(
            FunctionsX[I],
            FunctionsY[I],
            DerFunctionsX[I],
            DerFunctionsY[I],
            MinU, MaxU, StepU,
            LineRadius
        )
        translate -2*z
        translate -3*y
    }
#end // if

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7 ======= 8 ======= 9 ======= 10

background { color Gray20 }

cylinder {
    -10*x, +10*x, LineRadius
    pigment { color Red }
}
cylinder {
    -10*y, +10*y, LineRadius
    pigment { color Green }
}
cylinder {
    -10*z, +10*z, LineRadius
    pigment { color Blue }
}

/*
plane {
    y, 0
    pigment {
        checker
        color Black, color White
    }
}
*/

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7 ======= 8 ======= 9 ======= 10