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

#include "matrix.h"   // T.Nagano 2011 version

main()
{
    int i,j,steps,ks;
    MATRIX a,p,r;
    double x, pi=3.14159265359,e=1e-10;
    double g,sift,k;
    double U1;
    double Lat;
    FILE *fp1;

    
    SIZE=5;                              //行列のサイズ
    steps=10;                            //ヤコビ法のステップ

    printf("A :\n");  
    matprint(a);                         //行列の表示
    printf("\n");



    fp1=fopen("e-k.dat","w");            //計算結果を格納するファイル
   
    Lat=2.5;                            //格子定数


    g=2.0*pi/Lat;                           //逆格子
    U1=1.0;                             //ポテンシャルのフーリエ成分




    sift=(SIZE/2+1)-SIZE;               

//==========================================================================================
//     k点の繰り返し
//=========================
    for(k=-pi/Lat;k < pi/Lat; k=k+0.05){
      printf("k=  %lf\n",k);
      a=zeromat(a);


      for(i=0; i < SIZE; i++){
        a.ent[i][i]=(k+g*(i+sift))* (k+g*(i+sift)) ;         //対角要素
      }
      for(i=0; i < SIZE; i++){

//        for(j=i+1;j< SIZE;j++){
//          a.ent[i][j]=  U1;
//          a.ent[j][i]=a.ent[i][j];   
//        }
                                                             //非対角要素
          a.ent[i][i+1]=  U1;  
          a.ent[i+1][i]=a.ent[i][i+1];   

      }
      matprint(a);
      printf("\n");

//-------------------------------------------------------------------
//      ヤコビ法
//----------------------
      p=unitmat();
      for(ks=1;ks <= steps;ks++){
        for(i=0; i < SIZE; i++){
          for(j=i+1; j < SIZE; j++){
            // 回転角度設定
            if(fabs(a.ent[i][i]-a.ent[j][j])< e ){
                x=pi/4.0;
            }
            else{
                x=atan(2*a.ent[i][j]/(a.ent[i][i]-a.ent[j][j]))/2.0;
            }
            // 回転の行列
            r=rotmat(i,j,x);
            // 直交行列の更新
            p=matmul(p,r);
            // A の更新
            a=matmul(transmat(r),matmul(a,r));
          }
        }
      }
//-------------------------------------------------------------------



//-------------------------------------------------------------------
//    　結果出力
//-------------------------------------------------
      printf("%12.7f   ",k);
      fprintf(fp1,"%12.7f   ",k);
      for(i=0;i<SIZE;i++){
        printf("%12.7f   ",a.ent[i][i]);
        fprintf(fp1,"%12.7f   ",a.ent[i][i]);
      }
      printf("\n");
      fprintf(fp1,"\n");
      close(fp1);
//-------------------------------------------------------------------



      }
//==========================================================================================



}