HOME > natural science Laboratory > コンピュータ・シミュレーション講座 > 計算物理学

一様媒質中における一般解と具体例
1軸ガウシアンによる電子パルスの拡散

文責:遠藤 理平 (2011年10月 7日) カテゴリ:計算物理学(132)

【1-3】1軸ガウシアンによる光パルスと電子パルスの「電子パルスの波数空間の幅と実空間の幅の関係」では、実空間における幅が一定のまま伝搬する光パルスに対して、 電子パルスの実空間における幅は時間と共に拡散することを示しました。 その拡散の様子を見てみます。 前節で定義した1軸パルスの表式は

(1)

で与えられます。 中心波数 k_0 = 0 とすることで、パルスの中心は z=0 で留まるため、拡散の様子に着目することができます。 波数空間と実空間における半値全幅も前節で計算したとおり、

(1)

となります。下の図は t = 0 における波数空間と実空間における分布です。

また、波動関数ψの時間発展を計算した結果が次のとおりです。

電子パルスの拡散の様子

波動関数ψの時間発展を計算した結果が次のとおりです。
青:波動関数の絶対値
赤:波動関数の実部
緑:波動関数の虚部

図のx軸のスケールは ×10^{-11}[m] で、時間間隔は10^{-14}[s] とで計算しています。 時間と共にパルスが広がっている様子がわかります。

C++プログラムソース

/*
シュレディンガー方程式の1軸パルス
(公開:2011/10/07)
*/
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include <fstream>
#include <sstream>
#include <iostream>
#include <string>
#include <cstdio>
#include <iomanip>
#include <stdio.h>
#include <complex>
#if defined(_OPENMP)
  #include <omp.h>
#endif
#if defined(_MSC_VER)
  #include <direct.h>   // Windowsフォルダ作成用
#elif defined(__GNUC__)
  #include <sys/stat.h> //  UNIX系ディレクトリ作成用
#endif
using namespace std;

const double PI = acos(-1.0);
const double e = 2.7182818284590452354;
const complex<double> I = complex<double>(0.0,1.0);
const double c = 2.99792458E+8;
const double lambda0 = 500.0E-9;
const double mu0 = 4.0*PI*1.0E-7;
const double epsilon0 = 1.0/(4.0*PI*c*c)*1.0E+7;
const double h  = 6.6260896 * 1.0E-34;
double hbar = h/(2.0*PI);
const double me = 9.10938215 * 1.0E-31;
const double eV = 1.60217733 * 1.0E-19;

const int N  = 10000;//200 //波数の数
const int Lx=0, Lz=5000, Ld=10;
double dz  = 1.0E-11; //XY平面を描画する際には dz = 1.0E-10
double delta_z = 1.0E-8;
double sigma = 2.0*sqrt(2.0*log(2.0))/(delta_z);
double dk  = 30.0/(delta_z*double(N+1));
const int iz0 = Lz/2 , ix0 = Lx/2;

double E = 0;//1.0*eV;
double k0  = sqrt(2.0*me*E/pow(hbar,2)); 
double omega0  = hbar/(2.0*me) * pow(k0,2);

const int ts = 0, te = 200;                   //<------------ステップ数を指定
//double dt= 2.0*PI/omega0/double(te-ts+1);    //<------------時間間隔を指定
double dt  = 1.0 * 1.0E-14;                   //<------------時間間隔を指定
const double theta_s = 0.0, theta_e  = 0.0;  //<------------角度を指定
const double theta_k =20.0;
const double V0 = 0.0;
string folder = "Schrodinger_pulse1_2";//作成するフォルダ名

char str[200];
int main(){
  #if defined(_MSC_VER)
    _mkdir(folder.c_str());   // Windowsフォルダ作成
  #elif defined(__GNUC__)
    mkdir(folder.c_str(), S_IRWXU | S_IRWXG | S_IRWXO); // UNIX系のディレクトリ作成
  #endif
  #if defined(_OPENMP)
    omp_set_num_threads(8);
    cout << "OpenMPを利用します(最大スレッド数:"<< omp_get_max_threads() << ")" <<  endl;
  #endif

  //ガウス分布の確認
  ofstream ofs( "Gaussian.data" ); 
  for(int jz= 0; jz<=N ; jz++){
    double k = (k0 + dk * double(jz-N/2));
    double Gaussian = exp( - 1.0/2.0 * pow((k-k0)/sigma,2));
    ofs << k/sigma << " " << Gaussian << endl;
  }

  ////////////////////
  //角度ごとに計算
  ////////////////////
  for(double th=theta_s; th <=theta_e; th+=theta_k ){
    double Theta= PI/180.0 * th;     //弧度法へ変換
    char str[200];
    sprintf(str, "/%d", int(th));
    string folder_n = folder + str;
    #if defined(_MSC_VER)
      _mkdir(folder_n.c_str());  // Windowsフォルダ作成
    #elif defined(__GNUC__)
      mkdir(folder_n.c_str(), S_IRWXU | S_IRWXG | S_IRWXO); // UNIX系のディレクトリ作成
    #endif
    #pragma omp parallel
    {
      #pragma omp for //schedule(dynamic, 100)
      for(int tn=ts; tn<=te; tn++){
        double t_real = dt * double(tn);
        cout << th << " " << tn << endl;
        char str[200];
        string str1;
        ofstream fout_s;
        sprintf(str, "/%d/s%d.data",  int(th), tn); str1 = folder + str;
        fout_s.open(str1.c_str());

        for(int iz=0; iz<=Lz ;iz=iz+Ld){
          for(int ix=0 ; ix<=Lx ;ix=ix+Ld){
            double z = dz * double(iz-iz0);
            double x = dz * double(ix-ix0);
            complex<double> Psi   = complex<double>(0.0,0.0);
            complex<double> kz = k0 * cos(Theta);
            complex<double> kx = k0 * sin(Theta);


              for(int jz= 0; jz<=N ; jz++)
              {
                double k = (k0 + dk * double(jz-N/2));
                complex<double> kz = k * cos(Theta);
                complex<double> kx = k * sin(Theta);
                double omega = hbar/(2.0*me) * pow(k,2);
                
                Psi += exp( I*((kz)*(z) + (kx)*(x) -omega * t_real )) * exp( - 1.0/2.0 * pow((k-k0)/sigma,2));
              }
            fout_s  << z/dz  << " " ;
            if(Lx!=0)  fout_s << x/dz << " ";
            fout_s  << Psi.real()/sqrt(double(N)) << " " << Psi.imag()/sqrt(double(N)) << " " << abs(Psi)/sqrt(double(N)) << endl;

          }
          if(Lx!=0)  fout_s<<""<< endl;
        }
        fout_s.close();
      }
    }
  }
}

gnuplot テンプレート

1D.plt
1D-1.plt
1D-2.plt


【目次】シュレディンガー方程式とマクスウェル方程式



▲このページのトップNPO法人 natural science トップ

関連記事

計算物理学







▲このページのトップNPO法人 natural science トップ