/*----------------------------------------*/
/* ∞に高い壁に囲まれた粒子の波動関数     */
/*----------------------------------------*/

#include<stdio.h>
#include<math.h>

#define n      20      /*　grid数　            */ 
#define max_it 100000  /*　最大計算回数　　　  */
#define eps    1e-6   /*　誤差の閾値          */
#define Length 1.     /*  井戸の幅            */
void bsort( double u[], int index[]);
int main(void)
{
  int i,j,k,jj,kk,num;
  int jmax,kmax,imax,imax2,imax3;
  int index[n];
  double error, aa, theta, diag;
  double matrix_a[n][n];
  double matrix_a_new[n][n];
  double r[n][n], r_inv[n][n], r_stack[n][n], r_stack_new[n][n];
  double eigen_value[n];
  double delta_xi,dxi2;
  FILE *output_file, *output_file2, *output_file3;

  output_file=fopen("output.dat","w"); /* 出力ファイルオープン*/ 
  output_file2=fopen("output2.dat","w"); /* 出力ファイルオープン*/ 
  output_file3=fopen("output3.dat","w"); /* 出力ファイルオープン*/ 

  delta_xi = Length/(double)(n+1); dxi2 = delta_xi*delta_xi ;


  for( j = 0; j < n; j++){
    for( k = 0; k < n; k++){
      matrix_a[j][k]=0.;
    }
  }
  for( j = 0; j < n; j++){
      matrix_a[j][j]=2.;
      if( j != 0 ){ matrix_a[j][j-1] = -1; }
      if( j != n-1 ){ matrix_a[j][j+1] = -1; }
  }

  num=0; /* 繰り返し回数初期化 */

  for( j = 0; j < n; j++){
    for( k = 0; k < n; k++){
      if( j != k){
	r[j][k]=0.;  r_inv[j][k]=0.; r_stack[j][k]=0.;}
      else{
	r[j][k]=1.;  r_inv[j][k]=1.; r_stack[j][k]=1.;}
    }
  }

  for(i=0; i < max_it; i++){    /* i < max_it　である限りループ */

    num+=1; error=0;  /* numを1増やし、errorを初期化 */

    /* search for max component */
    aa=0.; jmax=1; kmax=0;
    for(j=0; j < n; j++){
      for(k=0; k < n; k++){
	if( k != j) {
	  if(aa < fabs(matrix_a[j][k])){ jmax=j; kmax=k; aa = fabs(matrix_a[j][k]);}
	}
      }
    }

    /* sum of the diagonal components */
    diag=0.;
    for(j=0; j < n; j++){
      diag += fabs(matrix_a[j][j]);
    }
    
    printf("error= %.16e \n", aa / diag);

    if( aa/diag < eps ) break; /* errorが十分小さければループを出る */



    /* set up kk,kl,lk,ll components of r */
    if( fabs( matrix_a[kmax][kmax]- matrix_a[jmax][jmax] ) < 1e-16 ){
      theta = atan(1.); }
    else{
      theta = 0.5 * atan( 2.* matrix_a[jmax][kmax]/(matrix_a[kmax][kmax]- matrix_a[jmax][jmax]));
    }

    for( j = 0; j < n; j++){
      for( k = 0; k < n; k++){
	if( j != k){
	  r[j][k]=0.;  r_inv[j][k]=0.;}
	else{
	  r[j][k]=1.;  r_inv[j][k]=1.;}
      }
    }

    r[jmax][jmax] =  cos(theta);
    r[jmax][kmax] =  sin(theta);
    r[kmax][jmax] = -sin(theta);
    r[kmax][kmax] =  cos(theta);

    r_inv[jmax][jmax] =  cos(theta);
    r_inv[jmax][kmax] = -sin(theta);
    r_inv[kmax][jmax] =  sin(theta);
    r_inv[kmax][kmax] =  cos(theta);



    /* calculate matrix_a_new */
    for(j=0; j < n; j++){
      for(k=0; k < n; k++){

	matrix_a_new[j][k]=0.;

	for(jj=0; jj < n; jj++){
	  for(kk=0; kk < n; kk++){
	    matrix_a_new[j][k] += r_inv[j][jj] * matrix_a[jj][kk] * r[kk][k] ;  
	  }
	}

      }
    }

    /* copy to matrix_a */
    for(j=0; j < n; j++){
      for(k=0; k < n; k++){
	    matrix_a[j][k] = matrix_a_new[j][k];  
      }
    }

   
    /* calculate r_stack_new */
    for(j=0; j < n; j++){
      for(k=0; k < n; k++){
	r_stack_new[j][k]=0.;
	for(jj=0; jj < n; jj++){
	  r_stack_new[j][k] += r_stack[j][jj] * r[jj][k];
	}
      }
    }
   
    
    /* copy to r_stack */
    for(j=0; j < n; j++){
       for(k=0; k < n; k++){
	 r_stack[j][k] = r_stack_new[j][k];
      }
    }

  }

  if( num == max_it){
  printf("収束しませんでした。条件をかえてください。\n");
  }
  else{
      printf("収束までにかかった繰り返し処理の数 = %d\n",num);
  }

  for(i=0; i < n ; i++){
    eigen_value[i]=matrix_a[i][i];
  }
  bsort(eigen_value,index);
  
  

  printf("e1=%f  e2=%f e3=%f\n",(float)eigen_value[index[0]]/dxi2,(float)eigen_value[index[1]]/dxi2,(float)eigen_value[index[2]]/dxi2);
  
  for(j=0; j < n ; j++){
    fprintf(output_file,"%f %f \n",(float)delta_xi*(j+1),(float)r_stack[j][index[0]]);
    fprintf(output_file2,"%f %f \n",(float)delta_xi*(j+1),(float)r_stack[j][index[1]]);
    fprintf(output_file3,"%f %f \n",(float)delta_xi*(j+1),(float)r_stack[j][index[2]]);

  }

  return (0);
}
void bsort( double u[], int index[]){
  int i,j,k;

  for(j=0; j < n ; j++){
    index[j]=j;
  }

  for(i=n-1; i > -1 ; i--){
    for(j=0; j < i ; j++){
      if(u[index[j]] > u[index[j+1]]){
	k=index[j+1]; index[j+1]=index[j]; index[j]=k;
      }
    }
  }    
  
  return;
}