src/qcsteadystate.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 - QCSteadyState.cpp
        1.0.0 - Original Implementation
*/

#include "qcsteadystate.h"
#include "qclinearalgebra.h"

QCSteadyState::QCSteadyState(void)
{
        tolerance = 1e-10;
}

QCSteadyState::~QCSteadyState(void)
{
        F.free();
}

void QCSteadyState::setConstant(QCComplexField<2>& field)
{
        complexType tmp(1.0,0.0);

        for(int k = 0; k < field.noOfColumns(); k ++)
                for(int j = 0; j < field.noOfRows(); j ++)
                        field.assignElement(j,k,tmp);
}

void QCSteadyState::usingTDGPFixedPoint(QCComplexField<2> &field, parameterType &deltaR, QCParameter &tdgpParameters)
{
        parameterType sum = 0.0, oldSum = 1.0, r = 0.0;
        //It is assumed that Small Omega Squared is 1st Parameter, Mu is 2nd Parameter and g the 3rd Parameter.
        parameterType smallOmegaSq = tdgpParameters.getParameter(0), mu = tdgpParameters.getParameter(1), g = tdgpParameters.getParameter(2);
        int radialExtent = field.noOfColumns()-1;

        //Adjust functional array
        F.free();
        F.resize(radialExtent+1);
        F = 1.0;

        //Use Fixed Point method to acquire boundary value problem solution
        while(abs(sum - oldSum) > tolerance && radialExtent > 0)
        {
                oldSum = sum;
                sum = 0.0;
                for(int j = 1; j < radialExtent; j ++)
                {
                        r = j*deltaR;
                        F(j) = 0.5*(F(j+1)+F(j-1)+0.5*deltaR*(F(j+1)-F(j-1))/j) - deltaR*deltaR*F(j)*(g*F(j)*F(j) + 0.5*smallOmegaSq*r*r - mu);
                        sum += F(j);
                }
                F(radialExtent) = 0.0;
                F(0) = 2*F(1) - F(2);
        }
        completedSolution = TDGPFixedPt;
}

void QCSteadyState::usingTDGPFixedPoint(QCComplexField<3> &field, parameterType &deltaR, QCParameter &tdgpParameters)
{
        parameterType sum = 0.0, oldSum = 1.0, r = 0.0;
        //It is assumed that Small Omega Squared is 1st Parameter, Mu is 2nd Parameter and g the 3rd Parameter.
        parameterType smallOmegaSq = tdgpParameters.getParameter(0), mu = tdgpParameters.getParameter(1), g = tdgpParameters.getParameter(2);
        int radialExtent = field.noOfColumns()-1;

        //Adjust functional array
        F.free();
        F.resize(radialExtent+1);
        F = 1.0;

        //Use Fixed Point method to acquire boundary value problem solution
        while(abs(sum - oldSum) > tolerance && radialExtent > 0)
        {
                oldSum = sum;
                sum = 0.0;
                for(int j = 1; j < radialExtent; j ++)
                {
                        r = j*deltaR;
                        F(j) = 0.5*(F(j+1)+F(j-1)+0.5*deltaR*(F(j+1)-F(j-1))/j) - deltaR*deltaR*F(j)*(g*F(j)*F(j) + 0.5*smallOmegaSq*r*r - mu);
                        sum += F(j);
                }
                F(radialExtent) = 0.0;
                F(0) = 2*F(1) - F(2);
        }
        completedSolution = TDGPFixedPt;
}

void QCSteadyState::usingTDGPVortexFixedPoint(QCComplexField<2> &field, parameterType &deltaR, QCParameter &tdgpParameters)
{
        parameterType sum = 0.0, oldSum = 1.0, r = 0.0, frac1, frac2, internal, extra, jSq, jFrac;
        //It is assumed that Small Omega Squared is 1st Parameter, Mu is 2nd Parameter, g the 3rd Parameter, Omega the 6th and n the 7th.
        parameterType smallOmegaSq = tdgpParameters.getParameter(0), mu = tdgpParameters.getParameter(1), g = tdgpParameters.getParameter(2), Omega = tdgpParameters.getParameter(5), n = tdgpParameters.getParameter(6);
        int radialExtent = field.noOfColumns()-1, count = 0;

        //Adjust functional array
        F.free();
        F.resizeAndPreserve(radialExtent+1);
        F = 1.0;
        //cerr << F << endl;
        if(radialExtent > 0)
                F(0) = 0.0;

        //Use Fixed Point method to acquire boundary value problem solution
        while(abs(sum - oldSum) > tolerance && radialExtent > 0)
        {
                oldSum = sum;
                sum = 0.0;
                for(int j = 1; j < radialExtent; j ++)
                {
                        r = j*deltaR;
                        internal = g*F(j)*F(j) + 0.5*smallOmegaSq*r*r - mu + n*Omega;
                        jFrac = (F(j+1)-F(j-1))/j;
                        frac1 = 0.5*deltaR*jFrac;
                        jSq = j*j;
                        frac2 = (n*n)/(jSq)*F(j);
                        extra = F(j+1)+F(j-1);
                        F(j) = 0.5*(extra+frac1-frac2) - deltaR*deltaR*F(j)*(internal);
                        //cerr << F(j) << endl;

                        //F(j) = 0.5*(F(j+1)+F(j-1)+0.5*deltaR*(F(j+1)-F(j-1))/j-(n*n)/(j*j)*F(j)) - deltaR*deltaR*F(j)*(g*F(j)*F(j) + 0.5*smallOmegaSq*r*r - mu + n*Omega);
                        sum += F(j);
                }
                //exit(0);
                //cerr << count << " ";
                //count ++;
                F(radialExtent) = 0.0;
        }
        completedSolution = TDGPFixedPt;
}

void QCSteadyState::usingTDGPVortexFixedPoint(QCComplexField<3> &field, parameterType xOffset, parameterType yOffset, parameterType zOffset, parameterType &deltaX, parameterType &deltaY, parameterType &deltaZ, parameterType &deltaR, QCParameter &tdgpParameters)
{
        parameterType sum = 0.0, oldSum = 1.0, r = 0.0, x = 0.0, y = 0.0, z = 0.0;
        //It is assumed that Small Omega Squared is 1st Parameter, Mu is 2nd Parameter, g the 3rd Parameter, Omega the 6th and n the 7th.
        parameterType smallOmegaSq = tdgpParameters.getParameter(0), mu = tdgpParameters.getParameter(1), g = tdgpParameters.getParameter(2), Omega = tdgpParameters.getParameter(5), n = tdgpParameters.getParameter(6);
        int radialExtent = field.noOfColumns()-1;
        parameterType m, rSquared;
        parameterType magnitude;
        int rInteger;

        //Adjust functional array
        F.free();
        F.resize(radialExtent+1);
        F = 1.0;
        if(radialExtent > 0)
                F(0) = 0.0;

        for(int l = 0; l < field.height(); l ++)
        {
                z = l*deltaZ - zOffset;
                //Use Fixed Point method to acquire boundary value problem solution
                while(abs(sum - oldSum) > tolerance && radialExtent > 0)
                {
                        oldSum = sum;
                        sum = 0.0;
                        for(int j = 1; j < radialExtent; j ++)
                        {
                                r = j*deltaR;
                                F(j) = 0.5*(F(j+1)+F(j-1)+0.5*deltaR*(F(j+1)-F(j-1))/j-(n*n)/(j*j)*F(j)) - deltaR*deltaR*F(j)*(g*F(j)*F(j) + 0.5*smallOmegaSq*(r*r + z*z) - mu + n*Omega);
                                sum += F(j);
                        }
                        F(radialExtent) = 0.0;
                }
                //Transfer Solution
                for(int k = 0; k < field.noOfColumns(); k++)
                {
                        y = (k)*deltaY - yOffset;
                        for(int j = 0; j < field.noOfRows(); j++)
                        {
                                x = (j)*deltaX - xOffset;
                                rSquared = x*x + y*y;
                                m = sqrt(rSquared)/deltaR;
                                rInteger = int(m);
                                magnitude = F(rInteger) + (m-rInteger)*( F(rInteger+1) - F(rInteger) );
                                field.assignElement(j,k,l,magnitude);
                        }
                }
        }
}

void QCSteadyState::usingTDGPGaussSeidel(QCComplexField<2>& currentField, parameterType xOffset, parameterType yOffset, parameterType &h_x, parameterType &h_y, QCParameter &tdgpParameters)
{
        //It is assumed that Small Omega Squared is 1st Parameter, Mu is 2nd Parameter, g the 3rd Parameter, hbar the 4th and m the 5th.
        /*parameterType x, y, omegaSq = tdgpParameters.getParameter(0), mu = tdgpParameters.getParameter(1), g = tdgpParameters.getParameter(2), hbar = tdgpParameters.getParameter(3), m = tdgpParameters.getParameter(4), inverseCoEff = 0.0;
        parameterType alpha = (hbar*hbar)/(8.0*m), dt = (h_x*h_x)/4.0;
        int xExtent = currentField.noOfRows()-1, yExtent = currentField.noOfColumns()-1;
        int maxIterations = static_cast<int> (0.25*10*(xExtent*xExtent));
        complexType Psi;

        setBoundaryVanish(currentField);

        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);
                                Psi =   Psi - alpha*(currentField.retrieveElement(j+1,k)+currentField.retrieveElement(j-1,k)+currentField.retrieveElement(j,k+1)+currentField.retrieveElement(j,k-1)-4.0*Psi)
                        + dt*Psi*(0.5*m*omegaSq*(x*x + y*y)+g*(Psi*Psi)-mu);
                        }
                }
        }*/
}

void QCSteadyState::usingTDGPFastSI(QCComplexField<2> &field, parameterType &h_r, QCParameter &tdgpParameters)
{
        //It is assumed that Small Omega Squared is 1st Parameter, Mu is 2nd Parameter, g the 3rd Parameter, Omega the 6th and n the 7th.
        /*parameterType h_t = 1e-4, omegaSq = tdgpParameters.getParameter(0), mu = tdgpParameters.getParameter(1), g = tdgpParameters.getParameter(2), hbar = tdgpParameters.getParameter(3), m = tdgpParameters.getParameter(4), n = tdgpParameters.getParameter(6);
        parameterType r = 0.0, alpha = (hbar*hbar*h_t)/(4.0*h_r*h_r*m), beta = (hbar*hbar*h_t)/(8.0*h_r*m), sum = 0.0, oldSum = 1.0;
        int radialExtent = field.noOfColumns()-1, maxIterations = 5;
        Array<parameterType,1> a(radialExtent),b(radialExtent),c(radialExtent),d(radialExtent),X(radialExtent);
        QCLinearAlgebra LA;
        parameterType Psi, oldPsi;
        Array<parameterType,1> G;

    //Adjust functional array
        F.free();
        F.resize(radialExtent+1);
        //for(int n = 0; n < radialExtent; n ++)
        //F(n) = 1.0;//(radialExtent - n)/100.0;
    F = 1.0;
        G.resize(radialExtent+1);

    //while(abs(sum - oldSum) > tolerance && radialExtent > 0)
        //{
    for(int n = 0; n < 5; n++)
    {
        //Boundary Conditions for Edges
        F(radialExtent-1) = 0.0;
            G = F;
            oldSum = sum;
            sum = 0.0;
        for(int m = 0; m < maxIterations; m ++)
        {
            for(int j = 1; j < radialExtent; j ++)
            {
                r = (j)*h_r;
                Psi = F(j);
                oldPsi = G(j);

                a(j) = alpha;
                b(j) = 1.0 - 2.0*alpha;
                c(j) = alpha;
                d(j) = G(j) - alpha*G(j+1) + 2.0*alpha*G(j) - alpha*G(j-1)
                        - (beta/r)*(F(j+1) - F(j-1))
                        - (beta/r)*(G(j+1) - G(j-1))
                        + 0.5*h_t*(0.5*m*omegaSq*(r*r) + g*Psi*Psi - mu)*Psi
                        + 0.5*h_t*(0.5*m*omegaSq*(r*r) + g*oldPsi*oldPsi - mu)*oldPsi;
            }
            //Free Boundary Conditions
            b(1) = b(1) + 2.0*a(1);
            c(1) = c(1) - a(1);
            a(radialExtent-1) = a(radialExtent-1) - c(radialExtent-1);
            b(radialExtent-1) = b(radialExtent-1) + 2.0*c(radialExtent-1);
            //Solve Tridiagonal
            LA.solveTridiagonal(a,b,c,X,d);
            for(int j = 1; j < radialExtent; j ++)
            {
                F(j) = X(j);
                sum += F(j);
            }
        }
    }*/
        //}
}

void QCSteadyState::usingTDGPFastSI(QCComplexField<2>& currentField, parameterType xOffset, parameterType yOffset, parameterType &h_x, parameterType &h_y, QCParameter &tdgpParameters)
{
    /*Array<complexType,2> oldField = *currentField.accessArray();
        int maxIterations = 5, 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, h_t, omegaSq = tdgpParameters.getParameter(0), g = tdgpParameters.getParameter(2), hbar = tdgpParameters.getParameter(3), m = tdgpParameters.getParameter(4);
    QCLinearAlgebra LA;

        //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)
                                                   + 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;
                        }
                        //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); //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));
                }
                //Boundary Conditions for Edges
                setBoundaryVanish(currentField);
        }*/
}

void QCSteadyState::transferSolutionTo(QCComplexField<2> &field, parameterType xOffset, parameterType yOffset, parameterType &deltaX, parameterType &deltaY, parameterType &deltaR)
{
        parameterType x, y;
        parameterType l, rSquared;
        parameterType magnitude;
        int rInteger;

        if(completedSolution == TDGPFixedPt)
        {
                for(int k = 0; k < field.noOfColumns(); k++)
                {
                        y = (k)*deltaY - yOffset;
                        for(int j = 0; j < field.noOfRows(); j++)
                        {
                                x = (j)*deltaX - xOffset;
                                rSquared = x*x + y*y;
                                l = sqrt(rSquared)/deltaR;
                                rInteger = int(l);
                                magnitude = F(rInteger) + (l-rInteger)*( F(rInteger+1) - F(rInteger) );
                                field.assignElement(j,k,magnitude);
                        }
                }
        }
        else if(completedSolution == TDGPCartesian)
        {
        }
}

void QCSteadyState::transferSolutionTo(QCComplexField<3> &field, parameterType xOffset, parameterType yOffset, parameterType zOffset, parameterType &deltaX, parameterType &deltaY, parameterType &deltaZ, parameterType &deltaR)
{
        parameterType x, y, z;
        parameterType m, rSquared;
        parameterType magnitude;
        int rInteger;

        if(completedSolution == TDGPFixedPt)
        {
                for(int l = 0; l < field.height(); l ++)
                {
                        z = (l)*deltaZ - zOffset;
                        for(int k = 0; k < field.noOfColumns(); k++)
                        {
                                y = (k)*deltaY - yOffset;
                                for(int j = 0; j < field.noOfRows(); j++)
                                {
                                        x = (j)*deltaX - xOffset;
                                        rSquared = x*x + y*y + z*z;
                                        m = sqrt(rSquared)/deltaR;
                                        rInteger = int(m);
                                        magnitude = F(rInteger) + (m-rInteger)*( F(rInteger+1) - F(rInteger) );
                                        field.assignElement(j,k,l,magnitude);
                                }
                        }
                }
        }
        else if(completedSolution == TDGPCartesian)
        {
        }
}

bool QCSteadyState::output1DSolution(const char* fileName, parameterType &deltaR)
{
        parameterType r;
        ofstream outFile(fileName,ios::out);

        if(outFile.fail())
        {
                cerr << "File Error: " << fileName << endl;
                return false;
        }
        else
        {
                //Plots half the function has its assumed to be symmetric
                for(int j = 0; j < F.extent(firstDim)/2; j++)
                {
                        r = j*deltaR;
                        outFile << r << "," << F(j) << endl;
                }
                outFile.close();
                return true;
        }
}

void QCSteadyState::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,complexZero);
                field.assignElement(j,yExtent,complexZero);
        }
        for(int k = 0; k <= yExtent; k ++)
        {
                field.assignElement(0,k,complexZero);
                field.assignElement(xExtent,k,complexZero);
        }
}

void QCSteadyState::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));
}

---