/*
2013.3.6 for a sample movie to show how to make movies from gnuplot outputs 
2013.3.6 values of wavefunctions of harmonic oscillator eigenstates are stored
2013.1.24-25 modifed by N.T.
extensions and modifications by R.D.
2012.12.17-19 modifed by N.T.
2012.12.17(mon)
>>>2013.2.6(wed) R.D.
調和振動子シュレディンガー方程式
dvm/dt = sigma(hmn/i*hbar)vn

*/

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

# define nmax 300     // 量子数nの最大値
# define N (nmax+1)   // 基底に含める（量子数n=0..nmaxをもつ）基底の個数
# define maxpower 6   // xのべきの最大値

const double hbar = 1.0;
const int ki = 2*maxpower;
double ol_b; // b (oscillator length of the basis)

int runge(double tmax,int itmax,int dit_obs, int dit_vec,int dit_wav,int cooling,double complex *v);
int set_hamiltonian_matrix(int option);
double deriv_xi_2 (int i, int j);
double xi_0 (int i, int j);
double xi_1 (int i, int j);
double xi_2 (int i, int j);
double xi_3 (int i, int j);
double xi_4 (int i, int j);
double xi_5 (int i, int j);
double xi_6 (int i, int j);
double complex timederivative( int i, double t, double complex *x);
double norm_square(double complex *x);
int normalize(double complex *x);
double harmonic_oscillator_wavefunction(int n, double x);
int write_vector(double time, double complex *v);
int write_wavefunction(double time, double complex *v);
int write_potential(double q0, double q1, double q2, double q3, double q4, double q5, double q6);
double deriv_xi_2 (int i, int j);
double xi_2 (int i, int j);
double energy_expectation_value( double complex *x);
double energy_2_expectation_value( double complex *x);
double x_expectation_value( double complex *x);
double x_2_expectation_value( double complex *x);

int main(){
   main2();
   return 0;
}

int main2(){
   double complex v[N]; // 状態ベクトル
   double tmax; // 時間発展を行う時間の長さ
   int itmax;  // 時間の刻みの総数（0番..itmax番の刻みがある)
   int dit_obs=0; // observablesを出力する時間ステップ番号の間隔
   int dit_vec=0; // 状態ベクトルを出力する時間ステップ番号の間隔
   int dit_wav=0; // 波動関数をファイル出力する時間ステップ番号の間隔
    // 注：dit_vec,dit_wavともにゼロ以下の値を設定すると出力が完全に抑止される。
   double omega = 0.5; // 注意：同一用途の同名変数の値を２箇所で定義している。
   double beta = 0.1; // 初期状態設定用パラメータ
                      // （温度の逆数、単位＝1/基底準位間隔）
   int cooling = 1; // =0で実時間発展、=1で虚時間発展
   int i; // 作業用変数(カウンター)

   /* 初期状態の設定 */

   for ( i=0; i<N; i++) v[i] = 0.0;
   v[0]=1.0;

   /* 虚時間発展 : 何らかのハミルトニアンの基底状態として初期状態を作る */
   /*
   set_hamiltonian_matrix(1);
   tmax = 10.0;
   itmax=10000000;
   dit_obs = 0;
   dit_vec = 0;
   dit_wav = 0;
   cooling = 1;
   runge(tmax, itmax, dit_obs, dit_vec, dit_wav, cooling, v);
   */
   /* 実時間発展 */

   set_hamiltonian_matrix(2) ;
   tmax =  8.0;
   itmax= 3000;
   dit_obs = 0;
   dit_vec = 0;
   dit_wav =10;
   cooling = 0;
   runge(tmax, itmax, dit_obs, dit_vec, dit_wav, cooling, v);

   return 0;
}

/* ルンゲ・クッタ法による時間発展 */

int runge(double tmax, int itmax, int dit_obs, int dit_vec, int dit_wav, int cooling, double complex *v){
   int nprint=6; // 出力する状態ベクトルの成分の個数（n=0..nprint-1を出力）

   double complex vtmp[N], k0[N], k1[N], k2[N], k3[N];
   double t,dt; // 時刻,時刻の刻み幅
   double complex tfct; // 実時間発展ではdt,虚時間発展では-I*dt
   double norm2; // 状態ベクトルのノルムの２乗
   double fct; // 作業変数（因子）
   double E; // エネルギー期待値E
   double E_2; // Eの2乗の期待値
   double X; // xの期待値
   double X_2; // xの2乗の期待値
   double delta_E; // ΔE
   double delta_x; // Δx
   int itime; // 時刻の刻みの番号=0番..itmax番
   int i; // 作業変数（カウンター)

   dt = tmax/itmax;
   if(nprint>N)nprint=N;

   /* t=0 */
   t = 0.0;              
   if(dit_obs > 0){  // t=0でのobservablesの出力
     E = energy_expectation_value(v);
     E_2 = energy_2_expectation_value(v);
     X = x_expectation_value(v);
     X_2 = x_2_expectation_value(v);
     delta_E = sqrt(fabs(E_2 - E*E));
     delta_x = sqrt(fabs(X_2 - X*X));
     printf("# t, E, delta_E, X, delta_x\n");
     printf("%f %.16f %.16f %.16f %.16f :obs\n", t, E, delta_E, X, delta_x);
   }

   if(dit_vec > 0) write_vector(t,v); // t=0でのvectorの出力

   if(dit_wav > 0) write_wavefunction(t,v);

   /* t>0 */
   if(cooling==0){ tfct=dt; } else { tfct=-dt*I; }
   for ( itime=1; itime<=itmax; itime++){

      t = dt*itime;

      for ( i=0; i<N; i++) k0[i] = timederivative( i, t, v)*tfct;

      for ( i=0; i<N; i++) vtmp[i]=v[i] + k0[i]/2.0;
      for ( i=0; i<N; i++) k1[i] = timederivative( i, t + dt/2.0, vtmp)*tfct;

      for ( i=0; i<N; i++) vtmp[i]=v[i] + k1[i]/2.0;
      for ( i=0; i<N; i++) k2[i] = timederivative( i, t + dt/2.0, vtmp)*tfct;

      for ( i=0; i<N; i++) vtmp[i]=v[i] + k2[i];
      for ( i=0; i<N; i++) k3[i] = timederivative( i, t + dt, vtmp)*tfct;
       
      for ( i=0; i<N; i++) v[i] += ( k0[i] + 2.0*k1[i] + 2.0*k2[i] + k3[i])/6.0;


      if(dit_obs > 0 && itime % dit_obs == 0){
        E = energy_expectation_value(v);
        E_2 = energy_2_expectation_value(v);
        X = x_expectation_value(v);
        X_2 = x_2_expectation_value(v);
        delta_E = sqrt(fabs(E_2 - E*E));
        delta_x = sqrt(fabs(X_2 - X*X));

        printf("%f %.16f %.16f %.16f %.16f :obs\n", t, E, delta_E, X, delta_x);
      }

      normalize(v); // 数値計算誤差によるノルムの１からのずれを修正する

      if(dit_vec > 0 && itime % dit_vec == 0) write_vector(t,v);

      if(dit_wav > 0 && itime % dit_wav == 0 ) write_wavefunction(t,v);

   }
   if(dit_wav < 0 ) write_wavefunction(t,v);
   return 0;
}

/* hamiltonianの計算 */

// double ham1[N][N];
double *ham1 = NULL;

int set_hamiltonian_matrix(int option){
   
   double mass_b = 1.0;
   double omega_b = 1.0;
   double mass = 1.0;
   double omega = 0.5;
   double x0 = 30.0;
   double x1 = 30.0;
   double hight1 = 1000.0;
   double hight2 = 2.0;
   double p, q1, q2, q3, q4, q5, q6, q0;
   int i, j;
   double x0_pot=7.0; // location of minimum
   double h_pot=5.0; // barrier height
   double d_pot=0.0; // 0.25; // splitting of the depth of two minima

   if(ham1 == NULL){
     ham1=malloc(N*(2*maxpower + 1)*sizeof(double));
     if(ham1 == NULL) {fprintf(stderr,"malloc:%d\n",N);exit(1);}
   }

   ol_b = sqrt(hbar/(mass_b*omega_b));

//　虚時間発展のハミルトニアンの各項の係数(coolingにより波束をつくる)
   if( option == 1 ){
      p = -hbar*hbar/(2.0*mass_b*ol_b*ol_b); // 2階微分の係数
      q1 = 0.0; // 1次の項の係数
      q2 = mass*omega*omega*ol_b*ol_b/2.0; // 2次の項の係数
      q3 = 0.0; // 3次の項の係数
      q4 = 0.0; // 4次の項の係数
      q5 = 0.0; // 5次の項の係数
      q6 = 0.0; // mass*omega_b*omega_b*ol_b*ol_b*ol_b*ol_b*ol_b*ol_b/2.0; // 6次の項の係数
      q0 = 0.0; // mass*omega_b*omega_b*x0*x0*x0*x0*x0*x0/2.0; // 定数項
   }
   else if( abs(option) == 2 ){
//　実時間発展のハミルトニアン(準備した波束を振動させる)
     double bxpow1=ol_b/x0_pot,bxpow=bxpow1;
      p = -hbar*hbar/(2.0*mass*ol_b*ol_b);
      q0 = 0.0;
      q1 = 0.0;
      q2 = 0.0;
      q3 = 0.0;
      q4 = 0.0;
      q5 = 0.0;
      q6 = h_pot/pow(2*x0_pot,6.0);
      if(option < 0){
//　虚数時間発展のハミルトニアン(波束を準備する)
        double vspp=h_pot/pow(2*x0_pot,6.0);
        q0 += x0_pot*x0_pot*x0_pot*x0_pot*x0_pot*x0_pot*vspp;
        q1 -= 6.0*x0_pot*x0_pot*x0_pot*x0_pot*x0_pot*vspp;
        q2 += 15.0*x0_pot*x0_pot*x0_pot*x0_pot*vspp;
        q3 -= 20.0*x0_pot*x0_pot*x0_pot*vspp;
        q4 += 15.0*x0_pot*x0_pot*vspp;
        q5 = -6.0*x0_pot*vspp;
        q6 = vspp;
      }
   }

   for( i=0; i<N; i++){

      int js = i-maxpower;  // jstart
      if (js < 0) js = 0;

      int je = i+maxpower;  // jend
      if (je > nmax) je = nmax;

      for ( j=js; j<=je; j++){
         int k = ki*i + j;
         ham1[k] = p*deriv_xi_2(i, j)+q0*xi_0(i, j)+q1*xi_1 (i, j)+q2*xi_2(i, j)
                   +q3*xi_3(i, j)+q4*xi_4(i, j)+q5*xi_5(i, j)+q6*xi_6(i, j);
      }
   }

   printf("# %f =mass\n",mass);
   printf("# %f =ol_b\n",ol_b);
   write_potential(q0, q1, q2, q3, q4, q5, q6);
   return 0;
}


// <i|(d^2/dxi^2)|j> 部分の計算
double deriv_xi_2 (int i, int j){

   if( i==(j-2) ){
      return sqrt(j*j-j)*0.5;
   }
   else if( i==j ){
      return -j-0.5;
   }
   else if( i==(j+2)){
      return sqrt((j+1.0)*(j+2.0))*0.5;
   }
   else {
      return 0;
   }
}

// <i|j> 部分の計算
double xi_0 (int i, int j){

   if( i==j ){
      return 1;
   }
   else {
      return 0;
   }
}

// <i|(xi^1)|j> 部分の計算
double xi_1 (int i, int j){

   if( i==(j-1) ){
      return sqrt(j*0.5);
   }
   else if( i==(j+1)){
      return sqrt((j+1.0)*0.5);
   }
   else {
      return 0;
   }
}

// <i|(xi^2)|j> 部分の計算
double xi_2 (int i, int j){

   if( i==(j-2) ){
      return sqrt(j*j-j)*0.5;
   }
   else if( i==j ){
      return j+0.5;
   }
   else if( i==(j+2)){
      return sqrt((j+1.0)*(j+2.0))*0.5;
   }
   else {
      return 0;
   }
}

// <i|(xi^3)|j> 部分の計算
double xi_3 (int i, int j){

   if( i==(j-3) ){
      return sqrt(j*(j-1.0)*(j-2.0)*0.5)*0.5;
   }
   else if( i==(j-1) ){
      return 3.0*j*0.5*sqrt(j*0.5);
   }
   else if( i==(j+1) ){
      return 3.0*(j+1.0)*0.5*sqrt((j+1.0)*0.5);
   }
   else if( i==(j+3) ){
      return 0.5*sqrt((j+1.0)*(j+2.0)*(j+3.0)*0.5);
   }
   else {
      return 0;
   }
}

// <i|(xi^4)|j> 部分の計算
double xi_4 (int i, int j){

   if( i==(j-4) ){
      return sqrt(j*(j-1.0)*(j-2.0)*(j-3.0))*0.25;
   }
   else if( i==(j-2) ){
      return (2.0*j-1.0)*0.5*sqrt(j*(j-1.0));
   }
   else if( i==j ){
      return 3.0*0.25*(j*j+(j+1.0)*(j+1.0));
   }
   else if( i==(j+2) ){
      return (2.0*j+3.0)*0.5*sqrt((j+1.0)*(j+2.0));
   }
   else if( i==(j+4) ){
      return sqrt((j+1.0)*(j+2.0)*(j+3.0)*(j+4.0))*0.25;
   }
   else {
      return 0;
   }
}

// <i|(xi^5)|j> 部分の計算
double xi_5 (int i, int j){

   if( i==(j-5) ){
      return sqrt(j*(j-1.0)*(j-2.0)*(j-3.0)*(j-4.0)/2.0)*0.25;
   }
   else if( i==(j-3) ){
      return (5.0*j-5.0)*sqrt(j*(j-1.0)*(j-2.0)/2.0)*0.25;
   }
   else if( i==(j-1) ){
      return (10.0*j*j+5.0)*sqrt(j/2.0)*0.25;
   }
   else if( i==(j+1) ){
      return (10.0*j*j+20.0*j+15.0)*sqrt((j+1.0)/2.0)*0.25;
   }
   else if( i==(j+3) ){
      return (5.0*j+10.0)*sqrt((j+1.0)*(j+2.0)*(j+3.0)/2.0)*0.25;
   }
   else if( i==(j+5) ){
      return sqrt((j+1.0)*(j+2.0)*(j+3.0)*(j+4.0)*(j+5.0)/2.0)*0.25;
   }
   else {
      return 0;
   }
}

// <i|(xi^6)|j> 部分の計算
double xi_6 (int i, int j){

   if( i==(j-6) ){
      return sqrt(j*(j-1.0)*(j-2.0)*(j-3.0)*(j-4.0)*(j-5.0))/8.0;
   }
   else if( i==(j-4) ){
      return (6.0*j-9.0)*sqrt(j*(j-1.0)*(j-2.0)*(j-3.0))/8.0;
   }
   else if( i==(j-2) ){
      return (15.0*j*j-15.0*j+15.0)*sqrt(j*(j-1.0))/8.0;
   }
   else if( i==j ){
      return (20.0*j*j*j+30.0*j*j+40.0*j+15.0)/8.0;
   }
   else if( i==(j+2) ){
      return (15.0*j*j+45.0*j+45.0)*sqrt((j+1.0)*(j+2.0))/8.0;
   }
   else if( i==(j+4) ){
      return (6.0*j+15.0)*sqrt((j+1.0)*(j+2.0)*(j+3.0)*(j+4.0))/8.0;
   }
   else if( i==(j+6) ){
      return sqrt((j+1.0)*(j+2.0)*(j+3.0)*(j+4.0)*(j+5.0)*(j+6.0))/8.0;
   }
   else {
      return 0;
   }
}

/* 右辺 = sigma(hmn/i*hbar)*vn の計算 */

double complex timederivative( int i, double t, double complex *x){

   double complex s;
   int j, js, je;

   js = i-maxpower;
   if (js < 0) js = 0;
   je = i+maxpower;
   if (je > nmax) je = nmax;
   s = 0;
   for ( j=js; j<=je; j++){

      s += ham1[ki*i+j]*x[j];

   }
   return -I/hbar*s;
}

double norm_square(double complex *x){ // 状態ベクトルのノルムの２乗を計算
   int i;
   double s;
   s = 0.0;
   for (i=0; i<N; i++){
     s += creal(x[i])*creal(x[i]) + cimag(x[i])*cimag(x[i]);
   }
   return s;
}

int normalize(double complex *x){ // 状態ベクトルを規格化
   int i;
   double fct;
   fct = 1.0/sqrt(norm_square(x));
   for ( i=0; i<N; i++){
      x[i] *= fct;
   }
   return 0;
}

double harmonic_oscillator_wavefunction(int n, double x){
  // y/m/d=13/1/17
  //    exp(-x^2/2) H_n(x)/sqrt(sqrt(pi) n^2 n!)
  // using H_n(x) = 2*x*H_{n-1}(x) -2*(n-1)*H_{n-2}(x)
  double a,b,c;
  double fnorm; // normalization factor
  double y; // retuned value psi_n(x) 
            // = value of the normalized n-node eigenwavefunction at x
  int k, k1=0, k2=0;
  double x2; // x*x
  const double fct=1.0e50,fctreci=1.0/fct;

  x2=x*x; 
  if(n>=2){
    b=1.0;
    c=2.0*x;
    for(k=1;k<n;k++){
      a=b;
      b=c;
      c=2*(x*b-k*a);
      if(fabs(c) > fct){
        k1++; a*=fctreci; b*=fctreci; c*=fctreci;
      }
    }
    fnorm=1;
    for(k=1;k<=n;k++){
      fnorm*=2*k;
      if(fabs(fnorm) > fct){
        k2++; fnorm*=fctreci;
      }
    }
    y=exp(-x2*0.5+(k1-k2*0.5)*log(fct))*c/sqrt(sqrt(M_PI)*fnorm);
  }
  else if(n==1) y=exp(-x2*0.5)*2*x/sqrt(sqrt(M_PI)*2.0);
  else if(n==0) y=exp(-x2*0.5)/sqrt(sqrt(M_PI));
  else { // if(n<0)
    printf("harmonic_oscillator_wavefunction:argument error %d %f\n",n,x);
    return 0.0;
  }
  return y;
}

int write_vector(double time, double complex *v){ 
// 状態ベクトル(v)を名前に番号が付いたファイルに出力する。
//   time: 時刻. 参考情報として出力ファイルに書き加える。
//   v   : 状態ベクトル（調和振動子基底での展開係数) 
  int i;
  static int filenumber = 0;
  double s;
  FILE *fo;
  char filename[256];

  if(filenumber > 999) return -1; // 出力ファイル数過大。ジョブ設定ミスの疑い。
  snprintf(filename,256,"vec%03d.dat",filenumber++);
  fo=fopen(filename,"w");
  fprintf(fo,"# time=%.5f\n",time);
  s=0.0;
  for(i=0;i<=nmax;i++){
    fprintf(fo,"%3d %.9f %.9f\n",i,creal(v[i]),cimag(v[i]));
    s+=creal(v[i])*creal(v[i])+cimag(v[i])*cimag(v[i]);
  }
  fprintf(fo,"# squared norm=%.10f\n",s);
  fclose(fo);
  return 0;
}

int write_wavefunction(double time, double complex *v){ 
// 波動関数(x)を名前に番号が付いたファイルに出力する。
//   time: 時刻. 参考情報として出力ファイルに書き加える。
//   v   : 状態ベクトル（調和振動子基底での展開係数) 
  const int npsi=150; // 500; // 波動関数の値を出力するx軸上の点の個数-1 (0番..npsi番の点)
  double xmin=-15.0; // 80.0; // 波動関数の値を出力する左端の点の座標
  double xmax= 15.0; // 80.0; // 波動関数の値を出力する右端の点の座標
  int i,j,n;
  static int filenumber = 0;
  double complex s;
  double norm2; // 状態のノルムの２乗（x表示で求めたもの)
  double x,dx;
  FILE *fo;
  char filename[256];
  static double *howf = NULL; // store values of wavefunctions of harmonic oscillator eigenstates

  if(filenumber > 999) return -1; // 出力ファイル数過大。ジョブ設定ミスの疑い。
  dx=(xmax-xmin)/npsi;
  if(howf == NULL){
    howf=(double *)malloc((nmax+1)*(npsi+1)*sizeof(double));
    for(j=i=0;i<=npsi;i++) {
      x=xmin+dx*i;
      for(n=0;n<=nmax;n++){
        howf[j]=harmonic_oscillator_wavefunction(n,x/ol_b);
        j++;
      }
    }
  }
  snprintf(filename,256,"wf%03d.dat",filenumber++);
  fo=fopen(filename,"w");
  fprintf(fo,"# time=%.5f\n",time);
  norm2=0.0;
  for(j=i=0;i<=npsi;i++) {
    x=xmin+dx*i;
    s=0.0;
    for(n=0;n<=nmax;n++){
      s+=v[n]*howf[j];
      j++;
    }
    s/=sqrt(ol_b);
    fprintf(fo,"%.3f %.6f %.6f\n",x,creal(s),cimag(s));
    norm2+=creal(s)*creal(s)+cimag(s)*cimag(s);
  }
  fprintf(fo,"# squared norm=%.10f\n",norm2*dx);
  fclose(fo);
  return 0;
}

int write_potential(double q0, double q1, double q2, double q3, double q4, double q5, double q6){ 
// potential energy を名前に番号が付いたファイルに出力する。
  const int npsi=500; // potentialの値を出力するx軸上の点の個数-1 (0番..npsi番の点)
  double xmin=-80.0; // potentialの値を出力する左端の点の座標
  double xmax= 80.0; // potentialの値を出力する右端の点の座標
  int i;
  static int filenumber = 0;
  double x,dx;
  FILE *fo;
  char filename[256];

  if(filenumber > 999) return -1; // 出力ファイル数過大。ジョブ設定ミスの疑い。
  snprintf(filename,256,"pot%03d.dat",filenumber++);
  fo=fopen(filename,"w");

  printf("# %f =q0\n",q0);
  printf("# %f =q1\n",q1);
  printf("# %f =q2\n",q2);
  printf("# %f =q3\n",q3);
  printf("# %f =q4\n",q4);
  printf("# %f =q5\n",q5);
  printf("# %f =q6\n",q6);

  dx=(xmax-xmin)/npsi;
  for(i=0;i<=npsi;i++) {
    x=xmin+dx*i;
    fprintf(fo,"%.3f %.6f\n",x,q0+x*(q1+x*(q2+x*(q3+x*(q4+x*(q5+x*q6))))));
  }
  fclose(fo);
  return 0;
}

// エネルギー期待値の計算
double energy_expectation_value( double complex *x){

   double complex s;
   int i, j, k, js, je;

   s = 0;
   for ( i=0; i<N; i++){   
      js = i-maxpower;
      if (js < 0) js = 0;
      je = i+maxpower;
      if (je > nmax) je = nmax;
      k = ki*i+j;
      for ( j=js; j<=je; j++){
          s += conj(x[i])*ham1[ki*i+j]*x[j];
      }
   }
   // sの虚部が0であることの確認
   if (fabs(cimag(s)) > 1.0e-20){
      printf("%f %f\n", creal(s), cimag(s));
   }

   return creal(s);
}

// エネルギーの2乗の期待値の計算
double energy_2_expectation_value( double complex *x){

   double complex s, s1;
   int i, j, k, js, je, ks, ke;

   s = 0.0;
   for ( i=0; i<N; i++){
       js = i-maxpower*2;
        if (js < 0) js = 0;
       je = i+maxpower*2;
        if (je > nmax) je = nmax;
      for ( j=js; j<=je; j++){
          if (i > j) ks = i-maxpower;
           else ks = j-maxpower;
          if (ks < 0) ks = 0;

          if (i > j) ke = j+maxpower;
           else ke = i+maxpower;
          if (ke > nmax) ke = nmax;

          s1 = 0.0;
         for (k=ks; k<=ke; k++){
            s1 += ham1[ki*i+k]*ham1[ki*k+j];
          }
         s += s1*conj(x[i])*x[j];
      }
   }
   // sの虚部が0であることの確認
   if (fabs(cimag(s)) > 1.0e-20){
      printf("%f %f\n", creal(s), cimag(s));
   }

   return creal(s);
}

// xの期待値の計算
double x_expectation_value( double complex *x){

   double complex s;
   int i, j;

   s = 0;
   for ( i=0; i<N; i++){   
      for ( j=0; j<N; j++){
          s += conj(x[i])*x[j]*ol_b*xi_1 (i, j);
      }
   }
   // sの虚部が0であることの確認
   if (fabs(cimag(s)) > 1.0e-20){
      printf("%f %f\n", creal(s), cimag(s));
   }

   return creal(s);
}

// xの2乗の期待値の計算
double x_2_expectation_value( double complex *x){

   double complex s;
   int i, j;

   s = 0;
   for ( i=0; i<N; i++){
      for ( j=0; j<N; j++){
          s += conj(x[i])*x[j]*ol_b*ol_b*xi_2 (i, j);
      }
   }
   // sの虚部が0であることの確認
   if (fabs(cimag(s)) > 1.0e-20){
      printf("%f %f\n", creal(s), cimag(s));
   }

   return creal(s);
}