home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l460 / 3.ddi / DEMOS.DI$ / KNOT.M < prev    next >
Encoding:
Text File  |  1993-03-07  |  3.2 KB  |  104 lines
  1. function knot();
  2. %KNOT    Tube surrounding a three-dimensional knot.
  3. % KNOT computes the parametric representation of tube-like
  4. % surfaces and displays the tube with SURF.
  5.  
  6. % Rouben Rostamian, March 1991.
  7. % Modified from Titan MATLAB to MATLAB V4.0, C. Moler, Sept. 91.
  8.  
  9. %    Copyright (c) 1984-93 by The MathWorks, Inc.
  10.  
  11. disp(' ');
  12. disp('Thanks to Rouben Rostamian, University of Maryland.')
  13.  
  14. % THE FOLLOWING PARAMETERS CAN BE CHANGED:
  15. m = 30;         % Number of grid points in each (circular) section of the tube.
  16. n = 60;         % Number of sections along the tube.
  17. R = .75;        % Radius of the tube.
  18. q = floor(n/3); % Symmetry index.  Try q=floor(n/3) (symetric) or q=floor(n/4)
  19.  
  20. % Do not change this!
  21. t = (0:n)/n;
  22.  
  23. % SPECIFY THE GENERATING FUNCTION HERE:
  24. % The generating function f0 must be 1-periodic.
  25. % f1 and f2 are the first and second derivatives of f0.
  26. a = 2; b = 3; c = 1.5;
  27. f0 =           sin(2*pi*t) +         a*sin(4*pi*t) -         b*cos(4*pi*t)/2 ...
  28.             +c*sin(6*pi*t); 
  29. f1 =    (2*pi)*cos(2*pi*t) +  a*(4*pi)*cos(4*pi*t) +  b*(4*pi)*sin(4*pi*t)/2 ...
  30.      +c*(6*pi)*cos(6*pi*t);
  31. f2 = -(2*pi)^2*sin(2*pi*t) -a*(4*pi)^2*sin(4*pi*t) +b*(4*pi)^2*cos(4*pi*t)/2 ...
  32.    -c*(6*pi)^2*sin(6*pi*t);
  33.  
  34. % DO NOT CHANGE ANYTHING BELOW THIS LINE.
  35.  
  36. % Extend f periodically to 2 period-intervals:
  37. f0 = [ f0(1:n) f0(1:n) ];
  38. f1 = [ f1(1:n) f1(1:n) ];
  39. f2 = [ f2(1:n) f2(1:n) ];
  40.  
  41. % [x10;x20;x30] is the parametric representation of the center-line of the tube:
  42. x10 = f0(1:n+1);
  43. x20 = f0(q+1:q+n+1);
  44. x30 = f0(2*q+1:2*q+n+1);
  45.  
  46. % [x11;x21;x31] is velocity (same as tangent) vector:
  47. x11 = f1(1:n+1);
  48. x21 = f1(q+1:q+n+1);
  49. x31 = f1(2*q+1:2*q+n+1);
  50.  
  51. % [x12;x22;x32] is acceleration vector:
  52. x12 = f2(1:n+1);
  53. x22 = f2(q+1:q+n+1);
  54. x32 = f2(2*q+1:2*q+n+1);
  55.  
  56. speed = sqrt(x11.^2 + x21.^2 + x31.^2);
  57.  
  58. % This is the dot-product of the velocity and acceleration vectors:
  59. velacc = x11.*x12 + x21.*x22 + x31.*x32;
  60.  
  61. % Here is the normal vector:
  62. nrml1 = speed.^2 .* x12 - velacc.*x11;
  63. nrml2 = speed.^2 .* x22 - velacc.*x21;
  64. nrml3 = speed.^2 .* x32 - velacc.*x31;
  65. normallength = sqrt(nrml1.^2 + nrml2.^2 + nrml3.^2);
  66.  
  67. % And here is the normalized normal vector:
  68. unitnormal1 = nrml1 ./ normallength;
  69. unitnormal2 = nrml2 ./ normallength;
  70. unitnormal3 = nrml3 ./ normallength;
  71.  
  72. % And the binormal vector ( B = T x N )
  73. binormal1 = (x21.*unitnormal3 - x31.*unitnormal2) ./ speed;
  74. binormal2 = (x31.*unitnormal1 - x11.*unitnormal3) ./ speed;
  75. binormal3 = (x11.*unitnormal2 - x21.*unitnormal1) ./ speed;
  76.  
  77. % s is the coordinate along the cicular cross-sections of the tube:
  78. s = (0:m)';
  79. s = (2*pi/m)*s;
  80.  
  81. % Finally, the parametric surface.  Each of x1, x2, x3 is an (m+1)x(n+1) matrix.
  82. % Rows represent coordinates along the tube.  Columns represent coordinates
  83. % in each (circular) cross-section of the tube.
  84.  
  85. xa1 = ones(m+1,1)*x10;
  86. xb1 = (cos(s)*unitnormal1 + sin(s)*binormal1);
  87. xa2 = ones(m+1,1)*x20;
  88. xb2 = (cos(s)*unitnormal2 + sin(s)*binormal2);
  89. xa3 = ones(m+1,1)*x30;
  90. xb3 = (cos(s)*unitnormal3 + sin(s)*binormal3);
  91. color = ones(m+1,1)*((0:n)*2/n-1);
  92.  
  93. x1 = xa1 + R*xb1;
  94. x2 = xa2 + R*xb2;
  95. x3 = xa3 + R*xb3;
  96.  
  97. clf
  98. colormap(hsv);
  99. surf(x1,x2,x3,color);
  100. axis([-4 4 -4 4 -4 4]);
  101. drawnow
  102.  
  103. disp(' ');
  104.