src/qcevolve.cpp

/***************************************************************************
 *   Quantum Construct (qC++)                                                                  *
 *   The Quantum Physics Computational Library                                         *
 *   Copyright (C) 2005 by Shekhar S. Chandra                                      *
 *   Shekhar.Chandra@sci.monash,edu.au                                     *
 *                                                                         *
 *   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.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/
/*
        Ver Info - QCEvolve.cpp
        1.0.0 - Original Implementation
*/

#include "qcevolve.h"

QCEvolve::QCEvolve(void)
{
        h_t = 1e-4;
}

QCEvolve::~QCEvolve(void)
{
}

void QCEvolve::evolveForwards(bool forwards)
{
        if(forwards && h_t < 0.0)
                h_t *= -1.0;
        else if(!forwards && h_t > 0.0)
                h_t *= -1.0;
}

bool QCEvolve::isForward()
{
        if(h_t < 0.0)
                return false;
        else
                return true;
}

void QCEvolve::usingTDGP(QCComplexField<2> &currentField, parameterType xOffset, parameterType yOffset, parameterType &h_x, parameterType &h_y, QCParameter &tdgpParameters)
{
        Array<complexType,2> oldField = *currentField.accessArray();
        int maxIterations = 4, xExtent = currentField.noOfRows()-1, yExtent = currentField.noOfColumns()-1;
        Array<complexType,1> a(xExtent),b(xExtent),c(xExtent),d(xExtent),X(xExtent);
        //It is assumed that Small Omega Squared is 1st Parameter, mu is the 2nd, g is 3rd Parameter and hbar the 4th Parameter and m the 5th in the QCParameter Class. All index referencing should be made C style.
    parameterType x, y, omegaSq = tdgpParameters.getParameter(0), g = tdgpParameters.getParameter(2), hbar = tdgpParameters.getParameter(3), m = tdgpParameters.getParameter(4);

        //Important Constants
        x = 4.0*h_x*h_x*m; //x is used as tmp storage
        y = 4.0*h_y*h_y*m; //y is used as tmp storage
        complexType     alpha_x(0.0,(h_t*hbar)/x),
        alpha_y(0.0,(h_t*hbar)/y),
        gamma(0.0,(m*omegaSq)/(2.0*hbar)),
        beta(0.0,g/hbar),
        Psi, oldPsi;

        for(int p = 0; p < maxIterations; p ++)
        {
                for(int k = 1; k < yExtent; k ++)
                {
                        y = (k)*h_y - yOffset;
                        for(int j = 1; j < xExtent; j ++)
                        {
                                x = (j)*h_x - xOffset;
                                Psi = currentField.retrieveElement(j,k);
                                oldPsi = oldField(j,k);

                                a(j) = -alpha_x;
                                b(j) = 1.0 + 2.0*alpha_x;
                                c(j) = -alpha_x;
                                d(j) = oldPsi
                                       + alpha_x*( oldField(j+1,k) + oldField(j-1,k)
                                                   + oldField(j,k+1) + oldField(j,k-1) - 4.0*oldPsi)
                                       - h_t*0.5*i/hbar*(0.5*m*omegaSq*(x*x + y*y) + g*norm(Psi))*Psi
                                       - h_t*0.5*i/hbar*(0.5*m*omegaSq*(x*x + y*y) + g*norm(oldPsi))*oldPsi;
                //cerr << abs(a(j)) << " " << abs(d(j)) << " " << abs(oldField(j,k)) << endl;
                        }
                        //exit(0);
                        //Free Boundary Conditions
                        b(1) = b(1) + 2.0*a(1);
                        c(1) = c(1) - a(1);
                        a(xExtent-1) = a(xExtent-1) - c(xExtent-1);
                        b(xExtent-1) = b(xExtent-1) + 2.0*c(xExtent-1);
                        //Solve Tridiagonal
                        LA.solveTridiagonal(a,b,c,X,d);
                        for(int j = 1; j < xExtent; j ++)
                        {
                                currentField.assignElement(j,k,X(j));
                                //if(p == 3) cerr << abs(X(j)) << ", " << arg(X(j)) << endl;
                        }
                        //if(p == 3) exit(0);
                }
        //cerr << "\n " << abs(a(1)) << ", " << abs(d(1)) << ", OldField: " << abs(oldField(1,1));
                for(int j = 1; j < xExtent; j ++)
                {
                        for(int k = 1; k < yExtent; k ++)
                        {
                                //a(k) = -alpha_y;
                                //b(k) = 1.0 + 2.0*alpha_y;
                                //c(k) = -alpha_y;
                                d(k) = currentField.retrieveElement(j,k);
                        }
                        //Free Boundary Conditions
                        //b(1) = b(1) + 2.0*a(1); //DO NOT USE CONTRACTED OPERATORS HERE (like =- etc.)
                        //c(1) = c(1) - a(1);
                        //a(yExtent-1) = a(yExtent-1) - c(yExtent-1);
                        //b(yExtent-1) = b(yExtent-1) + 2.0*c(yExtent-1);
                        //Solve Tridiagonal
                        LA.solveTridiagonal(a,b,c,X,d);
                        for(int k = 1; k < yExtent; k ++)
                        {
                                currentField.assignElement(j,k,X(k));
                                //cerr << abs(X(k)) << ", " << arg(X(k)) << endl;
                        }
                        //exit(0);
                }
                //Boundary Conditions for Edges
                setBoundaryGradZero(currentField,oldField);
        }
}

void QCEvolve::usingTDGPLz(QCComplexField<2> &currentField, parameterType xOffset, parameterType yOffset, parameterType &h_x, parameterType &h_y, QCParameter &tdgpParameters)
{
        int maxIterations = 4, xExtent = currentField.noOfRows()-1, yExtent = currentField.noOfColumns()-1;
        Array<complexType,1> a(xExtent),b(xExtent),c(xExtent),d(xExtent),X(xExtent);
        Array<complexType,2> oldField = *currentField.accessArray();
        //It is assumed that Small Omega Squared is 1st Parameter, mu is the 2nd, g is 3rd Parameter and hbar the 4th Parameter and m the 5th and Big Omega the 6th in the QCParameter Class. All index referencing should be made C style.
    parameterType x, y, omegaSq = tdgpParameters.getParameter(0), g = tdgpParameters.getParameter(2), hbar = tdgpParameters.getParameter(3), m = tdgpParameters.getParameter(4), bigOmega = tdgpParameters.getParameter(5);

        //Important Constants
        x = 4.0*h_x*h_x*m; //x is used as tmp storage
        y = 4.0*h_y*h_y*m; //y is used as tmp storage
        complexType     alpha_x(0.0,(h_t*hbar)/x),
        alpha_y(0.0,(h_t*hbar)/y),
        gamma(0.0,(m*omegaSq)/(2.0*hbar)),
        beta(0.0,g/hbar),
        Psi, oldPsi;

        for(int m = 0; m < maxIterations; m ++)
        {
                for(int k = 1; k < yExtent; k ++)
                {
                        y = (k)*h_y - yOffset;
                        for(int j = 1; j < xExtent; j ++)
                        {
                                x = (j)*h_x - xOffset;
                                Psi = currentField.retrieveElement(j,k);
                                oldPsi = oldField(j,k);

                                a(j) = -alpha_x;
                                b(j) = 1.0 + 2.0*alpha_x;
                                c(j) = -alpha_x;
                                d(j) = oldPsi
                                       + alpha_x*(oldField(j+1,k) + oldField(j-1,k) - 2.0*oldPsi)
                                       + alpha_y*(oldField(j,k+1) + oldField(j,k-1) - 2.0*oldPsi)
                                       - h_t*0.25*bigOmega*(y/h_x*(currentField.retrieveElement(j+1,k) - currentField.retrieveElement(j-1,k)) - x/h_y*(currentField.retrieveElement(j,k+1) - currentField.retrieveElement(j,k-1))
                                                            + y/h_x*(oldField(j+1,k) - oldField(j-1,k)) - x/h_y*(oldField(j,k+1) - oldField(j,k-1)))
                                       - h_t*0.5*i/hbar*(0.5*m*omegaSq*(x*x + y*y)
                                                         + g*norm(Psi))*Psi - h_t*0.5*i/hbar*(0.5*m*omegaSq*(x*x + y*y)
                                                                                              + g*norm(oldPsi))*oldPsi;
                        }
                        //Free Boundary Conditions
                        b(1) = b(1) + 2.0*a(1);
                        c(1) = c(1) - a(1);
                        a(xExtent-1) = a(xExtent-1) - c(xExtent-1);
                        b(xExtent-1) = b(xExtent-1) + 2.0*c(xExtent-1);
                        //Solve Tridiagonal
                        LA.solveTridiagonal(a,b,c,X,d);
                        for(int j = 1; j < xExtent; j ++)
                                currentField.assignElement(j,k,X(j));
                }

                for(int j = 1; j < xExtent; j ++)
                {
                        for(int k = 1; k < yExtent; k ++)
                        {
                                a(k) = -alpha_y;
                                b(k) = 1.0 + 2.0*alpha_y;
                                c(k) = -alpha_y;
                                d(k) = currentField.retrieveElement(j,k);
                        }
                        //Free Boundary Conditions
                        b(1) = b(1) + 2.0*a(1);
                        c(1) = c(1) - a(1);
                        a(yExtent-1) = a(yExtent-1) - c(yExtent-1);
                        b(yExtent-1) = b(yExtent-1) + 2.0*c(yExtent-1);
                        //Solve Tridiagonal
                        LA.solveTridiagonal(a,b,c,X,d);
                        for(int k = 1; k < yExtent; k ++)
                                currentField.assignElement(j,k,X(k));
                }
                //Boundary Conditions for Edges
                setBoundaryVanish(currentField);
        }
}

void QCEvolve::usingTDGP2Component(QCComplexField<2> &component1, QCComplexField<2> &component2, parameterType xOffset, parameterType yOffset, parameterType &h_x, parameterType &h_y, QCParameter &tdgpParameters)
{
        int maxIterations = 4, xExtent = component1.noOfRows()-1, yExtent = component1.noOfColumns()-1;
        Array<complexType,1> a1(xExtent),b1(xExtent),c1(xExtent),d1(xExtent),X1(xExtent);
        Array<complexType,1> a2(xExtent),b2(xExtent),c2(xExtent),d2(xExtent),X2(xExtent);
        Array<complexType,2> oldComponent1 = *component1.accessArray();
        Array<complexType,2> oldComponent2 = *component2.accessArray();
        parameterType x, y, omegaSq = tdgpParameters.getParameter(0), hbar = tdgpParameters.getParameter(3), m = tdgpParameters.getParameter(4), bigOmega = tdgpParameters.getParameter(5);
        parameterType g11 = tdgpParameters.getParameter(7), g22 = tdgpParameters.getParameter(8), g12 = tdgpParameters.getParameter(9), g21 = tdgpParameters.getParameter(10);

        //Important Constants
        x = 4.0*h_x*h_x*m; //x is used as tmp storage
        y = 4.0*h_y*h_y*m; //y is used as tmp storage
        complexType     alpha_x(0.0,(h_t*hbar)/x),
        alpha_y(0.0,(h_t*hbar)/y);

        for(int m = 0; m < maxIterations; m ++)
        {
                for(int k = 1; k < yExtent; k ++)
                {
                        y = k*h_y - yOffset;
                        for(int j = 1; j < xExtent; j ++)
                        {
                                x = j*h_x - xOffset;

                                a1(j) = a2(j) = -alpha_x;
                                b1(j) = b2(j) = 1.0 + 2.0*alpha_x;
                                c1(j) = c2(j) = -alpha_x;
                                d1(j) = oldComponent1(j,k)
                                        + alpha_x*(oldComponent1(j+1,k) + oldComponent1(j-1,k) - 2.0*oldComponent1(j,k))
                                        + alpha_y*(oldComponent1(j,k+1) + oldComponent1(j,k-1) - 2.0*oldComponent1(j,k))
                                        /*- h_t*0.25*bigOmega*(y/h_x*(component1.retrieveElement(j+1,k) - component1.retrieveElement(j-1,k)) - x/h_y*(component1.retrieveElement(j,k+1) - component1.retrieveElement(j,k-1))
                                                             + y/h_x*(oldComponent1(j+1,k) - oldComponent1(j-1,k)) - x/h_y*(oldComponent1(j,k+1) - oldComponent1(j,k-1)))*/
                                        - h_t*0.5*i/hbar*(0.5*m*omegaSq*(x*x + y*y)
                                                          + g11*norm(component1.retrieveElement(j,k))
                                                          + g12*norm(component2.retrieveElement(j,k)))*component1.retrieveElement(j,k)
                                        - h_t*0.5*i/hbar*(0.5*m*omegaSq*(x*x + y*y)
                                                          + g11*norm(oldComponent1(j,k))
                                                          + g12*norm(oldComponent2(j,k)))*oldComponent1(j,k);
                                d2(j) = oldComponent2(j,k)
                                        + alpha_x*(oldComponent2(j+1,k) + oldComponent2(j-1,k) - 2.0*oldComponent2(j,k))
                                        + alpha_y*(oldComponent2(j,k+1) + oldComponent2(j,k-1) - 2.0*oldComponent2(j,k))
                                        /*+ h_t*0.25*bigOmega*(y/h_x*(component2.retrieveElement(j+1,k) - component2.retrieveElement(j-1,k)) - x/h_y*(component2.retrieveElement(j,k+1) - component2.retrieveElement(j,k-1))
                                                             + y/h_x*(oldComponent2(j+1,k) - oldComponent2(j-1,k)) - x/h_y*(oldComponent2(j,k+1) - oldComponent2(j,k-1)))*/
                                        - h_t*0.5*i/hbar*(0.5*m*omegaSq*(x*x + y*y)
                                                          + g22*norm(component2.retrieveElement(j,k))
                                                          + g21*norm(component1.retrieveElement(j,k)))*component2.retrieveElement(j,k)
                                        - h_t*0.5*i/hbar*(0.5*m*omegaSq*(x*x + y*y)
                                                          + g22*norm(oldComponent2(j,k))
                                                          + g21*norm(oldComponent1(j,k)))*oldComponent2(j,k);
                        }
                        //Free Boundary Conditions
                        //Component 1
                        b1(1) = b1(1) + 2.0*a1(1);
                        c1(1) = c1(1) - a1(1);
                        a1(xExtent-1) = a1(xExtent-1) - c1(xExtent-1);
                        b1(xExtent-1) = b1(xExtent-1) + 2.0*c1(xExtent-1);
                        //Component 2
                        b2(1) = b2(1) + 2.0*a2(1);
                        c2(1) = c2(1) - a2(1);
                        a2(xExtent-1) = a2(xExtent-1) - c2(xExtent-1);
                        b2(xExtent-1) = b2(xExtent-1) + 2.0*c2(xExtent-1);
                        //Solve Tridiagonals
                        LA.solveTridiagonal(a1,b1,c1,X1,d1);
                        LA.solveTridiagonal(a2,b2,c2,X2,d2);
                        for(int j = 1; j < xExtent; j ++)
                        {
                                component1.assignElement(j,k,X1(j));
                                component2.assignElement(j,k,X2(j));
                        }
                }

                for(int j = 1; j < xExtent; j ++)
                {
                        for(int k = 1; k < yExtent; k ++)
                        {
                                a1(k) = a2(k) = -alpha_y;
                                b1(k) = b2(k) = 1.0 + 2.0*alpha_y;
                                c1(k) = c2(k) = -alpha_y;
                                d1(k) = component1.retrieveElement(j,k);
                                d2(k) = component2.retrieveElement(j,k);
                        }
                        //Free Boundary Conditions
                        //Component 1
                        b1(1) = b1(1) + 2.0*a1(1);
                        c1(1) = c1(1) - a1(1);
                        a1(yExtent-1) = a1(yExtent-1) - c1(yExtent-1);
                        b1(yExtent-1) = b1(yExtent-1) + 2.0*c1(yExtent-1);
                        //Component 2
                        b2(1) = b2(1) + 2.0*a2(1);
                        c2(1) = c2(1) - a2(1);
                        a2(yExtent-1) = a2(yExtent-1) - c2(yExtent-1);
                        b2(yExtent-1) = b2(yExtent-1) + 2.0*c2(yExtent-1);
                        //Solve Tridiagonals
                        LA.solveTridiagonal(a1,b1,c1,X1,d1);
                        LA.solveTridiagonal(a2,b2,c2,X2,d2);
                        for(int k = 1; k < yExtent; k ++)
                        {
                                component1.assignElement(j,k,X1(k));
                                component2.assignElement(j,k,X2(k));
                        }
                }
                //Boundary Conditions for Edges
                setBoundaryVanish(component1);
                setBoundaryVanish(component2);
        }
}

void QCEvolve::usingTDGP3D(QCComplexField<3> &currentField, parameterType xOffset, parameterType yOffset, parameterType zOffset, parameterType &h_x, parameterType &h_y, parameterType &h_z, QCParameter &tdgpParameters)
{
        Array<complexType,3> oldField = *currentField.accessArray();
        int maxIterations = 4, xExtent = currentField.noOfRows()-1, yExtent = currentField.noOfColumns()-1, zExtent = currentField.height()-1;
        Array<complexType,1> a(xExtent),b(xExtent),c(xExtent),d(xExtent),X(xExtent);
        //It is assumed that Small Omega Squared is 1st Parameter, mu is the 2nd, g is 3rd Parameter and hbar the 4th Parameter and m the 5th in the QCParameter Class. All index referencing should be made C style.
    parameterType x, y, z, omegaSq = tdgpParameters.getParameter(0), g = tdgpParameters.getParameter(2), hbar = tdgpParameters.getParameter(3), m = tdgpParameters.getParameter(4);

        //Important Constants
        x = 4.0*h_x*h_x*m; //x is used as tmp storage
        y = 4.0*h_y*h_y*m; //y is used as tmp storage
        z = 4.0*h_z*h_z*m; //z is used as tmp storage
        complexType     alpha_x(0.0,(h_t*hbar)/x),
        alpha_y(0.0,(h_t*hbar)/y),
        alpha_z(0.0,(h_t*hbar)/z),
        gamma(0.0,(m*omegaSq)/(2.0*hbar)),
        beta(0.0,g/hbar),
        Psi, oldPsi;

        for(int m = 0; m < maxIterations; m ++)
        {
                for(int l = 1; l < zExtent; l ++)
                {
                        z = (l)*h_z - zOffset;
                        for(int k = 1; k < yExtent; k ++)
                        {
                                y = (k)*h_y - yOffset;
                                for(int j = 1; j < xExtent; j ++)
                                {
                                        x = (j)*h_x - xOffset;
                                        Psi = currentField.retrieveElement(j,k,l);
                                        oldPsi = oldField(j,k,l);

                                        a(j) = -alpha_x;
                                        b(j) = 1.0 + 2.0*alpha_x;
                                        c(j) = -alpha_x;
                                        d(j) = oldPsi
                                                   + alpha_x*( oldField(j+1,k,l) + oldField(j-1,k,l)
                                                   + oldField(j,k+1,l) + oldField(j,k-1,l)
                                                   + oldField(j,k,l+1) + oldField(j,k,l-1) - 6.0*oldPsi)
                                                   - h_t*0.5*i/hbar*(0.5*m*omegaSq*(x*x + y*y + z*z) + g*norm(Psi))*Psi
                                                   - h_t*0.5*i/hbar*(0.5*m*omegaSq*(x*x + y*y + z*z) + g*norm(oldPsi))*oldPsi;
                                }
                                //Free Boundary Conditions
                                b(1) = b(1) + 2.0*a(1);
                                c(1) = c(1) - a(1);
                                a(xExtent-1) = a(xExtent-1) - c(xExtent-1);
                                b(xExtent-1) = b(xExtent-1) + 2.0*c(xExtent-1);
                                //Solve Tridiagonal
                                LA.solveTridiagonal(a,b,c,X,d);
                                for(int j = 1; j < xExtent; j ++)
                                        currentField.assignElement(j,k,l,X(j));
                        }
                }

                for(int l = 1; l < zExtent; l ++)
                {
                        for(int j = 1; j < xExtent; j ++)
                        {
                                for(int k = 1; k < yExtent; k ++)
                                {
                                        a(k) = -alpha_y;
                                        b(k) = 1.0 + 2.0*alpha_y;
                                        c(k) = -alpha_y;
                                        d(k) = currentField.retrieveElement(j,k,l);
                                }
                                //Free Boundary Conditions
                                b(1) = b(1) + 2.0*a(1); //DO NOT USE CONTRACTED OPERATORS HERE (like =- etc.)
                                c(1) = c(1) - a(1);
                                a(yExtent-1) = a(yExtent-1) - c(yExtent-1);
                                b(yExtent-1) = b(yExtent-1) + 2.0*c(yExtent-1);
                                //Solve Tridiagonal
                                LA.solveTridiagonal(a,b,c,X,d);
                                for(int k = 1; k < yExtent; k ++)
                                        currentField.assignElement(j,k,l,X(k));
                        }
                }

                for(int j = 1; j < xExtent; j ++)
                {
                        for(int k = 1; k < yExtent; k ++)
                        {
                                for(int l = 1; l < zExtent; l ++)
                                {
                                        a(l) = -alpha_z;
                                        b(l) = 1.0 + 2.0*alpha_z;
                                        c(l) = -alpha_z;
                                        d(l) = currentField.retrieveElement(j,k,l);
                                }
                                //Free Boundary Conditions
                                b(1) = b(1) + 2.0*a(1); //DO NOT USE CONTRACTED OPERATORS HERE (like =- etc.)
                                c(1) = c(1) - a(1);
                                a(zExtent-1) = a(zExtent-1) - c(zExtent-1);
                                b(zExtent-1) = b(zExtent-1) + 2.0*c(zExtent-1);
                                //Solve Tridiagonal
                                LA.solveTridiagonal(a,b,c,X,d);
                                for(int l = 1; l < zExtent; l ++)
                                        currentField.assignElement(j,k,l,X(l));
                        }
                }

                //Boundary Conditions for Edges
                setBoundaryVanish(currentField);
        }
}

void QCEvolve::setBoundaryGradZero(QCComplexField<2> &field, Array<complexType,2> &old)
{
        int xExtent = field.noOfRows()-1, yExtent = field.noOfColumns()-1;

        //Phase Gradient Zero at Boundary
        for(int j = 1; j < xExtent; j ++)
        {
                field.assignElement(j,0,old(j,0));
                field.assignElement(j,yExtent,old(j,yExtent));
        }
        for(int k = 1; k < yExtent; k ++)
        {
                field.assignElement(0,k,old(0,k));
                field.assignElement(xExtent,k,old(xExtent,k));
        }
        field.assignElement(0,0,old(0,0));
        field.assignElement(0,yExtent,old(0,yExtent));
        field.assignElement(xExtent,0,old(xExtent,0));
        field.assignElement(xExtent,yExtent,old(xExtent,yExtent));
}

void QCEvolve::setBoundaryFree(QCComplexField<2> &field)
{
        int xExtent = field.noOfRows()-1, yExtent = field.noOfColumns()-1;

        for(int k = 1; k < xExtent; k ++)
        {
                field.assignElement(k,0,2.0*field.retrieveElement(k,1) - field.retrieveElement(k,2));
                field.assignElement(k,yExtent,2.0*field.retrieveElement(k,yExtent-1) - field.retrieveElement(k,yExtent-2));
        }
        for(int l = 0; l <= yExtent; l ++)
        {
                field.assignElement(0,l,2.0*field.retrieveElement(1,l) - field.retrieveElement(2,l));
                field.assignElement(xExtent,l,2.0*field.retrieveElement(xExtent-1,l) - field.retrieveElement(xExtent-2,l));
        }
        field.assignElement(0,0,2.0*field.retrieveElement(1,1) - field.retrieveElement(2,2));
        field.assignElement(0,yExtent,2.0*field.retrieveElement(1,yExtent-1) - field.retrieveElement(2,yExtent-2));
        field.assignElement(xExtent,0,2.0*field.retrieveElement(xExtent-1,1) - field.retrieveElement(xExtent-2,2));
        field.assignElement(xExtent,yExtent,2.0*field.retrieveElement(xExtent-1,yExtent-1) - field.retrieveElement(xExtent-2,yExtent-2));
}

void QCEvolve::setBoundaryVanish(QCComplexField<2> &field)
{
        int xExtent = field.noOfRows()-1, yExtent = field.noOfColumns()-1;

        //Phase Gradient Zero at Boundary
        for(int j = 1; j < xExtent; j ++)
        {
                field.assignElement(j,0,0.0*exp(i*arg(field.retrieveElement(j,0))));
                field.assignElement(j,yExtent,0.0*exp(i*arg(field.retrieveElement(j,yExtent))));
        }
        for(int k = 0; k <= yExtent; k ++)
        {
                field.assignElement(0,k,0.0*exp(i*arg(field.retrieveElement(0,k))));
                field.assignElement(xExtent,k,0.0*exp(i*arg(field.retrieveElement(xExtent,k))));
        }
}

void QCEvolve::setBoundaryVanish(QCComplexField<3> &field)
{
        int xExtent = field.noOfRows()-1, yExtent = field.noOfColumns()-1, zExtent = field.height()-1;

        //Phase Gradient Zero at Boundary
        for(int l = 1; l < zExtent; l ++)
        {
                for(int j = 1; j < xExtent; j ++)
                {
                        field.assignElement(j,0,l,complexZero);
                        field.assignElement(j,yExtent,l,complexZero);
                }
                for(int k = 0; k <= yExtent; k ++)
                {
                        field.assignElement(0,k,l,complexZero);
                        field.assignElement(xExtent,k,l,complexZero);
                }
        }
        for(int j = 0; j <= xExtent; j ++)
                for(int k = 0; k <= yExtent; k ++)
                {
                        field.assignElement(j,k,0,complexZero);
                        field.assignElement(j,k,zExtent,complexZero);
                }
}

---