#*-------------------------------------------------------------------
* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
*
* This LIBRARY is free software; you can distribute it and/or modify
* it under the therms of the ALSOC FREE LICENSE as available at
* http://www.enq.ufrgs.br/alsoc.
*
* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
* from http://www.rps.eng.br Copyright (C) 2002-2004.
* All rights reserved.
*
* EMSO is distributed under the therms of the ALSOC LICENSE as
* available at http://www.enq.ufrgs.br/alsoc.
*
*--------------------------------------------------------------------
* The car axis problem
*
* The problem is a stiff DAE of index 3, consisting of 8
* differential and 2 algebraic equations.
*--------------------------------------------------------------------
* Author: Rafael de Pelegrini Soares
* $Id: sample_car.mso 249 2007-04-24 21:08:06Z rafael $
*--------------------------------------------------------------------*#
 
FlowSheet CarAxis

	PARAMETERS
	l as Real(Default=1);
	l0 as Real(Default=1/2);
	eps as Real(Default=1e-2);
	M as Real(Default=10);
	h as Real(Default=1/5);
	tau as Real(Default=3.1514/5);
	w as Real(Default=10);
	r as Real(Default=0.1);


	VARIABLES
	xl as Real(Default = 0, Lower=-1, Upper = 1);
	yl as Real(Default = 0.5, Lower=-1, Upper = 1);
	xr as Real(Default = 1, Lower=-2, Upper = 2);
	yr as Real(Default = 0.5, Lower=-1, Upper = 1);
	q(4) as Real(Default = -1);
	lambda(2) as Real(Default = 0);
	xb as Real(Default = 1);
	yb as Real(Default = 0);
	ll as Real(Default = 0.5, Brief="Left spring length");
	lr as Real(Default = 0.5, Brief="Right spring length");

	EQUATIONS
	diff([xl, yl, xr, yr])*'s' = q;
	
	eps^2*M/2 * diff(q(1))*'s' = (l0-ll)*xl/ll
		+lambda(1)*xb +2*lambda(1)*(xl-xr);
	
	eps^2*M/2 * diff(q(2))*'s' = (l0-ll)*yl/ll
		+lambda(1)*yb +2*lambda(2)*(yr-yl) - eps^2*M/2;

	
	eps^2*M/2 * diff(q(3))*'s' = (l0-lr)*(xr-xb)/lr
		-2*lambda(2)*(xl-xr);
	
	eps^2*M/2 * diff(q(4))*'s' = (l0-lr)*(yr-yb)/lr
		-2*lambda(2)*(yl-yr) - eps^2*M/2;
	
	xl*xb + yl*yb = 0;
	
	(xl-xr)^2 + (yl-yr)^2 = l^2;
	
	xb = sqrt(l^2-yb^2);
	
	yb = r*sin(w*time*'rad/s');
	
	ll = sqrt(xl^2 + yl^2);
	
	lr = sqrt((xr-xb)^2 + (yr-yb)^2);

	INITIAL
	yl = 0.5;
	xr = 1;
	q(3) = -0.5;
	q(4) = 0;
	
	OPTIONS
	TimeStart = 0;
	TimeStep = 0.01;
	TimeEnd = 3;
	Integration = "index0";# "index0"; # "original";
	NLASolver(
		File = "sundials",
		RelativeAccuracy = 1e-5,
		AbsoluteAccuracy = 1e-5
	);
	DAESolver(
		File = "dassl", # "mebdf",
		RelativeAccuracy = 1e-5,
		AbsoluteAccuracy = 1e-5
	);
end