CC BY
17ФИЗИКО-МАТЕМАТИЧЕСКИЕ НАУКИ
Уточненный расчет интегралов столкновений. Программа на фортране Попов В. Н.
Попов Валентин Николаевич /Popov Valentin Nikolaevich - кандидат технических наук,
старший научный сотрудник, Объединенный институт высоких температур РАН, г. Москва
Аннотация: расчеты интегралов столкновения (ИС) являются одним из основных элементов процедуры определения транспортных свойств. Использована модифицированная программа Тейлора для расчета классических (ИС) [1, 2]. Эта программа была введена в вычислительный комплекс BASEPA, программирование проведено в среде Builder 6 на языке C++[3-5]. В этой статье использован FORTRAN 77, транслятор Visual Fortran 6.5 compaq.
Ключевые слова: интеграл столкновений, потенциал, сечение рассеяния, квадратура.
УДК: 546.212 + 533.16
Введение
Нейтральные газы при относительно невысоких давлениях широко распространены в природе, являются рабочими веществами и теплоносителями многих технологических процессов. В рамках столкновений атомов и молекул остаются в основном состоянии с замкнутыми электронными оболочками. При расчете транспортных свойств чистых газов и их смесей с высокой точностью можно применять современный аппарат молекулярно-кинетической теории газов. Современные потенциалы позволяют проводить расчеты транспортных свойств газов с погрешностью порядка 1% [3-5]. Однако, как показал проверочный расчет по потенциалу Леннарда - Джонса (Л-Дж), отклонения по [1] могут достигать 1.8 %. Для исправления были сделаны уточняющие разработки Pr2 и Pr3, что дало отклонение 0.3 и 0.02%. При расчетах по потенциалу Морзе [4] выяснилось, что при ß>5 расчеты (ИС не идут. Для исправления в [4] введена уточняющая разработка Pr4. Эти уточняющие разработки введены в вычислительный комплекс BASEPA на языке C++ [3, 5]. В данной статье все разработки сделаны на языке фортран.
Алгоритм расчета интегралов столкновения
Используется для расчета классических интегралов столкновений Q
l =1,..4, s = 1,...n для сферически-симметричных потенциалов практически любого вида (отсутствует max на больших расстояниях). Алгоритм предполагает существование первой и второй производных потенциала по координате. Если производные имеют точки разрыва, то необходимо заменить используемый метод решения уравнений (метод Ньютона-Рафсона) на другой, который не требует непрерывности производных. Алгоритм основан на работе [1], см. также [2, 6, 7]. Здесь приводятся только те формулы, в которых произведены уточнения, подробности в [2]. Для расчета интегралов столкновений необходимо вычислить три интеграла [2, 6, 8].
Работа программы начинается с обращения к подпрограмме SCAN(R, R1, R2, R3, RC, EC). Задан потенциал, на выходе: R1 координата нуля потенциала ^, R2 координата минимума потенциала r2, EC критическая энергия E^ и др.
Первый интеграл: угол отклонения % (угол между векторами относительной скорости сталкивающихся частиц до и после столкновения).
Подпрограмма CHIGM(E,B,R,N, CHI), выходной параметр угол отклонения % (CHI).
Все расчеты традиционно проводятся для приведенных величин:
* r * ф * kT „* E
г = —; , ф = —; T =—; E = —,
Г S S S
где Т - абсолютная температура; k - постоянная Больцмана; S - значение потенциала в минимуме при r = r2
Второй интеграл: Приведенное сечение рассеяния.
Подпрограмма CROSS(E,EC,RC,R1, AL,N1,N2, UL,WL, Q) входной параметр E, выходной- 4 элемента массива Q. Для нахождения пределов интегрирования идет обращение к подпрограмме ORBIT (E,RE,ROP,RO,EC,RC,R1,EOLD) выходные параметры: RE нижний предел интегрирования первого интеграла (2) re, ROP верхний предел интегрирования первого интеграла (2) г0 , RO нижний предел интегрирования второго интеграла (2) r0.
Для сечения рассеяния, как универсальной функции приведенной энергии столкновения, прежде всего, определяется диапазон энергий. Он задается интервалом приведенных температур для расчета интегралов столкновений. Интеграл для приведенного сечения столкновения
Q ( ° '=ргг=^rn—JF (Гт) drm (1)
L *тв. сф. L *тв. сф.
где F(rm)=rm [2 (Е - фОт))-rm(p (rm)] • (l -cosl(x)), rm - расстояние
г «и \ 2(и-г) / наибольшего сближения, I О 1 ' I , = --- - величина сечения для модели
Лтв .сф.. 2
твердой сферы [6, 8]. Интеграл (1) для сечения рассеяния при Е < Екр в силу разрыва по ^ распадается на два интеграла; первый имеет нижний предел интегрирования гю , и верхний г00 - левая граница разрыва по второй интеграл берется от г0 (правая граница разрыва, значение ^ , которое при данной Е соответствует прицельному параметру закручивания) до да.
Для численного интегрирования сечения рассеяния используется (1) с переменной интегрирования г ^ Для энергий меньших Екр интеграл, как уже отмечалось, состоит из двух частей:
О® = Р (гт) йгт + (гт) йгт (2)
Для энергий E > Eкр г00 = г0 Интеграл от ге до ценивается с помощью квадратур Гаусса[5] после перехода к пределам [-1,1]
^ Р(Гт)йТт = 1 А? Р(х?) (3)
где х", Л" - узлы и веса квадратур для пределов интегрирования (-1,1). Значения
I"
■j' j
п
х", Л" - задаются из стандартных таблиц при выбранном числе шагов. Для интервала от г0 до да оценка интеграла с помощью квадратуры Гаусса-Лагерра [10]: ¡гу (Гт)йТт = аЦ= 1 В? ехр(у?)Р(у?) (4)
" ТЛ" п
где у у , В у - узлы и веса квадратур; Гт = Г0 + а у у , а - скалярный
фактор, который выбирается так, чтобы вклад в интеграл для ^ большего, чем
Тт = Г0 + ОС у", был несущественным. Оптимальное значение Г соответствует
минимальному числу точек в квадратурах (4).
Полученная таблица сечений столкновений затем интерполируется полинома Эйткена, описывающего зависимость
Третий интеграл: Расчет приведенного интеграла столкновений.
Интеграл считается в квадратурах Гаусса-Лагерра [10] по формуле:
>*) = (-F) -(F) S+1Q'00d(E /Г) =
^^"(yrrV0* (5)
y", B" - узлы и веса квадратур, Q(l) приведенные сечения столкновений. Диапазон энергий, для которых рассчитывается сечение столкновений, задается равенствами Е mm =Tminyj, Е max =Tmax y" , Tmm , Tmax - заданные приведенные температуры.
Программа на фортране по описанному алгоритму приведена ниже и обозначена как Pr1.
Основная программа COLL управляет ходом расчета, считывает данные, выводит результаты.
В программе проводится интерполяция сечений рассеяния полиномом Эйткена Программа COLL вычисляет приведенные интегралы столкновений, пригодные для расчета всех транспортных свойств при l,s = (1,1), (1,2), (1,3), (1,4), (1,5), (1,6), (1,7), (2,2), (2,3), (2,4),(2,5), (2,6), (3,3), (3,4), (3,5), (4,4). Этого набора интегралов достаточно для расчета транспортных свойств, включая коэффициенты термодиффузии газовых смесей, по крайней мере во втором приближении теории.
Подпрограмма POTENLJ вычисляет значение потенциала, первую и вторую производные.
ДОПОЛНЕНИЯ К ПРОГРАММЕ
Для повышения точности расчетов и возможности расчетов по потенциалу Morse при больших значениях ß были введены усовершенствование программы.
Программа Pr2 вместо интерполяции по Эйткену проводит прямой расчет по подпрограмме CROSS.
Программа Pr3 вместо (5) использует
äl's) *(Г) = ^С ехР ( - ехР (у) ) (ехР (у) ) s+2Q® ' (Ю d (у) (6)
где y=ln(E/T*). Интегрирование (6) проводится методом Гаусса [10] с переменным числом узлов (до 1000). Практически показано, что пределы интегрирования достаточно установить (-10, +5). Программа Pr4 введена для потенциала Morse, который может быть записан безразмерном виде: Ф*(г) = [ 1 -exp(ß(1 -(1 -ln(2)/ß) ß r))]2-1, (7)
где CT=rj= r2 (1-ln(2)/ß), а - диаметр столкновения (координата, где потенциальная функция проходит через ноль). r2 - равновесное состояние, где производная ф*(г) по r равна нулю, а ф*(г) = -1. Величина ß - параметр, характеризующий ширину потенциальной ямы. Значение ß может быть найдено из спектроскопических исследований [14]:
При малых значениях ß и больших T* нижний предел интегрирования Ге при
Ф*=E может быть не найден, так как линии не пересекаются. В работах [12, 13] эта проблема решается установкой жесткой стенки при малых r
Расчеты по потенциалу Морзе при ß большем 5 по описанной программе останавливают счет, причина скалярный фактор а (4). в работах [12, 13] в эту область
ß не зах°дили. В работе [11] расчеты проведены до ß=20.6931.
Для исправления программы интегрирования по (4) введена процедура:
Интеграл от r0 до r0 +DEL оценивается с помощью квадратур Гаусса[10], далее накапливается сумма и процесс повторяется от r0+DEL до r0 +2*DEL и.т.д. до тех пор, пока очередной интеграл не станет меньше EPS. Эта процедура введена в подпрограмму CROSS и названа Pr4 [4].
При подключении Pr4 расчеты по потенциалу Морзе проводились до ß=20.6931
ОПИСАНИЕ УПРАВЛЯЮЩИХ ЭЛЕМЕНТОВ IP2=0 Печать в ORBIT, IP2=1 нет печати IP3=0 Выдает все 16 интегралов столкновения,
Q*(1,1) Q*(1,2) Q*(1,3) Q*(1,4) Q*(1,5) Q*(1,6) Q*(1,7) Q*(2,2)
Q*(2,3) Q*(2,4) Q* (2,5) Q* (2,6) Q* (3,3) Q* (3,4) Q* (3,5) Q* (4,4)
IP3=1 Выдает только Q*(1,1) и Q* (2,2)
Вид потенциала: NN=0 (Л-Дж) 12-6, NN=1 Л-Дж ANN- AMM, NN =2 Morse
параметр р в программе BETMO
II=0 Ведет расчет по 286 значениям температуры, II=1 расчет ведется по заданным пользователем значениям температур их количество определяет NT. JJ2=0 интерполяция по методу Эйткена Pr1, JJ2=1 прямой расчет Pr2 JJ3=1 Расчет по Pr3, где число узлов NA=500 (до 1000), начало и конец интегрирования X1=-10., X2=5.
JJ4=1 Расчет по Pr4, где DEL=1.D0 шаг и EPS=0.0001 точность выхода из цикла , значения DEL и EPS могут быть заданы любые. ТЕСТИРОВАНИЕ ПРОГРАММЫ
Для контроля в табл. 1 приведены расчёты интегралов столкновения для потенциала (Л-Дж) 12-6.
Таблица 1. Сравнение расчётов интегралов столкновения
i • дм»
ill, Prl ifc, Pr2 CJj, Pi3 £2>, [16, 15] Yi Yi Y3
0.1 • 4.02461 • 4.09886 4.09950 4.09985 1,84 0,024 0,0085
0.3 • 2.80936 • 2.83550 2.84365 2.S4418 1,22 0,305 0,0186
0.7 • 1.87949 1.91758 1.92263 1.92299 2,26 0,281 0,0187
• 1.0 • 1.56197 • 1.59540 • 1.59318 • 1.59334 1,97 0,129 0,0100
• 2.0 • 1.16349 • 1.17639 • 1.17580 • 1.17592 1,06 0,040 0,0102
5.0 • 0.91991 • 0.92685 0.92681 0.92681 0,74 0,004 0
40 • 0.66765 • 0.67170 0.67170 0.67169 0,60 0,001 0,0015
Y1=100(Q1-Q*)/Q*, %
Таблица 2. Сравнение рассчитанных значений Q(1Д) * для потенциала (Л-Дж) (75,6) с
табличными данными.
T* Pr2 Pr3 [17] T* Pr2 Pr3 [17]
0.10 2.75326 2.75351 2.7535 4.00 1.01469 1.01488 1.0153
0.20 2.19242 2.19232 2.1916 10.00 0.95038 0.95043 0.9504
0.40 1.74900 1.74904 1.7485 40.00 0.90056 0.90055 0.9006
0.60 1.53270 1.53320 1.5331 60.00 0.88935 0.88935 0.8894
0.80 1.40168 1.40108 1.4011 100.0 0.87613 0.87613 0.8762
1.00 1.31167 1.31246 1.3129 200.0 0.85920 0.85920 0.8594
Таблица 3. Сравнение результатов расчета *(Т*) по потенциалу Морзе с литературными
данными, при /3=2.6931
T* PR2 PR3 [111 [12, 131
0.01 17.721 17.708 17.706 17.743
0.04 12.523 12.519 12.535 12.540
0.1 9.661 9.6976 9.7379 9.7135
0.4 5.036 5.0463 5.0400 5.0504
1 2.1462 2.1434 2.1403 2.1437
2 1.1801 1.1809 1.1794 1.1812
4 0.7367 0.7365 0.7362 0.7380
Таблица 4. Сравнение результатов расчета по потенциалу Морзе с литературными данными,
при /=20.6931
T* Q22 Pr3,Pr4 Q22[11] T* Q22 Pr3,Pr4 Q22[11]
0.004 1.82858 1.8187 2.000 1.10113 1.0821
0.010 1.70857 1.6986 4.000 1.02150 0.9783
0.020 1.62090 1.6105 10.000 0.95532 0.9359
0.040 1.53617 1.5246 20.000 0.92067 0.9174
0.100 1.42937 1.4171 40.000 0.88980 0.8892
0.200 1.35355 1.3410 100.000 0.85030 0.8502
0.400 1.28266 1.2692 200.000 0.82057 0.8205
1.000 1.18879 1.1741
ПРОГРАММА
PROGRAM COLL
COMMON/P/IP2;COMMON/P1/JJ4;COMMON/GR/DEL,EPS;COMMON/POT/NN
COMMON/POT1/BETMO,ANN,AMM;INTEGER NTMAX;PARAMETER
(NTMAX= 1000) REAL*8
Q(4),T(NTMAX),QL(NTMAX,4),EL(NTMAX),A(NTMAX),OMEGA(16),
* OMEGAS(NTMAX,16),DIF(4),AMIN,AMAX REAL*8 UL(50),WL(50),CC1,CC2,SS1,SS2,DD1,DD2,S1,S2
REAL*8 R,R1,R2,R3,RC,EC,RMIN,EINF, R34,PISS,X,CONST,EMIN,EMAX,
* DL,DU,E,EMINL,EMAXL,ECL, DE,EXL,BETMO,TT,ANN,AMM,DEL,EPS REAL*8 DT,TK, V2,DVDUMMY,DDVDUMMY,TS,XGG2(1001),WG2(1001) INTEGER N,N1,N2,I,IK,NT,J, K,KN,KX,K1,L, NUM,
* IT,ITN,JNT,IS,ISN,LNT,KKK,JIS,NIS,JS INTEGER II,JJ3,IP2,IP3,JJ2, JJ4, NN, MM,I1,I2,I3,NA,M,AL REAL*8 Z1,Z,XM,XG,PP,P3,P2,P1,S,X1,X2,SS,SA,SSA,DD,CC,KK,AI
OPEN (9,FILE='PF.DAT');R=1.0;N1=40;N=N1/2 NN=0;II=1;JJ2=1;JJ3=1;NA=500;X1=-10.;X2=5.;JJ4=1;IP2=1;IP3=1;DEL=1. EPS=0.0001;BETMD=2.6931;ANN=6.D0;AMM=75.D0;T(1)=0.1;T(2)=0.3 T(3)=0.7;T(4)=1.;T(5)=2.;T(6)=5;T(7)=40.;NT=7 IF(2*N .NE. N1) THEN
WRITE (*, '(/1X,''NUMBER OF GAUSS QUADRATURE POINTS '', & ''MUST BE EVEN'', /1X,'' N1='',I4)') N1
WRITE (9, '(/1X,''NUMBER OF GAUSS QUADRATURE POINTS '', & ''MUST BE EVEN'', /1X,'' N1='',I4)') N1 END IF
CALL SCAN(R, R1,R2,R3,RC,EC) RMIN=R2;EINF=1.E+05
IF (R1 .LE. 0D0) THEN
WRITE (*, '(/1X,''NON-ACCEPTABLE VALUE: R1='', E12.4)') R1 WRITE (9, '(/1X,''NON-ACCEPTABLE VALUE: R1='', E12.4)') R1 STOP;END IF IF (II .EQ. 1) GOTO 55 T(1)=0.1 ;DT=0. 05 ;TK=400. 0 DO I=1,NTMAX-1
DT=0.05;IF (T(I).GE. 5.0 -1E-3) DT=0.1 IF (T(I).GE. 10.0 -1E-3) DT=0.2;IF (T(I).GE. 20.0 -1E-3) DT=0.5 IF (T(I).GE. 50.0 -1E-3) DT=1.0;IF (T(I).GE.100.0 -1E-3) DT=2.0 IF (T(I).GE.200.0 -1E-3) DT=5.0;IF (T(I).GE.500.0 -1E-3) DT=10.0 T(I+1)=T(I)+DT
IF ( T(I+1) .GE. TK -1E-6) THEN IK=I+1 EXIT
END IF;END DO\ NT=286
55 UL(1)=0.0590198521D0;WL(1)=0.1428119733D0;UL(2)=0.3112391461D0 WL(2)=0.2587741075D0;UL(3)=0.7660969055D0;WL(3)=0.2588067072D0 UL(4)=1.4255975908D0;WL(4)=0.1833226889D0;UL(5)=2.2925620586D0 WL(5)=0.9816627262D-1;UL(6)=3.3707742642D0;WL(6)=0.4073247815D-1 UL(7)=4.6650837034D0;WL(7)=0.1322601940D-1;UL(8)=6.1815351187D0 WL(8)=0.3369349058D-2;UL(9)=7.9275392471D0;WL(9)=0.6721625640D-3 UL(10)=9.912098015D0;WL(10)=0.104461214D-3;UL(11)=12.146102711D0 WL(11)=0.125447219D-4;UL(12)=14.642732289D0;WL(12)=0.115131581D-5 UL(13)=17.417992646D0;WL(13)=0.796081295D-7;UL(14)=20.491460082D0 WL(14)=0.407285898D-8;UL(15)=23.887329848D0;WL(15)=0.150700822D-9 UL(16)=27.635937174D0;WL(16)=0.391773651D-11;UL(17)=31.776041352D0 WL(17)=0.689418105D-13;UL(18)=36.358405801D0; WL(18)=0.781980038D-15;WL(17)=0.689418105D-13; UL(18)=36.358405801D0;WL(18)=0.781980038D-15
UL(19)=41.451720484D0;WL(19)=0.535018881D-17;UL(20)=47.153106445D0 WL(20)=0.201051746D-19;UL(21)=53.608574544D0; WL(21)=0.360576586D-22;UL(22)=61.058531447D0; WL(22)=0.245181884D-25;UL(23)=69.962240035D0 WL(23)=0.408830159D-29;UL(24)=81.498279233D0; WL(24)=0.557534578D-33;N2=24; CALL POTENLJ (R2,V2,DVDUMMY,DDVDUMMY) WRITE (*, '(2X,''R1='',F10.6,3X,''R2='',F10.6,3X, & ''R3='',F10.6,3X, ''RC='',F10.6,
& /2X,''EC='',F10.6,3X,''V2=EMIN='',F10.6,3X,''EINF='',E12.5)') & R1,R2,R3,RC,EC,V2,EINF WRITE (9, '(2X,''R1='',F10.6,3X,''R2='',F10.6,3X, & ''R3='',F10.6,3X, ''RC='',F10.6,
& /2X,''EC='',F10.6,3X,''V2=EMIN='',F10.6,3X,''EINF='',E12.5)') & R1,R2,R3,RC,EC,V2,EINF KN=40 ;KX=2 *N2 ;R3 4=4. *R3 CALL POTENLJ (R34,PISS,X,X)
CONST=PISS*((4.*R3)**6)*1.D+04;CONST=-CONST;EMIN=T(1)*UL(1)/1.01 EMAX=T(NT)*UL(N2)* 1.01 ;DL=8.0/2.3026;DU=DL ;EMINL=DLOG(EMIN) EMAXL=DLOG(EMAX) ;ECL=DLOG(EC) ;K=1;K1=0 ;J=0 ;DE=DL*(ECL -EMINL)
IF (DE .LE. 0D0) THEN
ECL=EMINL-.69315;GOTO 3;END IF
K1=DE+.5;IF (K1 .LT. 1) K1=1;IF (K1 .GT. KN) K1=KN;DE=K1 DE=(ECL-EMINL)/DE;DE=DEXP(DE);E=EMIN 2 CONTINUE
IF(JJ3.EQ.1) GOTO 78 DO I=K,K1 ! DO 275 IF (E .GT. EINF) E=EINF-EMIN CALL CROSS(E,EC,RC,R1, CONST,N1,N2, UL,WL, Q) DO L=1,4 IF(Q(L) .LE. 0D0) THEN
QL(I,L)=0D0;ELSE;QL(I,L)=DLOG(Q(L));END IF;END DO ! (L=1,4) EL(I)=DLOG(E);E=DE*E;END DO ! (I=K,K1) L#275 IF (E .LE. 1.5 *EC) THEN
DE=((ECL+.69315)-EL(K1))/12.0;E=DEXP(EL(K1)+2.*DE);DE=DEXP(DE) K=K1+1;K1=K1+10;GOTO 2;END IF IF (J .GT. 0D0) GOTO 6 3 DE=DU*(EMAXL-(ECL+.69315)) IF (DE .LT. 0D0) GOTO 6
J=DE+.5;IF(J .GT. KX) J=KX;DE=J;DE=(EMAXL-(ECL+.69315))/DE DE=DEXP(DE);E=DEXP(ECL+.69315);K=K1+1;K1=K1+J+1;GOTO 2 6 CONTINUE DO I=1,NT ! DO 750 TS=T(I) DO L=1,16
OMEGA(L)=0.0;END DO DO K=1,N2 ! DO 700
E=TS*UL(K) ;NUM=5;ISN=NUM/2 ;ITN=(NUM-1 )/2 ;EXL=DLOG(E) DO JNT=1,K1 IF (EXL .LT. EL(JNT)) EXIT END DO
IS=JNT-ISN;IT=JNT+ITN;JNT=NUM DO WHILE (IS .LE. 0) IS=IS+1;JNT=JNT-1 ;END DO DO WHILE (IT .GT. K1) IT=IT-1;JNT=JNT-1 ;END DO LNT=JNT;JNT=JNT-1 DO L=1,4 DO J=1,LNT KKK=IS+J-1;A(J)=QL(KKK,L);END DO DO N=1,JNT DO J=N,JNT JIS=J+IS;NIS=JIS-N
A(J+1)=(A(J)*(EL(JIS)-EXL)-A(J+1)*(EL(NIS)-EXL))
& /(EL(JIS)-EL(NIS))
END DO;END DO
Q(L)=DEXP(A(JNT+1))
END DO ! (L=1,4)
IF(JJ2.EQ.1) CALL CROSS(E,EC,RC,R1, CONST,N1,N2, UL,WL, Q) JS=2
DO L=1,7
OMEGA(L)=OMEGA(L)+WL(K)*(UL(K)**JS)*Q(1);JS=JS+1 ;END DO JS=3
DO L=8,12
OMEGA(L)=OMEGA(L)+WL(K)*(UL(K)**JS)*Q(2);JS=JS+1;END DO JS=4
DO L=13,15
OMEGA(L)=OMEGA(L)+WL(K)*(UL(K)**JS)*Q(3);JS=JS+1 ;END DO OMEGA( 16)=OMEGA( 16)+WL(K)*(UL(K) **5)*Q(4) END DO ! (K=1,N2) L#700 JS=2
DO J=1,15 ! DO 725
N=1
DO L=2,JS N=L*N;END DO E=N;OMEGAS(I,J)=OMEGA(J)/E IF (J .EQ. 7) JS=2 IF (J .EQ. 12) JS=3 JS=JS+1
END DO ! (J=1,15) L#725 OMEGAS(1,16)=OMEGA( 16)/120.0 END DO ! (I=1,NT) L#750
WRITE (9,'(/1X,''TABLE OF DIMENSIONLESS OMEGA'' )' ) WRITE (9,'(1X,''FOR TEST CASE OF SIMPLE L.-J. (6-12)'',
* '' POTENTIAL 4*[(1/R)**12 - (1/R)**6]'' )' ) WRITE (*,'(1X,''N1='',I3,'' N2='',I3,'' NT='',I3 )' ) & N1,N2,NT
WRITE (9,'(1X,''N1='',I3,'' N2='',I3,'' NT='',I3 )' ) & N1,N2,NT IF(IP3.EQ.1) GOTO 57 WRITE (9, '(/3X,''T* '', \ )' ) DO L=1,4 DO JS=L,8-L
WRITE (9,'(3X,''('', I1, '', '', I1, '') '',1X, \ )' ) L,JS END DO END DO
WRITE (9, '( 1X )' ) GOTO 58
57 WRITE (9, '(10X,''T*'',10X,''(1,1)'',13X,''(2,2)'')' )
58 DO I=1,NT
IF(IP3.EQ.0)WRITE (9,'(1X,F8.4,1X,16F12.5)') T(I),
* (OMEGAS(I,L),L=1,16) IF(IP3.EQ.1)WRITE (9,'(1X,F8.4,1X,2F12.5)') T(I),
* OMEGAS(I,1), OMEGAS(I,8) END DO
GOTO 77 78 IF(IP3.EQ.1) GOTO 157 WRITE (9, '(/3X,''T* '', \ )' ) DO L=1,4 DO JS=L,8-L
WRITE (9,'(3X,''('', I1, '', '', I1, '') '',1X, \ )' ) L,JS END DO;END DO WRITE (9, '( 1X )' ) GOTO 158
157 WRITE (9, '(13X,''T*'',10X,''(1,1)'',13X,''(2,2)'')' )
158 A1=2;M=(NA+1)/2;XM=0.5*(X2+X1);XG=0.5*(X2-X1) DO I=1,M
Z=COS(3.141592654*(I-1.+0.75)/(NA+0.5)) DO KK=1,1100
P1=1.0;P2=0.0 DO J=1,NA
P3=P2;P2=P1;P1=((2.0*(J-1)+1)*Z*P2-(J-1)*P3)/(J) END DO
PP=NA*(Z*P1-P2)/(Z*Z-1.0);Z1=Z;Z=Z1-P1/PP IF (ABS(Z-Z1).GT.0.000000001) GOTO 44 END DO
44 XGG2(I)=XM-XG*Z;XGG2(NA-1 -I+2)=XM+XG*Z
WG2(I)=2.0*XG/((1.0-Z*Z)*PP*PP);WG2(NA-1-I+2)=WG2(I) END DO
DO K=1,NT DO L=1,16 OMEGA(L)=0.0 END DO DO I=1,NA-1 E=DEXP(XGG2(I))*T(K) CALL CROSS(E,EC,RC,R1,CONST,40,24, UL,WL, Q) JS=3
DD=DEXP(XGG2(I)) DO L=1,7
OMEGA(L)=OMEGA(L)+WG2(I)*(DD**JS)*Q(1)*DEXP(-DD);JS=JS+1;END
DO
JS=4
DO L=8,12
OMEGA(L)=OMEGA(L)+WG2(I)*(DD**JS)*Q(2)*DEXP(-DD);JS=JS+1;END
DO
JS=5
DO L=13,15
OMEGA(L)=OMEGA(L)+WG2(I)*(DD**JS)*Q(3)*DEXP(-DD);JS=JS+1 ;END DO
OMEGA( 16)=OMEGA( 16)+WG2(I)*(DD**6)*Q(4)*DEXP(-DD);END DO JS=2
DO J=1,15 N=1
DO L=2,JS N=L*N;END DO E=N;OMEGA(J)=OMEGA(J)/E IF (J .EQ. 7) JS=2;IF (J .EQ. 12) JS=3;JS=JS+1 ;END DO OMEGA( 16)=OMEGA( 16)/120. IF(IP3.EQ.0)WRITE (9,'(1X,F8.4,1X,16F12.5)') T(K),
* (OMEGA(L),L=1,16)
IF(IP3.EQ.1)WRITE (9,'(1X,F8.4,1X,2F12.5)') T(K),
* OMEGA(1), OMEGA(8) END DO
77 CONTINUE CLOSE(9) STOP END
SUBROUTINE POTENLJ(R,V,DV,DDV) COMMON/POT/NN;COMMON/POT1/BETMO,ANN,AMM REAL*8 R,V,DV,DDV,X,ANN,AMM,CC,QZ,QY,BETMO,AA
13
INTEGER::JJ2, JJ4, NN, MM X=R
IF(NN.EQ.0) THEN
V=4.0D0*((1.0D0/X)**12-(1.0D0/X)**6) DV=-24.0D0*(2*(1. 0D0/X) **13-(1.0D0/X) * *7) DDV=24.0D0*(26*(1.0D0/X)**14-7*(1.0D0/X)**8) END IF
IF(NN.EQ. 1) THEN
CC=(AMM/ANN)**(ANN/(AMM-ANN))*AMM/(AMM-ANN) V=CC*((1.0D0/X)**AMM-(1.0D0/X)**ANN)
DV=-CC*(AMM*(1.0D0/X)**(AMM+1.D0)-ANN*(1.0D0/X)**(ANN+1.D0)) DDV=CC*(AMM*(AMM+1.D0)*(1.0D0/X)**(AMM+2.D0)-ANN*(ANN+1.D0) **(1.0D0/X)**(ANN+2.D0)) END IF
IF(NN.EQ.2) THEN X=R*BETMO
CC=1.0D0-DLOG(2.0D0)/BETMO QZ=1.0D0- DEXP (BETMO -X*CC) QY=- DEXP ( BETMO-X*CC) V=(QZ*QZ-1.0D0) DV=-2.0D0*QZ*QY*CC*BETMO
DDV=2.0D0*CC*CC*QY*(1.0D0+2.0D0*QY)*BETMO*BETMO
END IF
RETURN
END
SUBROUTINE SCAN(R, R1,R2,R3,RC,EC) IMPLICIT NONE REAL*8 R1 ,R2,R3 ,RC,EC,R
REAL*8 PSI,DPSI,DDPSI, CORR,X,EMIN,U,DX, V,DV,DDV,Y, DR, & EINF, RINF, EMAX, RMAX INTEGER I IF (R .LE. 0D0) THEN WRITE (*, '(/1X, ''FUNCTION SCAN: '', & ''NON-ACCEPTABLE ARGUMENT R ='', E12.4)') R
STOP;END IF DO WHILE (.TRUE.) CALL POTENLJ(R,P SI,DPSI,DDP SI) IF (PSI .LT. 0D0) THEN IF (DPSI .LT. 0D0) THEN X=R-0.06*R;GOTO 16;END IF IF (DDPSI .LT. 0D0) THEN R=0.5*R;ELSE;R=0.75*R;END IF CYCLE END IF ! (PSI .LT. 0D0) IF (DDPSI .LE. 0D0) R=R+0.1*R
14 CORR=-PSI/DPSI
IF (DABS(CORR/R) .LE. 0.00001) EXIT
15 X=R+CORR
IF(X .LE. 0D0) X=R-0.06*R
16 CALL POTENLJ(X,PSI,DPSI,DDPSI) IF (DPSI .LE. 0D0) THEN
R=X;GOTO 14;ELSE;CORR=0. 5 * CORR;GOTO 15;END IF;END DO R1=R+CORR;R=R1;X=R;CORR=0.5*R
14
DO WHILE (.TRUE.) CALL POTENLJ(X,PSI,DPSI,DDPSI) IF(DDPSI .LE. 0D0) THEN CORR=0.5*CORR;GOTO 2;END IF R=X;CORR=-DPSI/DDPSI IF( DABS(CORR/R) .LE. 0.00001) EXIT 2 X=R+CORR
IF (X .LE. R1) X=R-0.06*(DABS (DPSI)/DPSI)*R END DO
R2=R;EMIN=P SI;U= 1.0D0 ;X=R2 ;DX=R2 -R1 DO WHILE (.TRUE.)
X=X+DX
CALL POTENLJ(X,V,DV,DDV) IF(U*DDV .EQ. 0D0) EXIT IF(U*DDV .GT. 0D0) CYCLE U=-U;DX=-0.5D0*DX IF(DABS(DX/X) .LE. 1.D-04) EXIT END DO
R3=X;U=1.D0;X=R3;DX=R3-R2 DO WHILE (.TRUE.)
X=X+DX
CALL POTENLJ(X,V,DV,DDV) Y=X*DDV+3.0*DV IF(U*Y .EQ. 0D0) EXIT IF(U*Y .GT. 0D0) CYCLE U=-U;DX=-.5*DX IF(DABS(DX/X) .LE. 1.D-05) EXIT END DO
RC=X;EC=V+0.5*RC*DV;X=-1.0D5*EMIN;DR=0.01*R1;R=R1 DO I=1,98
CALL POTENLJ(R,V,DV,DDV) IF (DDV .LT. 0D0) GOTO 5 IF (V .GE. X) EXIT R=R-DR;END DO
EINF=0.0;RINF=0.0;EMAX=V;RMAX=R+DR WRITE (*, '(1X, ''FUNCTION SCAN: '',
& ''POTENTIAL HAS NO SPURIOUS MAXIMUM FOR E LESS THAN'', & E14.7, 2X, ''AND R GREATER THAN'', F12.6)') EMAX,RMAX WRITE (9, '(1X, ''FUNCTION SCAN: '',
& ''POTENTIAL HAS NO SPURIOUS MAXIMUM FOR E LESS THAN'', & E14.7, 2X, ''AND R GREATER THAN'', F12.6)') EMAX,RMAX RETURN 5 CONTINUE Y=DDV;X=R;R=R+DR DO WHILE (.TRUE.) CALL POTENLJ(R,V,DV,DDV) CORR=DDV*(R-X)/(DDV-Y) IF (DABS(CORR/R) .LT. 0.0001) EXIT Y=DDV;X=R;R=R-CORR;END DO RINF=R;EINF=V;DR=0.1+DR;DV=-1.0 DO WHILE (DV .LT. 0D0) R=R+DR
CALL POTENLJ(R,V,DV,DDV)
END DO
DO WHILE (.TRUE.) CORR=DV/DDV
IF (DABS(CORR/R) .LE. 0.0001) EXIT R=R-CORR
CALL POTENLJ(R,V,DV,DDV) END DO
RMAX=R;EMAX=V
WRITE (*, '(1X, ''FUNCTION SCAN: '',
& ''POTENTIAL HAS SPURIOUS MAXIMUM'',
* /1X,''RINF='',F11.6,2X, ''EINF='',E14.7, 2X,
* ''RMAX='',F11.6,2X, ''EMAX='',E14.7)' ) RINF,EINF,RMAX,EMAX WRITE (9, '(1X, ''FUNCTION SCAN: '',
& ''POTENTIAL HAS SPURIOUS MAXIMUM'',
* /1X,''RINF='',F11.6,2X, ''EINF='',E14.7, 2X,
* ''RMAX='',F11.6,2X, ''EMAX='',E14.7)' ) RINF,EINF,RMAX,EMAX RETURN
END
SUBROUTINE ORBIT(E,EC,RC,R 1, RE,ROP,RO) COMMON/P/IP2
REAL*8 E,RE,ROP,RO,EC,RC,R1,AG,BETMO,ANN,AMM INTEGER J1,J2,J3,J4,K,AH
REAL*8 ER1,ER2,ER3, BIAS,R,DR, V,DV,DDV, CORR,X,DX,U,Y,DY,VO ER1=1.0D-06;ER2=1.0D-06;ER3=1.0D-06;BIAS=1.*ER1;J1=0 J2=0;J3=0 ;R=R1 ;DR=R1/100. DO J4=1,100 R=R-DR;CALL POTENLJ(R,V,DV,DDV) IF(E .LE. V) THEN IF (DDV .LE. 0D0) R=R+.5*DR EXIT;END IF IF(DDV .LE. 0D0) EXIT END DO RE=R
CALL POTENLJ(RE,V,DV,DDV) DO WHILE (J1 .LT. 30) CORR=(E-V)/DV;J 1=J 1+1 12 X=RE+CORR CALL POTENLJ(X,V,DV,DDV) IF (DV .GE. 0D0) THEN CORR=.5*CORR;GOTO 12 END IF ! (DV .GE. 0D0) IF (DABS(CORR/RE) .LE. ER1) GOTO 18 RE=X
END DO ! (J1 .LT. 30)
IF (DABS(CORR/RE) .LE. 100.0*ER1) GOTO 19
IF(IP2.EQ.0)WRITE(*,'(1X,''ORBIT:NO CONVERGENCE,RE='', E15.8)')RE IF(IP2.EQ.0)WRITE(9,'(1X,''ORBIT:NO CONVERGENCE,RE='', E15.8)')RE RE=X-CORR;BIAS=10. *ER 1 18 RE=RE+CORR IF (E .LT. 0.5*EC) THEN DO WHILE (.TRUE.) RE=RE+BIAS*RE;CALL POTENLJ(RE,V,DV,DDV) IF(E .GE. V) EXIT
END DO
END IF ! (E .LT. 0.5*EC)
IF( (RE.LE.0D0) .OR. (RE.GE.(R1+.0001*R1)) ) THEN IF(IP2.EQ.0) WRITE (*,'(1X, ''ORBIT: ERRORS IN INPUT DATA'')') IF(IP2.EQ.0) WRITE (9,'(1X, ''ORBIT: ERRORS IN INPUT DATA'')') RETURN END IF
IF (E .GE. EC) THEN
IF(IP2.EQ.0) WRITE (*,'(1X,''ORBIT: INTEGRAL CONVERGENCE: RE='',
E15.8,
& '' E='', E12.5)') RE, E IF(IP2.EQ.0) WRITE (9,'(1X,''ORBIT: INTEGRAL CONVERGENCE: RE='',
E15.8,
& '' E='', E12.5)') RE, E
RETURN END IF 19 CONTINUE K=0;U=1. ;X=RC;DX=RC DO WHILE (K .LT. 10) X=X+DX
CALL POTENLJ(X,V,DV,DDV)
Y=2.*(E-V)-X*DV
IF (U*Y .EQ. 0D0) EXIT
IF (U*Y .LT. 0D0) CYCLE
K=K+1;U=-U;DX=-.5*DX
END DO ! (K .LT. 10)
RO=X;R=RO
DO WHILE (J2 .LT. 30)
CALL POTENLJ(R,V,DV,DDV)
Y=2.*(E-V)-R*DV;DY=-3.0*DV-R*DDV;CORR=-Y/DY;R=R+CORR;lJ2=J2+1
IF (RC .GE. R) EXIT
IF (DABS(CORR/R) .LE. ER2) EXIT
END DO ! (J2 .LT. 30)
IF (DABS(CORR/R) .GT. 100.*ER2) THEN
IF(IP2.EQ.0) WRITE(*,'(1X,''ORBIT:NO CONVERGENCE:RO='',E 15.8)')RO IF(IP2.EQ.0) WRITE (9,'(1X,''ORBIT:NO CONVERGENCE:RO='',E 15.8)')RO END IF
K=0;U=1.0;X=R1;DX=(RC-R1)/100 DO WHILE (K .LT. 10)
X=X+DX
CALL POTENLJ(RO,V,DV,DDV)
VO=V;CALL POTENLJ(X,V,DV,DDV)
Y=(E-V) *X**2 -RO**2*(E-VO)
IF (U*Y .EQ. 0D0) EXIT
IF (U*Y .LT. 0D0) CYCLE
K=K+1;U=-U;DX=-.5*DX;END DO
ROP=X;R=ROP
DO WHILE (J3.LT.30)
CALL POTENLJ(R,V,DV,DDV)
Y=R**2*(E-V)-RO**2*(E-VO);DY=2*R*(E-V)-R**2*DV;CORR=-Y/DY
R=R+CORR;J3=J3+1
IF ((R.LT.R1) .OR. (RC.LT.R)) THEN
IF(IP2.EQ.0)WRITE (*,'(1X, ''ORBIT: ERRORS IN INPUT DATA'')')
17
IF(IP2.EQ.0)WRITE (9,'(1X, ''ORBIT: ERRORS IN INPUT DATA'')') RETURN END IF
ROP=R ! IT'S MY ADDITION (IS IT CORRECT?) ITS INFLUENCE IS NEGLIGIBLE
IF (DABS(CORR/R) .LE. ER3) RETURN END DO
IF (DABS(CORR/R) .LE. 100.*ER3) THEN
IF(IP2.EQ.0)WRITE (*,'(1X,''ORBIT:POOR
CONVERGENCE:ROP='',E15. 8)')ROP
IF(IP2.EQ.0)WRITE (9,'(1X,''ORBIT:POOR
CONVERGENCE:ROP='',E15. 8)')ROP ELSE
IF(IP2.EQ.0)WRITE (*,'(1X,''ORBIT:NO CONVERGENCE: ROP='',E15.8)')ROP
IF(IP2.EQ.0)WRITE (9,'(1X,''ORBIT:NO CONVERGENCE: ROP='',E15.8)')ROP
END IF
RETURN
END
SUBROUTINE CROSS(E,EC,RC,R1, AL,N1,N2, UL,WL, Q)
COMMON/P 1/JJ4;COMMON/GR/DEL,EPS
REAL*8 E,EC,RC,R1, AL,DEL,EPS,ZZ,UL(N2),WL(N2)
INTEGER N1,N2,KK,M(3,4),I,L,K, M1,M2, I1,I2
REAL*8 Q(4),QS(4,50),QT(4,50),Q1(4),Q2(4)
REAL*8 PI, RE,ROP,RO, A,RL, X, R, V,DV,DDV,F,DF,
& B, CHI, Y, BDB, Z, SCALE,U(50),W(50),AU
U(1)=0.9982377097D0;U(40)=-U(1);W(1)=0.0045212770D0;W(40)=W(1)
U(2)=0.9907262386D0;U(39)=-U(2);W(2)=0.0104982845D0;W(39)=W(2)
U(3)=0.9772599499D0;U(38)=-U(3);W(3)=0.0164210583D0;W(38)=W(3)
U(4)=0.9579168192D0;U(37)=-U(4);W(4)=0.0222458491D0;W(37)=W(4)
U(5)=0.9328128082D0;U(36)=-U(5);W(5)=0.0279370069D0;W(36)=W(5)
U(6)=0.9020988069D0;U(35)=-U(6);W(6)=0.0334601952D0;W(35)=W(6)
U(7)=0.8659595032D0;U(34)=-U(7);W(7)=0.0387821679D0;W(34)=W(7)
U(8)=0.8246223083D0;U(33)=-U(8);W(8)=0.0438709081D0;W(33)=W(8)
U(9)=0.7783056514D0;U(32)=-U(9);W(9)=0.0486958076D0;W(32)=W(9)
U(10)=0.7273182551D0;U(31)=-U(10);W(10)=0.0532278469D0;W(31)=W(10)
U(11)=0.6719566846D0;U(30)=-U(11);W(11)=0.0574397690D0;W(30)=W(11)
U(12)=0.6125538896D0;U(29)=-U(12);W(12)=0.0613062424D0;W(29)=W(12)
U(13)=0.5494671250D0;U(28)=-U(13);W(13)=0.0648040134D0;W(28)=W(13)
U(14)=0.4830758016D0;U(27)=-U(14);W(14)=0.0679120458D0;W(27)=W(14)
U(15)=0.4137792043D0;U(26)=-U(15);W(15)=0.0706116473D0;W(26)=W(15)
U(16)=0.3419940908D0;U(25)=-U(16);W(16)=0.0728865823D0;W(25)=W(16)
U(17)=0.2681521850D0;U(24)=-U(17);W(17)=0.0747231690D0;W(24)=W(17)
U(18)=0.1926975807D0;U(23)=-U(18);W(18)=0.0761103619D0;W(23)=W(18)
U(19)=0.1160840706D0;U(22)=-U(19);W(19)=0.0770398181D0;W(22)=W(19)
U(20)=0.0387724175D0;U(21)=-U(20);W(20)=0.0775059479D0;W(21)=W(20)
DO I1=1,3
DO I2=1,4
M(I1,I2)=25
END DO
END DO
M(1,2)=25;M( 1,3)=25;M( 1,4)=25;M(2,2)=25;M(2,3)=41 ;M(2,4)=71 M(3,2)=35;M(3,3)=45;M(3,4)=99;PI=DASIN(1.0D0)*2.0D0 DO L=1,4
Q1(L)=0.0;Q2(L)=0.0;END DO K=1;IF (E .LT. 1.25*EC) K=2;IF (E .LT. 0.1*EC) K=3 CALL ORBIT(E,EC,RC,R 1, RE,ROP,RO) IF (E .LE. 2.0*EC) THEN IF(E .GT. EC) THEN ROP=RC;RO=RC;END IF A=ROP-RE;RL=RO DO I=1, N1
X=0.5*(U(I)+1.0);L=2;IF (X .GT. 0.90) L=3;IF (X .GT. 0.95) L=4
M1 =M(K,L) ;R=A*X+RE;CALL POTENLJ (R,V,DV,DDV) F=R*DSQRT( 1.0D0 -V/E) ;DF=R*(2. 0D0*(E-V) -R*DV);B=F CALL CHIGM(E,B,R,M1 ,CHI) Y=DCOS(CHI);Z=Y;BDB=DF DO L=1,4
QS(L,I)=0.5D0*W(I)*BDB*(1.0D0-Y);Y=Y*Z;END DO END DO ! (I=1,N1) DO L=1,4 DO I=1,N1
Q1 (L)=Q 1 (L)+QS(L,I) END DO;END DO DO L=1,4 Q1(L)=A*Q1(L);END DO ELSE ! (E .GT. 2.0*EC) RL=RE
END IF ! (E .LE. 2.0*EC) SCALE=(((AL/E) **(1.0D0/6.0D0)) -RL)/UL(N2) IF(JJ4.EQ.0) THEN IF (SCALE .LE. 0D0) THEN WRITE (*, '(/1X, ''IN FUNCTION CROSS:'', /1X, & ''APPROXIMATION USED TO EVALUATE SCALE FACTOR IN '', & ''GAUSS-LAGUERRE QUADRATURE '', /1X, & ''FOR CROSS SECTIONS IS INVALID FOR ENERGY E ='',E12.4)') E
STOP END IF;END IF M2=10
IF (E .LE. 1.25*EC) THEN M2=40;END IF IF(JJ4.EQ.1) THEN A=DEL;ZZ=RL;ZZ=RL-DEL DO I=1,40 DO L=1,4 QT(L,I)=0.D0 END DO ;END DO DO KK=1,10 ZZ=ZZ+DEL; AU=0. DO I=1,40 X=0.5*(U(I)+1.);R=X*A+ZZ CALL POTENLJ(R,V,DV,DDV) IF(V/E.GT.1.) THEN F=0.1D-4 END IF
F=R*DSQRT( 1. -V/E) ;DF=R* (2.0D0*(E-V)-R*DV) ;B=F CALL CHIGM(E,B,R,M2,CHI)
Y=DCOS(CHI);Z=Y;BDB=DF;AU=AU+W(I)*BDB*(1.0D0-Y) DO L=1,4
QT(L,I)=W(I)*BDB*(1.0D0-Y)/2.D0;Y=Y*Z;END DO IF (CHI .GT. -PI) M2=15;IF (CHI .GT. -PI/2D0) M2=10;END DO ! (I=1,N2) DO L=1,4 DO I=1,40
Q2(L)=Q2(L)+QT(L,I) END DO;END DO IF(AU.LT.EPS) GOTO 34 END DO GOTO 34 END IF
33 DO L=1,4
Q2(L)=0.0;END DO DO I=1,N2 X=UL(I) ;R=X * S CALE+RL CALL POTENLJ(R,V,DV,DDV) F=R*DSQRT( 1 -V/E) ;DF=R*(2.0D0*(E-V) -R*DV);B=F CALL CHIGM(E,B,R,M2,CHI) Y=DCOS(CHI);Z=Y;BDB=SCALE*DF*DEXP(X) DO L=1,4
QT(L,I)=WL(I)*BDB *(1.0D0-Y);Y=Y*Z;END DO IF (CHI .GT. -PI) M2=15;IF (CHI .GT. -PI/2D0) M2=10 END DO ! (I=1,N2) DO L=1,4 DO I=1,N2
Q2(L)=Q2(L)+QT(L,I) END DO;END DO
34 DO L=1,4
Q(L)=Q1(L)+Q2(L);X=1.0D0-0.5D0*((1.0D0+(-1.0D0)**L)/(1.0D0+L)) Q(L)=Q(L)/(X*E) END DO RETURN END
SUBROUTINE CHIGM(E,B ,R,N, CHI) IMPLICIT NONE
REAL*8 PI,XL, SUM,AJ,AJN, U, V,DV,DDV,CHI,E,B,R INTEGER K,M,J,JC,N
PI=DASIN( 1.0D0)*2 ;K=N+1 ;M=K/2 ;XL=4 *N;SUM=0.0D0 ;JC=1 DO J=1,M
AJ=(2*J-1);AJ=DCOS(AJ*PI/XL);AJN=(2*(K-J)-1);AJN=DCOS(AJN*PI/XL) 1 CALL POTENLJ(R/AJ,V,DV,DDV) U=1.0D0-(V/E)-(B*B*AJ*AJ/(R*R)) IF (U.LE.0.0D0) THEN
U=0.0D0;ELSE;U=AJN/DSQRT(U);END IF;SUM=SUM+U IF ( (JC.EQ.0) .OR. (J.EQ.M .AND. K.EQ.2*M) ) THEN JC=1;ELSE;JC=0;U=AJN;AJN=AJ;AJ=U;GOTO 1 END IF
END DO ! (J=1,M)
U=N;CHI=PI*(1.0D0-B*SUM/(U*R))
RETURN
END
Заключение
В работе приведена программа на фортране, в которую введены уточняющие разработки Pr2, Pr3 и Pr4. Тестирование показало, что для потенциала (Л-Дж) отклонения не превысили 0.02 %, Табл. 1 Pr3. Введена программа Pr4 позволяющая рассчитывать (ИС) по потенциалу Морзе [8] при р>5.
Для практического использования с новым потенциалом надо ввести в подпрограмму POTENLJ расчет потенциала и его первой и второй производных.
Литература
1. Taylor W. Algorithms and FORTRAN programs to calculate classical collision integrals for realistic intermolecular potential. 29.11.1979. Mound Facility, USD of E.
2. Калиткин Н. Н., Славинская Н. А., Соколова И. А., «Прецизионный расчет интегралов столкновений от недостаточно гладких потенциалов», Матем. моделирование, 13:5, 2001. С. 97-109.
3. Попов В. Н. Теплофизические свойства ртути на основе модельного потенциала. Теплофизика высоких температур, 2012, том 50. № 5. С. 1-11.
4. Попов В. Н. Приближенный численный метод решения интеграло столкновений для потенциала Морзе в широком интервале параметров, Теплофизика высоких температур, 2013, том 51, № 1. С. 73-78
5. Мешков В. В., Попов В. Н., Фокин Л. Р. Диффузия в смеси газов ртуть-аргон. Журнал физической химии. 2014. Т. 88 №4. С. 586-5
6. Гиршфельдер Дж., Кертисс Ч., Берд Р. Молекулярная теория газов и жидкостей. М.: Изд-во иностр.лит., 1963. С. 933.
7. Славинская Н. А., Соколова И. А., Фокин Л. Р. Интегралы столкновений для потенциала Леннард - Джонса (6-6) // Мат. моделирование, 1998, № 10 (5). С. 3
8. Ферцигер Дж., Капер Г. Математическая теория процессов переноса в газах.М. :ИИЛ, 1976. С. 554.
9. Корн Г., Корн Е. Справочник по математике для научных работников и инженеров. М.: Наука, 1984.
10. Крылов В. И., Шульгина Л. Т. Справочная книга по численному интегрированию. М.: Наука, 1966.
11. Smith F. J. and Munn R. J. Automatic Calculation of the Transport Collision integrals Tables for the Morse Potential / J. Chem. Phys., 1964, 41, 3560-3568.
12. Самуилов Е. В., Цителаури Н. Н. Интегралы столкновения для потенциала Морзе // ТВТ. 1964. Т. 2. №4. С. 565-572.
13. Самуилов Е. В., Цителаури Н. Н. Интегралы столкновения для потенциала морзе // ТВТ. 1969. Т. 7. №1. С. 168-170.
14. Зацепин А. Ф., Благинина Л. А., Кухренко А. И., Пустоваров В. А., Чолах С. О. Нейтронно-индуцированный молекулярный дефект O"2 в ортогерманате бериллия // Физика твердого тела, 2007. Том 49. Вып. 5. С. 798-803.
15. Viehland L. A., Mason E. A., Morrison W. F., Flannery M. R. Tables of transport collision integrals for (n, 6, 4) ion-neutral potentials // Atomic Data/Nucl. Data Tables, 1975. V.16. P. 495-514.
16. Klein M., Hanley H. J. M., Smith F. J., Holland P. Tables of collision integrals and second virial coefficients for the (m, 6, 8) intermolecular potential function/ NSRDS-NBS 47. Wash.: GPO, 1974. P. 151.
17. Tables of collision integrals and second virial coefficients for the (m,6,8) intermolecular potential function /Klein H., Hanley H.I.M., Smith F.I, Holland P. Washington: Gov.Print.Off.,1974. (U.S.Dep. of Commerce.NBS. NSRD - NBS 47). 157
