Методические указания к курсу "Атомная спектроскопия". Шунина В.А - 15 стр.

                                 15


Z - заряд ядра; Nocc - кол-во электронов в 1s - состоянии;
[Rmin, Rmax] – отрезок, на котором вычисляется Pnℓ(r);
Pout и Pin - решения (2.15), полученное интегрированием "вперед" и
“назад” соответственно (энергия фиксирована);
Rm - точка сшивки Pout и Pin;
Pout вычисляется на отрезке [Rmin,Rm+del], а Pin – на [Rm-del,Rmax];
[emin,emax] - интервал, на котором ищется собственное значение;
e, e0 - любые значения до начала итераций, для которых |e-e0| > eps;
eps - точность самосогласования по энергии;
delta - шаг сканирования отрезка[emin, emax];
epsilon - точность вычисления интегралов;
InitD[...] - инициализация сетки точек по r и затравочной 1s
функции;
DerMatch[...] - функция для вычисления при фиксированной энергии
Pout,Pin и разности их производных в точке сшивки Rm;
Fmatch[...] - аналог DerMatch[...] с другими формальными параметрами
и типом результирующих данных;
CPotential[...] - функция для вычисления кулоновского члена Ф(r).

ClearAll[Z,Nocc,Rmin,Rmax,Rm,del,R,P,P0,Pout,Pin,x,r,Np,i,e,
e0,emin,emax,en0,en1,en2,energy,delta,eps,epsilon,iter,dash,
Rd,FiM,Fi,num,f1,f2,ren,V1,V2,norm,t0,t1];

InitD[Z_Real,Rmin_Real,Rmax_Real] :=
 Module[{ro1,h,step,num,R,Np,M1s,P0,i},
  (* генерация узлов сетки *)
  ro1 = -3.; h = 0.1; r = Exp[ro1]/Z; step = Exp[h]; num = 0;
  While[r < Rmax,r = r*step; num = num + 1];
  Np = num + 2; r = Exp[ro1]/Z; R = Table[0.,{Np}];
  For[i = 2, i <= Np, i++, R[[i]] = r; r = r*step];
  R[[1]] = Rmin; R[[Np]] = Rmax;
  If[Rmin >= R[[2]],
   Print["ОШИБКА: Rmin должно быть < R[[2]] !"]; Return[13]];
  (* P0 – нормировнная радиальная 1s-функция водорода *)
  M1s = Table[{R[[i]],(2*R[[i]])/E^R[[i]]},{i,Np}];
  P0 = Interpolation[M1s][x]; Return[{R,P0}]];

DerMatch[e_Real,Z_Real,Fi_,Rmin_Real,Rmax_Real,Rm_Real,del_Real] :=
 Module[{Eq,PoutM,Pout,PinM,Pin,h = 0.1,f},
  Eq = 0.5*P’’[x] + (Z/x - 0.5*Fi + e)*P[x] == 0;
  {{Pout}} = NDSolve[{Eq, P’[Rmin] == 1.,P[Rmin] == Rmin}, P[x],
     {x, Rmin, Rm + del}];
  PoutM = P[x] /. Pout /. x -> Rm;
  Pout = (P[x] /. Pout)/PoutM;
  {{Pin}} = NDSolve[{Eq, P’[Rmax] == -0.0001, P[Rmax] == 0.}, P[x],
     {x, Rmin, Rm + del}];
  PinM = P[x] /. Pin /. x -> Rm;
  Pin = (P[x] /. Pin)/PinM;
  f = (Pout /. x -> Rm - h) - (Pin /. x -> Rm - h);
  Return[{f,Pout,Pin}]];

Fmatch[energy_Real] :=