state.cc

使用量子轨道方法计算量子主方程的C++库

//   State.cc -- State algebra in Hilbert space.
//     
//   Copyright (C) 1995  Todd Brun and Ruediger Schack
//   
//   This program is free software; you can redistribute it and/or modify
//   it under the terms of the GNU General Public License as published by
//   the Free Software Foundation; either version 2 of the License, or
//   (at your option) any later version.
//   
//   This program is distributed in the hope that it will be useful,
//   but WITHOUT ANY WARRANTY; without even the implied warranty of
//   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//   GNU General Public License for more details.
//   
//   You should have received a copy of the GNU General Public License
//   along with this program; if not, write to the Free Software
//   Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
//
//   ----------------------------------------------------------------------
//   If you improve the code or make additions to it, or if you have
//   comments or suggestions, please contact us:
//
//   Dr. Todd Brun			        Tel    [phone]
//   Department of Physics                      FAX    [phone]
//   Queen Mary and Westfield College           email  [e-mail]
//   Mile End Road, London E1 4NS, UK
//
//   Dr. Ruediger Schack                        Tel    [phone]
//   Department of Mathematics                  FAX    [phone]
//   Royal Holloway, University of London       email  [e-mail]
//   Egham, Surrey TW20 0EX, UK
/////////////////////////////////////////////////////////////////////////////

#include <iostream.h>
#include <math.h>
#include <stdlib.h>

#include "State.h"
#include "PrimOp.h"

static const char rcsid[] = "$Id: State.cc,v 3.1 1996/11/19 10:05:08 rschack Exp $";

// These are global variables of file scope, used to avoid recalculating
// the square roots of the integers more often than necessary.

static double smallSize = 1.0e-8;
static int intSqrtMax = 0;
static double* intSqrtList = 0;

// This utility function is used to maintain the list of square roots.
// This function has file scope.

static void setIntSqrtMax(int maximum)
{
  int newsize,i;
  double* newPtr;

  if( maximum > intSqrtMax) {			// make sure of size
    newsize = 2*maximum;			// allow plenty of room
    newPtr = new double[newsize];		// allocate space
    for (i=0; i<intSqrtMax; i++)
      newPtr[i] = intSqrtList[i];		// re-use calculated values
    for (i=intSqrtMax; i<newsize; i++)
      newPtr[i] = sqrt(i);			// calculate new values
#ifndef NON_GNU_DELETE
    if(intSqrtList != 0) delete [] intSqrtList;	// recycle memory
#else // ---- NON GNU CODE ----
    if(intSqrtList != 0) delete [intSqrtMax] intSqrtList; // recycle memory
#endif // ---- NON GNU CODE ----
    intSqrtList = newPtr;			// point to new list	
    intSqrtMax = newsize;			// remember new size
  }
}

//  I/O operators for enum FreedomType; necessary of I/O of States

ostream& operator<<( ostream& s, FreedomType theType )
//
// Output a FreedomType to an output stream
{
  switch(theType) {
    case ALL:
      s << 0;
      break;
    case FIELD:
      s << 1;
      break;
    case SPIN:
      s << 2;
      break;
    case ATOM:
      s << 3;
      break;
    default:
      cerr << "Unknown FreedomType outputted in <<!" << endl;
      exit(1);
  }
  return s;
}

istream& operator>>( istream& s, FreedomType& theType )
//
// Input a FreedomType from an input stream
{
  int i;

  s >> i;
  switch(i) {
    case 0:
      theType = ALL;
      break;
    case 1:
      theType = FIELD;
      break;
    case 2:
      theType = SPIN;
      break;
    case 3:
      theType = ATOM;
      break;
    default:
      cerr << "Unknown FreedomType inputted in >>!" << endl;
      exit(1);
  }
  return s;
}

// copy and free utility functions used by Constructors, Destructors, and
// assignment operator.

void State::free() {
  if( data != 0) {			// if memory is allocated...
#ifndef NON_GNU_DELETE
    delete [] data;			// recycle storage
    if(nFreedoms > 1) {
      delete [] nDims;			// recycle parameter arrays
      delete [] sizes;
      delete [] nSkips;
      delete [] partDims;
      delete [] freedomTypes;
      delete [] coords;
    }
#else // ---- NON GNU CODE ----
    delete [totalDim] data;		// NON GNU recycle storage
    if( nFreedoms > 1 ) {
      delete [nFreedoms] nDims;		// NON GNU recycle parameter arrays
      delete [nFreedoms] sizes;
      delete [nFreedoms] nSkips;
      delete [nFreedoms] partDims;
      delete [nFreedoms] freedomTypes;
      delete [nFreedoms] coords;
    }
#endif // ---- NON GNU CODE ----
    nFreedoms = 0;			// zero variables
    totalDim = 0;
    maxSize = 0;
    myType = ALL;
    mySize = 0;
    nSkip = 0;
    coord = 0;
    myPointer = 0;			// reset pointers to zero
    data = 0;
    nDims = 0;
    sizes = 0;
    nSkips = 0;
    coords = 0;
    partDims = 0;
    freedomTypes = 0;
  }
}

void State::copy(const State& a)	// make a copy of an existing state
{
  int i;
  if( nFreedoms != a.nFreedoms) {	// if no memory allocated, allocate it
    nFreedoms = a.nFreedoms;
    if( nFreedoms > 1 ) {
      nDims = new int[nFreedoms];
      sizes = new int[nFreedoms];
      nSkips = new int[nFreedoms];
      partDims = new int[nFreedoms];
      freedomTypes = new FreedomType[nFreedoms];
      coords = new Complex[nFreedoms];
    }
  }
  if( totalDim != a.totalDim ) {	// if no memory allocated, allocate it
    totalDim = a.totalDim;
    data = new Complex[totalDim];
  }
  maxSize = a.maxSize;			// copy member data
  nSkip = a.nSkip;
  coord = a.coord;
  mySize = a.mySize;
  myType = a.myType;
  myPointer = data;
  for (i=0; i<maxSize; i++)
    data[i] = a.data[i];
  if( nFreedoms > 1 ) {
    for (i=0; i<nFreedoms; i++) {
      nDims[i] = a.nDims[i];
      sizes[i] = a.sizes[i];
      nSkips[i] = a.nSkips[i];
      partDims[i] = a.partDims[i];
      freedomTypes[i] = a.freedomTypes[i];
      coords[i] = a.coords[i];
    }
  }
}

void State::xerox(const State& a)
// Unlike copy, xerox makes a copy of an existing state including only
// the memory which is actually being used, rather than all the memory
// allocated.  This is useful as a more efficient way of creating
// temporary states than using the normal copy-constructor.
//
// Note that xerox should NOT be used to assign to a state which already
// has memory allocated.  This results in inefficient memory handling, and
// might cause errors.
{
  int i;
#ifndef OPTIMIZE_QSD
  if( data != 0 )
    error("State::xerox should not be used on a state with allocated memory!");
#endif
  nFreedoms = a.nFreedoms;
  if( nFreedoms > 1 ) {
    nDims = new int[nFreedoms];
    sizes = new int[nFreedoms];
    nSkips = new int[nFreedoms];
    partDims = new int[nFreedoms];
    freedomTypes = new FreedomType[nFreedoms];
    coords = new Complex[nFreedoms];
  }
  totalDim = a.maxSize;
  data = new Complex[totalDim];
  maxSize = a.maxSize;			// copy member data
  nSkip = a.nSkip;
  coord = a.coord;
  mySize = a.mySize;
  myType = a.myType;
  myPointer = data;
  for (i=0; i<maxSize; i++)
    data[i] = a.data[i];
  if( nFreedoms > 1 ) {
    for (i=0; i<nFreedoms; i++) {
      nDims[i] = a.sizes[i];
      sizes[i] = a.sizes[i];
      nSkips[i] = a.nSkips[i];
      partDims[i] = a.partDims[i];
      freedomTypes[i] = a.freedomTypes[i];
      coords[i] = a.coords[i];
    }
  }
}

State::State() {	// Constructor for uninitialized State
  nFreedoms = 0;
  totalDim = 0;
  maxSize = 0;
  nDims = 0;
  sizes = 0;
  nSkips = 0;
  partDims = 0;
  freedomTypes = 0;
  coords = 0;
  data = 0;
  myPointer = 0;
  nSkip = 0;
  myType = ALL;
  mySize = 0;
  coord = 0;
}

State::State(int n, FreedomType form)   	// Constructor for 1D
{						// ground state
#ifndef OPTIMIZE_QSD
  if (n<2) error("Invalid dimension size in State constructor!");
  if( (form==SPIN) && (n != 2) ) error("Spins must have dimension 2!");
#endif
  nFreedoms = 1;
  totalDim = maxSize = mySize = n;
  nSkip = 1;
  coord = 0;
  myType = form;
  data = new Complex[totalDim];
  myPointer = data;
  nDims = 0;
  sizes = 0;
  nSkips = 0;
  partDims = 0;
  freedomTypes = 0;
  coords = 0;
  data[0] = 1.0;				// put in the ground state
  for (int i=1; i<totalDim; i++) data[i] = 0;
}

State::State(int n, int* dimensions, FreedomType* forms) // Constructor for
{							 // multidimensional
  int i;						 // ground state
#ifndef OPTIMIZE_QSD
  if (n<1) error("Illegal number of freedoms in State constructor!");
  if (n=1) error("Do not call this constructor for 1-freedom states!");
#endif
  nFreedoms = n;
  nDims = new int[n];			// allocate memory
  sizes = new int[n];
  nSkips = new int[n];
  partDims = new int[n];
  freedomTypes = new FreedomType[n];
  coords = new Complex[n];
  totalDim = 1;
  for (i=0; i<n; i++) {
#ifndef OPTIMIZE_QSD
    if(dimensions[i] < 1) error("Illegal dimension in State constructor!");
    if( (forms[i] == SPIN) && (dimensions[i] != 2) )
      error("Spin states must have 2 dimensions!");
#endif
    nSkips[i] = totalDim;
    totalDim *= (nDims[i] = dimensions[i]);
    sizes[i] = nDims[i];
    partDims[i] = totalDim;
    freedomTypes[i] = forms[i];
    coords[i] = 0;
  }
  maxSize = totalDim;
  data = new Complex[totalDim];		// allocate memory
  myPointer = data;
  mySize = totalDim;
  nSkip = 1;
  coord = 0;
  myType = freedomTypes[0];
  data[0] = 1;
  for (i=1; i<totalDim; i++) data[i] = 0;
}

State::State(int n, int nstate, FreedomType form)	// Constructor for 
{							// Fock state
#ifndef OPTIMIZE_QSD
  if (n<2) error("Invalid dimension size in State constructor!");
  if( form != FIELD )
    error("Constructor only valid for FIELD type.");
#endif
  data = 0;
  fock(n,nstate);		// call Fock state function
}

State::State(int n, Complex alpha, FreedomType form)     // Constructor for
{							 // coherent state
#ifndef OPTIMIZE_QSD
  if (n<2) error("Invalid dimension size in State constructor!");
  if( form != FIELD )
    error("Constructor only valid for FIELD type.");
#endif
  data = 0;
  coherent(n,alpha);		// call coherent state function
}

State::State(int n, int nstate, Complex alpha, FreedomType form)
//
// Constructor for excited coherent state
{
#ifndef OPTIMIZE_QSD
  if (n<2) error("Invalid dimension size in State constructor!");
  if( form != FIELD )
    error("Constructor only valid for FIELD type.");
#endif
  data = 0;
  fock(n,nstate);	// call Fock state function
  coord = alpha;	// displace in phase space
}

State::State(int n, Complex* elements, FreedomType form) 
//
// Constructor for general 1 degree of freedom state
{
#ifndef OPTIMIZE_QSD
  if (n<2) error("Invalid dimension size in State constuctor!");
  if( (form==SPIN) && (n != 2) ) error("Spins must have dimension 2!");
#endif
  nFreedoms = 1;
  totalDim = maxSize = mySize = n;
  nSkip = 1;
  coord = 0;
  myType = form;
  data = new Complex[totalDim];
  myPointer = data;
  nDims = 0;
  sizes = 0;
  nSkips = 0;
  partDims = 0;
  freedomTypes = 0;
  coords = 0;
  for (int i=0; i<totalDim; i++) data[i] = elements[i];
}

State::State(const State& a) 			// Copy-initializer
{
  nFreedoms = 0;
  totalDim = 0;
  copy(a);
}

State::State(int n, State* stateList)		// Constructor for
{						// product state
#ifndef OPTIMIZE_QSD
  if (n<1) error("Illegal number of freedoms in State constructor!");
#endif
  data = 0;
  productState(n,stateList); 	// call product state function
}

State::~State() { free(); }		// Destructor

void State::fock(int n, int nstate)	// Create a Fock State
{
#ifndef OPTIMIZE_QSD
  if( (nstate >= n) || (n < 2) || (nstate < 0) )	// check dimension
    error("Illegal arguments to State::fock!");		// exit on error