--- a/main/mechanics/deprecated/ocframe_deprecated/SolveFrame.m	Fri Dec 09 12:29:53 2011 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,181 +0,0 @@
-## 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} {[@var{Reactions}, @var{Displacements}, @var{MemF}] = } SolveFrame(@var{joints}, @var{members}, @var{nodeloads}, @var{dist}, @var{point}) 
-##
-##
-## Solves a 2D frame with the matrix displacement method for
-## the following input parameters:
-##
-##    joints = [x , y, constraints ; ...]
-##
-##    constraints=[x , y, rotation] free=0, supported=1
-##
-##    members = [nodeN, nodeF, E, I, A; ...]
-##
-##    nodeloads = [node, Fx, Fy, Mz; ...]
-##
-##    loads on members:
-##
-##    dist = [membernum,FxN,FyN,FxF,FyF,a,b,local ; ...] for distributed loads
-##    where FxN and FyN are the loads on distance a from the near node
-##    (same with far node and distance b)
-##    local=1 if loads are on local axis, 0 if global
-##
-##    point = [membernum,Fx,Fy,a,local; ...]
-##    where Fx and Fy are the loads on distance a from the node near
-##    local=1 if loads are on local axis, 0 if global
-##
-##    Output is formated as follows (rownumber corresponds to
-##    node or member number):
-##
-##    Reactions = [Fx,Fy,Mz;...] where NaN if it was a non supported dof
-##
-##    Displacements = [x,y,rotation;...]
-##
-##    MemF = [FxN, FyN, MzN, FxF, FyF, MzF; ...] 
-## @end deftypefn
-
-function [Reactions,Displacements,MemF]=SolveFrame(joints,members,nodeloads,dist,point)
-	if (nargin < 5)
-		usage("[Reactions,Displacements,MemF]=SolveFrame(joints,members,nodeloads,dist,point) Use the help command for more information");
-	end
-	% calc info:
-	%	Elements Axis
-	%	y|
-	%	 |
-	%	 |___________________________ x
-	%      Near                         far
-	%	joints: [x, y, constraints; ...] 1= fixed
-	%	members [nodeN,nodeF,E,I,A; ...]
-	%	nodeloads [node,Fx,Fy,Mz; ...]
-	
-	P=D=zeros(rows(joints)*3,1);
-	%add nodal loads to P matrix
-	for load=1:rows(nodeloads)
-		c=node2code(nodeloads(load,1));
-		for i=0:2
-			P(c(1+i))=P(c(1+i))+nodeloads(load,2+i);
-		end
-	end
-	free=[]; %contains unconstrainted codenumbers
-	supported=[]; %contains constrainted codenumbers
-	for node=1:rows(joints)
-		c=node2code(node);
-		for i=3:5
-			if (joints(node,i)==0)
-				free=[free,c(i-2)];
-			else
-				supported=[supported,c(i-2)];
-			end
-		end
-	end
-	
-	%% global equation
-	%% { P_f }   [ K_ff   K_fs ]   { Delta_f }     { P_F_f }
-	%% {     } = [             ] . {         }  +  {       }
-	%% { P_s }   [ K_sf   K_ss ]   { Delta_s }     { P_F_s }
-	%% Solutions:
-	%% Delta_f = K_ff^-1 (P_f - P_F_f)
-	%% P_s = K_sf * Delta_f + P_F_s
-	
-	%stiffness matrix
-	[ K_ff, K_sf ] = GlobalStiffnessMatrixRegrouped (joints, members,free,supported); %K_ss, K_fs are not used and if not calculated script is faster
-
-	
-	%nodal forces and equivalent nodal ones
-	[P_F,MemFEF]=EquivalentJointLoads(joints,members,dist,point);
-	
-
-	%reduced matrices
-	P_f=P(free,:);
-	P_s=P(supported,:);
-	
-	P_F_f=P_F(free,:);
-	P_F_s=P_F(supported,:);
-	
-	%solution: find unknown displacements
-	
-	try
-		%A X = B => solve with cholesky => X = R\(R'\B)
-		r = chol (K_ff);
-		D_f=r\(r'\(P_f-P_F_f)); 
-	catch
-		error("System is unstable because K_ff is singular. Please check the support constraints!");
-	end_try_catch
-
-	%TODO: make use of iterative solver as an option. How????	
-
-	%[D_f, flag] = pcg(K_ff,P_f-P_F_f,1e-6,1000);
-	%if (flag==1)
-		%max iteration
-	%	printf('max iteration');
-	%	try
-	%		%A X = B => solve with cholesky => X = R\(R'\B)
-	%		r = chol (K_ff);
-	%		D_f=r\(r'\(P_f-P_F_f)); 
-	%	catch
-	%		error("System is unstable because K_ff is singular. Please check the support constraints!");
-	%	end_try_catch
-	%end
-	%if (flag==3)
-		%K_ff was found not positive definite
-	%	error("System is unstable because K_ff is singular. Please check the support constraints!");
-	%end
-	
-
-
-	
-	D(free,1)=D_f;
-	
-	%solution: find unknown (reaction) forces
-	P_s=K_sf*D_f+P_F_s;
-	P(supported,1)=P_s;
-	
-	
-	%solution: find unknown member forces
-	MemF=[];
-	for member=1:rows(members)
-		N=members(member,1);
-		F=members(member,2);
-		xN=joints(N,1);
-		yN=joints(N,2);
-		xF=joints(F,1);
-		yF=joints(F,2);
-		T=TransformationMatrix(xN,yN,xF,yF);
-		k=MemberStiffnessMatrix(sqrt((xN-xF)**2+(yN-yF)**2),members(member,3),members(member,4),members(member,5));
-		c=[node2code(N),node2code(F)];
-		MemF=[MemF;(k*T*D(c'))'];
-	end
-	MemF=MemF+MemFEF;%+fixed end forces
-
-	%format codenumbers to real output
-	%TODO: not efficient enough due to A=[A;newrow]
-	Displacements=[];
-	Reactions=[];
-	for i=0:rows(joints)-1
-		Displacements=[Displacements;D(1+3*i:3+3*i)'];
-		Reactions=[Reactions;P(1+3*i:3+3*i)'];
-	end
-	for i=1:rows(Reactions)
-		for j=1:columns(Reactions)
-			if (joints(i,j+2)==0)
-				Reactions(i,j)=NaN;
-			end
-		end
-	end
-end

--- a/main/mechanics/deprecated/ocframe_deprecated/Test.m	Fri Dec 09 12:29:53 2011 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,50 +0,0 @@
-## -*- texinfo -*-
-## @deftypefn {Function File} {} ocframe_ex3() 
-## Example of a planar frame.
-##
-## @end deftypefn
-function Test()
-	UNP120=[210e9,363.990e-8,17.000e-4];
-	RHS=[210e9,28.900e-8,8.730e-4];
-	joints=[0,0,1,1,0;
-	1.25,0,0,0,0;
-	1.575,0,0,0,0;
-	2.5,0,0,0,0;
-	3.425,0,0,0,0;
-	3.75,0,0,0,0;
-	5.,0,0,1,0;
-	0,1.5,0,0,0;
-	1.25,1.5,0,0,0;
-	2.5,1.5,0,0,0;
-	3.75,1.5,0,0,0;
-	5,1.5,0,0,0];
-	members=[1,2,UNP120;
-	2,3,UNP120;
-	3,4,UNP120;
-	4,5,UNP120;
-	5,6,UNP120;
-	6,7,UNP120;
-	8,9,UNP120;
-	9,10,UNP120;
-	10,11,UNP120;
-	11,12,UNP120;
-	1,8,RHS;
-	2,9,RHS;
-	4,10,RHS;
-	6,11,RHS;
-	7,12,RHS];
-
-	nodeloads=[3,0,-30e3,0;
-	5,0,-30e3,0;
-	12,30e3,0,0];
-	tic
-	[P1,D1,MemF1]=SolveFrame(joints,members,nodeloads,[],[]);
-	toc
-	tic
-	[P2,D2,MemF2]=SolveFrameOpt(joints,members,nodeloads,[],[]);
-	toc
-
-	P1-P2
-	D1-D2
-	MemF1-MemF2
-end

--- a/main/mechanics/deprecated/ocframe_deprecated/test_FixedEndForcesDist.m	Fri Dec 09 12:29:53 2011 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,43 +0,0 @@
-## 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/>.
-
-function [F1,F2,F3,F4,F5,F6]=test_FixedEndForcesDist(l,FxN,FyN,FxF,FyF,a,b)
-	% distributed load in local coords.
-	% a: distance from Near
-	% b: distance from Far
-	% FxN, FyN forces on start of load in local coordinates
-	% FxF, FyF forces on end of load in local coordinates
-	
-	x=linspace(a,l-b,100); #generate 100 intervals for riemann integration
-	dx=x(2)-x(1);
-	S=[];
-	for i=1:columns(x)
-		Fy = (FyF-FyN)/(l-b-a)*(x(i)-a)+FyN;
-		Fx = (FxF-FxN)/(l-b-a)*(x(i)-a)+FxN;
-		[ F1,F2,F3,F4,F5,F6 ] = test_FixedEndForcesPnt(l,x(i),Fx,Fy);
-		S=[S;[ F1,F2,F3,F4,F5,F6 ]];
-	end
-	for i=1:columns(S)
-		S(1,i)/=2;
-		S(rows(S),i)/=2;
-	end
-	F1 = sum(S(:,1))*dx;
-	F2 = sum(S(:,2))*dx;
-	F3 = sum(S(:,3))*dx;
-	F4 = sum(S(:,4))*dx;
-	F5 = sum(S(:,5))*dx;
-	F6 = sum(S(:,6))*dx;
-end

--- a/main/mechanics/deprecated/ocframe_deprecated/test_FixedEndForcesPnt.m	Fri Dec 09 12:29:53 2011 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,36 +0,0 @@
-## Copyright (C) 2011 Johan
-## 
-## This program 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 2 of the License, or
-## (at your option) any later version.
-## 
-## This program 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 Octave; see the file COPYING.  If not, see
-## <http://www.gnu.org/licenses/>.
-
-## test_FixedEndForcesPnt
-
-## Author: Johan <johan@johan-Satellite-L30>
-## Created: 2011-10-02
-
-function [ F1,F2,F3,F4,F5,F6 ] = test_FixedEndForcesPnt (l,a,Fx,Fy)
-	joints=[0,0,1,1,1;
-		a,0,0,0,0;
-		l,0,1,1,1];
-	nodeloads=[2,Fx,Fy,0];
-	members=[1,2,1,1,1;
-		2,3,1,1,1];
-	R=SolveFrame(joints,members,nodeloads,[],[]);
-	F1 = R(1,1);
-	F2 = R(1,2);
-	F3 = R(1,3);
-	F4 = R(3,1);
-	F5 = R(3,2);
-	F6 = R(3,3);
-endfunction

--- a/main/mechanics/inst/deprecated/area_poly2d.m	Fri Dec 09 12:29:53 2011 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,52 +0,0 @@
-%% Copyright (c) 2011 Juan Pablo Carbajal <carbajal@ifi.uzh.ch>
-%% 
-%%    This program 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
-%%    any later version.
-%%
-%%    This program 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 program. If not, see <http://www.gnu.org/licenses/>.
-
-%% -*- texinfo -*-
-%% @deftypefn {Function File} @var{A} = area_poly2d (@var{p})
-%% Calculates the Area of a 2D polygon.
-%%
-%% The polygon is described in @var{p}, where each row is a different vertex.
-%% The algorithm was adapted from P. Bourke web page
-%% @uref{http://local.wasp.uwa.edu.au/~pbourke/geometry/polyarea/}
-%%
-%% @seealso{inertia_moment_poly2d, center_mass_poly2d}
-%% @end deftypefn
-
-function A = area_poly2d(poly)
-
-  N = size(poly,1);
-  nxt = [2:N 1];
-  px = poly(:,1);
-  px_nxt = poly(nxt,1);
-  py = poly(:,2);
-  py_nxt = poly(nxt,2);
-
-  A = sum(px.*py_nxt - px_nxt.*py)/2;
-  
-end
-
-%!demo
-%! % A parametrized arbitrary triangle and its area
-%!
-%!  triangle = @(a,b,h) [0 0; b 0; a h];
-%!  h = linspace(0.1,1,10);
-%!  b = pi;
-%!  for i=1:length(h);
-%!    P = triangle(0.1,b,h(i));
-%!    area(i) = area_poly2d(P);
-%!  end
-%!
-%! % The area of the triangle is b*h/2
-%!  plot(h,area,'o',h,b*h/2); 

--- a/main/mechanics/inst/deprecated/center_mass_poly2d.m	Fri Dec 09 12:29:53 2011 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,60 +0,0 @@
-%% Copyright (c) 2011 Juan Pablo Carbajal <carbajal@ifi.uzh.ch>
-%% 
-%%    This program 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
-%%    any later version.
-%%
-%%    This program 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 program. If not, see <http://www.gnu.org/licenses/>.
-
-%% -*- texinfo -*-
-%% @deftypefn {Function File} @var{cm} = center_mass_poly2d (@var{p})
-%% Calculates the center of mass a 2D polygon.
-%%
-%% The polygon is described in @var{p}, where each row is a different vertex.
-%% The algorithm was adapted from P. Bourke web page
-%% @url{http://local.wasp.uwa.edu.au/~pbourke/geometry/polyarea/}
-%%
-%% @seealso{inertia_moment_poly2d, area_poly2d}
-%% @end deftypefn
-
-function cm = center_mass_poly2d(poly)
-
-  N = size(poly,1);
-  nxt = [2:N 1];
-  px = poly(:,1);
-  px_nxt = poly(nxt,1);
-  py = poly(:,2);
-  py_nxt = poly(nxt,2);
-  
-  cm = zeros(1,2);
-  cr_prod = (px.*py_nxt - px_nxt.*py);
-  
-  % Area
-  A = sum(cr_prod)/2;
-
-  % Center of mass  
-  cm(1) = sum( (px+px_nxt) .* cr_prod );
-  cm(2) = sum( (py+py_nxt) .* cr_prod );
-  cm = cm/(6*A);
-   
-end
-
-%!demo
-%! % The center of mass of this two triangles is the same
-%! % since both describe the same figure.
-%!
-%!  P = [0 0; 1 0; 0 1]; 
-%!  P2=[0 0; 0.1 0; 0.2 0; 0.25 0; 1 0; 0 1];
-%!  [center_mass_poly2d(P) center_mass_poly2d(P2)],
-%!
-%! % The centroid does not give the right answer
-%! 
-%!  [mean(P) mean(P2)],
-

--- a/main/mechanics/inst/deprecated/inertia_moment_ncpoly2d.m	Fri Dec 09 12:29:53 2011 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,118 +0,0 @@
-%% Copyright (c) 2011 Juan Pablo Carbajal <carbajal@ifi.uzh.ch>
-%% 
-%%    This program 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
-%%    any later version.
-%%
-%%    This program 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 program. If not, see <http://www.gnu.org/licenses/>.
-
-%% -*- texinfo -*-
-%% @deftypefn {Function File} @var{I} = inertia_moment_ncpoly2d (@var{p}, @var{m})
-%% Calculates the moment of inertia of a simple 2D polygon.
-%%
-%% The polygon is described in @var{p}, where each row is a different vertex.
-%% @var{m} is the total mass of the polygon, assumed uniformly distributed.
-%% The polygon is triangulated using delaunay algorithm and then the 
-%% Superposition Principle and the Parallel axis theorem are applied to each
-%% triangle.
-%%
-%% The position of the vertices is assumed to be given from the center of mass 
-%% of the polygon.
-%% To change a general polygon to this description you can use:
-%% @code{P = P - repmat(center_mass_poly2d(P),size(P,1))}.
-%%
-%% @seealso{inertia_moment_poly2d, center_mass_poly2d}
-%% @end deftypefn
-
-function I = inertia_moment_ncpoly2d (poly, M)
-  %% Get total area
-  A = area_poly2d (poly);
-
-  %% triangulate
-  T = delaunay (poly(:,1), poly(:,2));
-  nT = size(T,1);
-  
-% debug
-% triplot(T,poly(:,1),poly(:,2),'color',[0.8 0.8 0.8])
-% hold on
-% drawPolygon(poly,'k'); 
-% hold off;
-
-  I = 0;
-  for it = 1:nT
-    P = poly (T(it,:), :);
-    %% get centers of mass
-    cm =  center_mass_poly2d (P);
-    
-    % Check if triangle CoM is inside polygon
-    if ~inpolygon (cm(1), cm(2), poly(:,1), poly(:,2))
-      continue
-    end
-    
-    %% get the mass as fraction of total area
-    a = area_poly2d (P);
-    if a < 0
-      aux = P(1,:);
-      P(1,:) = P(2,:);
-      P(2,:) = aux;
-      a = -a;
-    end
-    m = M*a/A;
-
-% debug
-% patch(P(:,1),P(:,2),'facecolor','b','edgecolor','none');
-% pause
-
-
-
-    %% get moment of inertia
-    mo = inertia_moment_poly2d (P, m, cm);
-
-    %% Assemble: Superposition + parallel axis
-    I += mo + m * sumsq (cm);
-  end
-
-end
-
-%!demo
-%! % C shape polygon
-%!  poly = [0 0; 1 0; 1 0.25; 0.25 0.25; 0.25 0.75; 1 0.75; 1 1; 0 1];
-%! 
-%! % Take to center of mass
-%!  poly = poly - repmat(center_mass_poly2d(poly),size(poly,1),1);
-%!  A = area_poly2d(poly);
-%!
-%!  I = inertia_moment_ncpoly2d(poly,1)
-%!
-%!  % It should give (breaking C in rectangles)
-%!  r1 = [poly(1:3,:); poly(1,1) poly(3,2)]; a1 = area_poly2d(r1); m1 = abs(a1/A); 
-%!  c1 = center_mass_poly2d(r1); I1 = inertia_moment_poly2d(r1,m1,c1);
-%!  r2 = [r1(4,:); poly(4:5,:); poly(1,1) poly(5,2)]; a2 = area_poly2d(r2); m2 = abs(a2/A); 
-%!  c2 = center_mass_poly2d(r2); I2 = inertia_moment_poly2d(r2,m2,c2);
-%!  r3 = [poly(5:8,:); r2(4,:)]; a3 = area_poly2d(r3); m3 = abs(a3/A); 
-%!  c3 = center_mass_poly2d(r3); I3 = inertia_moment_poly2d(r3,m3,c3);
-%! 
-%! I1 + m1*sumsq(c1) + I2 + m2*sumsq(c2) + I3 + m3*sumsq(c3)
-
-%!test
-%! poly = [0 0; 1 0; 1 0.25; 0.25 0.25; 0.25 0.75; 1 0.75; 1 1; 0 1];
-%! poly = poly - repmat(center_mass_poly2d(poly),size(poly,1),1);
-%! A = area_poly2d(poly);
-%! I = inertia_moment_ncpoly2d(poly,1)
-%! % It should give (breaking C in rectangles)
-%! r1 = [poly(1:3,:); poly(1,1) poly(3,2)]; a1 = area_poly2d(r1); m1 = abs(a1/A); 
-%! c1 = center_mass_poly2d(r1); I1 = inertia_moment_poly2d(r1,m1,c1);
-%! r2 = [r1(4,:); poly(4:5,:); poly(1,1) poly(5,2)]; a2 = area_poly2d(r2); m2 = abs(a2/A); 
-%! c2 = center_mass_poly2d(r2); I2 = inertia_moment_poly2d(r2,m2,c2);
-%! r3 = [poly(5:8,:); r2(4,:)]; a3 = area_poly2d(r3); m3 = abs(a3/A); 
-%! c3 = center_mass_poly2d(r3); I3 = inertia_moment_poly2d(r3,m3,c3);
-%! I_shouldbe = I1 + m1*sumsq(c1) + I2 + m2*sumsq(c2) + I3 + m3*sumsq(c3);
-%! assert(I,I_shouldbe)
-

--- a/main/mechanics/inst/deprecated/inertia_moment_poly2d.m	Fri Dec 09 12:29:53 2011 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,69 +0,0 @@
-%% Copyright (c) 2011 Juan Pablo Carbajal <carbajal@ifi.uzh.ch>
-%% 
-%%    This program 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
-%%    any later version.
-%%
-%%    This program 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 program. If not, see <http://www.gnu.org/licenses/>.
-
-%% -*- texinfo -*-
-%% @deftypefn {Function File} @var{I} = inertia_moment_poly2d (@var{p}, @var{m}, @var{offset})
-%% Calculates the moment of inertia of a 2D star-shaped polygon.
-%%
-%% The polygon is described in @var{p}, where each row is a different vertex.
-%% @var{m} is the total mass of the polygon, assumed uniformly distributed.
-%% The optional argument @var{offset} is an origin translation vector. All vertex 
-%% are transformed to the reference frame with origin at @var{offset}. 
-%%
-%% This expression assumes that the polygon is star-shaped. The position of the
-%% vertices is assumed to be given from the center of mass of the polygon.
-%% To change a general polygon to this description you can use:
-%% @code{P = P - repmat(center_mass_poly2d(P),size(P,1))}.
-%% or call the function using the offset:
-%% @code{inertia_moment_poly2d (@var{p}, @var{m}, center_mass_poly2d(@var{p}))}
-%%
-%% @seealso{inertia_moment_ncpoly2d, center_mass_poly2d}
-%% @end deftypefn
-
-function I = inertia_moment_poly2d(poly,mass,offset=[])
-
-  numVerts = size(poly,1);
-  
-  if !isempty (offset)
-    V = bsxfun(@minus, poly, offset);
-    Vnext = bsxfun(@minus, poly([2:numVerts 1],:),offset);
-  else
-    V = poly; 
-    Vnext = poly([2:numVerts 1],:);
-  end
-  
-  %% Area of the parallelograms
-  A = sqrt(sumsq(cross([Vnext zeros(numVerts,1)],[V zeros(numVerts,1)],2),2));
-  %% Distance between points
-  B = dot(V,V,2) + dot(V,Vnext,2) + dot(Vnext,Vnext,2);
-  
-  C = sum(A.*B);
-  
-  I = mass*C/(6*sum(A));
-end
-
-%!demo
-%! 
-%! % The same triangle respect to one of its vertices described with two polygons
-%!  P = [0 0; 1 0; 0 1]; 
-%!  P2=[0 0; 0.1 0; 0.2 0; 0.25 0; 1 0; 0 1];
-%!
-%! % Now described from the center of mass of the triangle
-%! Pc = P - repmat(center_mass_poly2d(P), 3, 1);
-%! inertia_moment_poly2d(Pc,1)
-%!
-%! Pc = P2 -repmat(center_mass_poly2d(P2), size(P2,1), 1);
-%! inertia_moment_poly2d(Pc,1)
-

--- a/main/mechanics/inst/deprecated/second_moment_poly2d.m	Fri Dec 09 12:29:53 2011 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,70 +0,0 @@
-%% Copyright (c) 2011 Juan Pablo Carbajal <carbajal@ifi.uzh.ch>
-%% 
-%%    This program 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
-%%    any later version.
-%%
-%%    This program 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 program. If not, see <http://www.gnu.org/licenses/>.
-
-%% -*- texinfo -*-
-%% @deftypefn {Function File} @var{J} = second_moment_shape2d (@var{p})
-%% @deftypefnx{Function File} @var{J} = second_moment_shape2d (@var{p}),@var{type}, @var{matrix})
-%% Calculates the second moment of area of a 2D shape. 
-%%
-%% The polygon is described in @var{p}, where each row is a different vertex as seen
-%% from the center of mass of the shape (see @code{center_mass_shape2d}).
-%%
-%% The output @var{J} contains Ix, Ixy and Iy, in that order. It is assumed to be
-%% oriented counter clockwise, otherwise multiply output by -1.
-%%
-%% If @var{matrix} is @code{true} then @var{J} is the symmetric second moment of area
-%% matrix, instead of the vector @code{vech(J)}.
-%% @var{type} indicates how is the shape described. So far the only case handled
-%% is when @var{shape} is a 2D polygon, where each row of @var{shape} is a vertex.
-%%
-%% @seealso{inertia_moment_shape2d, center_mass_shape2d}
-%% @end deftypefn
-
-function J = second_moment_poly2d(shape, typep='polygon', matrix=false)
-
-  if strcmpi (typep, 'polygon')
-
-    [N dim] = size (shape);
-    nxt = [2:N 1];
-    px = shape(:,1);
-    px_nxt = shape(nxt,1);
-    py = shape(:,2);
-    py_nxt = shape(nxt,2);
-    
-    cr_prod = px.*py_nxt - px_nxt.*py;
-
-    J = zeros (3, 1);
-    J(1) = sum ( (py.^2 + py.*py_nxt + py_nxt.^2) .* cr_prod); % Jx
-    J(3) = sum ( (px.^2 + px.*px_nxt + px_nxt.^2) .* cr_prod); % Jy
-    J(2) = 0.5 * sum ( (px.*py_nxt + 2*(px.*py + px_nxt.*py_nxt) + px_nxt.*py) .* cr_prod);
-    J = J/12;
-
-  elseif strcmpi (typep, 'polyline')
-
-    error('second_moment_poly2d:Devel','Polyline, not implemented yet.');
-    #TODO
-  elseif strcmpi (typep, 'cbezier')
-
-    error('second_moment_poly2d:Devel', 'Cubic bezier, not implemented yet');
-    #TODO
-
-  end
-
-  if matrix
-    J = vech2mat (J);
-  end
-
-end
-

--- a/main/mechanics/inst/deprecated/vech2mat.m	Fri Dec 09 12:29:53 2011 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,55 +0,0 @@
-%% Copyright (c) 2010 Juan Pablo Carbajal <carbajal@ifi.uzh.ch>
-%% 
-%%    This program 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
-%%    any later version.
-%%
-%%    This program 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 program. If not, see <http://www.gnu.org/licenses/>.
-
-%% -*- texinfo -*-
-%% @deftypefn {Function File} {[ @var{r}, @var{c} ] = } vech2mat (@var{v}, @var{symm})
-%% Transforms a 1D array to a symmetric or antisymmetric matrix.
-%%
-%% Given a Nx1 array @var{v} describing the lower triangular part of a
-%% matrix (as obtained from @code{vech}), it returns the full matrix with the
-%% symmetry specified by @var{symm}. If @var{symm} is 1 the function returns
-%% a symmetric matrix. It returns an antisymmetric when @var{symm} is -1.
-%% If @var{symm} is omitted a symmetric matrix is returned.
-%% 
-%%
-%% @seealso{vech, ind2sub, sub2ind_tril}
-%% @end deftypefn
-
-function M = vech2mat(v, symm=1)
-
-    N = length (v);
-    dim = (sqrt ( 1 + 8*N ) - 1)/2;
-    [r c] = ind2sub_tril (dim, 1:N);
-    M = accumarray ([r; c].', v);
-    M += symm*tril (M, -1).';
-
-endfunction
-
-%!test %symmetric
-%! dim = 10;
-%! A = tril( floor ( 5*(2*rand(dim)-1) ) );
-%! A += A.';
-%! M = vech(A);
-%! M = vech2mat(M, 1);
-%! assert (A, M);
-
-%!test %antisymmetric
-%! dim = 10;
-%! A = tril( floor ( 5*(2*rand(dim)-1) ) );
-%! A -= A.';
-%! M = vech(A);
-%! M = vech2mat(M, -1);
-%! assert (A, M);
-

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/mechanics/inst/ocframe/deprecated/SolveFrame.m	Fri Dec 09 12:33:35 2011 +0000
@@ -0,0 +1,181 @@
+## 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} {[@var{Reactions}, @var{Displacements}, @var{MemF}] = } SolveFrame(@var{joints}, @var{members}, @var{nodeloads}, @var{dist}, @var{point}) 
+##
+##
+## Solves a 2D frame with the matrix displacement method for
+## the following input parameters:
+##
+##    joints = [x , y, constraints ; ...]
+##
+##    constraints=[x , y, rotation] free=0, supported=1
+##
+##    members = [nodeN, nodeF, E, I, A; ...]
+##
+##    nodeloads = [node, Fx, Fy, Mz; ...]
+##
+##    loads on members:
+##
+##    dist = [membernum,FxN,FyN,FxF,FyF,a,b,local ; ...] for distributed loads
+##    where FxN and FyN are the loads on distance a from the near node
+##    (same with far node and distance b)
+##    local=1 if loads are on local axis, 0 if global
+##
+##    point = [membernum,Fx,Fy,a,local; ...]
+##    where Fx and Fy are the loads on distance a from the node near
+##    local=1 if loads are on local axis, 0 if global
+##
+##    Output is formated as follows (rownumber corresponds to
+##    node or member number):
+##
+##    Reactions = [Fx,Fy,Mz;...] where NaN if it was a non supported dof
+##
+##    Displacements = [x,y,rotation;...]
+##
+##    MemF = [FxN, FyN, MzN, FxF, FyF, MzF; ...] 
+## @end deftypefn
+
+function [Reactions,Displacements,MemF]=SolveFrame(joints,members,nodeloads,dist,point)
+	if (nargin < 5)
+		usage("[Reactions,Displacements,MemF]=SolveFrame(joints,members,nodeloads,dist,point) Use the help command for more information");
+	end
+	% calc info:
+	%	Elements Axis
+	%	y|
+	%	 |
+	%	 |___________________________ x
+	%      Near                         far
+	%	joints: [x, y, constraints; ...] 1= fixed
+	%	members [nodeN,nodeF,E,I,A; ...]
+	%	nodeloads [node,Fx,Fy,Mz; ...]
+	
+	P=D=zeros(rows(joints)*3,1);
+	%add nodal loads to P matrix
+	for load=1:rows(nodeloads)
+		c=node2code(nodeloads(load,1));
+		for i=0:2
+			P(c(1+i))=P(c(1+i))+nodeloads(load,2+i);
+		end
+	end
+	free=[]; %contains unconstrainted codenumbers
+	supported=[]; %contains constrainted codenumbers
+	for node=1:rows(joints)
+		c=node2code(node);
+		for i=3:5
+			if (joints(node,i)==0)
+				free=[free,c(i-2)];
+			else
+				supported=[supported,c(i-2)];
+			end
+		end
+	end
+	
+	%% global equation
+	%% { P_f }   [ K_ff   K_fs ]   { Delta_f }     { P_F_f }
+	%% {     } = [             ] . {         }  +  {       }
+	%% { P_s }   [ K_sf   K_ss ]   { Delta_s }     { P_F_s }
+	%% Solutions:
+	%% Delta_f = K_ff^-1 (P_f - P_F_f)
+	%% P_s = K_sf * Delta_f + P_F_s
+	
+	%stiffness matrix
+	[ K_ff, K_sf ] = GlobalStiffnessMatrixRegrouped (joints, members,free,supported); %K_ss, K_fs are not used and if not calculated script is faster
+
+	
+	%nodal forces and equivalent nodal ones
+	[P_F,MemFEF]=EquivalentJointLoads(joints,members,dist,point);
+	
+
+	%reduced matrices
+	P_f=P(free,:);
+	P_s=P(supported,:);
+	
+	P_F_f=P_F(free,:);
+	P_F_s=P_F(supported,:);
+	
+	%solution: find unknown displacements
+	
+	try
+		%A X = B => solve with cholesky => X = R\(R'\B)
+		r = chol (K_ff);
+		D_f=r\(r'\(P_f-P_F_f)); 
+	catch
+		error("System is unstable because K_ff is singular. Please check the support constraints!");
+	end_try_catch
+
+	%TODO: make use of iterative solver as an option. How????	
+
+	%[D_f, flag] = pcg(K_ff,P_f-P_F_f,1e-6,1000);
+	%if (flag==1)
+		%max iteration
+	%	printf('max iteration');
+	%	try
+	%		%A X = B => solve with cholesky => X = R\(R'\B)
+	%		r = chol (K_ff);
+	%		D_f=r\(r'\(P_f-P_F_f)); 
+	%	catch
+	%		error("System is unstable because K_ff is singular. Please check the support constraints!");
+	%	end_try_catch
+	%end
+	%if (flag==3)
+		%K_ff was found not positive definite
+	%	error("System is unstable because K_ff is singular. Please check the support constraints!");
+	%end
+	
+
+
+	
+	D(free,1)=D_f;
+	
+	%solution: find unknown (reaction) forces
+	P_s=K_sf*D_f+P_F_s;
+	P(supported,1)=P_s;
+	
+	
+	%solution: find unknown member forces
+	MemF=[];
+	for member=1:rows(members)
+		N=members(member,1);
+		F=members(member,2);
+		xN=joints(N,1);
+		yN=joints(N,2);
+		xF=joints(F,1);
+		yF=joints(F,2);
+		T=TransformationMatrix(xN,yN,xF,yF);
+		k=MemberStiffnessMatrix(sqrt((xN-xF)**2+(yN-yF)**2),members(member,3),members(member,4),members(member,5));
+		c=[node2code(N),node2code(F)];
+		MemF=[MemF;(k*T*D(c'))'];
+	end
+	MemF=MemF+MemFEF;%+fixed end forces
+
+	%format codenumbers to real output
+	%TODO: not efficient enough due to A=[A;newrow]
+	Displacements=[];
+	Reactions=[];
+	for i=0:rows(joints)-1
+		Displacements=[Displacements;D(1+3*i:3+3*i)'];
+		Reactions=[Reactions;P(1+3*i:3+3*i)'];
+	end
+	for i=1:rows(Reactions)
+		for j=1:columns(Reactions)
+			if (joints(i,j+2)==0)
+				Reactions(i,j)=NaN;
+			end
+		end
+	end
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/mechanics/inst/ocframe/deprecated/Test.m	Fri Dec 09 12:33:35 2011 +0000
@@ -0,0 +1,50 @@
+## -*- texinfo -*-
+## @deftypefn {Function File} {} ocframe_ex3() 
+## Example of a planar frame.
+##
+## @end deftypefn
+function Test()
+	UNP120=[210e9,363.990e-8,17.000e-4];
+	RHS=[210e9,28.900e-8,8.730e-4];
+	joints=[0,0,1,1,0;
+	1.25,0,0,0,0;
+	1.575,0,0,0,0;
+	2.5,0,0,0,0;
+	3.425,0,0,0,0;
+	3.75,0,0,0,0;
+	5.,0,0,1,0;
+	0,1.5,0,0,0;
+	1.25,1.5,0,0,0;
+	2.5,1.5,0,0,0;
+	3.75,1.5,0,0,0;
+	5,1.5,0,0,0];
+	members=[1,2,UNP120;
+	2,3,UNP120;
+	3,4,UNP120;
+	4,5,UNP120;
+	5,6,UNP120;
+	6,7,UNP120;
+	8,9,UNP120;
+	9,10,UNP120;
+	10,11,UNP120;
+	11,12,UNP120;
+	1,8,RHS;
+	2,9,RHS;
+	4,10,RHS;
+	6,11,RHS;
+	7,12,RHS];
+
+	nodeloads=[3,0,-30e3,0;
+	5,0,-30e3,0;
+	12,30e3,0,0];
+	tic
+	[P1,D1,MemF1]=SolveFrame(joints,members,nodeloads,[],[]);
+	toc
+	tic
+	[P2,D2,MemF2]=SolveFrameOpt(joints,members,nodeloads,[],[]);
+	toc
+
+	P1-P2
+	D1-D2
+	MemF1-MemF2
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/mechanics/inst/ocframe/deprecated/test_FixedEndForcesDist.m	Fri Dec 09 12:33:35 2011 +0000
@@ -0,0 +1,43 @@
+## 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/>.
+
+function [F1,F2,F3,F4,F5,F6]=test_FixedEndForcesDist(l,FxN,FyN,FxF,FyF,a,b)
+	% distributed load in local coords.
+	% a: distance from Near
+	% b: distance from Far
+	% FxN, FyN forces on start of load in local coordinates
+	% FxF, FyF forces on end of load in local coordinates
+	
+	x=linspace(a,l-b,100); #generate 100 intervals for riemann integration
+	dx=x(2)-x(1);
+	S=[];
+	for i=1:columns(x)
+		Fy = (FyF-FyN)/(l-b-a)*(x(i)-a)+FyN;
+		Fx = (FxF-FxN)/(l-b-a)*(x(i)-a)+FxN;
+		[ F1,F2,F3,F4,F5,F6 ] = test_FixedEndForcesPnt(l,x(i),Fx,Fy);
+		S=[S;[ F1,F2,F3,F4,F5,F6 ]];
+	end
+	for i=1:columns(S)
+		S(1,i)/=2;
+		S(rows(S),i)/=2;
+	end
+	F1 = sum(S(:,1))*dx;
+	F2 = sum(S(:,2))*dx;
+	F3 = sum(S(:,3))*dx;
+	F4 = sum(S(:,4))*dx;
+	F5 = sum(S(:,5))*dx;
+	F6 = sum(S(:,6))*dx;
+end

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/mechanics/inst/ocframe/deprecated/test_FixedEndForcesPnt.m	Fri Dec 09 12:33:35 2011 +0000
@@ -0,0 +1,36 @@
+## Copyright (C) 2011 Johan
+## 
+## This program 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 2 of the License, or
+## (at your option) any later version.
+## 
+## This program 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 Octave; see the file COPYING.  If not, see
+## <http://www.gnu.org/licenses/>.
+
+## test_FixedEndForcesPnt
+
+## Author: Johan <johan@johan-Satellite-L30>
+## Created: 2011-10-02
+
+function [ F1,F2,F3,F4,F5,F6 ] = test_FixedEndForcesPnt (l,a,Fx,Fy)
+	joints=[0,0,1,1,1;
+		a,0,0,0,0;
+		l,0,1,1,1];
+	nodeloads=[2,Fx,Fy,0];
+	members=[1,2,1,1,1;
+		2,3,1,1,1];
+	R=SolveFrame(joints,members,nodeloads,[],[]);
+	F1 = R(1,1);
+	F2 = R(1,2);
+	F3 = R(1,3);
+	F4 = R(3,1);
+	F5 = R(3,2);
+	F6 = R(3,3);
+endfunction

--- a/main/mechanics/molecularDynamics/inst/demofunc1.m	Fri Dec 09 12:29:53 2011 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,26 +0,0 @@
-%% Copyright (c) 2011 Juan Pablo Carbajal <carbajal@ifi.uzh.ch>
-%% 
-%%    This program 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
-%%    any later version.
-%%
-%%    This program 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 program. If not, see <http://www.gnu.org/licenses/>.
-
-%% -*- texinfo -*-
-%% @deftypefn {Function File} {@var{f1},@var{f2}} = demofunc1 (@var{dx},@var{dv})
-%% function used in demos.
-%%
-%% @end deftypefn
-
-function [f1 f2] = demofunc1 (dx, dv)
-    f1 = lennardjones (dx, (0.15)^2, 1e-2);
-    f2 = -f1;
-endfunction
-

--- a/main/mechanics/molecularDynamics/inst/demofunc2.m	Fri Dec 09 12:29:53 2011 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,26 +0,0 @@
-%% Copyright (c) 2011 Juan Pablo Carbajal <carbajal@ifi.uzh.ch>
-%% 
-%%    This program 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
-%%    any later version.
-%%
-%%    This program 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 program. If not, see <http://www.gnu.org/licenses/>.
-
-%% -*- texinfo -*-
-%% @deftypefn {Function File} {@var{f1},@var{f2}} = demofunc2 (@var{x1},@var{x2},@var{v1},@var{v2})
-%% function used in demos.
-%%
-%% @end deftypefn
-
-function [f1 f2] = demofunc2 (x1,x2,v1,v2)
-   f1 = -15 *((x1-x2) + 0.2) -0.7*(v1-v2);
-   f2 = -f1;
-end
-

--- a/main/mechanics/molecularDynamics/inst/demofunc3.m	Fri Dec 09 12:29:53 2011 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,26 +0,0 @@
-%% Copyright (c) 2011 Juan Pablo Carbajal <carbajal@ifi.uzh.ch>
-%% 
-%%    This program 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
-%%    any later version.
-%%
-%%    This program 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 program. If not, see <http://www.gnu.org/licenses/>.
-
-%% -*- texinfo -*-
-%% @deftypefn {Function File} {@var{f1},@var{f2}} = demofunc3 (@var{dx},@var{dv})
-%% function used in demos.
-%%
-%% @end deftypefn
-
-function [f1 f2] = demofunc3 (dx, dv)
-    f1 = lennardjones (dx, (0.2)^2, 1);
-    f2 = -f1;
- end
- 

--- a/main/mechanics/molecularDynamics/inst/lennardjones.m	Fri Dec 09 12:29:53 2011 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,34 +0,0 @@
-%% Copyright (c) 2011 Juan Pablo Carbajal <carbajal@ifi.uzh.ch>
-%% 
-%%    This program 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
-%%    any later version.
-%%
-%%    This program 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 program. If not, see <http://www.gnu.org/licenses/>.
-
-%% -*- texinfo -*-
-%% @deftypefn {Function File} {@var{f},@var{v}} = lennardjones (@var{dx},@var{r02},@var{v0})
-%% Returns force and energy of the Lennard-Jonnes potential evaluated at @var{dx}.
-%%
-%% @end deftypefn
-
-function [f V] = lennardjones (dx,r02,V0)
-
-  r2 = sumsq(dx,2);
-  dr = dx ./ r2;
-
-  rr3 = (r02 ./ r2).^3;
-  rr6 = rr3.^2;
-
-  V = V0 * ( rr6 - 2*rr3 );
-
-  f = bsxfun(@times, 12 * V0 * ( rr6 - rr3 ), dr);
-
-end

--- a/main/mechanics/molecularDynamics/inst/mdsim.m	Fri Dec 09 12:29:53 2011 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,220 +0,0 @@
-%% Copyright (c) 2011 Juan Pablo Carbajal <carbajal@ifi.uzh.ch>
-%% 
-%%    This program 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
-%%    any later version.
-%%
-%%    This program 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 program. If not, see <http://www.gnu.org/licenses/>.
-
-%% -*- texinfo -*-
-%% @deftypefn {Function File} {@var{t},@var{pos},@var{vel}} = mdsim (@var{func},@var{tspan},@var{p0},@var{v0})
-%% Integrates a system of particles using velocity Verlet algorithm.
-%%
-%% @end deftypefn
-
-function [tspan y vy] = mdsim(func,tspan,x0,v0,varargin);
-
-  %% Parse options
-  options = {"mass","timestep","nonperiodic","substeps","boxlength"};
-  opt_given = {varargin{1:2:end}};
-  [tf idx] = ismember (options, opt_given);
-
-  % Parse masses
-  m = ones (size (x0,1), 1);
-  if tf(1)
-    m = varargin{idx(1)+1};
-  end
-
-  % Parse time step & substeps
-  nT = numel (tspan);
-  if nT == 1
-    tspan(2)=tspan;
-    tspan(1)=0;
-  end
-  
-  if tf(2) && tf(4)
-      error('mdsim:Error','You cannot define substeps and dt simultaneously.\n');
-  elseif tf(2) && ~tf(4)
-
-    dt = varargin{idx(2)+1};
-    if ~isscalar (dt) || any (dt < 0);
-      error('mdsim:Error','Timestep must be a positive scalar.\n');
-    end
-
-    substep = round( (tspan(2)-tspan(1)) / dt);
-
-  elseif ~tf(2) && tf(4)
-
-    substep = varargin{idx(4)+1};
-    if ~isscalar (substep) || any (substep < 2);
-      error('mdsim:Error','Substeps must be a natural number >= 2.\n');
-    end
-
-    dt = (tspan(2)-tspan(1)) / (substep-1);
-
-  else
-    substep = 1e1;
-    dt = (tspan(2)-tspan(1)) / (substep-1);
-  end
-
-  % Parse periodic bc or not
-  L = 1;
-  if tf(5)
-    L = varargin{idx(5)+1};
-    if ~isscalar (L) || any (L < 0);
-      error('mdsim:Error','Length of box side should be a positive scalar.\n');
-    end
-  end
-  if nargout(func) < 2 && nargin(func) < 4
-    error('mdsim:Error',['Interaction force must accept at least 4 arguments'...
-           'and return at least 2, when boundary conditions are periodic.\n']);
-  end
-  integrator = @(x_,v_,dt_) verletstep_boxed (x_, v_, m, dt_, func, L);
-
-  if tf(3)
-    if nargout(func) < 2 && nargin(func) < 2
-      error('mdsim:Error',['Interaction force must accept at least 2' ...
-             'arguments\nwhen boundary conditions are not periodic.\n']);
-    end
-    integrator = @(x_,v_,dt_) verletstep (x_, v_, m, dt_, func);
-  end
-
-  %% Allocate
-  [N dim] = size(x0);
-  y = zeros (N, dim, nT);
-
-  if nargout > 2
-    vy = y;
-    vy(:,:,1) = v0;  
-  end
-
-  %% Iterate ------------
-  if tf(3)
-    y(:,:,1) = x0;
-    auxX = x0;
-    auxV = v0;
-  else
-    auxX = x0;
-    ind = find(x0 < -L/2);
-    auxX(ind) = auxX(ind) + L * abs (round (auxX(ind)/L));
-    ind = find(x0 > -L/2);
-    auxX(ind) = auxX(ind) - L * abs (round (auxX(ind)/L));
-    y(:,:,1) = auxX;
-    auxV = v0;
-  end
-    
-  for it = 2:nT
-
-    for jt = 1:substep
-      [auxX auxV] = integrator (auxX, auxV, dt);
-    end
-
-    y(:,:,it) = auxX;
-    if nargout > 2
-     vy(:,:,it) = auxV;
-    end
-    
-  end
-  
-  if dim == 1
-    y = squeeze (y);
-    if nargout > 2
-      vy = squeeze (vy);
-    end
-  end
-  
-endfunction
-
-% Test arguments
-%!function [fx fy] = func4 (xq, x2, v1, v2)
-%!  fx = 0; fy = 0;
-%! end
-
-%!function [fx fy] = func2 (dx, dv)
-%!  fx = 0; fy = 0;
-%! end
-
-%!error (mdsim (@func2, [0 1], 0, 0, "nonperiodic", false))
-%!error (mdsim (@func4, [0 1], 0, 0, "nonperiodic", true))
-%!error (mdsim (@func4, [0 1], 0, 0))
-%!error (mdsim (@func2, [0 1], 0, 0,"timestep",-1))
-%!error (mdsim (@func2, [0 1], 0, 0,"timestep",[0 1]))
-%!error (mdsim (@func2, [0 1], 0, 0,"timestep",1,"substep",2))
-%!error (mdsim (@func2, [0 1], 0, 0,"substep",0))
-%!error (mdsim (@func2, [0 1], 0, 0,"substep",[0 1]))
-%!error (mdsim (@func2, [0 1], 0, 0,"boxlength",-1))
-%!error (mdsim (@func2, [0 1], 0, 0,"boxlength",[1 2]))
-%!error (mdsim ([], [0 1], 0, 0,"boxlength",[1 2]))
-
-%!demo
-%! N       = 6;
-%! P0      = linspace (-0.5+1/(N-1), 0.5-1/(N-1), N).';
-%! V0      = zeros (N, 1);
-%! nT      = 80;
-%! tspan   = linspace(0, 2, nT);
-%! [t P V] = mdsim ('demofunc1', tspan, P0, V0,'timestep',1e-3);
-%!
-%! figure (1)
-%! plot (P.',t,'.');
-%! xlabel ("Position")
-%! ylabel ("Time")
-%! axis tight
-%!
-%! disp("Initial values")
-%! disp ( sum ([V(:,1) P(:,1)], 1) )
-%! disp("Final values")
-%! disp ( sum ([V(:,end) P(:,end)], 1) )
-%!
-%! %-------------------------------------------------------------------
-%! % 1D particles with Lennard-Jones potential and periodic boundaries.
-%! % Velocity and position of the center of mass are conserved. 
-
-%!demo
-%! N     = 10;
-%! P0    = linspace (0,1,N).';
-%! V0    = zeros (N, 1);
-%! nT    = 80;
-%! tspan = linspace(0, 1, nT);
-%! [t P] = mdsim ('demofunc2', tspan, P0, V0,'nonperiodic', true);
-%!
-%! figure (1)
-%! plot (P.',t,'.');
-%! xlabel ("Position")
-%! ylabel ("Time")
-%!
-%! %-------------------------------------------------------------------
-%! % 1D array of springs with damping proportional to relative velocity and
-%! % nonzero rest length.
-
-
-%!demo
-%! 
-%! input("NOTE: It may take a while.\nPress Ctrl-C to cancel or <enter> to continue: ","s");
-%!
-%! N       = 4;
-%! [Px Py] = meshgrid (linspace (-0.5+0.5/(N-1), 0.5-0.5/(N-1), N));
-%! P0      = [Px(:) Py(:)];
-%! N       = size(P0,1);
-%! P0      = P0 + 0.1* 0.5/N *(2*rand (size (P0)) - 1);
-%! V0      = zeros (N, 2);
-%! nT      = 80;
-%! tspan   = linspace(0, 1, nT);
-%! [t P]   = mdsim ('demofunc3', tspan, P0, V0);
-%! x       = squeeze(P(:,1,:)); 
-%! y       = squeeze(P(:,2,:));
-%!
-%! figure (1)
-%! plot (x.',y.','.',x(:,end),y(:,end),'.k');
-%! xlabel ("X")
-%! ylabel ("Y")
-%! axis tight
-%!
-%! %-------------------------------------------------------------------
-%! % 2D particles with Lennard-Jones potential and periodic boundaries

--- a/main/mechanics/molecularDynamics/src/Makefile	Fri Dec 09 12:29:53 2011 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,7 +0,0 @@
-all: verletstep.oct verletstep_boxed.oct
-
-%.oct: %.cc
-	mkoctfile -Wall $<
-
-clean:
-	rm -f *.o octave-core core *.oct *~

--- a/main/mechanics/molecularDynamics/src/verletstep.cc	Fri Dec 09 12:29:53 2011 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,113 +0,0 @@
-#include <octave/oct.h>
-#include <octave/parse.h>
-
-DEFUN_DLD (verletstep, args, nargout, "Verlet velocity step")
-     {
-       int nargin = args.length();
-       octave_value_list retval;
-
-       unsigned int fcn_str = 0;
-       
-       if (nargin < 5)
-         print_usage ();
-       else
-         {
-            // Allocate inputs
-            octave_function *fcn;
-            if (args(4).is_function_handle () || args(4).is_inline_function ())
-              fcn_str = 0;
-            else if (args(4).is_string ())
-              fcn_str = 1; 
-            else
-              error ("verletstep: expected string,"," inline or function handle");
-
-            Matrix P = args(0).matrix_value ();
-            Matrix V = args(1).matrix_value ();
-            Matrix M = args(2).matrix_value ();
-            double dt = args(3).scalar_value ();
-            
-            dim_vector dv = P.dims();
-            octave_idx_type Nparticles = dv(0);
-            octave_idx_type Ndim = dv(1);
-
-            octave_value_list newargs;
-            Matrix A (dim_vector(Nparticles, Ndim),0);
-
-            // Evaluate interaction function
-            for (octave_idx_type ipart = 0; ipart < Nparticles; ipart++)
-            {
-              newargs(0) = P.row(ipart);
-              newargs(2) = V.row(ipart);
-
-              for (octave_idx_type jpart = ipart + 1; jpart < Nparticles; jpart++)
-              {
-                newargs(1) = P.row(jpart);
-                newargs(3) = V.row(jpart);
-
-                if (fcn_str)
-                  retval = feval (args (4).string_value (), newargs, nargout);
-                else
-                {
-                  fcn = args(4).function_value ();
-                  if (! error_state)
-                    retval = feval (fcn, newargs, nargout);
-                }
-
-                A.insert (retval(0).row_vector_value () + 
-                                      A.row (ipart), ipart, 0);
-                A.insert (retval(1).row_vector_value () + 
-                                      A.row (jpart), jpart, 0);
-                
-              }
-            }
-
-            for (octave_idx_type i = 0; i < A.rows (); i++) 
-              A.insert (A.row (i) / M(i), i, 0);
-
-            // Integrate by half time step velocity and full step position
-            V = V + A * dt/2;
-            P = P + V * dt;
-
-            A.fill(0);
-
-            // Evaluate interaction function
-            for (octave_idx_type ipart = 0; ipart < Nparticles; ipart++)
-            {
-              newargs(0) = P.row(ipart);
-              newargs(2) = V.row(ipart);
-
-              for (octave_idx_type jpart = ipart + 1; jpart < Nparticles; jpart++)
-              {
-                newargs(1) = P.row(jpart);
-                newargs(3) = V.row(jpart);
-
-                if (fcn_str)
-                  retval = feval (args (4).string_value (), newargs, nargout);
-                else
-                {
-                  fcn = args(4).function_value ();
-                  if (! error_state)
-                    retval = feval (fcn, newargs, nargout);
-                }
-
-                A.insert (retval(0).row_vector_value () + 
-                                                       A.row (ipart), ipart, 0);
-                A.insert (retval(1).row_vector_value () + 
-                                                       A.row (jpart), jpart, 0);
-                
-              }
-            }
-            for (octave_idx_type i = 0; i < A.rows (); i++) 
-              A.insert (A.row (i) / M(i), i, 0);
-
-            // Integrate velocity for the rest of the time step
-            V = V + A * dt/2;
-
-           retval(0) = P;
-           retval(1) = V;
-           retval(2) = A;
-       }
-       
-       return retval;
-     }
-

--- a/main/mechanics/molecularDynamics/src/verletstep_boxed.cc	Fri Dec 09 12:29:53 2011 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,135 +0,0 @@
-#include <octave/oct.h>
-#include <octave/parse.h>
-
-Matrix min_image (Matrix P, double L)
-{
-  // Returns the image of a particle inside a square box of side L centered in the origin
-  Matrix Pn = P;
-  for (octave_idx_type ip = 0; ip < P.rows (); ip++) 
-    for (octave_idx_type dim = 0; dim < P.columns (); dim++) 
-    {
-      if (P(ip,dim) > L/2)
-        Pn.fill(P(ip,dim) - L * labs (lrint (P(ip,dim) / L)), ip, dim, ip, dim);
-
-      else if (P(ip,dim) < -L/2)
-        Pn.fill(P(ip,dim) + L * labs (lrint (P(ip,dim) / L)), ip, dim, ip, dim);
-    }
-  return Pn;
-}
-
-/* When this version is used, the interaction function must depend on the
- relative position and velocities only. 
-*/
-DEFUN_DLD (verletstep_boxed, args, nargout, "Verlet velocity step in periodic space")
-     {
-       int nargin = args.length();
-       octave_value_list retval;
-
-       unsigned int fcn_str = 0;
-       
-       if (nargin < 6)
-         print_usage ();
-       else
-         {
-            // Allocate inputs
-            octave_function *fcn;
-            if (args(4).is_function_handle () || args(4).is_inline_function ())
-              fcn_str = 0;
-            else if (args(4).is_string ())
-              fcn_str = 1; 
-            else
-              error ("verletstep: expected string,"," inline or function handle");
-
-            Matrix P = args(0).matrix_value ();
-            Matrix V = args(1).matrix_value ();
-            Matrix M = args(2).matrix_value ();
-            double dt = args(3).scalar_value ();
-            
-            dim_vector dv = P.dims();
-            octave_idx_type Nparticles = dv(0);
-            octave_idx_type Ndim = dv(1);
-
-            octave_value_list newargs;
-            Matrix A (dim_vector(Nparticles, Ndim),0);
-
-            // Box containing the particles
-            double L = args(5).scalar_value ();
-
-            // Evaluate interaction function
-            for (octave_idx_type ipart = 0; ipart < Nparticles; ipart++)
-            {
-              
-              for (octave_idx_type jpart = ipart + 1; jpart < Nparticles; jpart++)
-              {
-                newargs(0) = min_image (P.row(ipart) - P.row(jpart), L);
-                newargs(1) = V.row(ipart) - V.row(jpart);
-
-                if (fcn_str)
-                  retval = feval (args (4).string_value (), newargs, nargout);
-                else
-                {
-                  fcn = args(4).function_value ();
-                  if (! error_state)
-                    retval = feval (fcn, newargs, nargout);
-                }
-
-                A.insert (retval(0).row_vector_value () + 
-                                      A.row (ipart), ipart, 0);
-                A.insert (retval(1).row_vector_value () + 
-                                      A.row (jpart), jpart, 0);
-                
-              }
-            }
-
-            for (octave_idx_type i = 0; i < A.rows (); i++) 
-              A.insert (A.row (i) / M(i), i, 0);
-
-            // Integrate by half time step velocity and full step position
-            V = V + A * dt/2;
-            P = P + V * dt;
-
-            A.fill(0);
-
-            // Evaluate interaction function
-            // Evaluate interaction function
-            for (octave_idx_type ipart = 0; ipart < Nparticles; ipart++)
-            {
-
-              for (octave_idx_type jpart = ipart + 1; jpart < Nparticles; jpart++)
-              {
-                newargs(0) = min_image (P.row(ipart) - P.row(jpart), L);
-                newargs(1) = V.row(ipart) - V.row(jpart);
-
-                if (fcn_str)
-                  retval = feval (args (4).string_value (), newargs, nargout);
-                else
-                {
-                  fcn = args(4).function_value ();
-                  if (! error_state)
-                    retval = feval (fcn, newargs, nargout);
-                }
-
-                A.insert (retval(0).row_vector_value () + 
-                                                       A.row (ipart), ipart, 0);
-                A.insert (retval(1).row_vector_value () + 
-                                                       A.row (jpart), jpart, 0);
-                
-              }
-            }
-            for (octave_idx_type i = 0; i < A.rows (); i++) 
-              A.insert (A.row (i) / M(i), i, 0);
-
-            // Integrate velocity for the rest of the time step
-            V = V + A * dt/2;
-
-            // Put all particles in a box of side L
-           P = min_image (P, L);
-
-           retval(0) = P;
-           retval(1) = V;
-           retval(2) = A;
-       }
-       
-       return retval;
-     }
-

--- a/main/mechanics/ocframe/inst/MSNForces.m	Fri Dec 09 12:29:53 2011 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,205 +0,0 @@
-## 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} {[@var{x}, @var{M}, @var{S}, @var{N}] =} MSNForces(@var{joints}, @var{members}, @var{dist}, @var{point}, @var{MemF}, @var{membernum}, @var{divisions}) 
-## 
-## This function returns the internal forces of a member for each position x. The member
-## is divided in 20 subelements if the argument is not given. The used sign convention is displayed in the help file.
-##
-## Input parameters are similar as with SolveFrame and PlotFrame with extra arguments:
-##    membernum = Number of the member to calculate
-##    divisions = Number of divisions for the member
-##
-## @end deftypefn
-
-function [x,M,S,N]=MSNForces(joints,members,dist,point,MemF,membernum,divisions)
-	if (nargin < 6)
-		usage("[x,M,S,N]=MSNForces(joints,members,dist,point,MemF,membernum,divisions) use the help command for more information");
-	end
-	if (nargin<7)
-		divisions=20;
-	end
-	Near=members(membernum,1);
-	Far=members(membernum,2);
-	L=sqrt((joints(Near,1)-joints(Far,1))^2+(joints(Near,2)-joints(Far,2))^2);
-	x=0:L/divisions:L;
-	M=zeros(1,columns(x));
-	S=zeros(1,columns(x));
-	N=zeros(1,columns(x));
-	FxN_end=MemF(membernum,1);
-	FyN_end=MemF(membernum,2);
-	MzN_end=MemF(membernum,3);
-	
-	%Moment forces
-	for i=1:columns(x)
-		%due to end forces
-		M(i)+=-MzN_end+FyN_end*x(i);
-		%due to dist
-		for j=1:rows(dist)
-			if (dist(j,1)==membernum)
-				FxN=dist(j,2);
-				FyN=dist(j,3);
-				FxF=dist(j,4);
-				FyF=dist(j,5);
-				a=dist(j,6);
-				b=dist(j,7);
-				local=dist(j,8);
-				if (local==0)
-					%convert FxN,... to local axes
-					temp=[FxN,FyN,0,FxF,FyF,0];
-					temp=(TransformationMatrix(joints(Near,1),joints(Near,2),joints(Far,1),joints(Far,2))')^-1*temp';
-					FxN=temp(1);
-					FyN=temp(2);
-					FxF=temp(4);
-					FyF=temp(5);
-				end
-				if (x(i)>L-b)
-					M(i)+= (FyN+FyF)/2*(L-a-b)* (x(i)-((FyN+2*FyF)*L+(2*a-b)*FyN+(a-2*b)*FyF)/(3*FyN+3*FyF));
-				end
-				if (x(i)>a && x(i)<=L-b)
-					%% bug :S zie bug.m
-					%M(i)+=(((3*FyN*x(i)+3*a*FyN)*L+(2*FyF-2*FyN)*x(i)^2+((-3*b-2*a)*FyN-a*FyF)*x(i)+(-3*a*b-2*a^2)*FyN-a^2*FyF)*((2*FyN*x(i)*L+(FyF-FyN)*x(i)^2+(-2*b*FyN-2*a*FyF)*x(i))/(2*L-2*b-2*a)-(2*a*FyN*L+(-2*a*b-a^2)*FyN-a^2*FyF)/(2*L-2*b-2*a)))/(6*FyN*L+(3*FyF-3*FyN)*x(i)+(-6*b-3*a)*FyN-3*a*FyF);
-					M(i)+=((3*FyN*x(i)^2-6*a*FyN*x(i)+3*a^2*FyN)*L+(FyF-FyN)*x(i)^3+(-3*b*FyN-3*a*FyF)*x(i)^2+((6*a*b+3*a^2)*FyN+3*a^2*FyF)*x(i)+(-3*a^2*b-2*a^3)*FyN-a^3*FyF)/(6*L-6*b-6*a);
-				end
-			end
-		end
-		%due to point load on member
-		for j=1:rows(point)
-			if (point(j,1)==membernum)
-				Fx=point(j,2);
-				Fy=point(j,3);
-				a=point(j,4);
-				local=point(j,5);
-				if (local==0)
-					%convert to global axes
-					temp=[Fx,Fy,0,0,0,0];
-					temp=(TransformationMatrix(joints(Near,1),joints(Near,2),joints(Far,1),joints(Far,2))')^-1*temp';
-					Fx=temp(1);
-					Fy=temp(2);
-				end
-				if (x(i)>a)
-					M(i)+=(x(i)-a)*Fy;
-				end
-			end
-		end
-		
-	end
-	M=M.*-1; %sign convention changed during development
-
-	%Shear forces
-	for i=1:columns(x)
-		%due to end forces
-		S(i)+=FyN_end;
-		%due to dist
-		for j=1:rows(dist)
-			if (dist(j,1)==membernum)
-				FxN=dist(j,2);
-				FyN=dist(j,3);
-				FxF=dist(j,4);
-				FyF=dist(j,5);
-				a=dist(j,6);
-				b=dist(j,7);
-				local=dist(j,8);
-				if (local==0)
-					%convert FxN,... to local axes
-					temp=[FxN,FyN,0,FxF,FyF,0];
-					temp=(TransformationMatrix(joints(Near,1),joints(Near,2),joints(Far,1),joints(Far,2))')^-1*temp';
-					FxN=temp(1);
-					FyN=temp(2);
-					FxF=temp(4);
-					FyF=temp(5);
-				end
-				if (x(i)>L-b)
-					S(i)+=((FyN+FyF)*L^2+((-2*b-2*a)*FyF-2*b*FyN)*L+b^2*FyN+(b^2+2*a*b)*FyF)/(2*L-2*b-2*a)-(2*a*FyN*L+(-2*a*b-a^2)*FyN-a^2*FyF)/(2*L-2*b-2*a);
-				end
-				if (x(i)>a && x(i)<=L-b)
-					S(i)+=(2*FyN*x(i)*L+(FyF-FyN)*x(i)^2+(-2*b*FyN-2*a*FyF)*x(i))/(2*L-2*b-2*a)-(2*a*FyN*L+(-2*a*b-a^2)*FyN-a^2*FyF)/(2*L-2*b-2*a);
-				end
-			end
-		end
-		%due to point load on member
-		for j=1:rows(point)
-			if (point(j,1)==membernum)
-				Fx=point(j,2);
-				Fy=point(j,3);
-				a=point(j,4);
-				local=point(j,5);
-				if (local==0)
-					%convert to global axes
-					temp=[Fx,Fy,0,0,0,0];
-					temp=(TransformationMatrix(joints(Near,1),joints(Near,2),joints(Far,1),joints(Far,2))')^-1*temp';
-					Fx=temp(1);
-					Fy=temp(2);
-				end
-				if (x(i)>a)
-					S(i)+=Fy;
-				end
-			end
-		end
-	end
-	%Normal forces
-	for i=1:columns(x)
-		%due to end forces
-		N(i)+=-FxN_end;
-		%due to dist
-		for j=1:rows(dist)
-			if (dist(j,1)==membernum)
-				FxN=dist(j,2);
-				FyN=dist(j,3);
-				FxF=dist(j,4);
-				FyF=dist(j,5);
-				a=dist(j,6);
-				b=dist(j,7);
-				local=dist(j,8);
-				if (local==0)
-					%convert FxN,... to local axes
-					temp=[FxN,FyN,0,FxF,FyF,0];
-					temp=(TransformationMatrix(joints(Near,1),joints(Near,2),joints(Far,1),joints(Far,2))')^-1*temp';
-					FxN=temp(1);
-					FyN=temp(2);
-					FxF=temp(4);
-					FyF=temp(5);
-				end
-				if (x(i)>L-b)
-					N(i)-=((FxN+FxF)*L^2+((-2*b-2*a)*FxF-2*b*FxN)*L+b^2*FxN+(b^2+2*a*b)*FxF)/(2*L-2*b-2*a)-(2*a*FxN*L+(-2*a*b-a^2)*FxN-a^2*FxF)/(2*L-2*b-2*a);
-				end
-				if (x(i)>a && x(i)<=L-b)
-					N(i)-=(2*FxN*x(i)*L+(FxF-FxN)*x(i)^2+(-2*b*FxN-2*a*FxF)*x(i))/(2*L-2*b-2*a)-(2*a*FxN*L+(-2*a*b-a^2)*FxN-a^2*FxF)/(2*L-2*b-2*a);
-				end
-			end
-		end
-		%due to point load on member
-		for j=1:rows(point)
-			if (point(j,1)==membernum)
-				Fx=point(j,2);
-				Fy=point(j,3);
-				a=point(j,4);
-				local=point(j,5);
-				if (local==0)
-					%convert to global axes
-					temp=[Fx,Fy,0,0,0,0];
-					temp=(TransformationMatrix(joints(Near,1),joints(Near,2),joints(Far,1),joints(Far,2))')^-1*temp';
-					Fx=temp(1);
-					Fy=temp(2);
-				end
-				if (x(i)>a)
-					N(i)-=Fx;
-				end
-			end
-		end
-	end
-end

--- a/main/mechanics/ocframe/inst/PlotDiagrams.m	Fri Dec 09 12:29:53 2011 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,198 +0,0 @@
-## 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} {} PlotDiagrams(@var{joints}, @var{members}, @var{dist}, @var{point}, @var{MemF}, @var{diagram}, @var{divisions}, @var{scale}) 
-##
-## This function plots the internal forces for all members. The force to be plotted can be selected with @var{diagram} 
-## which will be "M", "S" or "N" for the moment, shear or normal forces.
-##
-## Input parameters are similar as with SolveFrame and PlotFrame.
-## @end deftypefn
-
-function PlotDiagrams(joints,members,dist,point,MemF,diagram,divisions,scale)
-	%during development some errors occured on windows: http://wiki.octave.org/wiki.pl?OctaveForWindows see oct2mat and plot package
-	
-	%% scaling goes wrong when there al members have end moment 0
-	if (nargin < 6)
-		usage("PlotDiagrams(joints,members,dist,point,MemF,diagram,scale) use the help command for more information");
-	end
-	if (nargin<7)
-		divisions=20;
-	end
-	if (diagram=="M")
-		if (nargin<8)
-			scale=2/max(max(abs(MemF(:,[3,6]))))
-		end
-		%plot frame first
-		newplot();
-		clf();
-		hold off;
-		set (gca(), "dataaspectratio", [1,1,1]);
-		hold on;
-		
-		for i=1:rows(members)
-			Near=members(i,1); %near joint
-			Far=members(i,2);	%far joint
-			%make line from near to far joint
-			line([joints(Near,1),joints(Far,1)],[joints(Near,2),joints(Far,2)],"color","blue");
-			
-		end
-		%plot coord numbers
-		for i=1:rows(joints)
-			text(joints(i,1),joints(i,2),int2str(i)); %'color','blue'
-			hold on;
-		end
-		%plot member numbers
-		for i=1:rows(members)
-			text((joints(members(i,1),1)+joints(members(i,2),1))/2,(joints(members(i,1),2)+joints(members(i,2),2))/2,int2str(i),'color','blue');
-		end
-		
-		for i=1:rows(members)
-			[x,M]=MSNForces(joints,members,dist,point,MemF,i,divisions);
-			%finding maximum and minimum
-			[M_max,x_max]=max(M);
-			[M_min,x_min]=min(M);
-			%scale example scale factor: 2/max moment 
-			M=M.*scale; % moment curve consistent with deformation
-			Near=members(i,1);
-			Far=members(i,2);
-
-			%transformation to real member position!
-			L=sqrt((joints(Near,1)-joints(Far,1))^2+(joints(Near,2)-joints(Far,2))^2);
-			l=(joints(Far,1)-joints(Near,1))/L;
-			m=(joints(Far,2)-joints(Near,2))/L;
-			T=[l,m;-m,l]';
-			for i=1:columns(x)
-				temp=T*[x(i);M(i)]; %rotation
-				x(i)=temp(1)+joints(Near,1); %displacement x and y
-				M(i)=temp(2)+joints(Near,2);
-			end
-			plot(x,M,'r-');
-			%plot max and min
-			text(x(x_max),M(x_max),num2str(M_max,'%g'),'color','red');
-			text(x(x_min),M(x_min),num2str(M_min,'%g'),'color','red');
-		end
-	end
-	if (diagram=="S")
-		if (nargin<8)
-			scale=2/max(max(abs(MemF(:,[2,5]))));
-		end
-		%plot frame first
-		newplot();
-		clf();
-		hold off;
-		set (gca(), "dataaspectratio", [1,1,1]);
-		hold on;
-		
-		for i=1:rows(members)
-			Near=members(i,1); %near joint
-			Far=members(i,2);	%far joint
-			%make line from near to far joint
-			line([joints(Near,1),joints(Far,1)],[joints(Near,2),joints(Far,2)],"color","blue");
-			
-		end
-		%plot coord numbers
-		for i=1:rows(joints)
-			text(joints(i,1),joints(i,2),int2str(i)); %'color','blue'
-			hold on;
-		end
-		%plot member numbers
-		for i=1:rows(members)
-			text((joints(members(i,1),1)+joints(members(i,2),1))/2,(joints(members(i,1),2)+joints(members(i,2),2))/2,int2str(i),'color','blue');
-		end
-		
-		for i=1:rows(members)
-			[x,M,S]=MSNForces(joints,members,dist,point,MemF,i,divisions);
-			%finding maximum and minimum
-			[S_max,x_max]=max(S);
-			[S_min,x_min]=min(S);
-			%scale example scale factor: 2/max moment 
-			S=S.*scale;
-			Near=members(i,1);
-			Far=members(i,2);
-
-			%transformation to real member position!
-			L=sqrt((joints(Near,1)-joints(Far,1))^2+(joints(Near,2)-joints(Far,2))^2);
-			l=(joints(Far,1)-joints(Near,1))/L;
-			m=(joints(Far,2)-joints(Near,2))/L;
-			T=[l,m;-m,l]';
-			for i=1:columns(x)
-				temp=T*[x(i);S(i)];
-				x(i)=temp(1)+joints(Near,1);
-				S(i)=temp(2)+joints(Near,2);
-			end
-			plot(x,S,'r-');
-			%plot max and min
-			text(x(x_max),S(x_max),num2str(S_max,'%g'),'color','red');
-			text(x(x_min),S(x_min),num2str(S_min,'%g'),'color','red');
-		end
-	end
-	if (diagram=="N")
-		if (nargin<8)
-			scale=2/max(max(abs(MemF(:,[1,4]))));
-		end
-		%plot frame first
-		newplot();
-		clf();
-		hold off;
-		set (gca(), "dataaspectratio", [1,1,1]);
-		hold on;
-		
-		for i=1:rows(members)
-			Near=members(i,1); %near joint
-			Far=members(i,2);	%far joint
-			%make line from near to far joint
-			line([joints(Near,1),joints(Far,1)],[joints(Near,2),joints(Far,2)],"color","blue");
-			
-		end
-		%plot coord numbers
-		for i=1:rows(joints)
-			text(joints(i,1),joints(i,2),int2str(i)); %'color','blue'
-			hold on;
-		end
-		%plot member numbers
-		for i=1:rows(members)
-			text((joints(members(i,1),1)+joints(members(i,2),1))/2,(joints(members(i,1),2)+joints(members(i,2),2))/2,int2str(i),'color','blue');
-		end
-		
-		for i=1:rows(members)
-			[x,M,S,N]=MSNForces(joints,members,dist,point,MemF,i,divisions);
-			%finding maximum and minimum
-			[N_max,x_max]=max(N);
-			[N_min,x_min]=min(N);
-			%scale example scale factor: 2/max moment 
-			N=N.*scale;
-			Near=members(i,1);
-			Far=members(i,2);
-
-			%transformation to real member position!
-			L=sqrt((joints(Near,1)-joints(Far,1))^2+(joints(Near,2)-joints(Far,2))^2);
-			l=(joints(Far,1)-joints(Near,1))/L;
-			m=(joints(Far,2)-joints(Near,2))/L;
-			T=[l,m;-m,l]';
-			for i=1:columns(x)
-				temp=T*[x(i);N(i)];
-				x(i)=temp(1)+joints(Near,1);
-				N(i)=temp(2)+joints(Near,2);
-			end
-			plot(x,N,'r-');
-			%plot max and min
-			text(x(x_max),N(x_max),num2str(N_max,'%g'),'color','red');
-			text(x(x_min),N(x_min),num2str(N_min,'%g'),'color','red');
-		end
-	end
-end

--- a/main/mechanics/ocframe/inst/PlotFrame.m	Fri Dec 09 12:29:53 2011 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,78 +0,0 @@
-## 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} {} PlotFrame(@var{joints}, @var{members}, @var{D}, @var{factor}) 
-##
-##
-## Plots a 2D frame (with displacements if needed) using
-## the following input parameters:
-##
-##    joints = [x , y, constraints ; ...]
-##
-##    constraints=[x , y, rotation] free=0, supported=1
-##
-##    members = [nodeN, nodeF, E, I, A; ...]
-##
-##	  Optional arguments:
-##
-##    D = [x,y,rotation;...] Displacements as returned by SolveFrame
-##
-##    factor= Scaling factor for the discplacements (default: 10)
-##
-## @end deftypefn
-
-function PlotFrame(joints,members,D,factor)
-	if (nargin < 2)
-		usage("PlotFrame(joints,members,D,factor) D and factor ar optional. Use the help command for more information");
-	end
-	%during development some errors occured on windows: http://wiki.octave.org/wiki.pl?OctaveForWindows see oct2mat and plot package
-	newplot();
-	clf();
-	set (gca(), "dataaspectratio", [1,1,1]);
-	for i=1:rows(members)
-		N=members(i,1); %near joint
-		F=members(i,2);	%far joint
-		%make line from near to far joint
-		line([joints(N,1),joints(F,1)],[joints(N,2),joints(F,2)],"color","blue");
-		hold on;
-		if (nargin>=3)
-			%plot also displacements
-			if (nargin==3)
-				% default factor
-				factor=10;
-			end
-			cN=node2code(N);
-			cF=node2code(F);
-			line([joints(N,1)+D(N,1)*factor,joints(F,1)+D(F,1)*factor],[joints(N,2)+D(N,2)*factor,joints(F,2)+D(F,2)*factor],"color","red");
-			hold on;
-		end
-	end
-	%one meter extra as a border
-	set (gca (), "xlim", [min(joints(:,1))-1,max(joints(:,1))+1]);
-	set (gca (), "ylim", [min(joints(:,2))-1,max(joints(:,2))+1]);
-
-	%plot coord numbers
-	for i=1:rows(joints)
-		text(joints(i,1),joints(i,2),int2str(i)); %'color','blue'
-		hold on;
-	end
-	%plot member numbers
-	for i=1:rows(members)
-		text((joints(members(i,1),1)+joints(members(i,2),1))/2,(joints(members(i,1),2)+joints(members(i,2),2))/2,int2str(i),'color','blue');
-	end
-	hold off;
-end
\ No newline at end of file

--- a/main/mechanics/ocframe/inst/SolveFrame.m	Fri Dec 09 12:29:53 2011 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,184 +0,0 @@
-## 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} {[@var{Reactions}, @var{Displacements}, @var{MemF}] = } SolveFrame(@var{joints}, @var{members}, @var{nodeloads}, @var{dist}, @var{point}) 
-##
-##
-## Solves a 2D frame with the matrix displacement method for
-## the following input parameters:
-##
-##    joints = [x , y, constraints ; ...]
-##
-##    constraints=[x , y, rotation] free=0, supported=1
-##
-##    members = [nodeN, nodeF, E, I, A; ...]
-##
-##    nodeloads = [node, Fx, Fy, Mz; ...]
-##
-##    loads on members:
-##
-##    dist = [membernum,FxN,FyN,FxF,FyF,a,b,local ; ...] for distributed loads
-##    where FxN and FyN are the loads on distance a from the near node
-##    (same with far node and distance b)
-##    local=1 if loads are on local axis, 0 if global
-##
-##    point = [membernum,Fx,Fy,a,local; ...]
-##    where Fx and Fy are the loads on distance a from the node near
-##    local=1 if loads are on local axis, 0 if global
-##
-##    Output is formated as follows (rownumber corresponds to
-##    node or member number):
-##
-##    Reactions = [Fx,Fy,Mz;...] where NaN if it was a non supported dof
-##
-##    Displacements = [x,y,rotation;...]
-##
-##    MemF = [FxN, FyN, MzN, FxF, FyF, MzF; ...] 
-## @end deftypefn
-
-function [Reactions,Displacements,MemF]=SolveFrameOpt(joints,members,nodeloads,dist,point)
-	if (nargin < 5)
-		usage("[Reactions,Displacements,MemF]=SolveFrame(joints,members,nodeloads,dist,point) Use the help command for more information");
-	end
-	% calc info:
-	%	Elements Axis
-	%	y|
-	%	 |
-	%	 |___________________________ x
-	%      Near                         far
-	%	joints: [x, y, constraints; ...] 1= fixed
-	%	members [nodeN,nodeF,E,I,A; ...]
-	%	nodeloads [node,Fx,Fy,Mz; ...]
-	
-	P=D=zeros(rows(joints)*3,1);
-	%add nodal loads to P matrix
-	for load=1:rows(nodeloads)
-		c=node2code(nodeloads(load,1));
-		for i=0:2
-			P(c(1+i))=P(c(1+i))+nodeloads(load,2+i);
-		end
-	end
-	free=[]; %contains unconstrainted codenumbers
-	supported=[]; %contains constrainted codenumbers
-	for node=1:rows(joints)
-		c=node2code(node);
-		for i=3:5
-			if (joints(node,i)==0)
-				free=[free,c(i-2)];
-			else
-				supported=[supported,c(i-2)];
-			end
-		end
-	end
-	
-	%% global equation
-	%% { P_f }   [ K_ff   K_fs ]   { Delta_f }     { P_F_f }
-	%% {     } = [             ] . {         }  +  {       }
-	%% { P_s }   [ K_sf   K_ss ]   { Delta_s }     { P_F_s }
-	%% Solutions:
-	%% Delta_f = K_ff^-1 (P_f - P_F_f)
-	%% P_s = K_sf * Delta_f + P_F_s
-
-	%calculate transformation matrix and stiffnesmatrix for all members
-	k_array=MemberStiffnessMatrices(joints,members);
-	T_array=MemberTransformationMatrices(joints,members);
-
-	%stiffness matrix
-	[ K_ff, K_sf ] = GlobalStiffnessMatrixRegrouped (joints, members,free,supported,k_array,T_array); %K_ss, K_fs are not used and if not calculated script is faster
-
-	
-	%nodal forces and equivalent nodal ones
-	[P_F,MemFEF]=EquivalentJointLoads(joints,members,dist,point);
-	
-
-	%reduced matrices
-	P_f=P(free,:);
-	P_s=P(supported,:);
-	
-	P_F_f=P_F(free,:);
-	P_F_s=P_F(supported,:);
-	
-	%solution: find unknown displacements
-	
-	try
-		%A X = B => solve with cholesky => X = R\(R'\B)
-		r = chol (K_ff);
-		D_f=r\(r'\(P_f-P_F_f)); 
-		%D_f=cholinv(K_ff)*(P_f-P_F_f); %TODO: use this line but for testing purposes same method as old on
-	catch
-		error("System is unstable because K_ff is singular. Please check the support constraints!");
-	end_try_catch
-
-	%TODO: make use of iterative solver as an option. How???? (old code below endfunction)	
-
-	D(free,1)=D_f;
-	
-	%solution: find unknown (reaction) forces
-	P_s=K_sf*D_f+P_F_s;
-	P(supported,1)=P_s;
-	
-	
-	%solution: find unknown member forces
-	MemF=[];
-	for member=1:rows(members)
-		N=members(member,1);
-		F=members(member,2);
-		xN=joints(N,1);
-		yN=joints(N,2);
-		xF=joints(F,1);
-		yF=joints(F,2);
-		%T=TransformationMatrix(xN,yN,xF,yF);
-		%k=MemberStiffnessMatrix(sqrt((xN-xF)**2+(yN-yF)**2),members(member,3),members(member,4),members(member,5));
-		c=[node2code(N),node2code(F)];
-		MemF=[MemF;(k_array{member,1}*T_array{member,1}*D(c'))'];
-	end
-	MemF=MemF+MemFEF;%+fixed end forces
-
-	%format codenumbers to real output
-	%TODO: not efficient enough due to A=[A;newrow]
-	Displacements=[];
-	Reactions=[];
-	for i=0:rows(joints)-1
-		Displacements=[Displacements;D(1+3*i:3+3*i)'];
-		Reactions=[Reactions;P(1+3*i:3+3*i)'];
-	end
-	for i=1:rows(Reactions)
-		for j=1:columns(Reactions)
-			if (joints(i,j+2)==0)
-				Reactions(i,j)=NaN;
-			end
-		end
-	end
-end
-
-
-%[D_f, flag] = pcg(K_ff,P_f-P_F_f,1e-6,1000);
-%if (flag==1)
-	%max iteration
-%	printf('max iteration');
-%	try
-%		%A X = B => solve with cholesky => X = R\(R'\B)
-%		r = chol (K_ff);
-%		D_f=r\(r'\(P_f-P_F_f)); 
-%	catch
-%		error("System is unstable because K_ff is singular. Please check the support constraints!");
-%	end_try_catch
-%end
-%if (flag==3)
-	%K_ff was found not positive definite
-%	error("System is unstable because K_ff is singular. Please check the support constraints!");
-%end

--- a/main/mechanics/ocframe/inst/SolveFrameCases.m	Fri Dec 09 12:29:53 2011 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,158 +0,0 @@
-## 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} {[@var{results}] =} SolveFrameCases (@var{joints}, @var{members}, @var{loadcases}) 
-##
-##
-## Solves a 2D frame with the matrix displacement method for
-## the following input parameters:
-##
-##    joints = [x , y, constraints ; ...]
-##
-##    constraints=[x , y, rotation] free=0, supported=1
-##
-##    members = [nodeN, nodeF, E, I, A; ...]
-##
-##    loadcases is a struct array with for each loadcase the fields
-##
-##     - nodeloads = [node, Fx, Fy, Mz; ...]
-##
-##     - dist = [membernum,FxN,FyN,FxF,FyF,a,b,local ; ...]
-##
-##     - point = [membernum,Fx,Fy,a,local; ...]
-##
-##    input is as for the function SolveFrame.
-## 
-##    Output is a struct array with the fields: Displacements, Reactions and MemF  
-##
-##    (output formated as for the function SolveFrame.)
-## @end deftypefn
-
-function [results]=SolveFrameCases(joints,members,loadcases)
-	if (nargin != 3)
-		usage("[Reactions,Displacements,MemF]=SolveFrame(joints,members,loadcases) Use the help command for more information");
-	end
-	
-	free=[]; %contains unconstrainted codenumbers
-	supported=[]; %contains constrainted codenumbers
-	for node=1:rows(joints)
-		c=node2code(node);
-		for i=3:5
-			if (joints(node,i)==0)
-				free=[free,c(i-2)];
-			else
-				supported=[supported,c(i-2)];
-			end
-		end
-	end
-	
-	%stiffness matrix
-	[ K_ff, K_sf ] = GlobalStiffnessMatrixRegrouped (joints, members,free,supported); %K_ss, K_fs are not used and if not calculated script is faster
-
-	%from now  nodeloads,dist,point must depend on case: implement with ? !!!!!!!!
-	for lc=1:length(loadcases)
-
-		%convert loadcases to other variables
-		nodeloads=loadcases(lc).nodeloads;
-		dist=loadcases(lc).dist;
-		point=loadcases(lc).point;
-	
-		P=D=zeros(rows(joints)*3,1);
-		%add nodal loads to P matrix
-		for load=1:rows(nodeloads)
-			c=node2code(nodeloads(load,1));
-			for i=0:2
-				P(c(1+i))=P(c(1+i))+nodeloads(load,2+i);
-			end
-		end
-
-		%nodal forces and equivalent nodal ones
-		[P_F,MemFEF]=EquivalentJointLoads(joints,members,dist,point);
-		
-
-		%reduced matrices
-		P_f=P(free,:);
-		P_s=P(supported,:);
-		
-		P_F_f=P_F(free,:);
-		P_F_s=P_F(supported,:);
-		
-
-		%solution: find unknown displacements
-		%can be quicker for large systems with sparse matrix ! 
-		
-		[D_f, flag] = pcg(K_ff,P_f-P_F_f,1e-6,1000);
-		if (flag==1)
-			%max iteration
-			printf('max iteration');
-			try
-				%A X = B => solve with cholesky => X = R\(R'\B)
-				r = chol (K_ff);
-				D_f=r\(r'\(P_f-P_F_f)); 
-			catch
-				error("System is unstable because K_ff is singular. Please check the support constraints!");
-			end_try_catch
-		end
-		if (flag==3)
-			%K_ff was found not positive defnite
-			error("System is unstable because K_ff is singular. Please check the support constraints!");
-		end
-		
-		
-		
-		D(free,1)=D_f;
-		
-		%solution: find unknown (reaction) forces
-		P_s=K_sf*D_f+P_F_s;
-		P(supported,1)=P_s;
-		
-		
-		%find unknown member forces
-		MemF=[];
-		for member=1:rows(members)
-			N=members(member,1);
-			F=members(member,2);
-			xN=joints(N,1);
-			yN=joints(N,2);
-			xF=joints(F,1);
-			yF=joints(F,2);
-			T=TransformationMatrix(xN,yN,xF,yF);
-			k=MemberStiffnessMatrix(sqrt((xN-xF)**2+(yN-yF)**2),members(member,3),members(member,4),members(member,5));
-			c=[node2code(N),node2code(F)];
-			MemF=[MemF;(k*T*D(c'))'];
-		end
-		MemF=MemF+MemFEF;%+fixed end forces
-
-		%format codenumbers to real output
-		Displacements=[];
-		Reactions=[];
-		for i=0:rows(joints)-1
-			Displacements=[Displacements;D(1+3*i:3+3*i)'];
-			Reactions=[Reactions;P(1+3*i:3+3*i)'];
-		end
-		for i=1:rows(Reactions)
-			for j=1:columns(Reactions)
-				if (joints(i,j+2)==0)
-					Reactions(i,j)=NaN;
-				end
-			end
-		end
-		results(lc).Displacements=Displacements;
-		results(lc).Reactions=Reactions;
-		results(lc).MemF=MemF;
-	end
-end
\ No newline at end of file

--- a/main/mechanics/ocframe/inst/ocframe_ex1.m	Fri Dec 09 12:29:53 2011 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,24 +0,0 @@
-## -*- texinfo -*-
-## @deftypefn {Function File} {} ocframe_ex1() 
-## Example of a planar frame.
-##
-## @end deftypefn
-function [P,D,MemF]=ocframe_ex1()
-	joints=[0,0,1,1,1;
-	4,4,0,0,0;
-	8,4,1,1,0];
-
-	EIA=[210e9,23130*(10^-2)^4,84.5*(10^-2)^2];%IPE400
-
-	members=[1,2,EIA;2,3,EIA];
-
-	nodeloads=[];
-
-	dist=[1,0,-2e3,0,-2e3,0,0,0;2,0,-4e3,0,-4e3,0,0,1];
-	point=[];%1,0,-3e4,3,1
-
-	[P,D,MemF]=SolveFrame(joints,members,nodeloads,dist,point);
-	%PlotFrame(joints,members,D,10);
-	%plot moment diagram
-	PlotDiagrams(joints,members,dist,point,MemF,"S");
-end

--- a/main/mechanics/ocframe/inst/ocframe_ex2.m	Fri Dec 09 12:29:53 2011 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,36 +0,0 @@
-## -*- texinfo -*-
-## @deftypefn {Function File} {} ocframe_ex2() 
-## Example of a beam.
-##
-## @end deftypefn
-function [P,D,MemF]=ocframe_ex2()
-	joints=[0,0,1,1,0;
-	4,0,0,1,0;
-	8,0,0,1,0];
-	%Each 4 meter a node => 0,0 and 4,0 and 8,0 are the coordinates.
-	%The first node is a hinge and thus supported in the x and y direction => 1,1,0 for the constraints
-	%The following nodes are just rollers and thus supported in the y direction => 0,1,0 for the constraints 
-
-	EIA=[200e6,200e6*(10^-3)^4,6000*(10^-3)^2];
-	%EIA as a single vector to be used afterwards
-
-	members=[1,2,EIA;2,3,EIA];
-	%2 members, connection node 1 to node 2 and node 2 to node 3
-
-	nodeloads=[];
-	%there aren't nodal nodes in this example
-
-	dist=[1,0,-4e3,0,-4e3,0,0,1;
-	2,0,-4e3,0,-4e3,0,0,1];
-	%both members have a distributed load which takes the full length of the member. Notice the - sign caused
-	%by the axes conventions
-	%as we are working with newtons and meters the load is -4e3 N and not -4 kN
-	
-
-	[P,D,MemF]=SolveFrame(joints,members,nodeloads,dist,[])
-	%solve the frame with 
-	% P: reactions
-	% D: displacements
-	% MemF: the member end forces
-	
-end
\ No newline at end of file

--- a/main/mechanics/ocframe/inst/ocframe_ex3.m	Fri Dec 09 12:29:53 2011 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,47 +0,0 @@
-## -*- texinfo -*-
-## @deftypefn {Function File} {} ocframe_ex3() 
-## Example of a planar frame.
-##
-## @end deftypefn
-function [P,D,MemF]=ocframe_ex3()
-	UNP120=[210e9,363.990e-8,17.000e-4];
-	RHS=[210e9,28.900e-8,8.730e-4];
-	joints=[0,0,1,1,0;
-	1.25,0,0,0,0;
-	1.575,0,0,0,0;
-	2.5,0,0,0,0;
-	3.425,0,0,0,0;
-	3.75,0,0,0,0;
-	5.,0,0,1,0;
-	0,1.5,0,0,0;
-	1.25,1.5,0,0,0;
-	2.5,1.5,0,0,0;
-	3.75,1.5,0,0,0;
-	5,1.5,0,0,0];
-	members=[1,2,UNP120;
-	2,3,UNP120;
-	3,4,UNP120;
-	4,5,UNP120;
-	5,6,UNP120;
-	6,7,UNP120;
-	8,9,UNP120;
-	9,10,UNP120;
-	10,11,UNP120;
-	11,12,UNP120;
-	1,8,RHS;
-	2,9,RHS;
-	4,10,RHS;
-	6,11,RHS;
-	7,12,RHS];
-
-	nodeloads=[3,0,-30e3,0;
-	5,0,-30e3,0;
-	12,30e3,0,0];
-	
-	if (nargout>0)
-		[P,D,MemF]=SolveFrame(joints,members,nodeloads,[],[]);
-	else
-		[P,D,MemF]=SolveFrame(joints,members,nodeloads,[],[])
-		PlotFrame(joints,members,D,1);
-	end
-end

--- a/main/mechanics/ocframe/inst/ocframe_exLC.m	Fri Dec 09 12:29:53 2011 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,68 +0,0 @@
-## -*- texinfo -*-
-## @deftypefn {Function File} {} ocframe_exLC() 
-## Example of a beam with generation of eurocode ULS
-## load cases
-##
-## @end deftypefn
-
-function ocframe_exLC()
-
-	joints=[0,0,1,1,0;
-		4,0,0,1,0;
-		8,0,0,1,0;
-		12,0,0,1,0;
-		16,0,0,1,0];
-
-	EIA=[200e6,200e6*(10^-3)^4,6000*(10^-3)^2];
-	members=[1,2,EIA;
-	2,3,EIA;
-	3,4,EIA;
-	4,5,EIA];
-
-	%there is a variable load 4kN/m to be combined with 2kN/m permanent
-	%eurocode ULS
-	%Permanent actions: if favourable 1 if not 1,35
-	%Variable actions: if favourable 0 if not 1,5
-	%Category A: domestic phi0 = 0,7	phi1 = 0,5 phi2 = 0,3
-	% each case consists of permanent load * 1,35 plus one of the variable cases
-	% which makes 16 cases (try for each variable load 0 and 1,5 as the factor)
-	loadcases={};
-	for i=0:15
-		lc=toascii(dec2bin(i,4)).-48;
-		%no point loads and nodal loads
-		loadcases(i+1).nodeloads=[];
-		loadcases(i+1).point=[];
-		%dist load depends on case
-		dist=[];
-		for j=1:4
-			%add permanent
-			dist=[dist;[j,0,-2e3*1.35,0,-2e3*1.35,0,0,1]];
-			if (lc(j)==1)
-				dist=[dist;[j,0,-4e3*0.7*1.5,0,-4e3*0.7*1.5,0,0,1]];
-			end
-			loadcases(i+1).dist=dist;
-		end
-	end
-
-	[results]=SolveFrameCases(joints,members,loadcases);
-
-	%moment diagram envelope
-	moments=[];
-	for lc=1:16
-		x=[];
-		M=[];
-		for i=1:4
-			[u,Mo]=MSNForces(joints,members,loadcases(lc).dist,loadcases(lc).point,results(lc).MemF,i,40);
-			x=[x,u.+(4*(i-1))];
-			M=[M,Mo];
-		end
-		moments=[moments;M];
-	end
-
-	plot(x,max(moments),x,min(moments))
-
-end
-
-
-
-	
\ No newline at end of file

--- a/main/mechanics/ocframe/inst/ocframe_railwaybridge.m	Fri Dec 09 12:29:53 2011 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,46 +0,0 @@
-## -*- texinfo -*-
-## @deftypefn {Function File} {} ocframe_railwaybridge() 
-## Example taken from a real railwaybridge.
-##
-## @end deftypefn
-function [P,D,MemF]=ocframe_railwaybridge()
-	joints=[0:5.3636:59;zeros(4,12)]';
-	joints(1,3:4)=[1,1];
-	joints(12,3:4)=[1,1];
-
-	temp=[];
-	%parabola with height = 8 m => f(x)=(32*x)/59-(32*x^2)/3481
-	for i=2:11
-		temp=[temp;joints(i,1),(32*joints(i,1))/59-(32*joints(i,1)^2)/3481,0,0,0];
-	end
-
-	joints=[joints;temp];
-
-	%EIA
-	beam1=[200e9,135.7*(0.1)^4,3.204*(0.1)^2];
-	beam2=[200e9,80.5*(0.1)^4,2.011*(0.1)^2];
-	members=[];
-	for i=1:11
-		members=[members;i,i+1,beam1];
-	end
-	members=[members;1,13,beam1];
-	for i=13:21
-		members=[members;i,i+1,beam1];
-	end
-	members=[members;22,12,beam1];
-
-	for i=2:11
-		members=[members;i,i+11,beam2];
-	end
-
-	%own weight of beams neglected
-	
-	%some forces
-	nodeloads=[6,0,-50e3,0;7,0,-50e3,0];
-	
-	tic
-	[P,D,MemF]=SolveFrame(joints,members,nodeloads,[],[]);
-	toc
-	
-	%PlotFrame(joints,members,D,100);
-end
\ No newline at end of file

--- a/main/mechanics/ocframe/inst/private/EquivalentJointLoads.m	Fri Dec 09 12:29:53 2011 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,86 +0,0 @@
-## 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/>.
-
-function [P_F,MemFEF]=EquivalentJointLoads(joints,members,dist,point)
-	% joints: [x, y, constraints; ...] 1= fixed
-	% members [nodeN,nodeF,E,I,A; ...]
-	% dist containts distributed loads on members
-	%	[membernum,FxN,FyN,FxF,FyF,a,b,local ; ...]
-	% point contains the point loads on members
-	%   [membernum,Fx,Fy,a,local; ...]
-	% local=1 if in local coords
-
-	P_F=zeros(rows(joints)*3,1);
-	MemFEF=zeros(rows(members),6);
-	%distributed loads
-	for i=1:rows(dist)
-		membernum=dist(i,1);
-		FxN=dist(i,2);
-		FyN=dist(i,3);
-		FxF=dist(i,4);
-		FyF=dist(i,5);
-		a=dist(i,6);
-		b=dist(i,7);
-		local=dist(i,8);
-		N=members(membernum,1);
-		F=members(membernum,2);
-		if (local==0)
-			%convert FxN,... to local axes
-			temp=[FxN,FyN,0,FxF,FyF,0];
-			temp=(TransformationMatrix(joints(N,1),joints(N,2),joints(F,1),joints(F,2))')^-1*temp';
-			FxN=temp(1);
-			FyN=temp(2);
-			FxF=temp(4);
-			FyF=temp(5);
-		end
-		[F1,F2,F3,F4,F5,F6]=FixedEndForcesDist(sqrt((joints(N,1)-joints(F,1))^2+(joints(F,2)-joints(N,2))^2),FxN,FyN,FxF,FyF,a,b);
-		F_local=[F1,F2,F3,F4,F5,F6];
-		F_global=((TransformationMatrix(joints(N,1),joints(N,2),joints(F,1),joints(F,2))')*F_local')';
-		
-		c=[node2code(N),node2code(F)];
-		for j=1:6
-			P_F(c(j))+=F_global(j);
-			MemFEF(membernum,j)+=F_local(j);
-		end
-	end
-	for i=1:rows(point)
-		%[membernum,Fx,Fy,a,local; ...]
-		membernum=point(i,1);
-		Fx=point(i,2);
-		Fy=point(i,3);
-		a=point(i,4);
-		local=point(i,5);
-		N=members(membernum,1);
-		F=members(membernum,2);
-		if (local==0)
-			%convert to global axes
-			temp=[Fx,Fy,0,0,0,0];
-			temp=(TransformationMatrix(joints(N,1),joints(N,2),joints(F,1),joints(F,2))')^-1*temp';
-			Fx=temp(1);
-			Fy=temp(2);
-		end
-		%FixedEndForcesPnt(l,a,Fx,Fy)
-		[F1,F2,F3,F4,F5,F6]=FixedEndForcesPnt(sqrt((joints(N,1)-joints(F,1))^2+(joints(F,2)-joints(N,2))^2),a,Fx,Fy);
-		F_local=[F1,F2,F3,F4,F5,F6];
-		F_global=((TransformationMatrix(joints(N,1),joints(N,2),joints(F,1),joints(F,2))')*F_local')';
-		
-		c=[node2code(N),node2code(F)];
-		for j=1:6
-			P_F(c(j))+=F_global(j);
-			MemFEF(membernum,j)+=F_local(j);
-		end
-	end
-end

--- a/main/mechanics/ocframe/inst/private/FixedEndForcesDist.m	Fri Dec 09 12:29:53 2011 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,39 +0,0 @@
-## 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/>.
-
-function [F1,F2,F3,F4,F5,F6]=FixedEndForcesDist(l,FxN,FyN,FxF,FyF,a,b)
-	% distributed load in local coords.
-	% a: distance from Near
-	% b: distance from Far
-	% FxN, FyN forces on start of load in local coordinates
-	% FxF, FyF forces on end of load in local coordinates
-	
-	% formulas found by integrating the formulas for a point load with maxima maxima.sourceforge.net
-	F1=-((2*FxN+FxF)*l^2+((-b-4*a)*FxN+(b-2*a)*FxF)*l+(-b^2+a*b+2*a^2)*FxN+(-2*b^2-a*b+a^2)*FxF)/(6*l);
-	F2=-1*(((7*FyN+3*FyF)*l^4+((-3*b-13*a)*FyN+(3*b-7*a)*FyF)*l^3+((-3*b^2+4*a*b-3*a^2)*FyN+(3*b^2-4*a*b+3*a^2)*FyF)*l^2+
-((-3*b^3+a*b^2+a^2*b+17*a^3)*FyN+(-17*b^3-a*b^2-a^2*b+3*a^3)*FyF)*l+(2*b^4-2*a*b^3+2*a^2*b^2-2*a^3*b-8*a^4)*FyN+(8*b^4+2*a*b^3-2*a^2*b^2+2*a^3*b-2*a^4)*
-FyF)/(20*l^3));
-	F3=-1*(((3*FyN+2*FyF)*l^4+((3*a-2*b)*FyN+(2*b-3*a)*FyF)*l^3+((-2*b^2+a*b-27*a^2)*FyN+(2*b^2-a*b-3*a^2)*FyF)*l^2+
-((-2*b^3-a*b^2+4*a^2*b+33*a^3)*FyN+(-18*b^3+a*b^2-4*a^2*b+7*a^3)*FyF)*l+(3*b^4-3*a*b^3+3*a^2*b^2-3*a^3*b-12*a^4)*FyN+
-(12*b^4+3*a*b^3-3*a^2*b^2+3*a^3*b-3*a^4)*FyF)/(60*l^2));
-	F4=-((FxN+2*FxF)*l^2+((a-2*b)*FxN+(-4*b-a)*FxF)*l+(b^2-a*b-2*a^2)*FxN+(2*b^2+a*b-a^2)*FxF)/(6*l);
-	F5=-1*(((3*FyN+7*FyF)*l^4+((3*a-7*b)*FyN+(-13*b-3*a)*FyF)*l^3+((3*b^2-4*a*b+3*a^2)*FyN+(-3*b^2+4*a*b-3*a^2)*FyF)*l^2+
-((3*b^3-a*b^2-a^2*b-17*a^3)*FyN+(17*b^3+a*b^2+a^2*b-3*a^3)*FyF)*l+(-2*b^4+2*a*b^3-2*a^2*b^2+2*a^3*b+8*a^4)*FyN+(-8*b^4-2*a*b^3+2*a^2*b^2-2*a^3*b+2*a^4)*
-FyF)/(20*l^3));
-	F6=-1*(-((2*FyN+3*FyF)*l^4+((2*a-3*b)*FyN+(3*b-2*a)*FyF)*l^3+((-3*b^2-a*b+2*a^2)*FyN+(-27*b^2+a*b-2*a^2)*FyF)*l^2+
-((7*b^3-4*a*b^2+a^2*b-18*a^3)*FyN+(33*b^3+4*a*b^2-a^2*b-2*a^3)*FyF)*l+(-3*b^4+3*a*b^3-3*a^2*b^2+3*a^3*b+12*a^4)*FyN+
-(-12*b^4-3*a*b^3+3*a^2*b^2-3*a^3*b+3*a^4)*FyF)/(60*l^2));
-end
\ No newline at end of file

--- a/main/mechanics/ocframe/inst/private/FixedEndForcesPnt.m	Fri Dec 09 12:29:53 2011 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,33 +0,0 @@
-## 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/>.
-
-function [F1,F2,F3,F4,F5,F6]=FixedEndForcesPnt(l,a,Fx,Fy)
-	% point load in local coords.
-	% a: distance from Near
-	% l: length
-	% FxN, FyN forces on start of load in local coordinates
-	% FxF, FyF forces on end of load in local coordinates
-	
-	% b: distance from Far
-	b=l-a;
-	
-	F1=-b*Fx/l;
-	F2=-Fy*(b/l-a^2*b/l^3+a*b^2/l^3);
-	F3=-Fy*a*b^2/l^2;
-	F4=-a*Fx/l;
-	F5=-Fy*(a/l+a^2*b/l^3-a*b^2/l^3);
-	F6=Fy*a^2*b/l^2;
-end
\ No newline at end of file

--- a/main/mechanics/ocframe/inst/private/GlobalStiffnessMatrix.m	Fri Dec 09 12:29:53 2011 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,118 +0,0 @@
-## 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/>.
-
-function K=GlobalStiffnessMatrix(joints,members,memberstiffnessmatrices,membertransformationmatrices)
-	K=sparse(rows(joints)*3,rows(joints)*3);
-
-	# sparse matrix is even usefull for small frames:
-	#  example:
-	#  Attr Name           Size                     Bytes  Class
-	#  ==== ====           ====                     =====  ===== 
-	#       K_full        36x36                     10368  double
-	#       K_sparse      36x36                      2836  double
-
-	if (nargin==2)
-		#joints: [x, y, constraints; ...] 1= fixed
-		#members [nodeN,nodeF,E,I,A; ...]
-		for member=1:rows(members)
-			N=members(member,1);
-			F=members(member,2);
-			xN=joints(N,1);
-			yN=joints(N,2);
-			xF=joints(F,1);
-			yF=joints(F,2);
-			%T=TransformationMatrix(xN,yN,xF,yF);
-			%k=transpose(T)*MemberStiffnessMatrix(sqrt((xN-xF)**2+(yN-yF)**2),members(member,3),members(member,4),members(member,5))*T;
-			c=[node2code(N),node2code(F)];
-		
-			%TransformationMatrix stuff
-			L=sqrt((xN-xF)^2+(yN-yF)^2);
-			l=(xF-xN)/L;
-			m=(yF-yN)/L;
-			%MemberStiffnessMatrix stuff
-			E=members(member,3);
-			I=members(member,4);
-			A=members(member,5);
-		
-			%for i=1:rows(k)
-			%	for j=1:columns(k)
-			%		K(c(i),c(j))=K(c(i),c(j))+k(i,j);
-			%	end
-			%end
-		
-			%fill Global Stiffness Matrix
-			K(c(1),c(1))+=(l^2*A*E)/L+(12*m^2*E*I)/L^3;
-			K(c(1),c(2))+=(l*m*A*E)/L-(12*l*m*E*I)/L^3;
-			K(c(1),c(3))+=-(6*m*E*I)/L^2;
-			K(c(1),c(4))+=-(l^2*A*E)/L-(12*m^2*E*I)/L^3;
-			K(c(1),c(5))+=(12*l*m*E*I)/L^3-(l*m*A*E)/L;
-			K(c(1),c(6))+=-(6*m*E*I)/L^2;
-			K(c(2),c(1))+=(l*m*A*E)/L-(12*l*m*E*I)/L^3;
-			K(c(2),c(2))+=(m^2*A*E)/L+(12*l^2*E*I)/L^3;
-			K(c(2),c(3))+=(6*l*E*I)/L^2;
-			K(c(2),c(4))+=(12*l*m*E*I)/L^3-(l*m*A*E)/L;
-			K(c(2),c(5))+=-(m^2*A*E)/L-(12*l^2*E*I)/L^3;
-			K(c(2),c(6))+=(6*l*E*I)/L^2;
-			K(c(3),c(1))+=-(6*m*E*I)/L^2;
-			K(c(3),c(2))+=(6*l*E*I)/L^2;
-			K(c(3),c(3))+=(4*E*I)/L;
-			K(c(3),c(4))+=(6*m*E*I)/L^2;
-			K(c(3),c(5))+=-(6*l*E*I)/L^2;
-			K(c(3),c(6))+=(2*E*I)/L;
-			K(c(4),c(1))+=-(l^2*A*E)/L-(12*m^2*E*I)/L^3;
-			K(c(4),c(2))+=(12*l*m*E*I)/L^3-(l*m*A*E)/L;
-			K(c(4),c(3))+=(6*m*E*I)/L^2;
-			K(c(4),c(4))+=(l^2*A*E)/L+(12*m^2*E*I)/L^3;
-			K(c(4),c(5))+=(l*m*A*E)/L-(12*l*m*E*I)/L^3;
-			K(c(4),c(6))+=(6*m*E*I)/L^2;
-			K(c(5),c(1))+=(12*l*m*E*I)/L^3-(l*m*A*E)/L;
-			K(c(5),c(2))+=-(m^2*A*E)/L-(12*l^2*E*I)/L^3;
-			K(c(5),c(3))+=-(6*l*E*I)/L^2;
-			K(c(5),c(4))+=(l*m*A*E)/L-(12*l*m*E*I)/L^3;
-			K(c(5),c(5))+=(m^2*A*E)/L+(12*l^2*E*I)/L^3;
-			K(c(5),c(6))+=-(6*l*E*I)/L^2;
-			K(c(6),c(1))+=-(6*m*E*I)/L^2;
-			K(c(6),c(2))+=(6*l*E*I)/L^2;
-			K(c(6),c(3))+=(2*E*I)/L;
-			K(c(6),c(4))+=(6*m*E*I)/L^2;
-			K(c(6),c(5))+=-(6*l*E*I)/L^2;
-			K(c(6),c(6))+=(4*E*I)/L;
-		end
-	end
-	if (nargin==4)
-		#joints: [x, y, constraints; ...] 1= fixed
-		#members [nodeN,nodeF,E,I,A; ...]
-		for member=1:rows(members)
-			N=members(member,1);
-			F=members(member,2);
-			c=[node2code(N),node2code(F)];
-
-			xN=joints(N,1);
-			yN=joints(N,2);
-			xF=joints(F,1);
-			yF=joints(F,2);
-			%T=TransformationMatrix(xN,yN,xF,yF);
-			%k=transpose(T)*MemberStiffnessMatrix(sqrt((xN-xF)**2+(yN-yF)**2),members(member,3),members(member,4),members(member,5))*T;
-			k=transpose(membertransformationmatrices{member,1})*memberstiffnessmatrices{member,1}*membertransformationmatrices{member,1};
-		
-			for i=1:rows(k)
-				for j=1:columns(k)
-					K(c(i),c(j))=K(c(i),c(j))+k(i,j);
-				end
-			end
-		end
-	end
-end

--- a/main/mechanics/ocframe/inst/private/GlobalStiffnessMatrixModified.m	Fri Dec 09 12:29:53 2011 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,48 +0,0 @@
-## 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/>.
-
-function [ Kmod, Kcor ] = GlobalStiffnessMatrixModified(joints,members,supported)
-	%% \cite{Felippa2010} : "To apply support conditions without rearranging the equations we clear (set to zero) rows and columns
-	%% corresponding to prescribed zero displacements as well as the corresponding force components, and place
-	%% ones on the diagonal to maintain non-singularity. The resulting system is called the modified set of master
-	%% stiffness equations. "
-	
-	%% joints: [x, y, constraints; ...] 1= fixed
-	%% members [nodeN, nodeF, E, I, A; ...]
-	%% supported: vector with the supported (fixed) code numbers
-	
-	%% Kmod + Kcor = K
-	
-	%% Kcor is the correction matrix to make K
-	
-	Kmod=GlobalStiffnessMatrix(joints,members);
-	Kcor=sparse(rows(joints)*3,rows(joints)*3);
-	
-	%% Creat Kcor
-	for i=1:columns(supported)
-		Kcor(:,i)=Kmod(:,i);
-		Kcor(i,:)=Kmod(i,:);
-	end
-	%% delete rows and columns. This might be not the fastest method.
-	for i=1:columns(supported)
-		Kmod(:,i)=zeros(rows(joints)*3,1);
-		Kmod(i,:)=zeros(1,rows(joints)*3);
-		
-		%set 1 and -1 in Kmod and Kcor 
-		Kmod(i,i)+=1;
-		Kcor(i,i)-=1;
-	end
-end

--- a/main/mechanics/ocframe/inst/private/GlobalStiffnessMatrixRegrouped.m	Fri Dec 09 12:29:53 2011 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,57 +0,0 @@
-## 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/>.
-
-function [ K_ff, K_sf, K_ss, K_fs  ] = GlobalStiffnessMatrixRegrouped (joints, members,free,supported,memberstiffnessmatrices,membertransformationmatrices)
-	%% Returns the components of the global stiffness matrix regrouped as in the following equation:
-	%% { P_f }   [ K_ff   K_fs ]   { Delta_f }
-	%% {     } = [             ] . {         }
-	%% { P_s }   [ K_sf   K_ss ]   { Delta_s }
-	
-	%% joints: [x, y, constraints; ...] 1= fixed
-	%% members [nodeN, nodeF, E, I, A; ...]
-	%% free: vector with the free code numbers
-	%% supported: vector with the supported (fixed) code numbers
-	
-	%zeros(row,col)
-	
-	%number of free en supported joints
-	
-	%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-	%TODO: for large systems this function is still slow
-	%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-	if (nargin>4)
-		K=GlobalStiffnessMatrix(joints,members,memberstiffnessmatrices,membertransformationmatrices);
-	else
-		K=GlobalStiffnessMatrix(joints,members);
-	end
-	
-	%K_ff
-	K_ff=K(free,free);
-	
-	%K_sf
-	K_sf=K(supported,free);
-	
-	
-	%K_ss
-	if (nargout>=3)
-		K_ss=K(supported,supported);
-	end
-	
-	%K_fs
-	if (nargout==4)
-		K_fs=K(free,supported);
-	end
-end

--- a/main/mechanics/ocframe/inst/private/GlobalStiffnessMatrixRegrouped.m.old	Fri Dec 09 12:29:53 2011 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,98 +0,0 @@
-## 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/>.
-
-function [ K_ff, K_sf, K_ss, K_fs  ] = GlobalStiffnessMatrixRegrouped (joints, members,free,supported)
-	%% Returns the components of the global stiffness matrix regrouped as in the following equation:
-	%% { P_f }   [ K_ff   K_fs ]   { Delta_f }
-	%% {     } = [             ] . {         }
-	%% { P_s }   [ K_sf   K_ss ]   { Delta_s }
-	
-	%% joints: [x, y, constraints; ...] 1= fixed
-	%% members [nodeN, nodeF, E, I, A; ...]
-	%% free: vector with the free code numbers
-	%% supported: vector with the supported (fixed) code numbers
-	
-	%zeros(row,col)
-	
-	%number of free en supported joints
-	
-	%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-	%TODO: for large systems this function is still slow
-	% due to the indexof!
-	%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-	
-	f=columns(free);
-	s=columns(supported);
-	
-	K_ff=spalloc(f,f);
-	K_ss=spalloc(s,s);
-	K_sf=spalloc(s,f);
-	K_fs=spalloc(f,s);
-	for member=1:rows(members)
-		N=members(member,1);
-		F=members(member,2);
-		xN=joints(N,1);
-		yN=joints(N,2);
-		xF=joints(F,1);
-		yF=joints(F,2);
-		T=TransformationMatrix(xN,yN,xF,yF);
-		k=transpose(T)*MemberStiffnessMatrix(sqrt((xN-xF)**2+(yN-yF)**2),members(member,3),members(member,4),members(member,5))*T;%global coords.
-		c=[node2code(N),node2code(F)];%code numbers of the near and far nodes of the current member
-		
-		for i=1:rows(k)
-			for j=1:columns(k)
-				if (k(i,j)!=0)
-					c_row=c(i);
-					c_col=c(j);
-					
-					%K_ff
-					row=indexof(free,c_row);
-					col=indexof(free,c_col);
-					if (row!=-1 && col!=-1)
-						K_ff(row,col)=K_ff(row,col)+k(i,j);
-					end
-					
-					%K_sf
-					row=indexof(supported,c_row);
-					col=indexof(free,c_col);
-					if (row!=-1 && col!=-1)
-						K_sf(row,col)=K_sf(row,col)+k(i,j);
-					end
-					
-					
-					%K_ss
-					if (nargout>=3)
-						row=indexof(supported,c_row);
-						col=indexof(supported,c_col);
-						if (row!=-1 && col!=-1)
-							K_ss(row,col)=K_ss(row,col)+k(i,j);
-						end
-					end
-					
-					%K_fs
-					if (nargout==4)
-						row=indexof(free,c_row);
-						col=indexof(supported,c_col);
-						if (row!=-1 && col!=-1)
-							K_fs(row,col)=K_fs(row,col)+k(i,j);
-						end
-					end
-				end
-			end
-		end
-		
-	end
-end

--- a/main/mechanics/ocframe/inst/private/MemberStiffnessMatrices.m	Fri Dec 09 12:29:53 2011 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,32 +0,0 @@
-## 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/>.
-
-#MemberStiffnessMatrix(L,E,I,A)
-function k_array=MemberStiffnessMatrices(joints,members)
-
-	k_array = cell(rows(members),1); # is one column better than one row?
-	for member=1:rows(members)
-		N=members(member,1);
-		F=members(member,2);
-		xN=joints(N,1);
-		yN=joints(N,2);
-		xF=joints(F,1);
-		yF=joints(F,2);
-		
-		# calculate the member stiffness matrix in local coordinates and add to array
-		k_array{member,1}=MemberStiffnessMatrix(sqrt((xN-xF)**2+(yN-yF)**2),members(member,3),members(member,4),members(member,5));
-	end
-end

--- a/main/mechanics/ocframe/inst/private/MemberStiffnessMatrix.m	Fri Dec 09 12:29:53 2011 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,24 +0,0 @@
-## 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/>.
-
-function k=MemberStiffnessMatrix(L,E,I,A)
-	k=[A*E/L,0,0,-A*E/L,0,0;
-	0,12*E*I/L^3,6*E*I/L^2,0,-12*E*I/L^3,6*E*I/L^2;
-	0,6*E*I/L^2,4*E*I/L,0,-6*E*I/L^2,2*E*I/L;
-	-A*E/L,0,0,A*E/L,0,0;
-	0,-12*E*I/L^3,-6*E*I/L^2,0,12*E*I/L^3,-6*E*I/L^2;
-	0,6*E*I/L^2,2*E*I/L,0,-6*E*I/L^2,4*E*I/L];
-end
\ No newline at end of file

--- a/main/mechanics/ocframe/inst/private/MemberTransformationMatrices.m	Fri Dec 09 12:29:53 2011 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,30 +0,0 @@
-## 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/>.
-
-function T_array=MemberTransformationMatrices(joints,members)
-
-	T_array = cell(rows(members),1); # is one column better than one row?
-
-	for member=1:rows(members)
-		N=members(member,1);
-		F=members(member,2);
-		xN=joints(N,1);
-		yN=joints(N,2);
-		xF=joints(F,1);
-		yF=joints(F,2);
-		T_array{member,1}=TransformationMatrix(xN,yN,xF,yF);
-	end
-end

--- a/main/mechanics/ocframe/inst/private/TransformationMatrix.m	Fri Dec 09 12:29:53 2011 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,29 +0,0 @@
-## 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/>.
-
-function T=TransformationMatrix(x1,y1,x2,y2)
-	#Displacement transformation matrix = T
-	#Force transformation = T transposed
-	L=sqrt((x1-x2)^2+(y1-y2)^2);
-	l=(x2-x1)/L;
-	m=(y2-y1)/L;
-	T=[l,m,0,0,0,0;
-	-m,l,0,0,0,0;
-	0,0,1,0,0,0;
-	0,0,0,l,m,0;
-	0,0,0,-m,l,0;
-	0,0,0,0,0,1];
-end
\ No newline at end of file

--- a/main/mechanics/ocframe/inst/private/indexof.m	Fri Dec 09 12:29:53 2011 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,25 +0,0 @@
-## 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/>.
-
-function [index] = indexof(vector,number)
-	index=-1;
-	for i=1:columns(vector)
-		if (vector(i)==number)
-			index=i;
-			break
-		end
-	end
-end
\ No newline at end of file

--- a/main/mechanics/ocframe/inst/private/matrix2tex.m	Fri Dec 09 12:29:53 2011 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,33 +0,0 @@
-## 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/>.
-
-function matrix2tex(m)
-	printf("[tex]\\left(\\begin{array}{")
-	for c=1:columns(m)
-		printf("c")
-	end
-	printf("}\n")
-	for r=1:rows(m)
-		for c=1:columns(m)
-			if (c!=columns(m))
-				printf(" %3.3e &",m(r,c))
-			else
-				printf(" %3.3e \\\\\n",m(r,c))
-			end
-		end
-	end
-	printf("\\end{array}\\right)[/tex]\n\n")
-end

--- a/main/mechanics/ocframe/inst/private/node2code.m	Fri Dec 09 12:29:53 2011 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,19 +0,0 @@
-## 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/>.
-
-function c=node2code(node)
-	c=(node-1)*3+1:(node-1)*3+3;
-end
\ No newline at end of file

--- a/main/mechanics/ocframe/inst/private/test_FixedEndForcesPnt.m	Fri Dec 09 12:29:53 2011 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,36 +0,0 @@
-## Copyright (C) 2011 Johan
-## 
-## This program 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 2 of the License, or
-## (at your option) any later version.
-## 
-## This program 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 Octave; see the file COPYING.  If not, see
-## <http://www.gnu.org/licenses/>.
-
-## test_FixedEndForcesPnt
-
-## Author: Johan <johan@johan-Satellite-L30>
-## Created: 2011-10-02
-
-function [ F1,F2,F3,F4,F5,F6 ] = test_FixedEndForcesPnt (l,a,Fx,Fy)
-	joints=[0,0,1,1,1;
-		a,0,0,0,0;
-		l,0,1,1,1];
-	nodeloads=[2,Fx,Fy,0];
-	members=[1,2,1,1,1;
-		2,3,1,1,1];
-	R=SolveFrame(joints,members,nodeloads,[],[]);
-	F1 = R(1,1);
-	F2 = R(1,2);
-	F3 = R(1,3);
-	F4 = R(3,1);
-	F5 = R(3,2);
-	F6 = R(3,3);
-endfunction

--- a/main/mechanics/pre_install.m	Fri Dec 09 12:29:53 2011 +0000
+++ b/main/mechanics/pre_install.m	Fri Dec 09 12:33:35 2011 +0000
@@ -1,21 +1,27 @@
 function pre_install (desc)
-%% Prepares for installation a package that is organized in subfolders 
+%% Prepares for installation a package that is organized in subfolders
 
   %% List of subfolders
-  subfld = {"molecularDynamics","ocframe"};
+  subfld = {"molecularDynamics","ocframe", "core"};
 
   %% Create correct strings
-  subfld_ready = strcat({[pwd() filesep()]},
-                         subfld,[filesep() "*"]);
+  subfld_ready = strcat ({[pwd() filesep() "inst" filesep()]},
+                         subfld,[filesep() "src" filesep() "*"]);
 
   %% Destination folder
-  to_fld = strcat(pwd, filesep ());
+  to_fld = strcat (pwd (),[filesep() "src"]);
 
 
-  %% Copy files to package/ folder
+  %% Copy files to package/src folder
+  %% TODO handle merging of Makefiles
+  warning ("Copying subfolder src to package main dir, but multiple Makefiles are not handled")
+
+  if !exist("src","dir")
+    system(["mkdir " to_fld]);
+  end
+
   for from_fld = subfld_ready
-  %% TODO handle merging of Makefiles
-    system (["cp -Rf " from_fld{1} " " to_fld]);
+    system (["mv " from_fld{1} " " to_fld]);
     system (["rm -R " from_fld{1}(1:end-2)]);
   end