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