SHARE
TWEET

federgraph-special-sample-scilab-2D

federgraph Oct 16th, 2014 (edited) 196 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. //federgraph-special-sample-scilab-2D.sce
  2. //this file contains code for scilab 5.5.1
  3.  
  4. //the sample contains the special case, see
  5. //federgraph.blogspot.de/2014/10/federgraph-special-case-number-one.html
  6.  
  7. funcprot(0)
  8. function z = f(y)
  9. //the terms under the root    
  10. a1 = (x-x1)^2 + (y-y1)^2;
  11. a2 = (x-x2)^2 + (y-y2)^2;
  12. a3 = (x-x3)^2 + (y-y3)^2;
  13. //actual length of springs
  14. t1 = sqrt(a1);
  15. t2 = sqrt(a2);
  16. t3 = sqrt(a3);
  17. //expansion of springs
  18. f1 = (t1-l);
  19. f2 = (t2-l);
  20. f3 = (t3-l);
  21. //x-components of force
  22. u1 =f1 * (x-x1) / t1;
  23. u2 =f2 * (x-x2) / t2;
  24. u3 =f3 * (x-x3) / t3;
  25. //y-components of force
  26. v1 = f1 * (y-y1) / t1;
  27. v2 = f2 * (y-y2) / t2;
  28. v3 = f3 * (y-y3) / t3;
  29. //sum of components
  30. u = u1 + u2 + u3;
  31. v = v1 + v2 + v3;
  32. //absolute value of resultant force
  33. z = sqrt(u^2 + v^2);
  34. endfunction
  35.  
  36. function z = g(y)
  37.     //mirrored value
  38.     z = -f(y);
  39. endfunction
  40. funcprot(1)
  41.  
  42. //de.wikipedia.org/wiki/Gleichseitiges_Dreieck
  43. //setup of the constant equilateral triangle
  44. a = 100; //side length
  45. w3=sqrt(3); //Wurzel 3 (German!)
  46. h = w3/2 * a; //height
  47. ro = w3/3 * a; //radius of outer circle
  48. ri = w3/6; //radius of inner circle
  49.  
  50. //coordinates of the triangle, where the springs are attached
  51. x1 = -a/2;  
  52. x2 =  a/2;  
  53. x3 =0;
  54. y1 = -ri;
  55. y2 = -ri;
  56. y3 = ro;
  57.  
  58. clf()
  59. n = 128; //number of points in 2D plot
  60. x = 0; //a cross section cutting through the surface
  61. y=linspace(-100,100,n);
  62.  
  63. //l is the parameter, we explore special cases of the special sample
  64. //  for different values of l
  65. //l is the length of all springs, all springs have the same length l
  66. //l is special in the sample, as is the equilateral triangle, and k = 1
  67. //the spring constant k has already been eliminated from the equations
  68.  
  69. //a) -50 (all springs under tension)
  70. l = -a/2;
  71. plot(y, f, "k");
  72. plot(y, g, "k");
  73.  
  74. //b) 0 (force cone with peak at centroid)
  75. l = 0;
  76. plot(y, f, "k");
  77. plot(y, g, "k");
  78.  
  79. //c) 50 (single equilibrium at centroid)
  80. l = a/2;
  81. plot(y, f, "k");
  82. plot(y, g, "k");
  83.  
  84. //d) 57 (circles intersect at centroid)
  85. l = ro;
  86. plot(y, f, "c");
  87. plot(y, g, "c");
  88.  
  89. //e) 63 (touch down at corner C)
  90. l = w3 / (1 + w3) * a;
  91. plot(y, f, "m");
  92. plot(y, g, "m");
  93.  
  94. //f) 83 (bifurcation)
  95. l = 83;
  96. plot(y, f, "y");
  97. plot(y, g, "y");
  98.  
  99. //g) 86 (special case one, equilibrium at midpoint of side)
  100. l = h;
  101. plot(y, f, "r");
  102. plot(y, g, "r");
  103.  
  104. //h) 100 (no vertical edge / jump)
  105. l = a;
  106. plot(y, f, "g");
  107. plot(y, g, "g");
  108.  
  109. //i) 115 (centroid dome)
  110. l = 115;
  111. plot(y, f, "b");
  112. plot(y, g, "b");
  113.  
  114. //j) 237 (delta peak take off)
  115. l = w3/(w3-1);
  116. plot(y, f, "b");
  117. plot(y, g, "b");
RAW Paste Data
Top