scilab/modules/scicos_blocks/macros/Sources/CURVE_c.sci

   1 //  Scicos
   2 //
   3 //  Copyright (C) INRIA - Masoud Najafi <masoud.najafi@inria.fr>
   4 //                        Serge Steer <serge.steer@inria.fr>       1993
   5 //                        Habib Jreij                              1993
   6 //
   7 // This program is free software; you can redistribute it and/or modify
   8 // it under the terms of the GNU General Public License as published by
   9 // the Free Software Foundation; either version 2 of the License, or
  10 // (at your option) any later version.
  11 //
  12 // This program is distributed in the hope that it will be useful,
  13 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15 // GNU General Public License for more details.
  16 //
  17 // You should have received a copy of the GNU General Public License
  18 // along with this program; if not, write to the Free Software
  19 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  20 //
  21 // See the file ../license.txt
  22 //
  23
  24 function [x,y,typ]=CURVE_c(job,arg1,arg2)
  25     //** 07/01/2008 : Adapted fot Scilab 5.0 by Simone Mannori
  26     x=[];
  27     y=[];
  28     typ=[];
  29
  30     select job
  31
  32     case "set" then
  33
  34         x = arg1;
  35         model = arg1.model;
  36         graphics = arg1.graphics;
  37         exprs = graphics.exprs;
  38         ok = %f;
  39         SaveExit = %f;
  40
  41         while %t do
  42             Ask_again = %f;
  43             [ok,Method,xx,yy,PeriodicOption,graf,exprs] = scicos_getvalue("Spline data",["Spline"+...
  44             " Method (0..7)";"x";"y";"Periodic signal(y/n)?";"Launch"+...
  45             " graphic window(y/n)?"],list("vec",1,"vec",-1, ...
  46             "vec",-1,"str",1,"str",1),exprs)
  47             if  ~ok then
  48                 break;
  49             end
  50
  51             if PeriodicOption=="y" | PeriodicOption=="Y" then
  52                 PO=1;
  53             else
  54                 exprs(4)="n";
  55                 PO=0;
  56             end
  57             if graf=="y" | graf=="Y" then
  58                 gui=1;
  59             else
  60                 exprs(5)="n";
  61                 gui=0;
  62             end
  63
  64             mtd=int(Method);
  65             if mtd<0 then
  66                 mtd=0
  67             end;
  68
  69             if mtd>7 then
  70                 mtd=7;
  71             end
  72
  73             METHOD = getmethod(mtd);
  74
  75             if ~Ask_again then
  76                 xx=xx(:);
  77                 yy=yy(:);
  78                 [nx,mx] = size(xx);
  79                 [ny,my]=size(yy);
  80                 if ~((nx==ny)&(mx==my)) then
  81                     messagebox("Incompatible size of [x] and [y]","modal","error");
  82                     Ask_again = %t;
  83                 end
  84             end
  85
  86             if ~Ask_again then //+++++++++++++++++++++++++++++++++++++++
  87                 xy = [xx,yy];
  88                 [xy] = cleandata(xy); // just for sorting to be able to compare data before and after poke_point(.)
  89                 N= size(xy,"r");
  90                 exprs(5)="n";// exprs.graf='n'
  91
  92                 if gui then //_______Graphic editor___________
  93                     ipar=[N;mtd;PO];
  94                     rpar=[];
  95
  96                     if (winsid() == []) then
  97                         curwin = 0;
  98                     else
  99                         curwin = max(winsid())+1; //** prepare a brand new win
 100                     end
 101                     //** see below in this file; "poke_point" is very similar to "edit_curv"
 102                     [orpar,oipar,ok] = poke_point(xy,ipar,rpar); //** HERE WE ARE +++++++++++++++++++++++++++++++++++
 103                     if ~ok then
 104                         break;
 105                     end;//  exit without save
 106
 107                     // verifying the data change
 108                     N2=oipar(1);
 109                     xy2=[orpar(1:N2),orpar(N2+1:2*N2)];
 110                     New_methhod=oipar(2);
 111                     DChange=%f;
 112                     METHOD=getmethod(New_methhod);
 113                     if or(xy(:,1)<>xy2(:,1)) then,
 114                         DChange=%t;
 115                     end
 116                     if or(xy(1:N-1,2)<>xy2(1:N2-1,2)) then,
 117                         DChange=%t;
 118                     end
 119                     if (xy(N,2)<>xy2(N2,2) & (METHOD<>"periodic")) then,
 120                         DChange=%t;
 121                     end
 122                     if DChange then
 123                         exprs(2)=strcat(sci2exp(xy2(:,1)))
 124                         exprs(3)=strcat(sci2exp(xy2(:,2)))
 125                     end
 126
 127                     exprs(1)=sci2exp(New_methhod);
 128                     if oipar(3)==1 then,
 129                         perop="y";
 130                     else,
 131                         perop="n";
 132                     end
 133                     exprs(4)=perop;
 134                     SaveExit=%t
 135                 else//_____________________No graphics__________________________
 136                     [Xdummy,Ydummy,orpar]=Do_Spline(N,mtd,xy(:,1),xy(:,2));
 137                     if (METHOD=="periodic") then // periodic spline
 138                         xy(N,2)=xy(1,2);
 139                     end
 140                     if (METHOD=="order 2" | METHOD=="not_a_knot"|METHOD=="periodic" | METHOD=="monotone"| METHOD=="fast" | METHOD=="clamped") then
 141                         orpar=[xy(:,1);xy(:,2);orpar];
 142                     else
 143                         if (METHOD=="zero order"|METHOD=="linear")
 144                             orpar=[xy(:,1);xy(:,2);]
 145                         end
 146                     end
 147                     exprs(1)=sci2exp(mtd);// pour le cas methode>7 | method<0
 148                     oipar=[N;mtd;PO]
 149                     SaveExit=%t
 150                 end //___________________________________________________________
 151             end //++++++++++++++++++++++++++++++++++++++++++++++++++++++
 152
 153             if (SaveExit) then
 154                 xp=find(orpar(1:oipar(1))>=0);
 155                 if (xp<>[]) then
 156                     model.firing=orpar(xp(1)); //first positive event
 157                 else
 158                     model.firing=-1;
 159                 end
 160                 model.rpar=orpar
 161                 model.ipar=oipar
 162                 graphics.exprs=exprs;
 163                 x.model=model
 164                 x.graphics=graphics
 165                 break
 166             end
 167         end
 168     case "define" then
 169         model=scicos_model()
 170         xx=[0, 1, 2];
 171         yy=[10, 20, -30];
 172         N=3;
 173         Method=3;
 174         PeriodicOption="y";
 175         Graf="n"
 176         model.sim=list("curve_c",4)
 177         model.in=[]
 178         model.out=1
 179         model.rpar=[xx(:);yy(:)]
 180         model.ipar=[N;Method;1]
 181         model.blocktype="c"
 182         model.dep_ut=[%f %t]
 183         model.evtin=1
 184         model.evtout=1
 185         model.firing=0
 186         exprs=[sci2exp(Method);sci2exp(xx);sci2exp(yy);PeriodicOption;Graf]
 187
 188         gr_i=[]
 189         x=standard_define([2 2],model,exprs,gr_i)
 190     end
 191 endfunction
 192
 193
 194 function [rpar,ipar,ok] = poke_point(ixy,iparin,rparin)
 195
 196     [lhs,rhs]=argn(0)
 197
 198     //** get_click is already defined in "editi_curv"
 199     //in line definition of get_click
 200     deff("[btn,xc,yc,win,Cmenu]=get_click(flag)",[
 201     "if ~or(winsid() == curwin) then   Cmenu = ''Quit'';return,end,";
 202     "if argn(2) == 1 then";
 203     "  [btn, xc, yc, win, str] = xclick(flag);";
 204     "else";
 205     "  [btn, xc, yc, win, str] = xclick();";
 206     "end;";
 207     "if btn == -1000 then";
 208     "  if win == curwin then";
 209     "    Cmenu = ''Quit'';";
 210     "  else";
 211     "    Cmenu = ''Open/Set'';";
 212     "  end,";
 213     "  return,";
 214     "end";
 215     "if btn == -2 then";
 216     "  xc = 0;yc = 0;";
 217     "  try "    // added to handle unwanted menu actions in french version
 218     "    execstr(''Cmenu='' + part(str, 9:length(str) - 1));";
 219     "    execstr(''Cmenu='' + Cmenu);";
 220     "  catch"
 221     "    Cmenu=[]"
 222     "  end "
 223     "  return,";
 224     "end";
 225     "Cmenu=[]"])
 226
 227     ok = %f
 228     if rhs==0 then
 229         ixy=[];
 230     end;
 231
 232     if size(xy,"c")<2 then
 233         gcf().info_message = " No [y] is provided";
 234         return
 235     end
 236
 237     [xy] = cleandata(ixy)
 238
 239     N = size(xy,"r");
 240
 241     if rhs<=1 then
 242         NOrder = 1;
 243         PeridicOption = 0;
 244         ipar = [N;NOrder;PeridicOption]
 245         rpar = []
 246     else
 247         if rhs==2 then
 248             NOrder = iparin(2);
 249             PeridicOption = iparin(3);
 250             ipar = iparin;
 251             rpar = [];
 252         else
 253             if rhs==3 then
 254                 NOrder = iparin(2);
 255                 PeridicOption = iparin(3);
 256                 ipar = iparin;
 257                 rpar = rparin
 258             end
 259         end //** ???
 260     end //** ???
 261
 262     Amp=[];
 263     wp=[];
 264     phase=[];
 265     offset=[];
 266     np1=[];
 267     Sin_exprs = list(string(Amp),string(wp), string(phase),string(offset),string(np1));
 268     sAmp=[];
 269     sTp=[];
 270     sdelay=[];
 271     Sawt1_exprs = list(string(sAmp),string(sTp),string(sdelay));
 272     sAmp2=[];
 273     sTp2 = [];
 274     Sawt2_exprs = list(string(sAmp2),string(sTp2));
 275
 276     Amp3=[];
 277     Tp3=[];
 278     Pw3=[];
 279     Pd3=[];
 280     Bias3=[];
 281     Pulse_exprs=list(string(Amp3), string(Tp3),string(Pw3),string(Pd3),string(Bias3))
 282
 283     mean4=[];
 284     var4=[];
 285     seed4=[];
 286     sample4=[];
 287     np4=[];
 288     random_n_exprs=list(string(mean4),string(var4), string(seed4),string(sample4),string(np4))
 289
 290     min5=[];
 291     max5=[];
 292     seed5=[];
 293     sample5=[];
 294     np5=[];
 295     random_u_exprs=list(string(min5), string(max5), string(seed5),string(sample5),string(np5))
 296
 297     // bornes initiales du graphique
 298     xmx = max(xy(:,1));
 299     xmn=min(xy(:,1)),
 300     xmn=max(xmn,0);
 301     ymx = max(xy(:,2));
 302     ymn=min(xy(:,2));
 303     dx = xmx-xmn;
 304     dy = ymx-ymn
 305
 306     if dx==0 then
 307         dx=max(xmx/2,1),
 308     end;
 309     xmx = xmx+dx/50;
 310
 311     if dy==0 then
 312         dy=max(ymx/2,1),
 313     end;
 314     ymn = ymn-dy/50;
 315     ymx = ymx+dy/50;
 316
 317     rect = [xmn,ymn;xmx,ymx];
 318
 319     // initial draw
 320     f = scf(curwin);
 321
 322
 323     menu_r = [];
 324     menu_s = [];
 325     menu_o = ["zero order","linear","order 2","not_a_knot","periodic","monotone","fast","clamped"]
 326     menu_d = ["Clear","Data Bounds","Load from text file","Save to text file","Load from Excel","Periodic signal"]
 327     menu_t=["sine","sawtooth1","sawtooth2","pulse","random normal","random uniform"]
 328     menu_e=["Help","Exit without save","Save/Exit"]
 329     MENU=["Autoscale","Spline","Data","Standards","Exit"];
 330     menus = list(MENU,menu_s,menu_o,menu_d,menu_t,menu_e);
 331
 332     scam="menus(1)(1)"
 333     w="menus(3)(";
 334     r=")";
 335     Orderm=w(ones(menu_o))+string(1:size(menu_o,"*"))+r(ones(menu_o))
 336     w="menus(4)(";
 337     r=")";
 338     Datam=w(ones(menu_d))+string(1:size(menu_d,"*"))+r(ones(menu_d))
 339     w="menus(5)(";
 340     r=")";
 341     Standm=w(ones(menu_t))+string(1:size(menu_t,"*"))+r(ones(menu_t))
 342     w="menus(6)(";
 343     r=")";
 344     Exitm=w(ones(menu_e))+string(1:size(menu_e,"*"))+r(ones(menu_e))
 345
 346     execstr("Autoscale_"+string(curwin)+"=scam")
 347     execstr("Spline_"+string(curwin)+"=Orderm")
 348     execstr("Data_"+string(curwin)+"=Datam")
 349     execstr("Standards_"+string(curwin)+"=Standm")
 350     execstr("Exit_"+string(curwin)+"=Exitm")
 351
 352     addmenu(curwin,MENU(1))
 353     addmenu(curwin,MENU(2),menu_o)
 354     addmenu(curwin,MENU(3),menu_d)
 355     addmenu(curwin,MENU(4),menu_t)
 356     addmenu(curwin,MENU(5),menu_e)
 357     //===================================================================
 358
 359     drawlater();
 360     a = gca();
 361     a.data_bounds  = rect;
 362     a.axes_visible = "on";
 363     a.clip_state   = "on";
 364     xtitle( "", "time", "Output" ) ;
 365     a.title.font_size=2;
 366     a.title.font_style=4;
 367     a.title.foreground=2;
 368
 369     a.grid=[2 2];
 370     xpolys(xy(:,1),xy(:,2),[-1]);   //children(2)
 371     xpolys(xy(:,1),xy(:,2),[5]);    //children(1)
 372     splines = a.children(1).children
 373     points  = a.children(2).children
 374     //---------------------------------------
 375     [rpar,ipar]=AutoScale(a,xy,ipar,rpar)
 376     drawnow();
 377     // -- boucle principale
 378
 379     lines(0);
 380     while %t then //=================================================
 381         N = size(xy,"r");
 382         [btn,xc,yc,win,Cmenu] = get_click(); //** see
 383         if ((win>0) & (win<>curwin)) then
 384             Cmenu="Mouse click is Offside!";
 385         end
 386         if Cmenu==[] then
 387             Cmenu="edit",
 388         end
 389         if (Cmenu=="Exit") |(Cmenu=="Quit" ) then,
 390             ipar=[];
 391             rpar=[];
 392             ok=%f;
 393             return;
 394         end
 395         //-------------------------------------------------------------------
 396         if ((Cmenu=="zero order") | (Cmenu=="linear") | (Cmenu=="order 2")| ...
 397             (Cmenu=="not_a_knot")| (Cmenu=="periodic")| (Cmenu=="monotone")| ...
 398             (Cmenu=="fast")| (Cmenu=="clamped")) then
 399
 400             select  Cmenu
 401             case "zero order" then
 402                 NOrder=0;
 403             case "linear" then
 404                 NOrder=1;
 405             case "order 2" then
 406                 NOrder=2;
 407             case "not_a_knot" then
 408                 NOrder=3;
 409             case "periodic" then
 410                 NOrder=4;
 411             case "monotone" then
 412                 NOrder=5;
 413             case "fast" then
 414                 NOrder=6;
 415             case "clamped" then
 416                 NOrder=7;
 417             end
 418             ipar(2)=NOrder;
 419             [rpar,ipar]=AutoScale(a,xy,ipar,rpar)
 420         end
 421         //-------------------------------------------------------------------
 422         select Cmenu
 423         case "Data Bounds" then
 424             rectx=findrect(a);
 425             [mok, xmn1, xmx1, ymn1, ymx1] = scicos_getvalue("Enter new bounds",["xmin";"xmax"; "ymin";"ymax"], ..
 426             list("vec", 1,"vec", 1,"vec", 1,"vec", 1), string(rectx(:)))
 427             //drawlater();
 428             if mok then
 429                 if (xmn1 > xmx1 | ymn1 > ymx1) then
 430                     gcf().info_message = "Incorrect bounds"
 431                     mok=%f;
 432                 end
 433                 if xmn1<0 then
 434                     gcf().info_message = "X should be positive"
 435                     mok=%f;
 436                 end
 437                 if mok then
 438                     a.data_bounds=[xmn1, ymn1; xmx1, ymx1];
 439                 end
 440             end
 441             //drawnow();
 442             //-------------------------------------------------------------------
 443         case "Autoscale" then
 444             [rpar,ipar]=AutoScale(a,xy,ipar,rpar)
 445             //-------------------------------------------------------------------
 446         case "Periodic signal" then
 447             if PeridicOption==1 then,
 448                 ans0="y",
 449             else,
 450                 ans0="n",
 451             end;
 452             [mok,myans]=scicos_getvalue("Generating periodic signal",["y/n"],list("str",1),list(ans0));
 453             if ((myans=="y")|(myans=="Y")) then,
 454                 PeridicOption=1,
 455             else,
 456                 PeridicOption=0;
 457             end;
 458             ipar(3)=PeridicOption;
 459             [rpar,ipar]=AutoScale(a,xy,ipar,rpar)
 460             //-------------------------------------------------------------------
 461         case "sine" then
 462             [mok,Amp,wp,phase,offset,np1,Sin_exprs2]=scicos_getvalue(" Sine parameters", ...
 463             ["Amplitude";"Frequency(rad/sec)"; ...
 464             "Phase(rad)";"Bias";"number of points"],list("vec",1,"vec",1,"vec",1, ...
 465             "vec",1,"vec",1),Sin_exprs)
 466             if np1< 2 then
 467                 np1=2;
 468             end
 469             if mok & wp>0  then
 470                 NOrder=3;
 471                 ipar(2)=NOrder;
 472                 phase=atan(tan(phase));
 473                 xt=linspace(0,%pi*2/wp,np1)';
 474                 yt=Amp*sin(wp*xt+phase)+offset;
 475                 xy=[xt,yt];
 476                 [rpar,ipar]=AutoScale(a,xy,ipar,rpar)
 477                 Sin_exprs=Sin_exprs2
 478             end
 479             //-------------------------------------------------------------------
 480         case "sawtooth1" then
 481             [mok,sAmp,sTp,sdelay,Sawt1_exprs2]=scicos_getvalue("Sawtooth signal parameters", ...
 482             ["Amplitude";"Period";"delay"], ...
 483             list("vec",1,"vec",1,"vec",1),Sawt1_exprs)
 484             if mok & sTp>0 then
 485                 NOrder=1;
 486                 ipar(2)=NOrder;
 487                 if sdelay<sTp then
 488                     xt=[0;sdelay;sTp];
 489                     yt=[0;0;sAmp];
 490                 else
 491                     xt=[0];
 492                     yt=[0];
 493                 end
 494                 xy=[xt,yt];
 495                 [rpar,ipar]=AutoScale(a,xy,ipar,rpar);
 496                 Sawt1_exprs=Sawt1_exprs2
 497             end
 498             //-------------------------------------------------------------------
 499         case "sawtooth2" then
 500             [mok,sAmp2,sTp2,Sawt2_exprs2]=scicos_getvalue("Sawtooth signal parameters", ...
 501             ["Amplitude";"Period"],list("vec",1,"vec",1),Sawt2_exprs)
 502             if mok & sTp2>0 then
 503                 NOrder=1;
 504                 ipar(2)=NOrder;
 505                 xt=[0;sTp2];
 506                 yt=[sAmp2;-sAmp2];
 507                 xy=[xt,yt];
 508                 [rpar,ipar]=AutoScale(a,xy,ipar,rpar);
 509                 Sawt2_exprs=Sawt2_exprs2
 510             end
 511             //-------------------------------------------------------------------
 512         case "pulse" then
 513             [mok,Amp3,Tp3,Pw3,Pd3,Bias3,Pulse_exprs2] = scicos_getvalue("Square wave pulse signal", ...
 514             ["Amplitude";"Period (sec)";"Pulse width(% o"+...
 515             "f period)";"Phase delay (sec)";"Bias"],list("vec",1, ...
 516             "vec",1,"vec",1,"vec",1,"vec",1),Pulse_exprs);
 517             if mok & Tp3>0  then
 518                 NOrder=0;
 519                 ipar(2)=NOrder;
 520                 if (Pd3>0) then
 521                     xt=0;
 522                     yt=Bias3;
 523                 else
 524                     xt=[];
 525                     yt=[];
 526                 end
 527                 //otherwise there would be double points at 0
 528                 if Pd3<Tp3 then
 529                     if Pw3>0 then
 530                         xt=[xt;Pd3; Pw3*Tp3/100+Pd3;Tp3];
 531                         yt=[yt;Amp3+Bias3;Bias3;Bias3];
 532                     else
 533                         xt=[0;Tp3];yt=[Bias3;Bias3];
 534                     end
 535                 else
 536                     xt=[0;Tp3];yt=[Bias3;Bias3];
 537                 end
 538
 539                 xy=[xt,yt];
 540                 [rpar,ipar]=AutoScale(a,xy,ipar,rpar);
 541                 Pulse_exprs=Pulse_exprs2;
 542             end
 543             //-------------------------------------------------------------------
 544         case "random normal" then
 545             [mok,mean4,var4,seed4,sample4,np4,random_n_exprs2]=scicos_getvalue("Normal (Gaussian) random signal", ...
 546             ["Mean";"Variance";"Initial seed";"Sample time";"Number of points"],list("vec",1, ...
 547             "vec",1,"vec",1,"vec", ...
 548             1,"vec",1),random_n_exprs)
 549             if mok & sample4>0 then
 550                 NOrder=0;
 551                 ipar(2)=NOrder;
 552                 rand("normal");
 553                 rand("seed",seed4);
 554                 xt=0:sample4:sample4*(np4-1);
 555                 xt=xt(:);
 556                 yt=mean4+sqrt(var4)*rand(np4,1);
 557                 xy=[xt,yt];
 558                 [rpar,ipar]=AutoScale(a,xy,ipar,rpar);
 559                 random_n_exprs2=random_n_exprs;
 560             end
 561             //-------------------------------------------------------------------
 562         case "random uniform" then
 563             [mok,min5,max5,seed5,sample5,np5,random_u_exprs2]=scicos_getvalue("Uniform random signal", ...
 564             ["Minimum";"Maximum";"Initial seed";"Sample time";"Number of points"],list("vec",1, ...
 565             "vec",1,"vec",1,"vec", ...
 566             1,"vec",1),random_u_exprs)
 567             if mok & sample5>0 then
 568                 NOrder=0;
 569                 ipar(2)=NOrder;
 570                 rand("uniform");
 571                 rand("seed",seed5);
 572                 xt=0:sample5:sample5*(np5-1);
 573                 xt=xt(:);
 574                 yt=min5+(max5-min5)*rand(np5,1);
 575                 xy=[xt,yt];
 576                 [rpar,ipar]=AutoScale(a,xy,ipar,rpar);
 577                 random_u_exprs2=random_u_exprs;
 578
 579             end
 580             //-------------------------------------------------------------------
 581         case "Save/Exit" then
 582             NOrder=ipar(2);
 583             PeridicOption=ipar(3);
 584
 585             METHOD=getmethod(NOrder);
 586             if (METHOD=="periodic") then // periodic spline
 587                 xy(N,2)=xy(1,2);
 588             end
 589
 590             if (METHOD=="order 2" | METHOD=="not_a_knot"|METHOD=="periodic" | METHOD=="monotone"| METHOD=="fast" | METHOD=="clamped") then
 591                 rpar=[xy(:,1);xy(:,2);rpar];
 592             else
 593                 if (METHOD=="zero order"|METHOD=="linear")
 594                     rpar=[xy(:,1);xy(:,2);]
 595                 end
 596             end
 597
 598             ok=%t
 599             delete(f);
 600             return
 601             //-------------------------------------------------------------------
 602         case "Exit without save" then
 603             ipar=[];
 604             rpar=[];
 605             ok=%f
 606             delete(f);
 607             return
 608             //-------------------------------------------------------------------
 609         case "Clear" then
 610             xy=[0,0];
 611             NOrder=0;
 612             ipar(2)=NOrder;
 613             [rpar,ipar]=AutoScale(a,xy,ipar,rpar)
 614             //----------------------------------------------------------------
 615         case "Edit text data NOT IN USE" then
 616             //  editvar xy;
 617             [mok,xt,yt]=scicos_getvalue("Enter x and y data",["x";"y"],list("vec",-1,"vec",-1),list(strcat(sci2exp(xy(:,1))),strcat(sci2exp(xy(:,2)))));
 618             if mok then,
 619                 xy=[xt,yt];
 620                 [xy]=cleandata(xy),
 621                 [rpar,ipar]=AutoScale(a,xy,ipar,rpar)
 622             end
 623             //---------------------------------------------------------------
 624         case "Help" then
 625             t1="Mouse-left click: adding a new point"
 626             t2="Mouse-right click: remove a point"
 627             t3="Mouse-left double click: edit a point''s coordinates"
 628             t4="Mouse-left button press/drag/release: move a  point"
 629             t5="Change the window size: ''Data'' menu -> ''Databounds''"
 630             messagebox([t1;t2;t3;t4;t5],"modal","info");
 631             //---------------------------------------------------------------
 632         case "Load from Excel" then
 633             [tok,xytt]=ReadExcel()
 634             if tok then
 635                 xy=xytt;
 636                 NOrder=1
 637                 ipar(2)=NOrder;
 638                 [rpar,ipar]=AutoScale(a,xy,ipar,rpar)
 639             end
 640             //---------------------------------------------------------------
 641         case "Load from text file" then
 642             [tok,xytt]=ReadFromFile()
 643             if tok then
 644                 xy=xytt;
 645                 NOrder=1
 646                 ipar(2)=NOrder;
 647                 [rpar,ipar]=AutoScale(a,xy,ipar,rpar)
 648             end
 649             //---------------------------------------------------------------
 650         case "Save to text file" then
 651             [sok]=SaveToFile(xy)
 652             //---------------------------------------------------------------
 653         case "Replot" then
 654             if xy<>[] then
 655                 drawlater();
 656                 points.data=xy;
 657                 [rpar,ipar]=drawSplin(a,xy,ipar,rpar);
 658                 drawnow()
 659             end
 660             //----------------------------------------------------------
 661         case "edit" then
 662             HIT=%f
 663             if N<>0 then
 664                 xt=xy(:,1);
 665                 yt=xy(:,2);
 666                 dist=((xt-ones(N,1)*xc))^2+((yt-ones(N,1)*yc))^2
 667                 [dca,k]=min(dist);
 668                 rectx=a.data_bounds;
 669                 ex=abs(rectx(2,1)-rectx(1,1))/80;
 670                 ey=abs(rectx(2,2)-rectx(1,2))/80;
 671                 if (abs(xc-xt(k))<ex & abs(yc-yt(k))<ey) then
 672                     HIT=%t
 673                 end
 674             end
 675
 676             //_________________________
 677             //  if ~((NOrder==-1|NOrder==-2|NOrder==-3|NOrder==-4)) then
 678             if (~HIT)&(btn==0 | btn==3) then    // add point
 679                 if (xc>=0) then
 680                     if (xc==0) then
 681                         zz=find(x==0);
 682                         xy(zz,:)=[];
 683                     end
 684                     xy=[xy;xc,yc];
 685                     [xtt,k2]=gsort(xy(:,1),"r","i");xy=xy(k2,:)
 686                     drawlater();
 687                     points.data=xy;
 688                     [rpar,ipar]=drawSplin(a,xy,ipar,rpar);
 689                     drawnow();
 690                 end
 691             end
 692
 693             if (HIT)&(btn==2 | btn==5) then  //   remove point
 694                 if (xy(k,1)>0) |( xy(k,1)==0 & (size(find(xy(:,1)==0),"*")>1)) then
 695                     xy(k,:)=[];
 696                 end
 697                 drawlater();
 698                 points.data = xy;
 699                 [rpar,ipar] = drawSplin(a,xy,ipar,rpar);
 700                 drawnow();
 701             end
 702
 703             if (HIT)&(btn==0) then             // move point
 704                 [xy,rpar,ipar] = movept(a,xy,ipar,rpar,k)
 705             end
 706
 707             if (HIT)&(btn==10) then             // change data:: double click
 708                 [mok,xt,yt]=scicos_getvalue("Enter new x and y",["x";"y"],...
 709                 list("vec",1,"vec",1),list(sci2exp(xy(k,1)),sci2exp(xy(k,2))));
 710                 if mok then
 711                     xy(k,:) = [xt,yt];
 712                     [xy] = cleandata(xy)
 713                     drawlater();
 714                     points.data=xy;
 715                     [rpar,ipar]=AutoScale(a,xy,ipar,rpar)
 716                     drawnow()
 717                 end
 718             end
 719
 720             //  end
 721             //_________________________________
 722
 723         end
 724         //----------------------------------------------------------
 725     end
 726 endfunction
 727 //========================================================================
 728 function [orpar,oipar] = drawSplin(a,xy,iipar,irpar)
 729     N=size(xy,"r");// new size of xy
 730     x=xy(:,1);
 731     y=xy(:,2);
 732     points=a.children(2).children
 733     splines=a.children(1).children
 734     order=iipar(2);
 735     periodicoption=iipar(3);
 736     orpar=irpar;
 737
 738     METHOD=getmethod(order);
 739
 740     if periodicoption==1 then
 741         PERIODIC="periodic, T="+string(x(N)-x(1));
 742     else
 743         PERIODIC="aperiodic";
 744     end
 745     a.title.text=[string(N)+" points,  "+"Method: "+METHOD+",  "+PERIODIC];
 746
 747     if (N==0) then,
 748         return;
 749     end
 750     if (N==1) then,
 751         order=0;
 752     end
 753     //  NP=50;// number of intermediate points between two data points
 754     [X,Y,orpar]=Do_Spline(N,order,x,y);
 755     if (periodicoption==1) then
 756         X=[X;X($)];
 757         Y=[Y;Y(1)];
 758     else
 759         xmx=max(points.data(:,1));  xmn=min(points.data(:,1));
 760         XMX=max(0,xmx);
 761         XMN=max(0,xmn);
 762         xmx1=max(a.x_ticks.locations)
 763         XMX=max(XMX,xmx1);
 764         X=[X;XMX];
 765         Y=[Y;Y($)];
 766     end
 767     //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 768     splines.data=[X,Y];
 769     oipar=[N;iipar(2);periodicoption]
 770 endfunction
 771 //=============================================================
 772 function [xyt,orpar,oipar]=movept(a,xy,iipar,irpar,k)
 773     //on bouge un point existant
 774     points=a.children(2).children
 775     splines=a.children(1).children
 776     oipar=iipar
 777     orpar=irpar
 778     order=iipar(2);
 779     x=xy(:,1);
 780     y=xy(:,2);
 781
 782     if (x(k)==0) then
 783         zz=find(x==0);
 784         x(zz)=[];
 785         y(zz)=[];
 786         ZERO_POINT=%t
 787     else
 788         x(k)=[];
 789         y(k)=[];
 790         ZERO_POINT=%f
 791     end
 792
 793     btn=-1
 794
 795     while ~(btn==3 | btn==0| btn==10| btn==-5)
 796         rep=xgetmouse([%t %t]); xc=rep(1);yc=rep(2);btn=rep(3);
 797         if (ZERO_POINT) then
 798             xc=0;
 799         else
 800             if (xc<=0) then
 801                 zz=find(x==0);
 802                 x(zz)=[];
 803                 y(zz)=[];
 804                 ZERO_POINT=%t;
 805                 xc=0;
 806             end
 807         end
 808
 809         xt=[x;xc];
 810         yt=[y;yc];
 811         [xt,k2]=gsort(xt,"r","i");yt=yt(k2)
 812         xyt=[xt,yt];
 813
 814         drawlater();
 815         points.data=xyt;
 816         [orpar,oipar]=drawSplin(a,xyt,oipar,orpar);
 817         drawnow()
 818     end
 819
 820 endfunction
 821
 822 //==========================================================
 823 function rectx = findrect(a)
 824     splines=a.children(1).children
 825     points=a.children(2).children
 826
 827     if (points.data==[]) then
 828         rectx=a.data_bounds;
 829         return;
 830     end
 831
 832
 833     ymx1=max(splines.data(:,2));
 834     ymn1=min(splines.data(:,2))
 835
 836     xmx=max(points.data(:,1));
 837     xmn=min(points.data(:,1));
 838     ymx=max(points.data(:,2));
 839     ymn=min(points.data(:,2));
 840
 841
 842     XMX=max(0,xmx);
 843     XMN=max(0,xmn);
 844     YMX=max(ymx,ymx1);
 845     YMN=min(ymn,ymn1);
 846
 847     dx=XMX-XMN;
 848     dy=YMX-YMN
 849     if dx==0 then
 850         dx=max(XMX/2,1),
 851     end;
 852     XMX=XMX+dx/50
 853     if dy==0 then
 854         dy=max(YMX/2,1),
 855     end;
 856     YMN=YMN-dy/50;
 857     YMX=YMX+dy/50;
 858     rectx=[XMN,YMN;XMX,YMX];
 859 endfunction
 860
 861 //============================================================
 862 function [tok,xyo]=ReadExcel()
 863     TA=["A";"B";"C";"D";"E";"F";"G";"H";"I";"J";"K";"L";"M";"N";"O";"P"; ...
 864     "Q";"R";"S";"T";"U";"V";"W";"X";"Y";"Z";"a";"b";"c";"d";"e";"f"; ...
 865     "g";"h";"i";"j";"k";"l";"m";"n";"o";"p";"q";"r";"s";"t";"u";"v"; ...
 866     "w";"x";"y";"z"];
 867     TN=["0","1","2","3","4","5","6","7","8","9"];
 868     xyo=[];
 869     tok=%f;
 870     while %t
 871         [zok,filen,sheetN,xa,ya]=scicos_getvalue("Excel data file ",["Filename";"Sheet #"+...
 872         " ";"X[start:Stop]";"Y[start:stop]"],list("str",1, ...
 873         "vec",1,"str",1, ...
 874         "str",1), ...
 875         list(["Classeur1.xls"],["1"],["C5:C25"],["D5:D25"]));
 876         if ~zok then
 877             break,
 878         end
 879
 880         try
 881             [fd,SST,Sheetnames,Sheetpos] = xls_open(filen);
 882         catch
 883             gcf().info_message = "Scicos cannot find the excel file:"+filen;
 884             break;
 885         end
 886         try
 887             N=size(Sheetnames,"*");
 888             if ((sheetN<=N) &(sheetN>0)) then
 889                 [Value,TextInd] = xls_read(fd,Sheetpos(sheetN))
 890                 mclose(fd)
 891             end
 892             xa=strsubst(xa," ","");
 893             px=strindex(xa,":");
 894             ya=strsubst(ya," ","");
 895             py=strindex(ya,":");
 896             x1=part(xa,1:px-1);
 897             x2=part(xa,px+1:length(xa));
 898             y1=part(ya,1:py-1);
 899             y2=part(ya,py+1:length(ya));
 900
 901             x1p=min(strindex(x1,TN));
 902             if x1p==[] then,
 903                 gcf().info_message = "Bad address in X:"+x1;
 904                 break,
 905             end
 906             x11=part(x1,1:x1p-1);
 907             x12=part(x1,x1p:length(x1));
 908
 909             x2p=min(strindex(x2,TN));
 910             if x2p==[] then,
 911                 gcf().info_message = "Bad address in X:"+x2;
 912             break, end
 913             x21=par
 914             t(x2,1:x2p-1);
 915             x22=part(x2,x2p:length(x2));
 916
 917             y1p=min(strindex(y1,TN));
 918             if y1p==[] then,
 919                 gcf().info_message = "Bad address in Y:"+y1;
 920                 break,
 921             end
 922             y11=part(y1,1:y1p-1);
 923             y12=part(y1,y1p:length(y1));
 924
 925             y2p=min(strindex(y2,TN));
 926             if y2p==[] then,
 927                 gcf().info_message = "Bad address in Y:"+y2;
 928                 break,
 929             end
 930             y21=part(y2,1:y2p-1);
 931             y22=part(y2,y2p:length(y2));
 932
 933             // x11 x12: x21 x22
 934
 935             lx11=length(x11);
 936             lx21=length(x21);
 937             ly11=length(y11);
 938             ly21=length(y21)
 939             xstC=0;
 940             for i=1:lx11,
 941                 xstC=xstC+modulo(find(TA==part(x11,lx11-i+1)),26)*26^(i-1);
 942             end
 943             xenC=0;
 944             for i=1:lx21,
 945                 xenC=xenC+modulo(find(TA==part(x21,lx21-i+1)),26)*26^(i-1);
 946             end
 947             ystC=0;
 948             for i=1:ly11,
 949                 ystC=ystC+modulo(find(TA==part(y11,ly11-i+1)),26)*26^(i-1);
 950             end
 951             yenC=0;
 952             for i=1:ly11,
 953                 yenC=yenC+modulo(find(TA==part(y21,ly21-i+1)),26)*26^(i-1);
 954             end
 955
 956             xstR=evstr(x12);
 957             xenR=evstr(x22);
 958             ystR=evstr(y12);
 959             yenR=evstr(y22);
 960
 961             [mv,nv]=size(Value)
 962
 963             if ~(xstR<=mv & xstR>0 & xenR<=mv & xenR>0&ystR<=mv & ystR>0&yenR<=mv&yenR>0 ) then
 964                 gcf().info_message = "error in Row data addresses";
 965                 break
 966             end
 967             if ~(xstC<=nv & xstC>0 & xenC<=nv & xenC>0&ystC<=nv & ystC>0&yenC<=nv&yenC>0 ) then
 968                 gcf().info_message = "error in Column data addresses";
 969                 break
 970             end
 971
 972             xo=Value(min(xstR,xenR):max(xstR,xenR),min(xstC,xenC):max(xstC,xenC));
 973             yo=Value(min(ystR,yenR):max(ystR,yenR),min(ystC,yenC):max(ystC,yenC));
 974             [nx,mx]=size(xo);// adjusting the x and y size
 975             [ny,my]=size(yo);
 976             N=min(nx,ny);
 977             xo=xo(1:N,:);
 978             yo=yo(1:N,:);
 979
 980             xyo=[xo,yo];
 981             [xyo]=cleandata(xyo)
 982
 983             tok=%t;
 984             break,
 985         catch
 986             gcf().info_message = " Scicos cannot read your Excel file, please verify the parameters ";
 987             break
 988         end
 989     end
 990
 991 endfunction
 992 //---------------------------------------------------------------
 993 function [xyo]=cleandata(xye)
 994     xe=xye(:,1)
 995     ye=xye(:,2)
 996
 997     [nx,mx]=size(xe);// adjusting the x and y size
 998     [ny,my]=size(ye);
 999     N=min(nx,ny);
1000     xe=xe(1:N,:);
1001     ye=ye(1:N,:);
1002
1003     // checking for NULL data
1004     for i=1:N
1005         if (xe(i)<>xe(i)) then
1006             gcf().info_message = "x contains no data:x("+string(i)+")";
1007             return;
1008         end
1009         if (ye(i)<>ye(i)) then
1010             gcf().info_message = "Y contains no data:y("+string(i)+")";
1011             return;
1012         end
1013     end
1014     zz=find(xe<0);
1015     xe(zz)=[];
1016     ye(zz)=[]
1017     if (find(xe==0)==[]) then // add zero point
1018         xe($+1)=0;
1019         ye($+1)=0;
1020     end
1021
1022     [xo,k2]=gsort(xe,"r","i");
1023     yo=ye(k2)
1024
1025     xyo=[xo,yo];
1026 endfunction
1027 //---------------------------------------------------------------
1028 function  [orpar,oipar] = AutoScale(a,xy,inipar,inrpar)
1029     drawlater();
1030     oipar = inipar
1031     orpar = inrpar
1032     points = a.children(2).children
1033     splines = a.children(1).children
1034     points.data = xy;
1035     splines.data = xy;
1036     [orpar,oipar] = drawSplin(a,xy,oipar,orpar);
1037     rectx=findrect(a);
1038     a.data_bounds = rectx;
1039     drawnow()
1040 endfunction
1041 //============================
1042 function METHOD = getmethod(order)
1043     select order
1044     case 0 then,
1045         METHOD="zero order"
1046     case 1 then,
1047         METHOD="linear"
1048     case 2 then,
1049         METHOD="order 2"
1050     case 3 then,
1051         METHOD="not_a_knot"
1052     case 4 then,
1053         METHOD="periodic"
1054     case 5 then,
1055         METHOD="monotone"
1056     case 6 then,
1057         METHOD="fast"
1058     case 7 then,
1059         METHOD="clamped"
1060     end
1061 endfunction
1062 //=======================================
1063 function [sok,xye] = ReadFromFile()
1064     xye=[];sok=%f;
1065     while %t
1066         [sok,filen,Cformat,Cx,Cy]=scicos_getvalue("Text data file ",["Filename";"Reading [C] f"+...
1067         "ormat";"Abscissa column";"Output"+...
1068         " column"],list("str",1,"str",1,"vec",1,"vec",1), ...
1069         list(["mydatafile.dat"],["%g %g"],["1"],["2"]));
1070         if ~sok then
1071             break,
1072         end
1073         px=strindex(Cformat,"%");
1074         NC=size(px,"*");
1075         if NC==[] then,
1076             gcf().info_message = "Bad format in reading data file";
1077             sok=%f;
1078             break;
1079         end
1080         Lx=[];
1081         try
1082             fd=mopen(filen,"r");
1083             Lx=mfscanf(-1,fd,Cformat);
1084             mclose(fd);
1085         catch
1086             gcf().info_message = "Scicos cannot open the data file: " + filen;
1087             break;
1088         end
1089
1090         [nD,mD] = size(Lx);
1091         if ((mD==0) | (nD==0)) then,
1092             gcf().info_message = "No data read";
1093             sok=%f;
1094             break;
1095         end
1096         if (mD<>NC) then,
1097             gcf().info_message = "Bad format";
1098             sok=%f;
1099             break;
1100         end
1101
1102         xe=Lx(:,Cx);
1103         ye=Lx(:,Cy);
1104         xye=[xe,ye];
1105         [xye]=cleandata(xye)
1106         sok=%t;
1107         break,
1108     end
1109 endfunction
1110 //=======================================
1111 function [sok]=SaveToFile(xye)
1112     xe=xye(:,1)
1113     ye=xye(:,2)
1114     sok=%f;
1115     while %t
1116         [sok,filen,Cformat]=scicos_getvalue("Text data file ",["Filename";"Writing [C] f"+...
1117         "ormat"],list("str",1,"str",1), ...
1118         list(["mydatafile.dat"],["%g %g"]));
1119         if ~sok then
1120             break,
1121         end
1122         px=strindex(Cformat,"%");
1123         NC=size(px,"*");
1124         if NC<>2 then,
1125             gcf().info_message = "Bad format in writing data file";
1126             sok=%f;
1127             break;
1128         end
1129
1130         Cformat=Cformat+"\n";
1131
1132         try
1133             fd=mopen(filen,"w");
1134             mfprintf(fd,Cformat,xe,ye);
1135             mclose(fd);
1136         catch
1137             gcf().info_message = "Scicos cannot open the data file:" + filen;
1138             break;
1139         end
1140
1141         sok=%t;
1142         break,
1143     end
1144 endfunction
1145 //=========================================================
1146 function [X,Y,orpar]=Do_Spline(N,order,x,y)
1147     X=[];
1148     Y=[];
1149     orpar=[];
1150
1151     METHOD=getmethod(order);
1152
1153     if (METHOD=="zero order") then
1154         X=x(1);
1155         Y=y(1);
1156         for i=1:N-1
1157             X=[X;x(i);x(i+1);x(i+1)];
1158             Y=[Y;y(i);y(i);y(i+1)];
1159         end
1160         return
1161     end
1162     //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1163     if (METHOD=="linear") then
1164         X=[];
1165         for i=1:N
1166             X=[X;x(i)];
1167             Y=[Y;y(i)];
1168         end
1169         return
1170     end
1171     //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1172     if (N<25) then
1173         NP=10;
1174     else
1175         if (N<50) then
1176             NP=5;
1177         else
1178             if (N<100) then
1179                 NP=2;
1180             else
1181                 if (N<200) then
1182                     NP=1;
1183                 else
1184                     NP=0;
1185                 end;
1186             end;
1187         end;
1188     end
1189     for i=1:N-1
1190         X=[X;linspace(x(i),x(i+1),NP+2)']; // pour tous sauf "linear" et "zero order"
1191     end
1192     //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1193     if (N>2) & (METHOD=="order 2") then
1194         Z=ORDER2(x,y);
1195         A=Z(1:N-1);
1196         B=Z(N:2*N-2);
1197         C=Z(2*N-1:3*N-3);
1198
1199         for j=1:size(X,"*")
1200             for i=N-1:-1:1
1201                 if X(j)>=x(i) then,
1202                     break;
1203                 end
1204             end
1205             Y(j)=A(i)*(X(j)-x(i))^2+B(i)*(X(j)-x(i))+C(i);
1206         end
1207         orpar=matrix(Z,-1,1)
1208     end
1209     //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1210     if (METHOD=="not_a_knot") then
1211         try
1212             d = splin(x, y, METHOD);
1213             Y = interp(X, x, y, d);
1214             orpar=d(:);
1215         catch
1216         if gui then
1217         gcf().info_message = "ERROR in SPLINE: "+METHOD
1218         else
1219         messagebox(lasterror(), "ERROR in SPLINE: "+METHOD);
1220         end
1221         end
1222
1223     end
1224     //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1225     if (METHOD=="periodic") then
1226         if y(1)<>y(N) then
1227             y(N)=y(1)
1228         end
1229         try
1230             d = splin(x, y,METHOD);
1231             Y = interp(X, x, y, d);
1232             orpar=d(:);
1233         catch
1234         if gui then
1235         gcf().info_message = "ERROR in SPLINE: "+METHOD
1236         else
1237         messagebox(lasterror(), "ERROR in SPLINE: "+METHOD);
1238         end
1239         end
1240     end
1241     //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1242     if (METHOD=="monotone" ) then
1243         try
1244             d = splin(x, y, METHOD);
1245             Y = interp(X, x, y, d);
1246             orpar=d(:);
1247         catch
1248         if gui then
1249         gcf().info_message = "ERROR in SPLINE: "+METHOD
1250         else
1251         messagebox(lasterror(), "ERROR in SPLINE: "+METHOD);
1252         end
1253         end
1254
1255     end
1256     //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1257     if (METHOD=="fast") then
1258         try
1259             d = splin(x, y, METHOD);
1260             Y = interp(X, x, y, d);
1261             orpar=d(:);
1262         catch
1263         if gui then
1264         gcf().info_message = "ERROR in SPLINE: "+METHOD
1265         else
1266         messagebox(lasterror(), "ERROR in SPLINE: "+METHOD);
1267         end
1268         end
1269     end
1270     //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1271     if (METHOD=="clamped") then
1272         try
1273             d = splin(x, y, METHOD,[0;0]);
1274             Y = interp(X, x, y, d);
1275             orpar=d(:);
1276         catch
1277         if gui then
1278         gcf().info_message = "ERROR in SPLINE: "+METHOD
1279         else
1280         messagebox(lasterror(), "ERROR in SPLINE: "+METHOD);
1281         end
1282         end
1283     end
1284
1285 endfunction
1286 //=================================================
1287 function [Z]=ORDER2(x,y)
1288     N=size(x,"*")-1;
1289     A=zeros(3*N-1,N*3);
1290     B=zeros(3*N-1,1);
1291     for i=1:N
1292         j=3*(i-1)+1;
1293         A(j,i+2*N)=1;
1294         B(j)=y(i);
1295         A(j+1,i)=(x(i+1)-x(i))^2;
1296         A(j+1,i+N)=x(i+1)-x(i);
1297         A(j+1,i+2*N)=1;
1298         B(j+1)=y(i+1);
1299     end
1300
1301     for i=1:N-1
1302         j=3*(i-1)+1;
1303         A(j+2,i)=2*(x(i+1)-x(i));
1304         A(j+2,i+N)=1;
1305         A(j+2,i+N+1)=-1;
1306     end
1307
1308     Q=zeros(3*N,3*N);
1309     for i=1:N
1310         Q(i,i)=4*(x(i+1)-x(i))^2
1311         Q(i,i+N)=2*(x(i+1)-x(i))
1312         Q(i+N,i)=2*(x(i+1)-x(i))
1313         Q(i+N,i+N)=1;
1314     end
1315
1316     At=[Q,A';A,zeros(3*N-1,3*N-1)]
1317     Bt=[zeros(3*N,1);B]
1318     Zt=At\Bt;
1319     Z=Zt(1:3*N,1)
1320 endfunction
1321 //===================================================