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