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