"an electricity package"

"sinwave : generate a sinus signal"
"sqrwave : generate a square signal"
"triwave : generate a triangle signal"
"pd2mwave: generate a pd2 bridge signal (2 diodes and 2 rectifiers)"
"pd2fwave: generate a pd2 bridge signal (4 rectifiers)"
"scope   : plot an oscilloscope grid"
"plotZ   : plot of complex impedance (module and arg)"
"plotFFT : plot the Fourier Transform of a signal"
"bode    : plot the bode diagram"
j=1i;

function findYvrange(Vmin,Vmax)
	.. volt range
	cv=[100,50,20,10,5,2,1,0.5,0.2,0.1,0.05,0.02,0.01,0.005,0.002,0.001];
	.. find the right range
	cvi=1;
	for i=1 to length(cv)-1
		if (4*cv[i+1] < Vmax) || (-4*cv[i+1] > Vmin);
			cvi=i;
			break;
		endif;
	end;
	return cv[cvi];
endfunction

function findXvrange(Vmin,Vmax)
	.. volt range
	cv=[100,50,20,10,5,2,1,0.5,0.2,0.1,0.05,0.02,0.01,0.005,0.002,0.001];
	.. find the right range
	cvi=1;
	for i=1 to length(cv)-1
		if (5*cv[i+1] < Vmax) || (-5*cv[i+1] > Vmin);
			cvi=i;
			break;
		endif;
	end;
	return cv[cvi];
endfunction

function findtrange(Tmax)
	.. time range
	ct=[10,5,2,1,0.5,0.2,0.1,0.05,0.02,0.01,0.005,0.002,0.001,0.0005,0.0002,0.0001,0.00005,0.00002,0.00001,0.000005,0.000002,0.000001,0.0000005];
	.. find the right range
	cti=1;
	for i=1 to length(ct)-1
		if (10*ct[i+1]+epsilon < Tmax);
			cti=i;
			break;
		endif;
	end;
	return ct[cti];
endfunction

function drawgrid(dt,dv)
	ho=holding(1);
	lw=linewidth(1);
	c=color(1);
	.. draw vertical lines
	ls=linestyle(".");
	for i=1 to 9
		if i != 5
			plot([i*dt,i*dt],[-4*dv,4*dv]);
		endif;
	end;
	linestyle(ls);
	plot([5*dt,5*dt],[-4*dv,4*dv]);
	for j=1 to 49
		plot([j/5*dt,j/5*dt],[-0.1*dv,0.1*dv]);
	end;

	.. draw horizontal lines
	ls=linestyle(".");
	for i=-3 to 3
		if i != 0;
			plot([0,10*dt],[i*dv,i*dv]);
		endif;
	end;
	linestyle(ls);
	plot([0,10*dt],[0,0]);
	for j=-19 to 19
		plot([4.9*dt,5.1*dt],[j/5*dv,j/5*dv]);
	end;
	
	color(c);
	linewidth(lw);
	holding(ho);
	return 0;
endfunction;

function drawXYgrid(dv1,dv2)
	ho=holding(1);
	lw=linewidth(1);
	c=color(1);
	.. draw vertical lines
	ls=linestyle(".");
	for i=-4 to 4
		if i != 0
			plot([i*dv1,i*dv1],[-4*dv2,4*dv2]);
		endif;
	end;
	linestyle(ls);
	plot([0,0],[-4*dv2,4*dv2]);
	for j=-24 to 24
		plot([j/5*dv1,j/5*dv1],[-0.1*dv2,0.1*dv2]);
	end;

	.. draw horizontal lines
	ls=linestyle(".");
	for i=-3 to 3
		if i != 0;
			plot([-5*dv1,5*dv1],[i*dv2,i*dv2]);
		endif;
	end;
	linestyle(ls);
	plot([-5*dv1,5*dv1],[0,0]);
	for j=-19 to 19
		plot([-0.1*dv1,0.1*dv1],[j/5*dv2,j/5*dv2]);
	end;
	
	color(c);
	linewidth(lw);
	holding(ho);
	return 0;
endfunction;

function scope (x=0,y1=0,y2=0)
## o scope(x,y1,), scope(x,,y2) and scope(x,y1,y2) shows a scope grid and
##   the curves y1(x) or y2(x) or both.
## o scope(,y1,y2) shows the curve y2(y1) in XY mode.
## o scope() shows only the grid.
	if argn()==0;
		if !holding();clg; frame();endif;
		p=plot();
		dt=p[2]/10;
		dv1=p[4]/4;
		drawgrid(dt,dv1);
	elseif (argn()==1);
		error("not enough parameters (2 at least)");
	elseif (argn()==2) || ((argn()==3) && (y1==0));
		.. normal mode ; one curve
		if !holding();clg; frame();endif;
		dt=findtrange(max(x));
		if (y2==0);
			dv1=findYvrange(min(y1),max(y1));
		else
			dv1=findYvrange(min(y2),max(y2));
		endif;
		p=[0,10*dt,-4*dv1,4*dv1];setplot(p);
		drawgrid(dt,dv1);
		.. display ranges
		tp=toscreen([7.5*dt,3.8*dv1]);
		text("time : "|printf("%g s/div",dt),tp);
		tp[2]=tp[2]+textheight();
		if (y2==0);
			text("CH1  : "|printf("%g V/div",dv1),tp);
			ho=holding(1); plot(x,y1); holding(ho);
		else
			text("CH2  : "|printf("%g V/div",dv1),tp);
			ho=holding(1); plot(x,y2); holding(ho);
		endif;
	elseif (argn()==3) && (x==0);
		.. XY mode
		if !holding();clg; frame();endif;
		dv1=findXvrange(min(y1),max(y1));
		dv2=findYvrange(min(y2),max(y2));
		p=[-5*dv1,5*dv1,-4*dv2,4*dv2];setplot(p);
		drawXYgrid(dv1,dv2);
		.. display ranges
		tp=toscreen([2.5*dv1,3.8*dv2]);
		text("CH1 (X) : "|printf("%g V/div",dv1),tp);
		tp[2]=tp[2]+textheight();
		text("CH2 (Y) : "|printf("%g V/div",dv2),tp);
		ho=holding(1); plot(y1,y2); holding(ho);
	else
		.. normal display ; two curves
		if !holding();clg; frame();endif;
		dt=findtrange(max(x));
		dv1=findYvrange(min(y1),max(y1));
		dv2=findYvrange(min(y2),max(y2));
		p=[0,10*dt,-4*dv1,4*dv1];setplot(p);
		drawgrid(dt,dv1);

		.. display ranges
		tp=toscreen([7.5*dt,3.8*dv1]);
		text("time : "|printf("%g s/div",dt),tp);
		tp[2]=tp[2]+textheight();
		text("CH1  : "|printf("%g V/div",dv1),tp);
		tp[2]=tp[2]+textheight();
		text("CH2  : "|printf("%g V/div",dv2),tp);

		.. display a CH1 / CH2 mark
		tp0=toscreen([x[1],y1[1]]);
		tp0[1]=tp0[1]-5;tp0[2]=tp0[2]-textheight()/2;
		tp1=toscreen([x[1],dv1/dv2*y2[1]]);
		tp1[1]=tp1[1]-5;tp1[2]=tp1[2]-textheight()/2;
		for i=1 to 10;
			if abs(tp1[2]-tp0[2])>textheight(); break; endif;
			tp2=toscreen([x[i*dt/x[2]+1],y1[i*dt/x[2]+1]]);
			tp2[1]=tp2[1]-5;
			if tp2[2]<tp0[2]; tp2[2]=tp2[2]-textheight(); endif;
			tp3=toscreen([x[i*dt/x[2]+1],dv1/dv2*y2[i*dt/x[2]+1]]);
			tp3[1]=tp3[1]-5;
			if tp3[2]<tp1[2]; tp3[2]=tp3[2]+textheight(); endif;
			tp0=tp2;
			tp1=tp3;
		end;
		rtext("CH1>",tp0);
		rtext("CH2>",tp1);
		ho=holding(1); plot(x,y1); plot(x,dv1/dv2*y2); holding(ho);
	endif;
	return p;
endfunction

function sinwave(t,f,E=1,phi=0)
## returns a vector containing the square signal
##
## t     : time
## f     : frequency
## E     : RMS value
## phi   : origin phase (rad)
##
	useglobal;
	return E*sqrt(2)*sin(2*pi*f*t+phi);
endfunction;

function sqrwave(t,f,E,alpha=0.5)
## returns a vector containing the square signal
##
## t     : time
## f     : frequency
## E     : amplitude (from 0 to max)
## alpha : period ratio
##
	m=t;
	T=1/f;
	i=1;k=0;
	repeat;
		if t[i]-k*T<alpha*T;
			m[i]=E;
			i=i+1;
		elseif (t[i]-k*T>=alpha*T) && (t[i]-k*T<T);
			m[i]=-E;
			i=i+1;
		else
			k=k+1;
		endif;
		if i>length(t); break; endif;
	end;
	return m;
endfunction;

function triwave(t,f,E,alpha=0.5)
## returns a vector containing the triangle signal
##
## t     : time
## f     : frequency
## E     : amplitude (from 0 to max)
## alpha : period ratio
##
	m=t;
	T=1/f;
	i=1;k=0;
	a=2*E/(alpha*T);
	b=-2*E/((1-alpha)*T);
	b0=E*(1+alpha)/(1-alpha);
	repeat;
		if t[i]-k*T<alpha*T;
			m[i]=a*(t[i]-k*T)-E;
			i=i+1;
		elseif (t[i]-k*T>=alpha*T) && (t[i]-k*T<T);
			m[i]=b*(t[i]-k*T)+b0;
			i=i+1;
		else
			k=k+1;
		endif;
		if i>length(t); break; endif;
	end;
	return m;
endfunction;

function pd2mwave(t,f,V,alpha=0)
## returns a vector containing the values of the voltage across an RL rectifier bridge load.
## t     : time array
## f     : frequency
## V     : RMS value of the sine wave
## alpha : angle (°)
  m = t;
  T=1/f;
  i=1;k=0;
  repeat;
    if (t[i]-k*T<alpha/360*T) || ((t[i]-k*T>=T/2) && (t[i]-k*T<T/2+alpha/360*T));
      m[i]=0;i=i+1;
    elseif (t[i]-k*T>=alpha/360*T) && (t[i]-k*T<T/2);
      m[i]=V*sqrt(2)*sin(2*pi*f*(t[i]-k*T));i=i+1;
    elseif ((t[i]-k*T>=T/2+alpha/360*T) && (t[i]-k*T<T));
      m[i]=-V*sqrt(2)*sin(2*pi*f*(t[i]-k*T));i=i+1;
    else k=k+1;
    endif;
    if i>length(t);break;endif;
  end;
  return m;
endfunction

function pd2fwave(t,f,V,alpha=0)
## returns a vector the values of the voltage across an RL rectifier bridge load.
## t     : time array
## f     : frequency
## V     : RMS value of the sine wave
## alpha : angle (°)
  m = t;
  T=1/f;
  i=1;k=0;
  repeat;
    if (t[i]-k*T>=alpha/360*T) && (t[i]-k*T<T/2+alpha/360*T);
      m[i]=V*sqrt(2)*sin(2*pi*f*(t[i]-k*T));i=i+1;
    elseif (t[i]-k*T<alpha/360*T) || ((t[i]-k*T>=T/2+alpha/360*T) && (t[i]-k*T<T));
      m[i]=-V*sqrt(2)*sin(2*pi*f*(t[i]-k*T));i=i+1;
    else
      k=k+1;
    endif;
    if i>length(t);break;endif;
  end;
  return m;
endfunction


function plotZ(func,fmin,fmax,nb=200,fmes=-1,Zmes=-1,phmes=0)
## plots the amplitude and phase of a complex impedance
## from fmin Hz to fmax Hz
##
## func              : a function describing the impedance
## fmin, fmax        : frequency range
## nb                : number of points
## fmes, Zmes, phmes : experimental data
##
	h=textheight();b=textwidth();dh=floor((1024-7.5*h)/2);

	f=fmin:(fmax-fmin)/(nb-1):fmax;

	if isfunction(func);
		z=func(f,args());
	else
		z=expreval(func,f);
	endif;

	m=abs(z);
	ph=180/pi*arg(z);

	.. affichage de la courbe en fonction du temps
	window([8*b,1.5*h,1024-2*b,1.5*h+dh]);
	{t,form}=ticks(0,max(m));
	if max(m)>t[length(t)];
		ymax=2*t[length(t)]-t[length(t)-1];
	else
		ymax=t[length(t)];
	endif

	setplot([fmin,fmax,0,ymax]);
	xplot(f,m);
	rtext("f (Hz)",[1024-2*b,2.5*h+dh]);
	vcutext("Z (ohm)",[5,1.5*h+.5*dh]);
	hd=holding(1);
	if length(fmes)>0 && fmes!=-1;
		c=color(10);
		m=markerstyle("+");
		mark(fmes,Zmes);
		color(c);
	endif;

	window([8*b,4.5*h+dh,1024-2*b,1024-3*h]);
	xplot(f,ph);
	rtext("f (Hz)",[1024-2*b,1024-2*h]);
	vcutext("phase (°)",[5,4.5*h+1.5*dh]);
	if length(fmes)>0 && fmes!=-1;
		c=color(10);
		mark(fmes,phmes);
		markerstyle(m);
		color(c);
	endif;
	holding(hd);
	return z;
endfunction;

function plotFFT(func,T,nbe=256)
## plots the spectra of a signal defined by the function parameter
## over T seconds.
##
## func : the signal (function)
## T    : time
## nbe  : sample number (power of 2)
##
	h=textheight();b=textwidth();dh=floor((1024-7.5*h)/2);

	t=0:T/(nbe-1):T;
	if isfunction(func);
		y=func(t,args());
	else
		y=expreval(func,t);
	endif;
	.. affichage de la courbe en fonction du temps
	window([8*b,1.5*h,1024-2*b,1.5*h+dh]);xplot(t,y);
	rtext("t (s)",[1024-2*b,2.5*h+dh]);
	vcutext("signal",[5,1.5*h+.5*dh]);

	.. période et fréquence d'échantillonnage
	Te=T/(nbe-1);Fe=1/Te;
	.. calcul de la fft
	c=fft(y);
	f=0:Fe/(nbe-1):Fe/2;
	cn=matrix([1,nbe/2],0);
	cn[1]=abs(c[1])/nbe;
	for i=2 to nbe/2;
		cn[i]=(abs(c[i])+abs(c[nbe+2-i]))/nbe;
	end;
	hd=holding(1);
	window([8*b,4.5*h+dh,1024-2*b,1024-3*h]);xplot(f,cn);
	rtext("f (Hz)",[1024-2*b,1024-2*h]);
	vcutext("spectre",[5,4.5*h+1.5*dh]);
	holding(hd);
	return Fe;
endfunction;

..load logplot;

function bode(tf,fmin=0.01,fmax=100,pt=100,what=0,fmes=-1,Gmes=0,phmes=0)
## draws the bode's diagram of the transfert function tf(x)
##
## Example:
##
##   bode("1/(1+j*2*pi*x/1000)",10,10e3);
##
## draws the diagram from 10 and 1000 Hz for the specified
## transfert function.
##
## tf   = transfert function (variable -> x)
## fmin = start frequency
## fmax = end frequency
## pt   = samples
## what = 0 magnitude and phase both on the same graphic.
##        1 magnitude only
##        2 phase only
	expi = floor(log10(fmin));
	expf = ceil(log10(fmax));
	f   = (10)^linspace(expi,expf,pt);
	if isfunction(tf);
		T=tf(f,args());
	else
		T=expreval(tf,f);
	endif;
	gdb = 20*log10(abs(T));
	ph=arg(T);

	h=textheight();b=textwidth();

	if what==0;
		dh=floor((1023-7.5*h)/2);
		window([8*b,1.5*h,1023-2*b,1.5*h+dh]);
		xlogplot(f,gdb,10);
		hd=holding(1);
		if length(fmes)>0 && fmes!=-1;
			c=color(10);
			m=markerstyle("+");
			mark(log10(fmes),Gmes);
			color(c);
		endif;
		ctext("f (Hz)",[8*b+(1023-10*b)/2,3*h+dh]);
		vcutext("Gain (dB)",[5,1.5*h+dh/2]);
		window([8*b,4.5*h+dh,1023-2*b,1023-3*h]);
		xlogplot(f,ph/pi*180,10);
		if length(fmes)>0 && fmes!=-1;
			c=color(10);
			mark(log10(fmes),phmes);
			markerstyle(m);
			color(c);
		endif;
		ctext("f (Hz)",[8*b+(1023-10*b)/2,1023-1.5*h]);
		vcutext("Phase (°)",[5,4.5*h+1.5*dh]);
		holding(hd);
	elseif what==1;
		shrinkwindow();
		xlogplot(f,gdb,10);
		if length(fmes)>0 && fmes!=-1;
			hd=holding(1);
			c=color(10);
			m=markerstyle("+");
			mark(log10(fmes),Gmes);
			markerstyle(m);
			color(c);
			holding(hd);
		endif;
		ctext("f (Hz)",[8*b+(1023-10*b)/2,1023-h]);
		vcutext("Gain (dB)",[5,512]);
	else
		shrinkwindow();
		xlogplot(f,ph/pi*180,10);
		if length(fmes)>0 && fmes!=-1;
			hd=holding(1);
			c=color(10);
			m=markerstyle("+");
			mark(log10(fmes),phmes);
			markerstyle(m);
			color(c);
			holding(hd);
		endif;
		ctext("f (Hz)",[8*b+(1023-10*b)/2,1023-h]);
		vcutext("Phase (°)",[5,512]);
	endif;
	return 0;
endfunction;