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