/*				 19.MAY.2005			 fittacurve.cc
							   by Giovanni Moruzzi


*/
#include <X11/Xlib.h>
#include <X11/Xutil.h>

#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <unistd.h>
#include <string.h>
#include <time.h>

#include "corso.h"

extern XVideoData xvd;

main()
{
    int i,iComm,j,l0,nc,nComm,nData,nn,nPar,nShape;
    int *ParFix;
    double a,b,bias,c,center,delta,dx,height,hwhm,path,slope,tolerance,xx;
    double *eps,*par,*ParMax,*ParMin,*shape,*x,*y,*ycal;
    char **comandi,**LimName,**ParName,**ParVal,**titolo;
    char answ[16][25];
    
    // ....................................................... fit tolerance
    tolerance=1.0e-6;
    // .................................................. apri la finestra X
    StartXWindow(800,600,0,0,"Fittacurve");
    // .......................................................... dimensioni
    nData=760;
    nn=1024;
    nPar=7;
    nShape=512;
    // .....................................................................
    ParFix=new int [nPar];
    eps=new double [nData];
    par=new double [nPar];
    ParMax=new double [nPar];
    ParMin=new double [nPar];
    shape=new double [nShape];
    x=new double [nData];
    y=new double [nData];
    ycal=new double [nData];
    // ............................................................ x values
    for (i=0;i<nData;i++) x[i]=double(i);
    // ............................................................ y values
    for (i=0;i<nData;i++) y[i]=ycal[i]=0.0;
    // .............................................................. errors
    for (i=0;i<nData;i++) eps[i]=2.0e-4;
    // .................................................... fixed parameters
    for (i=0;i<nPar;i++) ParFix[i]=FALSE;
    // ParFix[2]=TRUE;
    // .................................................. parametri iniziali
    ParName=new (char *) [nPar];
    ParVal=new (char *) [2*nPar];
    for (i=0;i<(2*nPar);i++) ParVal[i]=new char [24];
    ParName[0]="Gauss width";
    ParName[1]="Lorentz width";
    ParName[2]="Path length";
    ParName[3]="bias";
    ParName[4]="slope";
    ParName[5]="center";
    ParName[6]="height";
    par[0]=32.0;     //   Gaussian width
    par[1]=30.0;     //   Lorentzian width
    par[2]=30.0;     //   mirror path
    par[3]=0.0;      //   bias
    par[4]=0.0;      //   bias slope
    par[5]=300.0;    //   line center
    par[6]=300.0;    //   line height
    for (i=0;i<nPar;i++) sprintf(ParVal[i],"%.5lf",par[i]);
    // ..................................................... nomi dei limiti
    LimName=new (char *) [2*nPar];
    for (i=0;i<(2*nPar);i++) LimName[i]=new char [24];
    for (i=0;i<nPar;i++)
    {
	sprintf(LimName[2*i],"minimum %s",ParName[i]);
	sprintf(LimName[1+2*i],"maximum %s",ParName[i]);
    }
    // ............................................................. comandi
    nComm=3;
    comandi=new (char *) [nComm];
    comandi[0]="accept";
    comandi[1]="Exit";
    // ............................................................. titoli
    titolo=new (char *) [2];
    titolo[0]="Scegli la curva da fittare";
    titolo[1]="Scegli l'approssimazione iniziale";
    // .............................................. loop curva da fittare
    for (;;)
    { 
	XSetForeground(xvd.display,xvd.gc,xvd.white.pixel);
	XFillRectangle(xvd.display,xvd.backpix,xvd.gc,0,0,
		       xvd.width,xvd.height);
	// .................... ................... calcola la convoluzione
	SincVoigt(shape,nShape,&hwhm,par[0],par[1],par[2]);
	// ............................................. calcola il profilo
	for (i=0;i<nData;i++)
	{ 
	    y[i]=par[3]+x[i]*par[4];
	    xx=fabs(par[5]-x[i]);
	    if (xx>double(nShape-1)) continue;
	    l0=int(xx);
	    dx=xx-double(l0);
	    if (dx>0.5) {l0++;dx-=1.0;}
	    c=shape[l0];
	    a=(shape[abs(l0-1)]-2*c+shape[l0+1])/2.0;
	    b=(shape[l0+1]-shape[abs(l0-1)])/2.0;
	    delta=c+dx*(a*dx+b);
	    y[i]+=delta*par[6];
	}
	XPlot2Arrays(20,10,760,400,1,x,y,ycal,nData);
	// .......................................................... menu
	iComm=XMenu(0,titolo[0],ParName,ParVal,nPar,20,comandi,2,30,420);
	if (iComm==1) break;
	for (i=0;i<nPar;i++) par[i]=atof(ParVal[i]);
    }
    // ======================================================== loop di fit
    for (;;)
    { // ............................. scegli l'approssimazione di partenza
	comandi[0]="accept";
	for (i=0;i<nPar;i++) sprintf(ParVal[i],"%.5lf",par[i]);
	for (;;)
	{ 
	    XSetForeground(xvd.display,xvd.gc,xvd.white.pixel);
	    XFillRectangle(xvd.display,xvd.backpix,xvd.gc,0,0,
			   xvd.width,xvd.height);
	    // .................... ............... calcola la convoluzione
	    SincVoigt(shape,nShape,&hwhm,par[0],par[1],par[2]);
	    // ......................................... calcola il profilo
	    for (i=0;i<nData;i++)
	    { 
		ycal[i]=par[3]+x[i]*par[4];
		xx=fabs(par[5]-x[i]);
		if (xx>double(nShape-1)) continue;
		l0=int(xx);
		dx=xx-double(l0);
		if (dx>0.5) {l0++;dx-=1.0;}
		c=shape[l0];
		a=(shape[abs(l0-1)]-2*c+shape[l0+1])/2.0;
		b=(shape[l0+1]-shape[abs(l0-1)])/2.0;
		delta=c+dx*(a*dx+b);
		ycal[i]+=delta*par[6];
	    }
	    XPlot2Arrays(20,10,760,400,1,x,y,ycal,nData);
	    // ....................................................... menu
	    iComm=XMenu(0,titolo[1],ParName,ParVal,nPar,20,comandi,
			2,30,420);
	    if (iComm==1) break;
	    for (i=0;i<nPar;i++) par[i]=atof(ParVal[i]);
	}
	// ................................................. fissa i limiti
	XSetForeground(xvd.display,xvd.gc,xvd.white.pixel);
	XFillRectangle(xvd.display,xvd.backpix,xvd.gc,0,0,
		       xvd.width,xvd.height);
	for (i=0;i<nPar;i++)
	{
	    sprintf(ParVal[2*i],"%.5lf",0.6*par[i]);
	    sprintf(ParVal[1+2*i],"%.5lf",1.4*par[i]);
	}
	XMenu(1,"Fissa i limiti",LimName,ParVal,2*nPar,20,comandi,2,30,30);
	for (i=0;i<nPar;i++)
	{
	    ParMin[i]=atof(ParVal[2*i]);
	    ParMax[i]=atof(ParVal[1+2*i]);
	    if (ParMin[i]==ParMax[i]) ParFix[i]=TRUE;
	}
	// ............................................... fit di Marquardt
	MarqLoop(x,y,ycal,eps,nData,par,ParMin,ParMax,ParFix,nPar,
		 ShapeDer,tolerance);
	// ................................................................
	XPlot2Arrays(20,10,760,400,1,x,y,ycal,nData);
	for (i=0;i<nPar;i++) XPrintf(20+400*(i/4),420+20*(i%4),1,
				     "%15s%15.5lf",ParName[i],par[i]);
	XSync(xvd.display,TRUE);
	comandi[0]="restart";
	iComm=XMenu(1,0,0,0,0,0,comandi,2,280,520);
	if (iComm==1) break;
    }
    // ......................................................... free memory
    delete [] ycal;
    delete [] y;
    delete [] x;
    delete [] shape;
    delete [] ParMin;
    delete [] ParMax;
    delete [] par;
    delete [] eps;
    delete [] ParFix;
    CloseXWindow();
}