3f09d19df273477daf1f0a9bd23133bebe5edffa
[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., 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 //===================================================