federgraph-special-sample-scilab-2D

//federgraph-special-sample-scilab-2D.sce
//this file contains code for scilab 5.5.1

//the sample contains the special case, see
//federgraph.blogspot.de/2014/10/federgraph-special-case-number-one.html

funcprot(0)
function z = f(y)
//the terms under the root
a1 = (x-x1)^2 + (y-y1)^2;
a2 = (x-x2)^2 + (y-y2)^2;
a3 = (x-x3)^2 + (y-y3)^2;
//actual length of springs
t1 = sqrt(a1);
t2 = sqrt(a2);
t3 = sqrt(a3);
//expansion of springs
f1 = (t1-l);
f2 = (t2-l);
f3 = (t3-l);
//x-components of force
u1 =f1 * (x-x1) / t1;
u2 =f2 * (x-x2) / t2;
u3 =f3 * (x-x3) / t3;
//y-components of force
v1 = f1 * (y-y1) / t1;
v2 = f2 * (y-y2) / t2;
v3 = f3 * (y-y3) / t3;
//sum of components
u = u1 + u2 + u3;
v = v1 + v2 + v3;
//absolute value of resultant force
z = sqrt(u^2 + v^2);
endfunction

function z = g(y)
    //mirrored value
    z = -f(y);
endfunction
funcprot(1)

//de.wikipedia.org/wiki/Gleichseitiges_Dreieck
//setup of the constant equilateral triangle
a = 100; //side length
w3=sqrt(3); //Wurzel 3 (German!)
h = w3/2 * a; //height
ro = w3/3 * a; //radius of outer circle
ri = w3/6; //radius of inner circle

//coordinates of the triangle, where the springs are attached
x1 = -a/2;
x2 =  a/2;
x3 =0;
y1 = -ri;
y2 = -ri;
y3 = ro;

clf()
n = 128; //number of points in 2D plot
x = 0; //a cross section cutting through the surface
y=linspace(-100,100,n);

//l is the parameter, we explore special cases of the special sample
//  for different values of l
//l is the length of all springs, all springs have the same length l
//l is special in the sample, as is the equilateral triangle, and k = 1
//the spring constant k has already been eliminated from the equations

//a) -50 (all springs under tension)
l = -a/2;
plot(y, f, "k");
plot(y, g, "k");

//b) 0 (force cone with peak at centroid)
l = 0;
plot(y, f, "k");
plot(y, g, "k");

//c) 50 (single equilibrium at centroid)
l = a/2;
plot(y, f, "k");
plot(y, g, "k");

//d) 57 (circles intersect at centroid)
l = ro;
plot(y, f, "c");
plot(y, g, "c");

//e) 63 (touch down at corner C)
l = w3 / (1 + w3) * a;
plot(y, f, "m");
plot(y, g, "m");

//f) 83 (bifurcation)
l = 83;
plot(y, f, "y");
plot(y, g, "y");

//g) 86 (special case one, equilibrium at midpoint of side)
l = h;
plot(y, f, "r");
plot(y, g, "r");

//h) 100 (no vertical edge / jump)
l = a;
plot(y, f, "g");
plot(y, g, "g");

//i) 115 (centroid dome)
l = 115;
plot(y, f, "b");
plot(y, g, "b");

//j) 237 (delta peak take off)
l = w3/(w3-1);
plot(y, f, "b");
plot(y, g, "b");