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