 # federgraph-special-sample-scilab-2D

Oct 16th, 2014
337
0
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");