numericalcomputmethods.cpp

Schrodinger Equation 数值计算中的方程

#include "stdAfx.h"
#include "nrutil.h"
#include "nrtypes.h"
#include "NumericalComputMethods.h"
//#include <iostream>
#include<cmath>
//using namespace std;
NumericalComputMethods::NumericalComputMethods()
{

}
NumericalComputMethods::~NumericalComputMethods()
{


}
void NumericalComputMethods::flmoon(const int n,const int nph, int &jd,DP &frac)
{
	const DP RAD=3.14;
	int i;
	DP am,as,c, t,t2,xtra;
	c=n+nph/4.0;
	t=c/1236.85;
	t2=t*t;
	as=359.22+29.10*t2;
	jd=2415020+28*n+7*nph;
	xtra=0.75933+1.53*c*((1.178e-4)-(1.55e-7)*t)*t2;
	if(nph==0||nph==2)
		xtra+=(0.1734-3.93e-4*t)*sin(RAD*as)-0.6280*sin(RAD*am);
	//else NR::nrerror("nph is unknownin flmoon");
	i=int(xtra>=0.0?floor(xtra):ceil(xtra-1.0));
	jd+=i;
	frac=xtra-i;
}
void NumericalComputMethods::rot(Mat_IO_DP &a,const DP s, const DP tau,const int i,
		const int j,const int k,const int l)
{
	DP g,h;
	g=a[i][j];
	h=a[k][l];
	a[i][j]=g-s*(h+g*tau);
	a[k][l]=h+s*(g-h*tau);
}
void NumericalComputMethods::determineIndexMiniElem(int& i1,int& i2,int& i3,int& i4,Vec_DP v)
{
	int n=v.size();
	/////////////////////////////////////////////////////////////////////
	/////////////////////////////////////////////////////////////////////
	////////////////determine the first mini-value of v//////////////////
	double mini_1=v[0];
	double max=v[0];
	for(int i=0;i<n;i++)
	{
		max=MAX(max,v[i]);
	}
    for(int i=0;i<n;i++)
	{
      mini_1=MIN(mini_1,v[i]);

	}
	for(int i=0;i<n;i++)
	{
		if(v[i]==mini_1)
			i1=i;
	}
	/////////////////////////////////////////////////////////////////////
	/////////////////////////////////////////////////////////////////////
	////////////////determine the second mini-value of v//////////////////
    double mini_2=max;
    for(int i=0;i<n;i++)
	{
		if(i!=i1)
			mini_2=MIN(mini_2,v[i]);		

	}
	for(int i=0;i<n;i++)
	{
		if(mini_2==v[i])
			i2=i;
	}
	/////////////////////////////////////////////////////////////////////
	/////////////////////////////////////////////////////////////////////
	////////////////determine the third mini-value of v//////////////////
    double mini_3=max;
    for(int i=0;i<n;i++)
	{
		if(i!=i1&&i!=i2)
            mini_3=MIN(mini_3,v[i]);

	}
	for(int i=0;i<n;i++)
	{
		if(mini_3==v[i])
			i3=i;
	}
	/////////////////////////////////////////////////////////////////////
	/////////////////////////////////////////////////////////////////////
	////////////////determine the fourth mini-value of v//////////////////
    double mini_4=max;
    for(int i=0;i<n;i++)
	{
		if(i!=i1&&i!=i2&&i!=i3)
            mini_4=MIN(mini_4,v[i]);

	}
	for(int i=0;i<n;i++)
	{
		if(mini_4==v[i])
			i4=i;
	}


}
void NumericalComputMethods::ProduceMatrix(Mat_IO_DP& M)
{
    double m1=0.067;//the effect mass in the well;
	double m2=0.092;//the effect mass outside well;
	double v_0=0.23*1.6e-19;//the offset of heterojunction.
	double g=0.28e-9;//the grid of the z axis;
	double v_=2*v_0*9.1e-31*g*g/SQR(1.35e-34);
	for(int i=0;i<61;i++)
	{
		for(int j=0;j<61;j++)
		{
            if(j==i)
                M[i][j]=2.0/m2+v_;
			else if(j==(i+1))
				M[i][j]=-1/m2;
			else if(j==(i-1))
				M[i][j]=-1/m2;
			else 
				M[i][j]=0.0;
		}
	}
	M[20][20]=1/m2+1/m1+v_;
	M[20][21]=-1/m1;
	M[21][20]=-1/m1;
	M[19][20]=-1/m2;
	M[20][19]=-1/m2;
	for(int i=21;i<40;i++)
	{
		for(int j=21;j<40;j++)
		{
            if(j==i)
                M[i][j]=2.0/m1;
			else if(j==(i+1))
				M[i][j]=-1/m1;
			else if(j==(i-1))
				M[i][j]=-1/m1;
		}
	}
	M[40][40]=1/m2+1/m1+v_;
	M[40][41]=-1/m2;
	M[41][40]=-1/m2;
	M[39][40]=-1/m1;
	M[40][39]=-1/m1;

	for(int jj=0;jj<61;jj++)
	{
		for(int jjj=0;jjj<61;jjj++)
		{
			cout<<M[jj][jjj]<<' ';
		}
		cout<<endl;
	}
	

}
void NumericalComputMethods::jacobitransformation(Mat_IO_DP &a,Vec_O_DP &d, Mat_O_DP &v,
		int &nrot)
//computes all eigenvalues and eigenvectors of a real symmetric matrix
// a[0..n-1][0..n-1]. On output,elements of a above the diagonal are 
//destroyed. d[0..n-1] returns the eigenvalues of a. v[0..n-1][0..n-1]
//is a matrix whose columns contain, on output, the normalized eigenvectors
//of a. nrot returns the number of Jacobi rotations that were required.

{

	int i, j, ip,iq;
	DP tresh, theta, tau,t, sm, s, h, g, c;
	int n=d.size();
	Vec_DP b(n),z(n);
	for(ip=0;ip<n;ip++)
	{
		for(iq=0;iq<n;iq++)
		{
			v[ip][iq]=0.0;
		}
		v[ip][ip]=1.0;
	}
	for(ip=0;ip<n;ip++)
	{
		b[ip]=d[ip]=a[ip][ip];
		z[ip]=0.0;
	}
	nrot=0;
	for(i=1;i<=50;i++)
	{
		sm=0.0;
		for(ip=0;ip<n-1;ip++)
		{
			for(iq=ip+1;iq<n;iq++)
				sm +=fabs(a[ip][iq]);
		}
		if(sm==0)
		{
			return;//
		}
		if(i<4)
			tresh=0.2*sm/(n*n);
		else
			tresh=0.0;
		for(ip=0;ip<n-1;ip++)
		{
			for(iq=ip+1;iq<n;iq++)
			{
				g=100.0*fabs(a[ip][iq]);
				if(i>4&&(fabs(d[ip])+g)==fabs(d[ip])
					&&(fabs(d[iq])+g)==fabs(d[iq]))
					a[ip][iq]=0.0;
				else if(fabs(a[ip][iq])>tresh)
				{
					h=d[iq]-d[ip];
					if((fabs(h)+g)==fabs(h))
						t=(a[ip][iq])/h;
					else
					{
						theta=0.5*h/(a[ip][iq]);
						t=1.0/(fabs(theta)+sqrt(1.0+theta*theta));
						if(theta<0.0)t=-t;
					}
					c=1.0/sqrt(1+t*t);
					s=t*c;
					tau=s/(1.0+c);
					h=t*a[ip][iq];
					z[ip]-=h;
					z[iq]+=h;
					d[ip]-=h;
					d[iq]+=h;
					a[ip][iq]=0.0;
					for(j=0;j<ip;j++)
						rot(a,s,tau,j,ip,j,iq);
					for(j=ip+1;j<iq;j++)
						rot(a,s,tau,ip,j,j,iq);
					for(j=iq+1;j<n;j++)
						rot(a,s,tau,ip,j,iq,j);
					for(j=0;j<n;j++)
						rot(v,s,tau,j,ip,j,iq);
					++nrot;
				}
			}
		}
		for(ip=0;ip<n;ip++)
		{
			b[ip]+=z[ip];
			d[ip]=b[ip];			
			z[ip]=0.0;

		}

	}
	//NR::nrerror("need more iterations in routine jacobi");

}