/* ===================================================================== 
   井戸型ポテンシャルの固有値を計算しグラフィック出力 Jan 12 2003
   ===================================================================== */
#include <stdio.h>
#include <math.h>
#include "cpgplot.h"



/* ---------------------------------------------------------------------
   potential function V(x)
   --------------------------------------------------------------------- */
double V(double x, double V0, double a)
{
    double potential;
    if( fabs(x) > a )
        potential = V0;
    else
        potential = 0.0;
    return ( potential );
}



/* ---------------------------------------------------------------------
   parity check
   --------------------------------------------------------------------- */
void ParityCheck( int parity, double *phi, double *dphi )
{
  if(parity == -1)
    {
      *phi = 0;   // initial values at x = 0
      *dphi = 1;  // first derivative
    }
  else
    {
      *phi = 1;
      *dphi = 0;
    }
}



/* ---------------------------------------------------------------------
   Euler algorithm
   --------------------------------------------------------------------- */
void Euler( double V0, double a, double x, 
	    double dx ,double E ,double *dphi, double *phi )
{
  double d2phi;
  /* \frac{2m}{hbar^2} = 1 となる無次元の単位系で計算 */
  d2phi = (V( x, V0, a) - E) * (*phi);
  /* オイラー-クロマー・アルゴリズム */
  *dphi += d2phi * dx;         
  *phi  += *dphi  * dx;
}



/* ---------------------------------------------------------------------
   井戸型ポテンシャルの描画 (青線)
   --------------------------------------------------------------------- */
void DrawPotential( double a, double V0, double xmax )
{
  cpgsci(4); cpgslw(1);
  cpgmove( -xmax, V0 );
  cpgdraw( -a, V0 );
  cpgdraw( -a, 0 );
  cpgdraw( a, 0 );
  cpgdraw( a, V0 );
  cpgdraw( xmax, V0 );
}



/* ---------------------------------------------------------------------
   main
   --------------------------------------------------------------------- */
int main(void)
{
  double V0, a, dx;
  double E, deltaE,  d2phi, dphi, phi, x;
  double x_old, phi_old, xmax, bottom, top;
  int    parity, i;

  /* 物理的設定 */
  V0   = 4.0;
  a    = 4.0;
  printf( "Parity = ?\n" ); scanf( "%d", &parity );

  /* 計算設定 */
  dx   = 0.001;
  deltaE = 0.0002;

  /* 固有値の推定値 */
  E = 0.1;
  printf( "E = ?\n" ); scanf( "%lf", &E );

  /* パリティにより phi, dphi の初期値を決定 */
  ParityCheck( parity, &phi, &dphi );

  for(;;){
    x = 0;
    /* 大きな x までオイラー法を繰り返す */
    while( x <= 6 ){
      Euler( V0, a, x, dx, E, &dphi, &phi);
      x += dx;
    }
    /* 大きな x で phi がほぼ 0 となるようにEを決める */
    if( fabs(phi) < 0.02 ) break;
    /* Eの値を deltaE だけ変え, 波動関数を初期値に戻し, 計算を繰り返す */
    E += deltaE;
    ParityCheck( parity, &phi, &dphi );
  }

  printf("Eigen value = %lf\n", E);

  /* グラフ描画設定 */
  xmax = 6.0;
  bottom = 3.0;
  top = 5.0;

  cpgopen( "/xserv" );
  cpgpap( 5.0, 1.0 );
  cpgenv( -xmax, xmax, -bottom, top, 0, 0 );

  /* 井戸型ポテンシャルの描画 */
  DrawPotential( a, V0, xmax);

  /* 固有値の描画 */
  cpgsci(2); cpgmove( -a, E ); cpgdraw( a, E);

  /* 再計算しながらグラフ描画 */
  x = 0; ParityCheck( parity, &phi, &dphi );
  cpgsci(9); cpgslw(1);
  while( x<=xmax ){
    x_old = x;
    phi_old = phi;
    x += dx;
    Euler( V0, a, x, dx, E , &dphi, &phi);
    cpgmove( x_old, phi_old ); cpgdraw( x, phi );
    cpgmove( -x_old, phi_old * parity ); cpgdraw( -x, phi * parity );
  }
  
  cpgclos();
  return 0;
}