/* ELQ.c */
/* Eigenvalues for a tilted washboard potential */
/* by the Jacobi method */ 
/* S_Jacobi() by K. Minemura http://www-in.aut.ac.jp/~minemura/pub/Csimu/C/Jacobi.html */
/* 2012.2.23, 2020.6.13 by M. Suzuki */
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>

#define	TOL	1.0e-10   
#define	N	201
#define	MAX	500000
#define	xL	-0.6
#define	xH	1.0
#define DISPLAYMAX	5
#define EIGEN_OUT_NUMBER 10


void SetMatrix(int n, double a[n][n], double Ic, double C, double alpha)
{
	int i, j;
	double hbar=1.05457266e-34;		// [Js]
	double e=1.60217733e-19;		// [C]
	double EJ, EC, x, h;
	double potential();
	char buffer[100];
	
	EJ=(hbar/2/e)*Ic;	// [J]
	EJ/=2*M_PI*hbar;	// [Hz]
	EC=2*e*e/C;			// [J]
	EC/=2*M_PI*hbar;	// [Hz]
	h=(xH-xL)/(n-1);
	
	for(i=0;i<n;i++)
	{
		for(j=0;j<n;j++)
		{
			a[i][j]=0;
		}
	}
	for(i=1;i<n-1;i++)
	{
		x=xL+h*i;
		a[i][i]=2.0*EC/(h*h)+potential(x, EJ, alpha);
		a[i][i+1]=-EC/(h*h);
		a[i][i-1]=-EC/(h*h);
	}
	a[0][0]=2.0*EC/(h*h)+potential(xL, EJ, alpha);
	a[0][1]=-EC/(h*h);
	a[n-1][n-1]=2.0*EC/(h*h)+potential(xH, EJ, alpha);
	a[n-1][n-2]=-EC/(h*h);
}


double potential(double x, double EJ, double alpha)
{
				
	double y, z, U, U0;
	z=asin(alpha);
	U0=-EJ*(alpha*z+cos(z));
	y=x+z;
	U=-EJ*(alpha*y+cos(y))-U0;
	return U;
}


int S_Jacobi(int n, double a[n][n], double x[n][n])
{
	char buffer[100];
    int    i, j, k, m, count, status;
    double amax, amax0, theta, co, si, co2, si2, cosi, pi=M_PI;
    double aii, aij, ajj, aik, ajk;
    /* set initial values */
    for(i=0; i<n; i++){
		for(j=0; j<n; j++){
            if( i == j )  x[i][j] = 1.0; else  x[i][j] = 0.0;
        }
    }

    /* repeat Jacobi similarity transformation */
	count=0;
	status=9;
	while(count <= MAX)
	{
		/* seek for the maximum off-diagonal element */
        amax=0.0;
		for(k=0; k<n-1; k++)
		{
			for(m=k+1; m<n; m++)
			{
                amax0=fabs(a[k][m]);
				if(amax0 > amax)
				{
					i=k;
					j=m;
					amax=amax0;
				}
            }
		}
        /* judge the convergence */
		if(amax <= TOL)
		{
			status=0;
			break;
		}
		else
		{
			aii=a[i][i];
			aij=a[i][j];
			ajj=a[j][j];
			
            /* calculate the angle of rotation */
			if(fabs(aii-ajj) < TOL)
			{
                theta=0.25*pi*aij/fabs(aij);
			}
			else
			{
                theta=0.5*atan(2.0*aij/(aii-ajj));
            }
			co=cos(theta);
			si=sin(theta);
			co2=co*co;
			si2=si*si;
			cosi=co*si;

            /* similarity transformation matrix */
            a[i][i]=co2*aii+2.0*cosi*aij+si2*ajj;
            a[j][j]=si2*aii-2.0*cosi*aij+co2*ajj;
			a[i][j]=0.0;
			a[j][i]=0.0;
			for(k=0; k<n; k++)
			{
				if(k!=i && k!=j)
				{
					aik=a[k][i];
					ajk=a[k][j];
					a[k][i]=co*aik+si*ajk;
					a[i][k]=a[k][i];
					a[k][j]=-si*aik+co*ajk;
					a[j][k]=a[k][j];
                }
            }

            /* eigen vectors */
			for(k=0; k<n; k++)
			{
				aik=x[k][i];
				ajk=x[k][j];
                x[k][i]=co*aik+si*ajk;
                x[k][j]=-si*aik+co*ajk;
            }
            count++;
		}
		if(count % 100 ==0)
		{
			printf("\r S_Jacobi > iteration=%d",count);
			fflush(stdout);
 		}
    }
	printf("\n");
    return status;
}


void S_OutMat(int m, int n, double a[m][n])
{
	int i,j, k;
	if(n>DISPLAYMAX) k=DISPLAYMAX;
	else k=n;
	for(i=0;i<k;i++)
	{
		for(j=0; j<k; j++)
		{
			printf("  %10.6le", a[i][j]);
		}
		printf("\n");
   }
}

int main(int argc, char *argv[])
{
	FILE *fp;
	char filenameout[200];
	char buffer[100];
	int    i, j, status, k, kk, mn;
	double hbar=1.05457266e-34;		// [Js]
	double e=1.60217733e-19;		// [C]
	double a[N][N], x[N][N], d[N], dx, p;   
	double Ic, C, alpha, EJ;


	if(argc<4)
	{
		printf("usage ELQ 'Ic[uA]' 'C[fF]' 'alpha'\n");
		exit(0);
	}
	
	if(N>DISPLAYMAX) kk=DISPLAYMAX;
	else kk=N;
	
	dx=(xH-xL)/N;
	Ic=atof(argv[1]);
	Ic*=1e-6;
	printf("Ic %le\n",Ic);
	C=atof(argv[2]);
	C*=1e-15;
	printf("C %le\n",C);
	alpha=atof(argv[3]);
	printf("alpha %lf\n",alpha);

	EJ=(0.5*hbar/e)*Ic;	// [J]
	EJ/=2*M_PI*hbar;	// [Hz]
	
	strcpy(filenameout, "ELQ_Ic");
	strcat(filenameout, argv[1]);
	strcat(filenameout, "C");
	strcat(filenameout, argv[2]);
	strcat(filenameout, "a");
	strcat(filenameout, argv[3]);
	strcat(filenameout, "d");
	sprintf(buffer, "%d", N);
	strcat(filenameout, buffer);
	strcat(filenameout, ".txt");
	printf("Output filename  %s\n", filenameout);
	SetMatrix(N, a, Ic, C, alpha);
	
/*		if(N>DISPLAYMAX) k=DISPLAYMAX;
		else k=N;
	for(i=0;i<N;i++)
	{
		printf("%d ",i);
		for(j=i-1;j<i+2;j++)
		{
			printf("%le ",a[i][j]);
		}
			printf("\n");
	}*/
	
	
    printf("matix to be diagonalized\n");
    S_OutMat(N, N, a);
	
	printf("carry out the Jacobi method for engenvalue problems\n");
	status=S_Jacobi(N, a, x);
	
	
	/* sort in order of increasing eigenvalue */
	for(i=0;i<N;i++)
	{
		d[i]=a[i][i];	/* diagonal components */
	}
	
	/* exchanges eigenvalues, together with eigenvectors, i.e. columns of x,  */
	mn=0;	/* number of negative eigenvalues */
	for(i=0;i<N-mn;i++)
	{
		if(d[i]<0 && i!=N-mn-1)
		{
			p=d[k=i];
			while(d[N-mn-1]<0) mn++;
			d[k]=d[N-mn-1];
			d[N-mn-1]=p;
			for(j=0;j<N;j++)
			{
				p=x[j][k];
				x[j][k]=x[j][N-mn-1];
				x[j][N-mn-1]=p;
			}
			mn++;
		}
	}
	for(i=0;i<N-mn;i++)	/* positive eigenvalues */
	{
		p=d[k=i];
		for(j=i+1;j<N-mn;j++)
		{
			if(d[j]<=p) p=d[k=j];	/* seeking the smallest component */
		}
		if(k!=i)
		{
			d[k]=d[i];
			d[i]=p;
			for(j=0;j<N;j++)
			{
				p=x[j][i];
				x[j][i]=x[j][k];
				x[j][k]=p;
			}
		}
	}
	for(i=N-mn;i<N;i++)	/* negative eigenvalues */
	{
		p=d[k=i];
		for(j=i+1;j<N;j++)
		{
			if(d[j]<=p) p=d[k=j];	/* seeking the smallest component */
		}
		if(k!=i)
		{
			d[k]=d[i];
			d[i]=p;
			for(j=0;j<N;j++)
			{
				p=x[j][i];
				x[j][i]=x[j][k];
				x[j][k]=p;
			}
		}
	}
	
	fp=fopen(filenameout, "w");
	for(i=0;i<N;i++)
	{
		fprintf(fp, "%lf\t%le", xL+i*dx, d[i]);
		for(j=0;j<EIGEN_OUT_NUMBER;j++) fprintf(fp, "\t%lf", x[i][j]);
		fprintf(fp, "\n");
	}
	fclose(fp);
	
	if(status==0)
	{
		printf("eigenvalues \n");
		for(i=0; i<kk; i++) printf("%10.6le ", d[i]); 
		printf("\n");
		printf("eigenvectors\n");
		S_OutMat(N, N, x);
	}
}