FEATool Multiphysics  v1.16.5
Finite Element Analysis Toolbox
ex_schrodinger1.m File Reference

Description

EX_SCHRODINGER1 Schrodinger equation for the Hydrogen atom.

[ FEA, OUT ] = EX_SCHRODINGER1( VARARGIN ) Computation of energy levels and electron orbits for a 1-particle Hydrogen atom system using the Schrodinger equation in cylindrical coordiantes. Accepts the following property/value pairs.

Input       Value/{Default}        Description
-----------------------------------------------------------------------------------
hmax        scalar {0.1}           Grid cell size (in meter scale)
sfun        string {sflag2}        Finite element shape function
iphys       scalar 0/{1}           Use physics mode to define problem    (=1)
                                   or directly define fea.eqn/bdr fields (=0)
                                   or use core assembly functions        (<0)
iplot       scalar 0/{1}           Plot solution (=1)
                                                                                  .
Output      Value/(Size)           Description
-----------------------------------------------------------------------------------
fea         struct                 Problem definition struct
out         struct                 Output struct

Code listing

 cOptDef = { 'hmax',     0.05;
             'sfun',     'sflag2';
             'iplot',    true;
             'tol',      0.02;
             'fid',      1 };
 [got,opt] = parseopt(cOptDef,varargin{:});
 fid       = opt.fid;


% Geometry definition (in meters).
 fea.sdim = {'r', 'z'};
 fea.geom.objects = { gobj_circle([0,0], 3, 'C1'), ...
                      gobj_circle([0,0], 0.5, 'C2'), ...
                      gobj_rectangle(-3, 0, -3, 3, 'R1'), ...
                      gobj_rectangle(-3, 0, -3, 3, 'R2')};
 fea = geom_apply_formula( fea, 'C1-R1' );
 fea = geom_apply_formula( fea, 'C2-R2' );


% Grid generation.
 fea.grid = gridgen( fea, 'hmax', opt.hmax, 'fid', fid );

% Scale geometry and grid to nanometer scale.
 s = 1e-9;
 fea = geom_apply_transformation( fea, fea.geom.tags, 0, s );
 fea.grid.p = fea.grid.p * s;


% Constants and expressions.
 me = 9.10939e-31;     % Electron mass.
 ec = 1.6021773e-19;   % Electron charge.
 hp = 6.626076e-34;    % Planck's constant.
 e0 = 8.8541878e-12;   % Electric capacitivity in vacuum.
 m  = 0;               % Ground state energy.
 fea.expr = { 'me', me;
              'e',  ec;
              'hp', hp;
              'e0', e0;
              'C',  '8*me*pi^2/hp^2';
              'm',  m};


% Problem definition.
 fea = addphys(fea,@poisson);          % Add Poisson equation physics mode.
 fea.phys.poi.sfun = { opt.sfun };     % Set FEM shape function.

 fea.phys.poi.eqn.seqn = ...
 'C*r*u'' - r*(ur_r+uz_z) = 1+(C*r*e^2/(4*pi*e0*sqrt(r^2+z^2))-m^2/r)*u_t';
% fea.phys.poi.eqn.coef{1,end} = { 'C*r' };   % Set time coefficient.
% fea.phys.poi.eqn.coef{2,end} = { 'r' };     % Set diffusion coefficient.
% fea.phys.poi.eqn.coef{3,end} = { '1+(C*r*e^2/(4*pi*e0*sqrt(r^2+z^2))-m^2/r)*u' };       % Define source term coefficient.


% Define boundary conditions (Homogenous Neumann on symmetry axis, and Dirichlet on sphere boundary.)
 fea.phys.poi.bdr.sel = [1 1 2 2 2];


% Parse and solve problem.
 fea = parsephys(fea);
 fea = parseprob(fea);
 [fea.sol.u, fea.sol.l] = solveeig( fea, 'sigma', -2.18e-18, ...
                                    'neigs', 10, 'fid', fid );


% Error checking.
 lam = [-2.18e-18; -5.45e-19; -5.45e-19; -2.422e-19; -2.422e-19];
 out.err = abs((lam - fea.sol.l(1:5))./lam);
 out.pass = all(out.err < opt.tol);


% Postprocessing.
 if( opt.iplot>0 )
   isol = 5;
   postplot( fea, 'surfexpr', 'u', 'solnum', isol )
   title(sprintf('u(lambda(%d)) = %g', isol, fea.sol.l(isol)))
 end

 if( nargout==0 )
   clear fea out
 end