// Scicos
//
// Copyright (C) INRIA - Masoud Najafi <masoud.najafi@inria.fr>
// Serge Steer <serge.steer@inria.fr> 1993
// Habib Jreij 1993
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program; if not, write to the Free Software
// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
//
// See the file ../license.txt
//
function [x,y,typ]=CURVE_c(job,arg1,arg2)
//** 07/01/2008 : Adapted fot Scilab 5.0 by Simone Mannori
x=[];
y=[];
typ=[];
select job
case "set" then
x = arg1;
model = arg1.model;
graphics = arg1.graphics;
exprs = graphics.exprs;
ok = %f;
SaveExit = %f;
while %t do
Ask_again = %f;
[ok,Method,xx,yy,PeriodicOption,graf,exprs] = scicos_getvalue("Spline data",["Spline"+...
" Method (0..7)";"x";"y";"Periodic signal(y/n)?";"Launch"+...
" graphic window(y/n)?"],list("vec",1,"vec",-1, ...
"vec",-1,"str",1,"str",1),exprs)
if ~ok then
break;
end
if PeriodicOption=="y" | PeriodicOption=="Y" then
PO=1;
else
exprs(4)="n";
PO=0;
end
mtd=int(Method);
if mtd<0 then
mtd=0
end;
if mtd>7 then
mtd=7;
end
METHOD = getmethod(mtd);
if ~Ask_again then
xx=xx(:);
yy=yy(:);
[nx,mx] = size(xx);
[ny,my]=size(yy);
if ~((nx==ny)&(mx==my)) then
messagebox("Incompatible size of [x] and [y]","modal","error");
Ask_again = %t;
end
end
if ~Ask_again then //+++++++++++++++++++++++++++++++++++++++
xy = [xx,yy];
[xy] = cleandata(xy); // just for sorting to be able to compare data before and after poke_point(.)
N= size(xy,"r");
exprs(5)="n";// exprs.graf='n'
if graf=="y" | graf=="Y" then //_______Graphic editor___________
ipar=[N;mtd;PO];
rpar=[];
if (winsid() == []) then
curwin = 0;
else
curwin = max(winsid())+1; //** prepare a brand new win
end
//** see below in this file; "poke_point" is very similar to "edit_curv"
[orpar,oipar,ok] = poke_point(xy,ipar,rpar); //** HERE WE ARE +++++++++++++++++++++++++++++++++++
if ~ok then
break;
end;// exit without save
// verifying the data change
N2=oipar(1);
xy2=[orpar(1:N2),orpar(N2+1:2*N2)];
New_methhod=oipar(2);
DChange=%f;
METHOD=getmethod(New_methhod);
if or(xy(:,1)<>xy2(:,1)) then,
DChange=%t;
end
if or(xy(1:N-1,2)<>xy2(1:N2-1,2)) then,
DChange=%t;
end
if (xy(N,2)<>xy2(N2,2) & (METHOD<>"periodic")) then,
DChange=%t;
end
if DChange then
exprs(2)=strcat(sci2exp(xy2(:,1)))
exprs(3)=strcat(sci2exp(xy2(:,2)))
end
exprs(1)=sci2exp(New_methhod);
if oipar(3)==1 then,
perop="y";
else,
perop="n";
end
exprs(4)=perop;
SaveExit=%t
else//_____________________No graphics__________________________
[Xdummy,Ydummy,orpar]=Do_Spline(N,mtd,xy(:,1),xy(:,2));
if (METHOD=="periodic") then // periodic spline
xy(N,2)=xy(1,2);
end
if (METHOD=="order 2" | METHOD=="not_a_knot"|METHOD=="periodic" | METHOD=="monotone"| METHOD=="fast" | METHOD=="clamped") then
orpar=[xy(:,1);xy(:,2);orpar];
else
if (METHOD=="zero order"|METHOD=="linear")
orpar=[xy(:,1);xy(:,2);]
end
end
exprs(1)=sci2exp(mtd);// pour le cas methode>7 | method<0
oipar=[N;mtd;PO]
SaveExit=%t
end //___________________________________________________________
end //++++++++++++++++++++++++++++++++++++++++++++++++++++++
if (SaveExit) then
xp=find(orpar(1:oipar(1))>=0);
if (xp<>[]) then
model.firing=orpar(xp(1)); //first positive event
else
model.firing=-1;
end
model.rpar=orpar
model.ipar=oipar
graphics.exprs=exprs;
x.model=model
x.graphics=graphics
break
end
end
case "define" then
model=scicos_model()
xx=[0, 1, 2];
yy=[10, 20, -30];
N=3;
Method=3;
PeriodicOption="y";
Graf="n"
model.sim=list("curve_c",4)
model.in=[]
model.out=1
model.rpar=[xx(:);yy(:)]
model.ipar=[N;Method;1]
model.blocktype="c"
model.dep_ut=[%f %t]
model.evtin=1
model.evtout=1
model.firing=0
exprs=[sci2exp(Method);sci2exp(xx);sci2exp(yy);PeriodicOption;Graf]
gr_i=[]
x=standard_define([2 2],model,exprs,gr_i)
end
endfunction
function [rpar,ipar,ok] = poke_point(ixy,iparin,rparin)
[lhs,rhs]=argn(0)
//** get_click is already defined in "editi_curv"
//in line definition of get_click
deff("[btn,xc,yc,win,Cmenu]=get_click(flag)",[
"if ~or(winsid() == curwin) then Cmenu = ''Quit'';return,end,";
"if argn(2) == 1 then";
" [btn, xc, yc, win, str] = xclick(flag);";
"else";
" [btn, xc, yc, win, str] = xclick();";
"end;";
"if btn == -1000 then";
" if win == curwin then";
" Cmenu = ''Quit'';";
" else";
" Cmenu = ''Open/Set'';";
" end,";
" return,";
"end";
"if btn == -2 then";
" xc = 0;yc = 0;";
" try " // added to handle unwanted menu actions in french version
" execstr(''Cmenu='' + part(str, 9:length(str) - 1));";
" execstr(''Cmenu='' + Cmenu);";
" catch"
" Cmenu=[]"
" end "
" return,";
"end";
"Cmenu=[]"])
ok = %f
if rhs==0 then
ixy=[];
end;
if size(xy,"c")<2 then
xinfo(" No [y] is provided");
return
end
[xy] = cleandata(ixy)
N = size(xy,"r");
if rhs<=1 then
NOrder = 1;
PeridicOption = 0;
ipar = [N;NOrder;PeridicOption]
rpar = []
else
if rhs==2 then
NOrder = iparin(2);
PeridicOption = iparin(3);
ipar = iparin;
rpar = [];
else
if rhs==3 then
NOrder = iparin(2);
PeridicOption = iparin(3);
ipar = iparin;
rpar = rparin
end
end //** ???
end //** ???
Amp=[];
wp=[];
phase=[];
offset=[];
np1=[];
Sin_exprs = list(string(Amp),string(wp), string(phase),string(offset),string(np1));
sAmp=[];
sTp=[];
sdelay=[];
Sawt1_exprs = list(string(sAmp),string(sTp),string(sdelay));
sAmp2=[];
sTp2 = [];
Sawt2_exprs = list(string(sAmp2),string(sTp2));
Amp3=[];
Tp3=[];
Pw3=[];
Pd3=[];
Bias3=[];
Pulse_exprs=list(string(Amp3), string(Tp3),string(Pw3),string(Pd3),string(Bias3))
mean4=[];
var4=[];
seed4=[];
sample4=[];
np4=[];
random_n_exprs=list(string(mean4),string(var4), string(seed4),string(sample4),string(np4))
min5=[];
max5=[];
seed5=[];
sample5=[];
np5=[];
random_u_exprs=list(string(min5), string(max5), string(seed5),string(sample5),string(np5))
// bornes initiales du graphique
xmx = max(xy(:,1));
xmn=min(xy(:,1)),
xmn=max(xmn,0);
ymx = max(xy(:,2));
ymn=min(xy(:,2));
dx = xmx-xmn;
dy = ymx-ymn
if dx==0 then
dx=max(xmx/2,1),
end;
xmx = xmx+dx/50;
if dy==0 then
dy=max(ymx/2,1),
end;
ymn = ymn-dy/50;
ymx = ymx+dy/50;
rect = [xmn,ymn;xmx,ymx];
// initial draw
f = scf(curwin);
menu_r = [];
menu_s = [];
menu_o = ["zero order","linear","order 2","not_a_knot","periodic","monotone","fast","clamped"]
menu_d = ["Clear","Data Bounds","Load from text file","Save to text file","Load from Excel","Periodic signal"]
menu_t=["sine","sawtooth1","sawtooth2","pulse","random normal","random uniform"]
menu_e=["Help","Exit without save","Save/Exit"]
MENU=["Autoscale","Spline","Data","Standards","Exit"];
menus = list(MENU,menu_s,menu_o,menu_d,menu_t,menu_e);
scam="menus(1)(1)"
w="menus(3)(";
r=")";
Orderm=w(ones(menu_o))+string(1:size(menu_o,"*"))+r(ones(menu_o))
w="menus(4)(";
r=")";
Datam=w(ones(menu_d))+string(1:size(menu_d,"*"))+r(ones(menu_d))
w="menus(5)(";
r=")";
Standm=w(ones(menu_t))+string(1:size(menu_t,"*"))+r(ones(menu_t))
w="menus(6)(";
r=")";
Exitm=w(ones(menu_e))+string(1:size(menu_e,"*"))+r(ones(menu_e))
execstr("Autoscale_"+string(curwin)+"=scam")
execstr("Spline_"+string(curwin)+"=Orderm")
execstr("Data_"+string(curwin)+"=Datam")
execstr("Standards_"+string(curwin)+"=Standm")
execstr("Exit_"+string(curwin)+"=Exitm")
addmenu(curwin,MENU(1))
addmenu(curwin,MENU(2),menu_o)
addmenu(curwin,MENU(3),menu_d)
addmenu(curwin,MENU(4),menu_t)
addmenu(curwin,MENU(5),menu_e)
//===================================================================
drawlater();
a = gca();
a.data_bounds = rect;
a.axes_visible = "on";
a.clip_state = "on";
xtitle( "", "time", "Output" ) ;
a.title.font_size=2;
a.title.font_style=4;
a.title.foreground=2;
a.grid=[2 2];
xpolys(xy(:,1),xy(:,2),[-1]); //children(2)
xpolys(xy(:,1),xy(:,2),[5]); //children(1)
splines = a.children(1).children
points = a.children(2).children
//---------------------------------------
[rpar,ipar]=AutoScale(a,xy,ipar,rpar)
drawnow();
// -- boucle principale
lines(0);
while %t then //=================================================
N = size(xy,"r");
[btn,xc,yc,win,Cmenu] = get_click(); //** see
if ((win>0) & (win<>curwin)) then
Cmenu="Mouse click is Offside!";
end
if Cmenu==[] then
Cmenu="edit",
end
if (Cmenu=="Exit") |(Cmenu=="Quit" ) then,
ipar=[];
rpar=[];
ok=%f;
return;
end
//-------------------------------------------------------------------
if ((Cmenu=="zero order") | (Cmenu=="linear") | (Cmenu=="order 2")| ...
(Cmenu=="not_a_knot")| (Cmenu=="periodic")| (Cmenu=="monotone")| ...
(Cmenu=="fast")| (Cmenu=="clamped")) then
select Cmenu
case "zero order" then
NOrder=0;
case "linear" then
NOrder=1;
case "order 2" then
NOrder=2;
case "not_a_knot" then
NOrder=3;
case "periodic" then
NOrder=4;
case "monotone" then
NOrder=5;
case "fast" then
NOrder=6;
case "clamped" then
NOrder=7;
end
ipar(2)=NOrder;
[rpar,ipar]=AutoScale(a,xy,ipar,rpar)
end
//-------------------------------------------------------------------
select Cmenu
case "Data Bounds" then
rectx=findrect(a);
[mok, xmn1, xmx1, ymn1, ymx1] = scicos_getvalue("Enter new bounds",["xmin";"xmax"; "ymin";"ymax"], ..
list("vec", 1,"vec", 1,"vec", 1,"vec", 1), string(rectx(:)))
//drawlater();
if mok then
if (xmn1 > xmx1 | ymn1 > ymx1) then
xinfo("Incorrect bounds")
mok=%f;
end
if xmn1<0 then
xinfo("X should be positive")
mok=%f;
end
if mok then
a.data_bounds=[xmn1, ymn1; xmx1, ymx1];
end
end
//drawnow();
//-------------------------------------------------------------------
case "Autoscale" then
[rpar,ipar]=AutoScale(a,xy,ipar,rpar)
//-------------------------------------------------------------------
case "Periodic signal" then
if PeridicOption==1 then,
ans0="y",
else,
ans0="n",
end;
[mok,myans]=scicos_getvalue("Generating periodic signal",["y/n"],list("str",1),list(ans0));
if ((myans=="y")|(myans=="Y")) then,
PeridicOption=1,
else,
PeridicOption=0;
end;
ipar(3)=PeridicOption;
[rpar,ipar]=AutoScale(a,xy,ipar,rpar)
//-------------------------------------------------------------------
case "sine" then
[mok,Amp,wp,phase,offset,np1,Sin_exprs2]=scicos_getvalue(" Sine parameters", ...
["Amplitude";"Frequency(rad/sec)"; ...
"Phase(rad)";"Bias";"number of points"],list("vec",1,"vec",1,"vec",1, ...
"vec",1,"vec",1),Sin_exprs)
if np1< 2 then
np1=2;
end
if mok & wp>0 then
NOrder=3;
ipar(2)=NOrder;
phase=atan(tan(phase));
xt=linspace(0,%pi*2/wp,np1)';
yt=Amp*sin(wp*xt+phase)+offset;
xy=[xt,yt];
[rpar,ipar]=AutoScale(a,xy,ipar,rpar)
Sin_exprs=Sin_exprs2
end
//-------------------------------------------------------------------
case "sawtooth1" then
[mok,sAmp,sTp,sdelay,Sawt1_exprs2]=scicos_getvalue("Sawtooth signal parameters", ...
["Amplitude";"Period";"delay"], ...
list("vec",1,"vec",1,"vec",1),Sawt1_exprs)
if mok & sTp>0 then
NOrder=1;
ipar(2)=NOrder;
if sdelay<sTp then
xt=[0;sdelay;sTp];
yt=[0;0;sAmp];
else
xt=[0];
yt=[0];
end
xy=[xt,yt];
[rpar,ipar]=AutoScale(a,xy,ipar,rpar);
Sawt1_exprs=Sawt1_exprs2
end
//-------------------------------------------------------------------
case "sawtooth2" then
[mok,sAmp2,sTp2,Sawt2_exprs2]=scicos_getvalue("Sawtooth signal parameters", ...
["Amplitude";"Period"],list("vec",1,"vec",1),Sawt2_exprs)
if mok & sTp2>0 then
NOrder=1;
ipar(2)=NOrder;
xt=[0;sTp2];
yt=[sAmp2;-sAmp2];
xy=[xt,yt];
[rpar,ipar]=AutoScale(a,xy,ipar,rpar);
Sawt2_exprs=Sawt2_exprs2
end
//-------------------------------------------------------------------
case "pulse" then
[mok,Amp3,Tp3,Pw3,Pd3,Bias3,Pulse_exprs2] = scicos_getvalue("Square wave pulse signal", ...
["Amplitude";"Period (sec)";"Pulse width(% o"+...
"f period)";"Phase delay (sec)";"Bias"],list("vec",1, ...
"vec",1,"vec",1,"vec",1,"vec",1),Pulse_exprs);
if mok & Tp3>0 then
NOrder=0;
ipar(2)=NOrder;
if (Pd3>0) then
xt=0;
yt=Bias3;
else
xt=[];
yt=[];
end
//otherwise there would be double points at 0
if Pd3<Tp3 then
if Pw3>0 then
xt=[xt;Pd3; Pw3*Tp3/100+Pd3;Tp3];
yt=[yt;Amp3+Bias3;Bias3;Bias3];
else
xt=[0;Tp3];yt=[Bias3;Bias3];
end
else
xt=[0;Tp3];yt=[Bias3;Bias3];
end
xy=[xt,yt];
[rpar,ipar]=AutoScale(a,xy,ipar,rpar);
Pulse_exprs=Pulse_exprs2;
end
//-------------------------------------------------------------------
case "random normal" then
[mok,mean4,var4,seed4,sample4,np4,random_n_exprs2]=scicos_getvalue("Normal (Gaussian) random signal", ...
["Mean";"Variance";"Initial seed";"Sample time";"Number of points"],list("vec",1, ...
"vec",1,"vec",1,"vec", ...
1,"vec",1),random_n_exprs)
if mok & sample4>0 then
NOrder=0;
ipar(2)=NOrder;
rand("normal");
rand("seed",seed4);
xt=0:sample4:sample4*(np4-1);
xt=xt(:);
yt=mean4+sqrt(var4)*rand(np4,1);
xy=[xt,yt];
[rpar,ipar]=AutoScale(a,xy,ipar,rpar);
random_n_exprs2=random_n_exprs;
end
//-------------------------------------------------------------------
case "random uniform" then
[mok,min5,max5,seed5,sample5,np5,random_u_exprs2]=scicos_getvalue("Uniform random signal", ...
["Minimum";"Maximum";"Initial seed";"Sample time";"Number of points"],list("vec",1, ...
"vec",1,"vec",1,"vec", ...
1,"vec",1),random_u_exprs)
if mok & sample5>0 then
NOrder=0;
ipar(2)=NOrder;
rand("uniform");
rand("seed",seed5);
xt=0:sample5:sample5*(np5-1);
xt=xt(:);
yt=min5+(max5-min5)*rand(np5,1);
xy=[xt,yt];
[rpar,ipar]=AutoScale(a,xy,ipar,rpar);
random_u_exprs2=random_u_exprs;
end
//-------------------------------------------------------------------
case "Save/Exit" then
NOrder=ipar(2);
PeridicOption=ipar(3);
METHOD=getmethod(NOrder);
if (METHOD=="periodic") then // periodic spline
xy(N,2)=xy(1,2);
end
if (METHOD=="order 2" | METHOD=="not_a_knot"|METHOD=="periodic" | METHOD=="monotone"| METHOD=="fast" | METHOD=="clamped") then
rpar=[xy(:,1);xy(:,2);rpar];
else
if (METHOD=="zero order"|METHOD=="linear")
rpar=[xy(:,1);xy(:,2);]
end
end
ok=%t
delete(f);
return
//-------------------------------------------------------------------
case "Exit without save" then
ipar=[];
rpar=[];
ok=%f
delete(f);
return
//-------------------------------------------------------------------
case "Clear" then
xy=[0,0];
NOrder=0;
ipar(2)=NOrder;
[rpar,ipar]=AutoScale(a,xy,ipar,rpar)
//----------------------------------------------------------------
case "Edit text data NOT IN USE" then
// editvar xy;
[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)))));
if mok then,
xy=[xt,yt];
[xy]=cleandata(xy),
[rpar,ipar]=AutoScale(a,xy,ipar,rpar)
end
//---------------------------------------------------------------
case "Help" then
t1="Mouse-left click: adding a new point"
t2="Mouse-right click: remove a point"
t3="Mouse-left double click: edit a point''s coordinates"
t4="Mouse-left button press/drag/release: move a point"
t5="Change the window size: ''Data'' menu -> ''Databounds''"
messagebox([t1;t2;t3;t4;t5],"modal","info");
//---------------------------------------------------------------
case "Load from Excel" then
[tok,xytt]=ReadExcel()
if tok then
xy=xytt;
NOrder=1
ipar(2)=NOrder;
[rpar,ipar]=AutoScale(a,xy,ipar,rpar)
end
//---------------------------------------------------------------
case "Load from text file" then
[tok,xytt]=ReadFromFile()
if tok then
xy=xytt;
NOrder=1
ipar(2)=NOrder;
[rpar,ipar]=AutoScale(a,xy,ipar,rpar)
end
//---------------------------------------------------------------
case "Save to text file" then
[sok]=SaveToFile(xy)
//---------------------------------------------------------------
case "Replot" then
if xy<>[] then
drawlater();
points.data=xy;
[rpar,ipar]=drawSplin(a,xy,ipar,rpar);
drawnow()
end
//----------------------------------------------------------
case "edit" then
HIT=%f
if N<>0 then
xt=xy(:,1);
yt=xy(:,2);
dist=((xt-ones(N,1)*xc))^2+((yt-ones(N,1)*yc))^2
[dca,k]=min(dist);
rectx=a.data_bounds;
ex=abs(rectx(2,1)-rectx(1,1))/80;
ey=abs(rectx(2,2)-rectx(1,2))/80;
if (abs(xc-xt(k))<ex & abs(yc-yt(k))<ey) then
HIT=%t
end
end
//_________________________
// if ~((NOrder==-1|NOrder==-2|NOrder==-3|NOrder==-4)) then
if (~HIT)&(btn==0 | btn==3) then // add point
if (xc>=0) then
if (xc==0) then
zz=find(x==0);
xy(zz,:)=[];
end
xy=[xy;xc,yc];
[xtt,k2]=gsort(xy(:,1),"r","i");xy=xy(k2,:)
drawlater();
points.data=xy;
[rpar,ipar]=drawSplin(a,xy,ipar,rpar);
drawnow();
end
end
if (HIT)&(btn==2 | btn==5) then // remove point
if (xy(k,1)>0) |( xy(k,1)==0 & (size(find(xy(:,1)==0),"*")>1)) then
xy(k,:)=[];
end
drawlater();
points.data = xy;
[rpar,ipar] = drawSplin(a,xy,ipar,rpar);
drawnow();
end
if (HIT)&(btn==0) then // move point
[xy,rpar,ipar] = movept(a,xy,ipar,rpar,k)
end
if (HIT)&(btn==10) then // change data:: double click
[mok,xt,yt]=scicos_getvalue("Enter new x and y",["x";"y"],...
list("vec",1,"vec",1),list(sci2exp(xy(k,1)),sci2exp(xy(k,2))));
if mok then
xy(k,:) = [xt,yt];
[xy] = cleandata(xy)
drawlater();
points.data=xy;
[rpar,ipar]=AutoScale(a,xy,ipar,rpar)
drawnow()
end
end
// end
//_________________________________
end
//----------------------------------------------------------
end
endfunction
//========================================================================
function [orpar,oipar] = drawSplin(a,xy,iipar,irpar)
N=size(xy,"r");// new size of xy
x=xy(:,1);
y=xy(:,2);
points=a.children(2).children
splines=a.children(1).children
order=iipar(2);
periodicoption=iipar(3);
orpar=irpar;
METHOD=getmethod(order);
if periodicoption==1 then
PERIODIC="periodic, T="+string(x(N)-x(1));
else
PERIODIC="aperiodic";
end
a.title.text=[string(N)+" points, "+"Method: "+METHOD+", "+PERIODIC];
if (N==0) then,
return;
end
if (N==1) then,
order=0;
end
// NP=50;// number of intermediate points between two data points
[X,Y,orpar]=Do_Spline(N,order,x,y);
if (periodicoption==1) then
X=[X;X($)];
Y=[Y;Y(1)];
else
xmx=max(points.data(:,1)); xmn=min(points.data(:,1));
XMX=max(0,xmx);
XMN=max(0,xmn);
xmx1=max(a.x_ticks.locations)
XMX=max(XMX,xmx1);
X=[X;XMX];
Y=[Y;Y($)];
end
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
splines.data=[X,Y];
oipar=[N;iipar(2);periodicoption]
endfunction
//=============================================================
function [xyt,orpar,oipar]=movept(a,xy,iipar,irpar,k)
//on bouge un point existant
points=a.children(2).children
splines=a.children(1).children
oipar=iipar
orpar=irpar
order=iipar(2);
x=xy(:,1);
y=xy(:,2);
if (x(k)==0) then
zz=find(x==0);
x(zz)=[];
y(zz)=[];
ZERO_POINT=%t
else
x(k)=[];
y(k)=[];
ZERO_POINT=%f
end
btn=-1
while ~(btn==3 | btn==0| btn==10| btn==-5)
rep=xgetmouse([%t %t]); xc=rep(1);yc=rep(2);btn=rep(3);
if (ZERO_POINT) then
xc=0;
else
if (xc<=0) then
zz=find(x==0);
x(zz)=[];
y(zz)=[];
ZERO_POINT=%t;
xc=0;
end
end
xt=[x;xc];
yt=[y;yc];
[xt,k2]=gsort(xt,"r","i");yt=yt(k2)
xyt=[xt,yt];
drawlater();
points.data=xyt;
[orpar,oipar]=drawSplin(a,xyt,oipar,orpar);
drawnow()
end
endfunction
//==========================================================
function rectx = findrect(a)
splines=a.children(1).children
points=a.children(2).children
if (points.data==[]) then
rectx=a.data_bounds;
return;
end
ymx1=max(splines.data(:,2));
ymn1=min(splines.data(:,2))
xmx=max(points.data(:,1));
xmn=min(points.data(:,1));
ymx=max(points.data(:,2));
ymn=min(points.data(:,2));
XMX=max(0,xmx);
XMN=max(0,xmn);
YMX=max(ymx,ymx1);
YMN=min(ymn,ymn1);
dx=XMX-XMN;
dy=YMX-YMN
if dx==0 then
dx=max(XMX/2,1),
end;
XMX=XMX+dx/50
if dy==0 then
dy=max(YMX/2,1),
end;
YMN=YMN-dy/50;
YMX=YMX+dy/50;
rectx=[XMN,YMN;XMX,YMX];
endfunction
//============================================================
function [tok,xyo]=ReadExcel()
TA=["A";"B";"C";"D";"E";"F";"G";"H";"I";"J";"K";"L";"M";"N";"O";"P"; ...
"Q";"R";"S";"T";"U";"V";"W";"X";"Y";"Z";"a";"b";"c";"d";"e";"f"; ...
"g";"h";"i";"j";"k";"l";"m";"n";"o";"p";"q";"r";"s";"t";"u";"v"; ...
"w";"x";"y";"z"];
TN=["0","1","2","3","4","5","6","7","8","9"];
xyo=[];
tok=%f;
while %t
[zok,filen,sheetN,xa,ya]=scicos_getvalue("Excel data file ",["Filename";"Sheet #"+...
" ";"X[start:Stop]";"Y[start:stop]"],list("str",1, ...
"vec",1,"str",1, ...
"str",1), ...
list(["Classeur1.xls"],["1"],["C5:C25"],["D5:D25"]));
if ~zok then
break,
end
try
[fd,SST,Sheetnames,Sheetpos] = xls_open(filen);
catch
xinfo("Scicos cannot find the excel file:"+filen);
break;
end
try
N=size(Sheetnames,"*");
if ((sheetN<=N) &(sheetN>0)) then
[Value,TextInd] = xls_read(fd,Sheetpos(sheetN))
mclose(fd)
end
xa=strsubst(xa," ","");
px=strindex(xa,":");
ya=strsubst(ya," ","");
py=strindex(ya,":");
x1=part(xa,1:px-1);
x2=part(xa,px+1:length(xa));
y1=part(ya,1:py-1);
y2=part(ya,py+1:length(ya));
x1p=min(strindex(x1,TN));
if x1p==[] then,
xinfo("Bad address in X:"+x1);
break,
end
x11=part(x1,1:x1p-1);
x12=part(x1,x1p:length(x1));
x2p=min(strindex(x2,TN));
if x2p==[] then,
xinfo("Bad address in X:"+x2);
break, end
x21=part(x2,1:x2p-1);
x22=part(x2,x2p:length(x2));
y1p=min(strindex(y1,TN));
if y1p==[] then,
xinfo("Bad address in Y:"+y1);
break,
end
y11=part(y1,1:y1p-1);
y12=part(y1,y1p:length(y1));
y2p=min(strindex(y2,TN));
if y2p==[] then,
xinfo("Bad address in Y:"+y2);
break,
end
y21=part(y2,1:y2p-1);
y22=part(y2,y2p:length(y2));
// x11 x12: x21 x22
lx11=length(x11);
lx21=length(x21);
ly11=length(y11);
ly21=length(y21)
xstC=0;
for i=1:lx11,
xstC=xstC+modulo(find(TA==part(x11,lx11-i+1)),26)*26^(i-1);
end
xenC=0;
for i=1:lx21,
xenC=xenC+modulo(find(TA==part(x21,lx21-i+1)),26)*26^(i-1);
end
ystC=0;
for i=1:ly11,
ystC=ystC+modulo(find(TA==part(y11,ly11-i+1)),26)*26^(i-1);
end
yenC=0;
for i=1:ly11,
yenC=yenC+modulo(find(TA==part(y21,ly21-i+1)),26)*26^(i-1);
end
xstR=evstr(x12);
xenR=evstr(x22);
ystR=evstr(y12);
yenR=evstr(y22);
[mv,nv]=size(Value)
if ~(xstR<=mv & xstR>0 & xenR<=mv & xenR>0&ystR<=mv & ystR>0¥R<=mv¥R>0 ) then
xinfo("error in Row data addresses");
break
end
if ~(xstC<=nv & xstC>0 & xenC<=nv & xenC>0&ystC<=nv & ystC>0¥C<=nv¥C>0 ) then
xinfo("error in Column data addresses");
break
end
xo=Value(min(xstR,xenR):max(xstR,xenR),min(xstC,xenC):max(xstC,xenC));
yo=Value(min(ystR,yenR):max(ystR,yenR),min(ystC,yenC):max(ystC,yenC));
[nx,mx]=size(xo);// adjusting the x and y size
[ny,my]=size(yo);
N=min(nx,ny);
xo=xo(1:N,:);
yo=yo(1:N,:);
xyo=[xo,yo];
[xyo]=cleandata(xyo)
tok=%t;
break,
catch
xinfo(" Scicos cannot read your Excel file, please verify"+...
" the parameters ");
break
end
end
endfunction
//---------------------------------------------------------------
function [xyo]=cleandata(xye)
xe=xye(:,1)
ye=xye(:,2)
[nx,mx]=size(xe);// adjusting the x and y size
[ny,my]=size(ye);
N=min(nx,ny);
xe=xe(1:N,:);
ye=ye(1:N,:);
// checking for NULL data
for i=1:N
if (xe(i)<>xe(i)) then
xinfo("x contains no data:x("+string(i)+")");
return;
end
if (ye(i)<>ye(i)) then
xinfo("Y contains no data:y("+string(i)+")");
return;
end
end
zz=find(xe<0);
xe(zz)=[];
ye(zz)=[]
if (find(xe==0)==[]) then // add zero point
xe($+1)=0;
ye($+1)=0;
end
[xo,k2]=gsort(xe,"r","i");
yo=ye(k2)
xyo=[xo,yo];
endfunction
//---------------------------------------------------------------
function [orpar,oipar] = AutoScale(a,xy,inipar,inrpar)
drawlater();
oipar = inipar
orpar = inrpar
points = a.children(2).children
splines = a.children(1).children
points.data = xy;
splines.data = xy;
[orpar,oipar] = drawSplin(a,xy,oipar,orpar);
rectx=findrect(a);
a.data_bounds = rectx;
drawnow()
endfunction
//============================
function METHOD = getmethod(order)
select order
case 0 then,
METHOD="zero order"
case 1 then,
METHOD="linear"
case 2 then,
METHOD="order 2"
case 3 then,
METHOD="not_a_knot"
case 4 then,
METHOD="periodic"
case 5 then,
METHOD="monotone"
case 6 then,
METHOD="fast"
case 7 then,
METHOD="clamped"
end
endfunction
//=======================================
function [sok,xye] = ReadFromFile()
xye=[];sok=%f;
while %t
[sok,filen,Cformat,Cx,Cy]=scicos_getvalue("Text data file ",["Filename";"Reading [C] f"+...
"ormat";"Abscissa column";"Output"+...
" column"],list("str",1,"str",1,"vec",1,"vec",1), ...
list(["mydatafile.dat"],["%g %g"],["1"],["2"]));
if ~sok then
break,
end
px=strindex(Cformat,"%");
NC=size(px,"*");
if NC==[] then,
xinfo("Bad format in reading data file");
sok=%f;
break;
end
Lx=[];
try
fd=mopen(filen,"r");
Lx=mfscanf(-1,fd,Cformat);
mclose(fd);
catch
xinfo("Scicos cannot open the data file: " + filen);
break;
end
[nD,mD] = size(Lx);
if ((mD==0) | (nD==0)) then,
xinfo("No data read");
sok=%f;
break;
end
if (mD<>NC) then,
xinfo("Bad format");
sok=%f;
break;
end
xe=Lx(:,Cx);
ye=Lx(:,Cy);
xye=[xe,ye];
[xye]=cleandata(xye)
sok=%t;
break,
end
endfunction
//=======================================
function [sok]=SaveToFile(xye)
xe=xye(:,1)
ye=xye(:,2)
sok=%f;
while %t
[sok,filen,Cformat]=scicos_getvalue("Text data file ",["Filename";"Writing [C] f"+...
"ormat"],list("str",1,"str",1), ...
list(["mydatafile.dat"],["%g %g"]));
if ~sok then
break,
end
px=strindex(Cformat,"%");
NC=size(px,"*");
if NC<>2 then,
xinfo("Bad format in writing data file");
sok=%f;
break;
end
Cformat=Cformat+"\n";
try
fd=mopen(filen,"w");
mfprintf(fd,Cformat,xe,ye);
mclose(fd);
catch
xinfo("Scicos cannot open the data file:"+filen);
break;
end
sok=%t;
break,
end
endfunction
//=========================================================
function [X,Y,orpar]=Do_Spline(N,order,x,y)
X=[];
Y=[];
orpar=[];
METHOD=getmethod(order);
if (METHOD=="zero order") then
X=x(1);
Y=y(1);
for i=1:N-1
X=[X;x(i);x(i+1);x(i+1)];
Y=[Y;y(i);y(i);y(i+1)];
end
return
end
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
if (METHOD=="linear") then
X=[];
for i=1:N
X=[X;x(i)];
Y=[Y;y(i)];
end
return
end
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
if (N<25) then
NP=10;
else
if (N<50) then
NP=5;
else
if (N<100) then
NP=2;
else
if (N<200) then
NP=1;
else
NP=0;
end;
end;
end;
end
for i=1:N-1
X=[X;linspace(x(i),x(i+1),NP+2)']; // pour tous sauf "linear" et "zero order"
end
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
if (N>2) & (METHOD=="order 2") then
Z=ORDER2(x,y);
A=Z(1:N-1);
B=Z(N:2*N-2);
C=Z(2*N-1:3*N-3);
for j=1:size(X,"*")
for i=N-1:-1:1
if X(j)>=x(i) then,
break;
end
end
Y(j)=A(i)*(X(j)-x(i))^2+B(i)*(X(j)-x(i))+C(i);
end
orpar=matrix(Z,-1,1)
end
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
if (METHOD=="not_a_knot") then
try
d = splin(x, y, METHOD);
Y = interp(X, x, y, d);
orpar=d(:);
catch
xinfo("ERROR in SPLINE: "+METHOD)
end
end
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
if (METHOD=="periodic") then
if y(1)<>y(N) then
y(N)=y(1)
end
try
d = splin(x, y,METHOD);
Y = interp(X, x, y, d);
orpar=d(:);
catch
xinfo("ERROR in SPLINE: "+METHOD)
end
end
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
if (METHOD=="monotone" ) then
try
d = splin(x, y, METHOD);
Y = interp(X, x, y, d);
orpar=d(:);
catch
xinfo("ERROR in SPLINE: "+METHOD)
end
end
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
if (METHOD=="fast") then
try
d = splin(x, y, METHOD);
Y = interp(X, x, y, d);
orpar=d(:);
catch
xinfo("ERROR in SPLINE: "+METHOD)
end
end
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
if (METHOD=="clamped") then
try
d = splin(x, y, METHOD,[0;0]);
Y = interp(X, x, y, d);
orpar=d(:);
catch
xinfo("ERROR in SPLINE: "+METHOD)
end
end
endfunction
//=================================================
function [Z]=ORDER2(x,y)
N=size(x,"*")-1;
A=zeros(3*N-1,N*3);
B=zeros(3*N-1,1);
for i=1:N
j=3*(i-1)+1;
A(j,i+2*N)=1;
B(j)=y(i);
A(j+1,i)=(x(i+1)-x(i))^2;
A(j+1,i+N)=x(i+1)-x(i);
A(j+1,i+2*N)=1;
B(j+1)=y(i+1);
end
for i=1:N-1
j=3*(i-1)+1;
A(j+2,i)=2*(x(i+1)-x(i));
A(j+2,i+N)=1;
A(j+2,i+N+1)=-1;
end
Q=zeros(3*N,3*N);
for i=1:N
Q(i,i)=4*(x(i+1)-x(i))^2
Q(i,i+N)=2*(x(i+1)-x(i))
Q(i+N,i)=2*(x(i+1)-x(i))
Q(i+N,i+N)=1;
end
At=[Q,A';A,zeros(3*N-1,3*N-1)]
Bt=[zeros(3*N,1);B]
Zt=At\Bt;
Z=Zt(1:3*N,1)
endfunction
//===================================================