home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: Graphics / Graphics.zip / povsrc31.zip / torus.c < prev    next >
C/C++ Source or Header  |  1999-10-20  |  17KB  |  933 lines
  1. /****************************************************************************
  2. *                   torus.c
  3. *
  4. *  This module implements functions that manipulate torii.
  5. *
  6. *  This module was written by Dieter Bayer [DB].
  7. *
  8. *  from Persistence of Vision(tm) Ray Tracer
  9. *  Copyright 1996,1999 Persistence of Vision Team
  10. *---------------------------------------------------------------------------
  11. *  NOTICE: This source code file is provided so that users may experiment
  12. *  with enhancements to POV-Ray and to port the software to platforms other
  13. *  than those supported by the POV-Ray Team.  There are strict rules under
  14. *  which you are permitted to use this file.  The rules are in the file
  15. *  named POVLEGAL.DOC which should be distributed with this file.
  16. *  If POVLEGAL.DOC is not available or for more info please contact the POV-Ray
  17. *  Team Coordinator by email to team-coord@povray.org or visit us on the web at
  18. *  http://www.povray.org. The latest version of POV-Ray may be found at this site.
  19. *
  20. * This program is based on the popular DKB raytracer version 2.12.
  21. * DKBTrace was originally written by David K. Buck.
  22. * DKBTrace Ver 2.0-2.12 were written by David K. Buck & Aaron A. Collins.
  23. *
  24. *****************************************************************************/
  25.  
  26. /****************************************************************************
  27. *
  28. *  Explanation:
  29. *
  30. *  ---
  31. *
  32. *  June 1994 : Creation. [DB]
  33. *
  34. *****************************************************************************/
  35.  
  36. #include "frame.h"
  37. #include "povray.h"
  38. #include "vector.h"
  39. #include "povproto.h"
  40. #include "bbox.h"
  41. #include "polysolv.h"
  42. #include "matrices.h"
  43. #include "objects.h"
  44. #include "torus.h"
  45.  
  46.  
  47.  
  48. /*****************************************************************************
  49. * Local preprocessor defines
  50. ******************************************************************************/
  51.  
  52. /* Minimal depth for a valid intersection. */
  53.  
  54. #define DEPTH_TOLERANCE 1.0e-4
  55.  
  56. /* Tolerance used for order reduction during root finding. */
  57.  
  58. #define ROOT_TOLERANCE 1.0e-4
  59.  
  60.  
  61.  
  62. /*****************************************************************************
  63. * Static functions
  64. ******************************************************************************/
  65.  
  66. static int intersect_torus (RAY *Ray, TORUS *Torus, DBL *Depth);
  67. static int   All_Torus_Intersections (OBJECT *Object, RAY *Ray, ISTACK *Depth_Stack);
  68. static int   Inside_Torus (VECTOR point, OBJECT *Object);
  69. static void  Torus_Normal (VECTOR Result, OBJECT *Object, INTERSECTION *Inter);
  70. static TORUS *Copy_Torus (OBJECT *Object);
  71. static void  Translate_Torus (OBJECT *Object, VECTOR Vector, TRANSFORM *Trans);
  72. static void  Rotate_Torus (OBJECT *Object, VECTOR Vector, TRANSFORM *Trans);
  73. static void  Scale_Torus (OBJECT *Object, VECTOR Vector, TRANSFORM *Trans);
  74. static void  Transform_Torus (OBJECT *Object, TRANSFORM *Trans);
  75. static void  Invert_Torus (OBJECT *Object);
  76. static void  Destroy_Torus (OBJECT *Object);
  77.  
  78.  
  79.  
  80. /*****************************************************************************
  81. * Local variables
  82. ******************************************************************************/
  83.  
  84. static METHODS Torus_Methods =
  85. {
  86.   All_Torus_Intersections, Inside_Torus, Torus_Normal, Default_UVCoord,
  87.   (COPY_METHOD)Copy_Torus, Translate_Torus, Rotate_Torus,
  88.   Scale_Torus, Transform_Torus, Invert_Torus, Destroy_Torus
  89. };
  90.  
  91.  
  92. /*****************************************************************************
  93. *
  94. * FUNCTION
  95. *
  96. *   All_Torus_Intersections
  97. *
  98. * INPUT
  99. *
  100. *   Object      - Object
  101. *   Ray         - Ray
  102. *   Depth_Stack - Intersection stack
  103. *   
  104. * OUTPUT
  105. *
  106. *   Depth_Stack
  107. *   
  108. * RETURNS
  109. *
  110. *   int - TRUE, if an intersection was found
  111. *   
  112. * AUTHOR
  113. *
  114. *   Dieter Bayer
  115. *   
  116. * DESCRIPTION
  117. *
  118. *   Determine ray/torus intersection and clip intersection found.
  119. *
  120. * CHANGES
  121. *
  122. *   Jun 1994 : Creation.
  123. *
  124. ******************************************************************************/
  125.  
  126. static int All_Torus_Intersections(OBJECT *Object, RAY *Ray, ISTACK *Depth_Stack)
  127. {
  128.   int i, max_i, Found;
  129.   DBL Depth[4];
  130.   VECTOR IPoint;
  131.  
  132.   Found = FALSE;
  133.  
  134.   if ((max_i = intersect_torus(Ray, (TORUS *)Object, Depth)) > 0)
  135.   {
  136.     for (i = 0; i < max_i; i++)
  137.     {
  138.       if ((Depth[i] > DEPTH_TOLERANCE) && (Depth[i] < Max_Distance))
  139.       {
  140.         VEvaluateRay(IPoint, Ray->Initial, Depth[i], Ray->Direction);
  141.  
  142.         if (Point_In_Clip(IPoint, Object->Clip))
  143.         {
  144.           push_entry(Depth[i], IPoint, Object, Depth_Stack);
  145.  
  146.           Found = TRUE;
  147.         }
  148.       }
  149.     }
  150.   }
  151.  
  152.   return(Found);
  153. }
  154.  
  155.  
  156.  
  157. /*****************************************************************************
  158. *
  159. * FUNCTION
  160. *
  161. *   intersect_torus
  162. *
  163. * INPUT
  164. *
  165. *   Ray   - Ray
  166. *   Torus - Torus
  167. *   Depth - Intersections found
  168. *   
  169. * OUTPUT
  170. *
  171. *   Depth
  172. *   
  173. * RETURNS
  174. *
  175. *   int - Number of intersections found
  176. *   
  177. * AUTHOR
  178. *
  179. *   Dieter Bayer
  180. *   
  181. * DESCRIPTION
  182. *
  183. *   Determine ray/torus intersection.
  184. *
  185. *   Note that the torus is rotated about the y-axis!
  186. *
  187. * CHANGES
  188. *
  189. *   Jun 1994 : Creation.
  190. *
  191. ******************************************************************************/
  192.  
  193. static int intersect_torus(RAY *Ray, TORUS *Torus, DBL *Depth)
  194. {
  195.   int i, n;
  196.   DBL len, R2, Py2, Dy2, PDy2, k1, k2;
  197.   DBL y1, y2, r1, r2;
  198.   DBL c[5], r[4];
  199.   VECTOR P, D;
  200.  
  201.   Increase_Counter(stats[Ray_Torus_Tests]);
  202.  
  203.   /* Transform the ray into the torus space. */
  204.  
  205.   MInvTransPoint(P, Ray->Initial, Torus->Trans);
  206.  
  207.   MInvTransDirection(D, Ray->Direction, Torus->Trans);
  208.  
  209.   VLength(len, D);
  210.  
  211.   VInverseScaleEq(D, len);
  212.  
  213.   i = 0;
  214.  
  215.   y1 = -Torus->r;
  216.   y2 =  Torus->r;
  217.   r1 = Sqr(Torus->R - Torus->r);
  218.   r2 = Sqr(Torus->R + Torus->r);
  219.  
  220. #ifdef TORUS_EXTRA_STATS
  221.   Increase_Counter(stats[Torus_Bound_Tests]);
  222. #endif
  223.  
  224.   if (Test_Thick_Cylinder(P, D, y1, y2, r1, r2))
  225.   {
  226. #ifdef TORUS_EXTRA_STATS
  227.     Increase_Counter(stats[Torus_Bound_Tests_Succeeded]);
  228. #endif
  229.  
  230.     R2   = Sqr(Torus->R);
  231.     r2   = Sqr(Torus->r);
  232.  
  233.     Py2  = P[Y] * P[Y];
  234.     Dy2  = D[Y] * D[Y];
  235.     PDy2 = P[Y] * D[Y];
  236.  
  237.     k1   = P[X] * P[X] + P[Z] * P[Z] + Py2 - R2 - r2;
  238.     k2   = P[X] * D[X] + P[Z] * D[Z] + PDy2;
  239.  
  240.     c[0] = 1.0;
  241.  
  242.     c[1] = 4.0 * k2;
  243.  
  244.     c[2] = 2.0 * (k1 + 2.0 * (k2 * k2 + R2 * Dy2));
  245.  
  246.     c[3] = 4.0 * (k2 * k1 + 2.0 * R2 * PDy2);
  247.  
  248.     c[4] = k1 * k1 + 4.0 * R2 * (Py2 - r2);
  249.  
  250.     n = Solve_Polynomial(4, c, r, Test_Flag(Torus, STURM_FLAG), ROOT_TOLERANCE);
  251.  
  252.     while(n--)
  253.     {
  254.       Depth[i++] = r[n] / len;
  255.     }
  256.   }
  257.  
  258.   if (i)
  259.   {
  260.     Increase_Counter(stats[Ray_Torus_Tests_Succeeded]);
  261.   }
  262.  
  263.   return(i);
  264. }
  265.  
  266.  
  267.  
  268. /*****************************************************************************
  269. *
  270. * FUNCTION
  271. *
  272. *   Inside_Torus
  273. *
  274. * INPUT
  275. *
  276. *   IPoint - Intersection point
  277. *   Object - Object
  278. *   
  279. * OUTPUT
  280. *   
  281. * RETURNS
  282. *
  283. *   int - TRUE if inside
  284. *   
  285. * AUTHOR
  286. *
  287. *   Dieter Bayer
  288. *   
  289. * DESCRIPTION
  290. *
  291. *   Test if a point lies inside the torus.
  292. *
  293. * CHANGES
  294. *
  295. *   Jun 1994 : Creation.
  296. *
  297. ******************************************************************************/
  298.  
  299. static int Inside_Torus(VECTOR IPoint, OBJECT *Object)
  300. {
  301.   DBL r, r2;
  302.   VECTOR P;
  303.   TORUS *Torus = (TORUS *)Object;
  304.  
  305.   /* Transform the point into the torus space. */
  306.  
  307.   MInvTransPoint(P, IPoint, Torus->Trans);
  308.  
  309.   r  = sqrt(Sqr(P[X]) + Sqr(P[Z]));
  310.  
  311.   r2 = Sqr(P[Y]) + Sqr(r - Torus->R);
  312.  
  313.   if (r2 <= Sqr(Torus->r))
  314.   {
  315.     return(!Test_Flag(Torus, INVERTED_FLAG));
  316.   }
  317.   else
  318.   {
  319.     return(Test_Flag(Torus, INVERTED_FLAG));
  320.   }
  321. }
  322.  
  323.  
  324.  
  325. /*****************************************************************************
  326. *
  327. * FUNCTION
  328. *
  329. *   Torus_Normal
  330. *
  331. * INPUT
  332. *
  333. *   Result - Normal vector
  334. *   Object - Object
  335. *   Inter  - Intersection found
  336. *   
  337. * OUTPUT
  338. *
  339. *   Result
  340. *   
  341. * RETURNS
  342. *   
  343. * AUTHOR
  344. *
  345. *   Dieter Bayer
  346. *   
  347. * DESCRIPTION
  348. *
  349. *   Calculate the normal of the torus in a given point.
  350. *
  351. * CHANGES
  352. *
  353. *   Jun 1994 : Creation.
  354. *
  355. ******************************************************************************/
  356.  
  357. static void Torus_Normal(VECTOR Result, OBJECT *Object, INTERSECTION *Inter)
  358. {
  359.   DBL dist;
  360.   VECTOR P, N, M;
  361.   TORUS *Torus = (TORUS *)Object;
  362.  
  363.   /* Transform the point into the torus space. */
  364.  
  365.   MInvTransPoint(P, Inter->IPoint, Torus->Trans);
  366.  
  367.   /* Get normal from derivatives. */
  368.  
  369.   dist = sqrt(P[X] * P[X] + P[Z] * P[Z]);
  370.  
  371.   if (dist > EPSILON)
  372.   {
  373.     M[X] = Torus->R * P[X] / dist;
  374.     M[Y] = 0.0;
  375.     M[Z] = Torus->R * P[Z] / dist;
  376.   }
  377.   else
  378.   {
  379.     Make_Vector(M, 0.0, 0.0, 0.0);
  380.   }
  381.  
  382.   VSub(N, P, M);
  383.  
  384.   /* Transform the normalt out of the torus space. */
  385.  
  386.   MTransNormal(Result, N, Torus->Trans);
  387.  
  388.   VNormalize(Result, Result);
  389. }
  390.  
  391.  
  392.  
  393. /*****************************************************************************
  394. *
  395. * FUNCTION
  396. *
  397. *   Translate_Torus
  398. *
  399. * INPUT
  400. *
  401. *   Object - Object
  402. *   Vector - Translation vector
  403. *   
  404. * OUTPUT
  405. *
  406. *   Object
  407. *   
  408. * RETURNS
  409. *   
  410. * AUTHOR
  411. *
  412. *   Dieter Bayer
  413. *
  414. * DESCRIPTION
  415. *
  416. *   Translate a torus.
  417. *
  418. * CHANGES
  419. *
  420. *   Jun 1994 : Creation.
  421. *
  422. ******************************************************************************/
  423.  
  424. static void Translate_Torus(OBJECT *Object, VECTOR Vector, TRANSFORM *Trans)
  425. {
  426.   Transform_Torus(Object, Trans);
  427. }
  428.  
  429.  
  430.  
  431. /*****************************************************************************
  432. *
  433. * FUNCTION
  434. *
  435. *   Rotate_Torus
  436. *
  437. * INPUT
  438. *
  439. *   Object - Object
  440. *   Vector - Rotation vector
  441. *   
  442. * OUTPUT
  443. *
  444. *   Object
  445. *   
  446. * RETURNS
  447. *   
  448. * AUTHOR
  449. *
  450. *   Dieter Bayer
  451. *   
  452. * DESCRIPTION
  453. *
  454. *   Rotate a torus.
  455. *
  456. * CHANGES
  457. *
  458. *   Jun 1994 : Creation.
  459. *
  460. ******************************************************************************/
  461.  
  462. static void Rotate_Torus(OBJECT *Object, VECTOR Vector, TRANSFORM *Trans)
  463. {
  464.   Transform_Torus(Object, Trans);
  465. }
  466.  
  467.  
  468.  
  469. /*****************************************************************************
  470. *
  471. * FUNCTION
  472. *
  473. *   Scale_Torus
  474. *
  475. * INPUT
  476. *
  477. *   Object - Object
  478. *   Vector - Scaling vector
  479. *   
  480. * OUTPUT
  481. *
  482. *   Object
  483. *   
  484. * RETURNS
  485. *   
  486. * AUTHOR
  487. *
  488. *   Dieter Bayer
  489. *   
  490. * DESCRIPTION
  491. *
  492. *   Scale a torus.
  493. *
  494. * CHANGES
  495. *
  496. *   Jun 1994 : Creation.
  497. *
  498. ******************************************************************************/
  499.  
  500. static void Scale_Torus(OBJECT *Object, VECTOR Vector, TRANSFORM *Trans)
  501. {
  502.   Transform_Torus(Object, Trans);
  503. }
  504.  
  505.  
  506.  
  507. /*****************************************************************************
  508. *
  509. * FUNCTION
  510. *
  511. *   Transform_Torus
  512. *
  513. * INPUT
  514. *
  515. *   Object - Object
  516. *   Trans  - Transformation to apply
  517. *   
  518. * OUTPUT
  519. *
  520. *   Object
  521. *   
  522. * RETURNS
  523. *   
  524. * AUTHOR
  525. *
  526. *   Dieter Bayer
  527. *   
  528. * DESCRIPTION
  529. *
  530. *   Transform a torus and recalculate its bounding box.
  531. *
  532. * CHANGES
  533. *
  534. *   Jun 1994 : Creation.
  535. *
  536. ******************************************************************************/
  537.  
  538. static void Transform_Torus(OBJECT *Object, TRANSFORM *Trans)
  539. {
  540.   Compose_Transforms(((TORUS *)Object)->Trans, Trans);
  541.  
  542.   Compute_Torus_BBox((TORUS *)Object);
  543. }
  544.  
  545.  
  546.  
  547. /*****************************************************************************
  548. *
  549. * FUNCTION
  550. *
  551. *   Invert_Torus
  552. *
  553. * INPUT
  554. *
  555. *   Object - Object
  556. *   
  557. * OUTPUT
  558. *
  559. *   Object
  560. *   
  561. * RETURNS
  562. *   
  563. * AUTHOR
  564. *
  565. *   Dieter Bayer
  566. *   
  567. * DESCRIPTION
  568. *
  569. *   Invert a torus.
  570. *
  571. * CHANGES
  572. *
  573. *   Jun 1994 : Creation.
  574. *
  575. ******************************************************************************/
  576.  
  577. static void Invert_Torus(OBJECT *Object)
  578. {
  579.   Invert_Flag(Object, INVERTED_FLAG);
  580. }
  581.  
  582.  
  583.  
  584. /*****************************************************************************
  585. *
  586. * FUNCTION
  587. *
  588. *   Create_Torus
  589. *
  590. * INPUT
  591. *   
  592. * OUTPUT
  593. *   
  594. * RETURNS
  595. *
  596. *   TORUS * - new torus
  597. *   
  598. * AUTHOR
  599. *
  600. *   Dieter Bayer
  601. *   
  602. * DESCRIPTION
  603. *
  604. *   Create a new torus.
  605. *
  606. * CHANGES
  607. *
  608. *   Jun 1994 : Creation.
  609. *
  610. ******************************************************************************/
  611.  
  612. TORUS *Create_Torus()
  613. {
  614.   TORUS *New;
  615.  
  616.   New = (TORUS *)POV_MALLOC(sizeof(TORUS), "torus");
  617.  
  618.   INIT_OBJECT_FIELDS(New,TORUS_OBJECT,&Torus_Methods)
  619.  
  620.   New->Trans = Create_Transform();
  621.  
  622.   New->R =
  623.   New->r = 0.0;
  624.  
  625.   return (New);
  626. }
  627.  
  628.  
  629.  
  630. /*****************************************************************************
  631. *
  632. * FUNCTION
  633. *
  634. *   Copy_Torus
  635. *
  636. * INPUT
  637. *
  638. *   Object - Object
  639. *   
  640. * OUTPUT
  641. *   
  642. * RETURNS
  643. *
  644. *   void * - New torus
  645. *   
  646. * AUTHOR
  647. *
  648. *   Dieter Bayer
  649. *   
  650. * DESCRIPTION
  651. *
  652. *   Copy a torus.
  653. *
  654. * CHANGES
  655. *
  656. *   Jun 1994 : Creation.
  657. *
  658. *   Sep 1994 : fixed memory leakage [DB]
  659. *
  660. ******************************************************************************/
  661.  
  662. static TORUS *Copy_Torus(OBJECT *Object)
  663. {
  664.   TORUS *New, *Torus = (TORUS *)Object;
  665.  
  666.   New = Create_Torus();
  667.  
  668.   /* Get rid of the transformation created in Create_Torus(). */
  669.  
  670.   Destroy_Transform(New->Trans);
  671.  
  672.   /* Copy torus. */
  673.  
  674.   *New = *Torus;
  675.  
  676.   New->Trans = Copy_Transform(Torus->Trans);
  677.  
  678.   return (New);
  679. }
  680.  
  681.  
  682.  
  683. /*****************************************************************************
  684. *
  685. * FUNCTION
  686. *
  687. *   Destroy_Torus
  688. *
  689. * INPUT
  690. *
  691. *   Object - Object
  692. *   
  693. * OUTPUT
  694. *
  695. *   Object
  696. *   
  697. * RETURNS
  698. *   
  699. * AUTHOR
  700. *
  701. *   Dieter Bayer
  702. *   
  703. * DESCRIPTION
  704. *
  705. *   Destroy a torus.
  706. *
  707. * CHANGES
  708. *
  709. *   Jun 1994 : Creation.
  710. *
  711. ******************************************************************************/
  712.  
  713. static void Destroy_Torus (OBJECT *Object)
  714. {
  715.   Destroy_Transform(((TORUS *)Object)->Trans);
  716.  
  717.   POV_FREE (Object);
  718. }
  719.  
  720.  
  721.  
  722. /*****************************************************************************
  723. *
  724. * FUNCTION
  725. *
  726. *   Compute_Torus_BBox
  727. *
  728. * INPUT
  729. *
  730. *   Torus - Torus
  731. *   
  732. * OUTPUT
  733. *
  734. *   Torus
  735. *   
  736. * RETURNS
  737. *   
  738. * AUTHOR
  739. *
  740. *   Dieter Bayer
  741. *   
  742. * DESCRIPTION
  743. *
  744. *   Calculate the bounding box of a torus.
  745. *
  746. * CHANGES
  747. *
  748. *   Jun 1994 : Creation.
  749. *
  750. ******************************************************************************/
  751.  
  752. void Compute_Torus_BBox(TORUS *Torus)
  753. {
  754.   DBL r1, r2;
  755.  
  756.   r1 = Torus->r;
  757.   r2 = Torus->R + Torus->r;
  758.  
  759.   Make_BBox(Torus->BBox, -r2, -r1, -r2, 2.0 * r2, 2.0 * r1, 2.0 * r2);
  760.  
  761.   Recompute_BBox(&Torus->BBox, Torus->Trans);
  762. }
  763.  
  764.  
  765.  
  766. /*****************************************************************************
  767. *
  768. * FUNCTION
  769. *
  770. *   Test_Thick_Cylinder
  771. *
  772. * INPUT
  773. *
  774. *   P  - Ray initial
  775. *   D  - Ray direction
  776. *   h1 - Height 1
  777. *   h2 - Height 2
  778. *   r1 - Square of inner radius
  779. *   r2 - Square of outer radius
  780. *   
  781. * OUTPUT
  782. *   
  783. * RETURNS
  784. *
  785. *   int - TRUE, if hit
  786. *   
  787. * AUTHOR
  788. *
  789. *   Dieter Bayer
  790. *   
  791. * DESCRIPTION
  792. *
  793. *   Test if a given ray defined in the lathe's coordinate system
  794. *   intersects a "thick" cylinder (rotated about y-axis).
  795. *
  796. * CHANGES
  797. *
  798. *   Jun 1994 : Creation.
  799. *
  800. ******************************************************************************/
  801.  
  802. int Test_Thick_Cylinder(VECTOR P, VECTOR  D, DBL h1, DBL  h2, DBL  r1, DBL  r2)
  803. {
  804.   DBL a, b, c, d;
  805.   DBL u, v, k, r, h;
  806.  
  807.   if (fabs(D[Y]) < EPSILON)
  808.   {
  809.     if ((P[Y] < h1) || (P[Y] > h2))
  810.     {
  811.       return(FALSE);
  812.     }
  813.   }
  814.   else
  815.   {
  816.     /* Intersect ray with the cap-plane. */
  817.  
  818.     k = (h2 - P[Y]) / D[Y];
  819.  
  820.     u = P[X] + k * D[X];
  821.     v = P[Z] + k * D[Z];
  822.  
  823.     if ((k > EPSILON) && (k < Max_Distance))
  824.     {
  825.       r = u * u + v * v;
  826.  
  827.       if ((r >= r1) && (r <= r2))
  828.       {
  829.         return(TRUE);
  830.       }
  831.     }
  832.  
  833.     /* Intersectionersect ray with the base-plane. */
  834.  
  835.     k = (h1 - P[Y]) / D[Y];
  836.  
  837.     u = P[X] + k * D[X];
  838.     v = P[Z] + k * D[Z];
  839.  
  840.     if ((k > EPSILON) && (k < Max_Distance))
  841.     {
  842.       r = u * u + v * v;
  843.  
  844.       if ((r >= r1) && (r <= r2))
  845.       {
  846.         return(TRUE);
  847.       }
  848.     }
  849.   }
  850.  
  851.   a = D[X] * D[X] + D[Z] * D[Z];
  852.  
  853.   if (a > EPSILON)
  854.   {
  855.     /* Intersect with outer cylinder. */
  856.  
  857.     b = P[X] * D[X] + P[Z] * D[Z];
  858.  
  859.     c = P[X] * P[X] + P[Z] * P[Z] - r2;
  860.  
  861.     d = b * b - a * c;
  862.  
  863.     if (d >= 0.0)
  864.     {
  865.       d = sqrt(d);
  866.  
  867.       k = (-b + d) / a;
  868.  
  869.       if ((k > EPSILON) && (k < Max_Distance))
  870.       {
  871.         h = P[Y] + k * D[Y];
  872.  
  873.         if ((h >= h1) && (h <= h2))
  874.         {
  875.           return(TRUE);
  876.         }
  877.       }
  878.  
  879.       k = (-b - d) / a;
  880.  
  881.       if ((k > EPSILON) && (k < Max_Distance))
  882.       {
  883.         h = P[Y] + k * D[Y];
  884.  
  885.         if ((h >= h1) && (h <= h2))
  886.         {
  887.           return(TRUE);
  888.         }
  889.       }
  890.     }
  891.  
  892.     /* Intersect with inner cylinder. */
  893.  
  894.     c = P[X] * P[X] + P[Z] * P[Z] - r1;
  895.  
  896.     d = b * b - a * c;
  897.  
  898.     if (d >= 0.0)
  899.     {
  900.       d = sqrt(d);
  901.  
  902.       k = (-b + d) / a;
  903.  
  904.       if ((k > EPSILON) && (k < Max_Distance))
  905.       {
  906.         h = P[Y] + k * D[Y];
  907.  
  908.         if ((h >= h1) && (h <= h2))
  909.         {
  910.           return(TRUE);
  911.         }
  912.       }
  913.  
  914.       k = (-b - d) / a;
  915.  
  916.       if ((k > EPSILON) && (k < Max_Distance))
  917.       {
  918.         h = P[Y] + k * D[Y];
  919.  
  920.         if ((h >= h1) && (h <= h2))
  921.         {
  922.           return(TRUE);
  923.         }
  924.       }
  925.     }
  926.   }
  927.  
  928.   return(FALSE);
  929. }
  930.  
  931.  
  932.  
  933.