f2444410d0a58c0fe18d72f88320f168980df1ea
[scilab.git] / scilab / modules / differential_equations / demos / dae / dae2 / pendg.sci
1 //
2 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 // Copyright (C) INRIA
4 //
5 // This file is distributed under the same license as the Scilab package.
6 //
7
8 funcprot(0);
9
10 function [res,ires]=pendg(t,y,ydot)
11     x          = y(1:3);
12     u          = y(4:6);
13     lambda     = y(7);
14     xp         = ydot(1:3);
15     up         = ydot(4:6);
16     res(1:3)   = xp-u;
17     res(4)     = (M+m)*up(1)+m*l*cos(x(3))*up(3)-m*l*sin(x(3))*u(3)^2+lambda*fx(x(1),x(2))+k*u(1);
18     res(5)     = (M+m)*up(2)+m*l*sin(x(3))*up(3)+m*l*cos(x(3))*u(3)^2+(M+m)*g+lambda*fy(x(1),x(2))+k*u(2);
19     res(6)     = m*l*cos(x(3))*up(1)+m*l*sin(x(3))*up(2)+m*l^2*up(3)+m*g*sin(x(3));
20     res(7)     = -(fx(x(1),x(2))*u(1)+fy(x(1),x(2))*u(2));
21     ires       = 0;
22 endfunction
23
24
25
26 function  H=build_sliding_pendulum ()
27
28     //build the sliding pendulum figure and graphic objects,
29     //return the handle on moving parts
30
31     f = scf(100001);
32     clf(f,"reset");
33
34     f             = gcf();
35     f.pixmap      = "on";
36     a             = gca();
37     drawlater();
38     f.axes_size = [610,676] ;
39
40     xmin = -1.5;
41     xmax =  1.5;
42     ymin = -1.1;
43     ymax =  2.25;
44     a.data_bounds=[xmin ymin;xmax ymax]
45     a.isoview = "on";
46
47     //the framework
48     xrect(xmin,ymax,xmax-xmin,ymax-ymin)
49
50     //the curve
51     vx=[xmin:0.1:xmax]';
52     vy=vx^2;
53     xpoly(vx,vy,"lines")
54     c=gce();
55     c.foreground=5;
56
57     //the pendulum segment
58     x           = 0;
59     y           = 0;
60     teta        = 0;
61     xp          = x+l*sin(teta); yp=y-l*cos(teta);
62     r           = 0.05 // the bullet half diameter
63     xp1         = x+(l-r)*sin(teta);
64     yp1         = y-(l-r)*cos(teta);
65     xpoly([x;xp1],[y;yp1],"lines")
66     p           = gce();
67     p.thickness = 2;
68     //the pendulum bullet
69     xfarc(xp-r,yp+r,2*r,2*r,0,360*64)
70     b           = gce()
71     H           = glue([p,b]) //return the handle on segment and bullet
72 endfunction
73
74
75 function  draw_sliding_pendulum (H,state)
76     //draw the pendulum in a given state
77     x=state(1); y=state(2); teta=state(3);
78     // bullet half diameter
79     b = H.children(1);r=b.data(3)/2
80
81     xp=x+l*sin(teta); yp=y-l*cos(teta);
82     xp1=x+(l-r)*sin(teta); yp1=y-(l-r)*cos(teta);
83     drawlater()
84     p = H.children(2);p.data=[x, y; xp1, yp1];
85     b = H.children(1); b.data=[xp-r,yp+r,2*r,2*r,0,360*64];
86     drawnow();
87     a=gca();
88     a.title.text="sliding pendulum, parabola";
89     a.title.font_size=3;
90 endfunction
91
92 funcprot(1);