CURVE_c
function [x,y,typ]=CURVE_c(job,arg1,arg2)
// Masoud Najafi 07/2007 --------
// origine: serge Steer, Habib Jreij INRIA 1993
// Copyright INRIA
x=[];y=[];typ=[];
select job
case 'plot' then
standard_draw(arg1)
case 'getinputs' then
[x,y,typ]=standard_inputs(arg1)
case 'getoutputs' then
[x,y,typ]=standard_outputs(arg1)
case 'getorigin' then
[x,y]=standard_origin(arg1)
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]=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, x_message('incompatible size of x and y'); 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 ~exists('curwin') then
gh=gcf();
curwin=gh.figure_id
end
save_curwin=curwin;
curwin=max(winsid())+1;
[orpar,oipar,ok]=poke_point(xy,ipar,rpar);
curwin=save_curwin;
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=['rpar=arg1.model.rpar;n=model.ipar(1);order=model.ipar(2);';
'xx=rpar(1:n);yy=rpar(n+1:2*n);';
'[XX,YY,rpardummy]=Do_Spline(n,order,xx,yy)';
'xmx=maxi(XX);xmn=mini(XX);';
'ymx=maxi(YY);ymn=mini(YY);';
'dx=xmx-xmn;if dx==0 then dx=maxi(xmx/2,1);end';
'xmn=xmn-dx/20;xmx=xmx+dx/20;';
'dy=ymx-ymn;if dy==0 then dy=maxi(ymx/2,1);end;';
'ymn=ymn-dy/20;ymx=ymx+dy/20;';
'xx2=orig(1)+sz(1)*((XX-xmn)/(xmx-xmn));';
'yy2=orig(2)+sz(2)*((YY-ymn)/(ymx-ymn));';
'xset(''color'',2)';
'xpoly(xx2,yy2,''lines'');']
x=standard_define([2 2],model,exprs,gr_i)
end
endfunction
function [rpar,ipar,ok]=poke_point(ixy,iparin,rparin)
[lhs,rhs]=argn(0)
//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 == -100 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 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=maxi(xy(:,1));xmn=mini(xy(:,1)),xmn=max(xmn,0);
ymx=maxi(xy(:,2));ymn=mini(xy(:,2));
dx=xmx-xmn;dy=ymx-ymn
if dx==0 then dx=maxi(xmx/2,1),end;
xmx=xmx+dx/50;
if dy==0 then dy=maxi(ymx/2,1),end;
ymn=ymn-dy/50;ymx=ymx+dy/50;
rect=[xmn,ymn;xmx,ymx];
//===================================================================
f=scf();
if ~MSDOS then
delmenu(curwin,'3D Rot.')
delmenu(curwin,'Edit')
delmenu(curwin,'File')
delmenu(curwin,'Insert')
else
hidetoolbar(curwin)
// French
delmenu(curwin,'&Fichier')
delmenu(curwin,'&Editer')
delmenu(curwin,'&Outils')
delmenu(curwin,'&Inserer')
// English
delmenu(curwin,'&File')
delmenu(curwin,'&Edit')
delmenu(curwin,'&Tools')
delmenu(curwin,'&Insert')
end
//menuss=menus;menuss(1)=menus(1)(2:$);menubar(curwin,menuss)
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 f"+...
"ile','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)
//===================================================================
//initial draw
f.pixmap='off';
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();
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]=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 soululd be positive')
mok=%f;
end
if mok then
a.data_bounds=[xmn1,ymn1;xmx1,ymx1];
end
end
//drawnow();//show_pixmap();
//-------------------------------------------------------------------
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]=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]=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]=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]=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]=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]=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]=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]=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'''
x_message([t1;t2;t3;t4;t5]);
//---------------------------------------------------------------
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]=mini(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,:)
f.pixmap='on';
drawlater();
points.data=xy;
[rpar,ipar]=drawSplin(a,xy,ipar,rpar);
show_pixmap();
drawnow()
f.pixmap='off';
end
end
if (HIT)&(btn==2 | btn==5) then // remove point
f.pixmap='on';
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);
show_pixmap();
drawnow()
f.pixmap='off';
end
if (HIT)&(btn==0) then // move point
f.pixmap='on';
[xy,rpar,ipar]=movept(a,xy,ipar,rpar,k)
f.pixmap='off';
end
if (HIT)&(btn==10) then // change data:: double click
[mok,xt,yt]=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)
f.pixmap='on';
drawlater();
points.data=xy;
[rpar,ipar]=AutoScale(a,xy,ipar,rpar)
show_pixmap();
drawnow()
f.pixmap='off';
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=maxi(points.data(:,1)); xmn=mini(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);
show_pixmap();
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=maxi(points.data(:,1));xmn=mini(points.data(:,1));
ymx=maxi(points.data(:,2));ymn=mini(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=maxi(XMX/2,1),end;
XMX=XMX+dx/50
if dy==0 then dy=maxi(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]=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 canot 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;
show_pixmap();
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]=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 canot 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]=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 canot 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
//===================================================