--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/mechanics/inst/ocframe/private/beamEB.m	Tue May 28 15:47:06 2013 +0000
@@ -0,0 +1,111 @@
+## Copyright (C) 2010 Johan Beke
+##
+## This software is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 3 of the License, or (at
+## your option) any later version.
+##
+## This software is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with this software; see the file COPYING.  If not, see
+## <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {} beamEB() 
+## Beam on an elastic bed. (Winkler)
+## NOT READY YET!
+## @end deftypefn
+
+function [v,M,V] = beamEB(L, b, div, pointloads, distributedloads, E, I, k)
+	#number of nodes
+	NN = div + 1 + 4;
+	#length of one division
+	delta = L / div;
+	#vector with the load (delta/2 at both size of the point)
+	q = zeros(1, NN-1);
+
+	#load the vector q
+	x = delta/2;
+	small= 10e-15;
+	for i=1:NN-2
+		#point loads are spread over the interval 
+		for index=1:rows(pointloads)
+			if (pointloads(index,1) < x + small && pointloads(index,1) > x - delta - small)
+				q(i) = pointloads(index,2)/delta;
+			endif
+		endfor
+		#dist loads
+		for index=1:rows(distributedloads)
+			#TODO: 
+		endfor
+		x+=delta;
+	endfor
+	# create and solve the system of finite differences
+	A = zeros(NN,NN);
+	B = zeros(NN,1);
+	
+	for node=3:NN-2
+		#first 2 and last 2 node are dummy nodes for constraints
+		if (node==3) # first node: M = 0 and V = 0
+			# M = 0 = v" 
+			A(node-1,node-1) = 1;
+			A(node-1,node)   =-2;
+			A(node-1,node+1) = 1;
+			# V = 0 = v'''  
+			# (using −1/2 	1 	0 	−1 	1/2 	: http://en.wikipedia.org/wiki/Finite_difference_coefficient)
+			A(node-2,node-2) = -1/2;
+			A(node-2,node-1) = 1;
+			A(node-2,node) = 0;
+			A(node-2,node+1) = -1;
+			A(node-2,node+2) = 1/2;
+		endif
+		if (node == NN-2) # last node: M = 0 and V = 0
+			# M = 0 = v" constraint on last row
+			A(node+1,node-1) = 1;
+			A(node+1,node)   =-2;
+			A(node+1,node+1) = 1;
+			# V = 0 = v''' constraint on row "node"  
+			# (using −1/2 	1 	0 	−1 	1/2 	: http://en.wikipedia.org/wiki/Finite_difference_coefficient)
+			A(node+2,node-2) = -1/2;
+			A(node+2,node-1) = 1;
+			A(node+2,node) = 0;
+			A(node+2,node+1) = -1;
+			A(node+2,node+2) = 1/2;
+		endif
+		# normal node:  EI v'''' + b k v = q
+		# pattern for v'''' around node: 1 -4 6 -4 1
+		# source: Abramowitz and Stegun
+		A(node, node-2) = 1/delta^4;
+		A(node, node-1) = -4/delta^4;
+		A(node, node) = 6/delta^4 + b*k/(E*I);
+		A(node, node+1) = -4/delta^4;
+		A(node, node+2) = 1/delta^4;
+		
+		B(node,1) = q(node-2)/(E*I);
+	endfor
+
+	#solve with LU decomposition
+	[L,U,P]=lu(A);
+	v=inv(U)*inv(L)*B;
+
+	#calculate moments and shear
+	M = zeros(NN);
+	V = zeros(NN);
+	for node=3:NN-2
+		# 1	-2	1
+		M(node)=-E*I/(delta^2)*(v(node-1)-2*v(node)+v(node+1));
+		# −1/2 	1 	0 	−1 	1/2
+		V(node)=-E*I/(delta^3)*(-0.5*v(node-2)+v(node-1)-v(node+1)+0.5*v(node+2));
+	endfor
+	#delete dummy node at output
+	v = v(3:NN-2);
+	M = M(3:NN-2);
+	V = V(3:NN-2);
+endfunction
+
+#[v,M,V]=beamEB(15,1,500,[3,500;6.50,800],[],30.0e6, 1*1.077^3/12, 20.0e-3/0.01^3);
+