//=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=
//    Program huckel
//=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=

#include "cip.h"
#include "nxgraph.h"


//---- declaration of class Array2d

class Array2d
{
  int dim;         // dimension of this array
  double*  head;   // pointer to assign head of this array
  double** array;  // pointers to assign head of each rows
public:
       Array2d( int _dim );
      ~Array2d();
  double* operator [] ( int i );
  int  Dim( void );
  void SetZero( void );
  void SetUnit( void );
  void Print( void );
};

//---- definition of class Array2d
//---- constructor, allocates memory and set pointers
Array2d::Array2d( int _dim )
{
  dim   = _dim;
  head  = new (double)  [dim*dim];
  array = new (double*) [dim];
  for( int i=0 ; i<dim ; i++ ) array[i] = &head[i*dim];
}

//---- destructor, free the memory
Array2d::~Array2d()
{
  delete [] head;
  delete [] array;
}

//---- represents head address of each rows
inline double* Array2d::operator [] ( int i )
{
  return array[i];
}
//---- represent dimension of this array
inline int Array2d::Dim( void )
{
  return dim;
}
//---- sets zero to all elements
void Array2d::SetZero( void )
{
  for( int i=0 ; i<dim ; i++ ){
    for( int j=0 ; j<dim ; j++ ){
      array[i][j] = 0.0;
    }
  }
}
//---- sets unit matrix
void Array2d::SetUnit( void )
{
  for( int i=0 ; i<dim ; i++ ){
    for( int j=i+1 ; j<dim ; j++ ){
      array[i][j] = array[j][i] = 0.0;
    }
    array[i][i] = 1.0;
  }
}

//---- prints all elements
void Array2d::Print( void )
{
  for( int i=0 ; i<dim ; i++ ){
    for( int j=0 ; j<dim ; j++ ){
      printf("%+6.4f ", array[i][j] );
    }
    printf("\n");
  }
  printf("\n");
}

//---- prototypes of functions
void Jacobi( Array2d& A, Array2d& V );
void DrawPz( int x, int y, double c );
void Draw( Array2d& A, Array2d& V );

//---- main function
int main( void )
{
  const int ORBITS = 4;

  Array2d A(ORBITS), V(ORBITS);

  A.SetZero();

  A[0][1] = A[1][0] = 1.0;
  A[1][2] = A[2][1] = 1.0;
  A[2][3] = A[3][2] = 1.0;

  Jacobi( A, V );

  A.Print();
  V.Print();

  const int WIN_WIDTH  = 24 + (ORBITS-1)*28 + 16 + 76;
  const int WIN_HEIGHT = 32 + (ORBITS-1)*64 + 32;

  NXOpenWindow( "Huckel approximation", WIN_WIDTH, WIN_HEIGHT );

  Draw( A, V );

  NXPause();
  NXCloseWindow();
  return 0;
}


//---- sign of given argument
inline int sgn( double theta )
{
  if( theta >= 0.0 ) return +1;
  return -1;
}

//---- rotates Vp and Vq by give cos and sin
inline void Rotate( double c, double s, double& Vp, double& Vq )
{
  double _Vp;
  Vp = c*(_Vp=Vp) - s*Vq;
  Vq = s*_Vp + c*Vq;
}

//---- calculates eigen values and eigen vectors
void Jacobi( Array2d& A, Array2d& V )
{
  double theta, s, c, t;
  int p, q, r, nonzero;

  V.SetUnit();
  do{
    nonzero = 0;
    for( p=0 ; p<A.Dim() ; p++ ){
      for( q=p+1 ; q<A.Dim() ; q++ ){
	if( A[p][q] == 0.0 ) continue;
	nonzero++;

	theta = (A[q][q]-A[p][p])/(2*A[p][q]);

	if( -1e8<theta && theta<+1e8 ){
	  t = 1.0/(theta+sgn(theta)*sqrt(theta*theta+1.0));
	  c = 1.0/sqrt(t*t+1);
	  s = t*c;
	}else{
	  t = 0.5/theta;
	  c = 1.0;
	  s = t;
	}

	A[p][p] -= A[p][q]*t;
	A[q][q] += A[p][q]*t;
	A[p][q]  = A[q][p] = 0.0;

	for( r=0 ; r<A.Dim() ; r++ ){
	  Rotate( c, s, V[r][p], V[r][q] );

	  if( r==q || r==p ) continue;

	  Rotate( c, s, A[r][p], A[r][q] );
	  A[p][r] = A[r][p];
	  A[q][r] = A[r][q];
	}
      }
    }
  }while(nonzero);
}

//---- draws pz orbital with given magnitude at given position
void DrawPz( int x, int y, double c )
{
  int r;
  if( c>0.0 ){
    r = int(c*24.0);
    NXSetColor( NX_PINK );    NXFillCircle( x, y-r, r );
    NXSetColor( NX_CYAN );    NXFillCircle( x, y+r, r );
  }else{
    r = int(-c*24.0);
    NXSetColor( NX_CYAN );    NXFillCircle( x, y-r, r );
    NXSetColor( NX_PINK );    NXFillCircle( x, y+r, r );
  }
}

//---- draws orbitals
void Draw( Array2d& A, Array2d& V )
{
  int i,j;

  NXClearWindow();
  for( j=0 ; j<V.Dim() ; j++ ){
    for( i=0 ; i<V.Dim() ; i++ ){
      DrawPz( i*28+24, j*64+32, V[i][j] );
    }
    NXSetColor( NX_BLACK );
    NXDrawPrintf( i*28+16, j*64+36, "X = %+6.4f", A[j][j] );
  }
  NXFlush();
}