// ............................................................ convoluzione.cc
#include <curses.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include "corso.h"

extern XVideoData xvd;

double gauss(double x,double w);
double lorentz(double x,double w);

main()
{
	int biPoints,halfPoints,i,ii,lg2n,nComm,nPar,nPoints,NuoviPunti,OldPoints,
	ridisegna,x0m,x0p,x1m,x1p,y0,y1;
	clock_t musec1,musec2;
	timespec ts1,ts2;
	double cOffs,cScale,gOffs,gScale,lOffs,lScale,oldWg,oldWl,wg,wl,x,yMax,yMin;
	double *convol,*gfun,*gtransf,*lfun,*ltransf,*convtransf;
	char answ[122],buff[40];
	char **comandi,**ParName,**ParVal;
	FILE *hnd;
	
	// ........................................................... tempo di sonno
	ts1.tv_sec=0;
	ts1.tv_nsec=10000000;
	// .......................................................... numero di punti
	//nPoints=512;oldWg=wg;
	oldWl=wl;
	nPoints=1024;
	biPoints=nPoints*2;
	halfPoints=nPoints/2;
	// ............................................................. modo grafico
	StartXWindow(800,600,0,0,"Convoluzione",2);
	// ................................................................... comandi
	nComm=1;
	comandi=new char * [nComm];
	comandi[0]="Fine";
	// ................................................................ parametri
	nPar=3;
	ParName=new char * [nPar];
	ParName[0]="Larghezza Gaussiana";
	ParName[1]="Larghezza Lorenziana";
	ParName[2]="Numero di punti";
	ParVal=new char * [nPar];
	for (i=0;i<nPar;i++) ParVal[i]=new char [28];
	// .......................................................... valori iniziali
	wg=10.0;
	wl=10.0;
	sprintf(ParVal[0],"%.5lf",wg);
	sprintf(ParVal[1],"%.5lf",wl);
	nPoints=1024;
	sprintf(ParVal[2],"%d",nPoints);
	OldPoints=0;
	// ..................................................................... menu
	XSync(xvd.display,TRUE);
	XWinPermMenu(1,"Convoluzione",ParName,ParVal,nPar,28,comandi,nComm,0);
	// ..........................................................................
	oldWg=-1.0;
	oldWl=-1.0;
	ridisegna=TRUE;
	NuoviPunti=TRUE;
	// .................................................................... ciclo
	for (;;)
	{
		nanosleep(&ts1,&ts2);
		if (NuoviPunti)
		{ // ....................................... il numero di punti e' cambiato
			if (OldPoints>0)
			{ // ................................................. pulisci la memoria
				delete [] convtransf;
				delete [] ltransf;
				delete [] gtransf;
				delete [] lfun;
				delete [] gfun;
				delete [] convol;
			}
			// ...................................................... nuove dimensioni
			biPoints=nPoints*2;
			halfPoints=nPoints/2;
			// ................................................... rialloca la memoria
			convol=new double [biPoints];
			gfun=new double [biPoints];
			lfun=new double [biPoints];
			gtransf=new double [biPoints];
			ltransf=new double [biPoints];
			convtransf=new double[biPoints];
		}
		// ............................................................... gaussiana
		yMin=1.0e100;
		yMax=1.0e-100;oldWg=wg;
		oldWl=wl;
		for (i=0;i<halfPoints;i++)
		{
			x=double(i)/4.0;
			gfun[2*i]=gfun[biPoints-2-2*i]=gauss(x,wg);//     parte reale
			gfun[2*i+1]=gfun[biPoints-1-2*i]=0.0;//           parte immaginaria
			if (gfun[2*i]<yMin) yMin=gfun[2*i];
			if (gfun[2*i]>yMax) yMax=gfun[2*i];
		}
		gScale=400.0/(yMax-yMin);
		gOffs=yMin;
		// ............................................................ lorentziana
		yMin=1.0e100;
		yMax=1.0e-100;
		for (i=0;i<halfPoints;i++)
		{
			x=double(i)/4.0;
			lfun[2*i]=lfun[biPoints-2-2*i]=lorentz(x,wl);//     parte reale
			lfun[2*i+1]=lfun[biPoints-1-2*i]=0.0;//           parte immaginaria
			if (lfun[2*i]<yMin) yMin=lfun[2*i];
			if (lfun[2*i]>yMax) yMax=lfun[2*i];
		}
		lScale=400.0/(yMax-yMin);
		lOffs=yMin;
		// .................................................. trasformata gaussiana
		for (i=0;i<biPoints;i++) gtransf[i]=gfun[i];
		fft(gtransf,nPoints,+1);
		// ................................................ trasformata lorentziana
		for (i=0;i<biPoints;i++) ltransf[i]=lfun[i];
		fft(ltransf,nPoints,+1);
		// ............................................... trasformata convoluzione
		for (i=0;i<nPoints;i++)
		{
			convtransf[2*i]=gtransf[2*i]*ltransf[2*i]-gtransf[1+2*i]*ltransf[1+2*i];
			convtransf[1+2*i]=gtransf[2*i]*ltransf[1+2*i]+gtransf[1+2*i]*ltransf[2*i];
		}
		// ........................................... antitrasformata convoluzione
		for (i=0;i<biPoints;i++) convol[i]=convtransf[i];
		fft(convol,nPoints,-1);
		// ........................................... normalizzazione convoluzione
		yMin=1.0e100;
		yMax=1.0e-100;
		for (i=0;i<halfPoints;i++)
		{
			if (convol[2*i]<yMin) yMin=convol[2*i];
			if (convol[2*i]>yMax) yMax=convol[2*i];
		}
		cScale=400.0/(yMax-yMin);
		cOffs=yMin;
		// ................................................................ disegna
		if (ridisegna)
		{ // .............................. seleziona e pulisci il contesto grafico
			XSetForeground(xvd.display,xvd.gc,xvd.white);
			XFillRectangle(xvd.display,xvd.backpix[0],xvd.gc,0,0,800,600);
			// ................................................. disegna la gaussiana
			XSetForeground(xvd.display,xvd.gc,xvd.red[1]);
			y0=int((499.0-gScale*(gfun[0]-gOffs))+0.5);
			x0m=x0p=400;
			for (i=0;i<halfPoints;i++)
			{
				y1=int((499.0-gScale*(gfun[2*i]-gOffs))+0.5);
				x1m=x0m-1;
				x1p=x0p+1;
				XDrawLine(xvd.display,xvd.backpix[0],xvd.gc,x1m,y1,x0m,y0);
				XDrawLine(xvd.display,xvd.backpix[0],xvd.gc,x1p,y1,x0p,y0);
				y0=y1;
				x0m=x1m;
				x0p=x1p;
			}
			XPrintf(0,10,540,1,TRUE,"Gaussiana");
			// ............................................... disegna la lorentziana
			XSetForeground(xvd.display,xvd.gc,xvd.blue[1]);
			y0=int((499.0-lScale*(lfun[0]-lOffs))+0.5);
			x0m=x0p=400;
			for (i=0;i<halfPoints;i++)
			{
				y1=int((499.0-lScale*(lfun[2*i]-lOffs))+0.5);
				x1m=x0m-1;
				x1p=x0p+1;
				XDrawLine(xvd.display,xvd.backpix[0],xvd.gc,x1m,y1,x0m,y0);
				XDrawLine(xvd.display,xvd.backpix[0],xvd.gc,x1p,y1,x0p,y0);
				y0=y1;
				x0m=x1m;
				x0p=x1p;
			}
			XPrintf(0,10,560,1,TRUE,"Lorentziana");
			// .............................................. disegna la convoluzione
			XSetForeground(xvd.display,xvd.gc,xvd.black);
			y0=int((499.0-cScale*(convol[0]-cOffs))+0.5);
			x0m=x0p=400;
			for (i=0;i<halfPoints;i++)
			{
				y1=int((499.0-cScale*(convol[2*i]-cOffs))+0.5);
				x1m=x0m-1;
				x1p=x0p+1;
				XDrawLine(xvd.display,xvd.backpix[0],xvd.gc,x1m,y1,x0m,y0);
				XDrawLine(xvd.display,xvd.backpix[0],xvd.gc,x1p,y1,x0p,y0);
				y0=y1;
				x0m=x1m;
				x0p=x1p;
				if (x0m<0 || x0p>xvd.width[0]) break;
			}
			XPrintf(0,10,580,1,TRUE,"Convoluzione (Voigt)");
			XRedraw();
		}
		// ....................................................... chiama il menu
		for (;;)
		{
			nanosleep(&ts1,&ts2);
			XSync(xvd.display,FALSE);
			i=XWinPermMenu(1,"Convoluzione",ParName,ParVal,nPar,28,comandi,nComm,1);
			if (i==0) break;
			nPoints=atoi(ParVal[2]);
			for(lg2n=1;;lg2n++) if ((nPoints>>lg2n)==1) break;
			if ((1<<lg2n)!=nPoints) snprintf(ParVal[2],40,"%d non va bene!",nPoints);
			else break;
			XRedraw();
		}
		if (i==0) break;
		wg=atof(ParVal[0]);
		wl=atof(ParVal[1]);
		ridisegna=(wg != oldWg || wl !=oldWl || nPoints!=OldPoints);
		NuoviPunti=(nPoints!=OldPoints);
		OldPoints=nPoints;
		oldWg=wg;
		oldWl=wl;
		XSync(xvd.display,TRUE);
		XRedraw();
	}
	// ..........................................................................
	XWinPermMenu(1,"Convoluzione",ParName,ParVal,nPar,28,comandi,nComm,2);
	// ................................................................. clean-up
	for (i=0;i<nPar;i++) delete ParVal[i];
	delete [] ParVal;
	delete [] ParName;
	delete [] convtransf;
	delete [] ltransf;
	delete [] gtransf;
	delete [] lfun;
	delete [] gfun;
	delete [] convol;
	delete [] comandi;
	CloseXWindow();
}

// ------------------------------------------------------------------------

double gauss(double x,double w)
{
	double xx;
	
	xx=x/w;
	return exp(-xx*xx);
}
double lorentz(double x,double w)
{
	return w*w/(x*x+w*w);
}