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);