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