#include <math.h>
#include <stdio.h>
#include <iostream.h>
#include "sgeev.cc"
#include "dgemm.cc"
#define DBG 0
/* CODING: 
 need to: 
   0. set up chebychev spacing grid Te[i]
   1. set up the matrixes for De=>DD=>A.
   2. solve for eigenvalues, and eigenvectors.
									*/
double De[130][130],Dee[130][130],DD[130][130];
float A[130][130];
double *Te,te[521];

double so=1.;


void cheby(double a,double b,int n, double ti[])
{
  int i,j;
  for (i=0;i<=n;i++) ti[i]=a+(1+cos(M_PI*(i*1./n-1.)))*(b-a)/2;
}							


void make_Te(int dim) 
{ 
  cheby(-so,so,dim+1,te);
  Te=&(te[1]);
}


void make_D(int);    	/*  makes the spectral deriv matrix  */ 

void make_DD(int dim) 	/*  makes the spec 2nd deriv matrix  */
{
  int i;
  make_D(dim);
  dgemmSimple dgemmy(dim);
  dgemmy.mul(De,Dee,DD);
}

/* a is the chopped off (by 2 dim's) DD matrix */
void make_A(int dim) 
{
  int i,j;
  make_DD(dim+2);
  for (i=0;i<dim;i++)
    for (j=0;j<dim;j++) 
      A[i][j]=-DD[i+1][j+1]/2;
}

void make_D(int dim)	/* makes the two idendical matrixes De=Dee */
{			
 int i,j,k;	
 double x;     
 for (i=0;i<dim;i++)
   {for (j=0;j<dim;j++) 
     {if(i!=j)
      { 
        x=1.;
	for (k=0;k<dim;k++)  
	  {
	    if (k!=j)
	      {
	        if (k!=i) x*=(te[i]-te[k]);
	        x/=(te[j]-te[k]);
	      }
	  }
	De[j][i]=x;   Dee[j][i]=x; 
      }
   }}
 for(i=0;i<dim;i++) 
  {
     for (x=0.,k=0;k<dim;k++) if (k!=i) x+=1./(te[i]-te[k]);
     De[i][i]=x; Dee[i][i]=x;
  }
}
/* note: fortran stores the matrixes in reverse notation, relative to C 
   notation. 
*/


void setup(int dim) 
{
 make_Te(dim);
 make_A(dim); 
}


void modify_A(int dim)   /* modifies A from D^2 to D^2 -x x Y, confined QHO */ 
{
  int i;
  for (i=0;i<dim;i++) A[i][i]+=Te[i]*Te[i]/2.;
}


main() 
{
 float f;
 FILE * F;
 FILE * file[3];
 int i,j,k,l,dim;
 double a,b,c;
 sgeevClass sgeevy(0);
 F=fopen("dat/es.dat","w");
 file[0]=fopen("dat/0.dat","w");
 file[1]=fopen("dat/1.dat","w");
 file[2]=fopen("dat/2.dat","w");

// PART I. COMPUTE EIGENFNS FOR USER.

 puts("s_well=a_well/a_classcalturn?\n");
 scanf("%lf",&so);
 
   puts("dim?  dim=N-1");
   scanf("%d",&dim); 
   
   setup(dim); 
   modify_A(dim);
   sgeevy.set(dim);
   sgeevy.get_eigensys(A);
   sgeevy.get_order();
   for (j=0;j<3;j++)   /* print out wavefn info into files: */
   {
   sgeevy.select(j);
   printf("\neigenvalue %22.20f    (+i %22.20f)\n",sgeevy.sel_val.re,sgeevy.sel_val.im); 
   for (i=0;i<dim;i++)
       fprintf(file[j],"%22.20g    %22.20f\n",Te[i], sgeevy.sel_vec[i].re);
   }
   fclose(file[0]); fclose(file[1]); fclose(file[2]);
   
// PART II. PLOT E VS INF.WELL.POS, LET so RANGE FROM .5 TO 3   

   dim=16;
   for (i=0;i<50;i++)
   {
     so=.5+2*i*2.5/100;
     setup(dim); 
     modify_A(dim);
     sgeevy.set(dim);
     sgeevy.get_eigensys(A,0);
     sgeevy.get_order();
     sgeevy.select(0); 
     fprintf(F,"%22.20g    %22.20g\n",  so, sgeevy.sel_val.re);
   }
   
   return 1; 
}