ПРИЛОЖЕНИЕ Ж. РАСПЕЧАТКИ ПРОГРАММ НА ЯЗЫКЕ ФОРТРАН: ЧАСТЬ 2

Ниже приводятся распечатки перевода восьми программ, приведенных в ч. 2 текста, с языка True BASIC на Фортран-77. Программы пропускались иа машине VAX 11/750. Поскольку язык Фортран разрабатывался в те годы, когда машинная графика еще не получила широкого распространения, то стандартный Фортран не содержит никаких графических инструкций и подпрограмм. В качестве примера использования графического пакета в Фортране в наших программах проиллюстрировано применение PLOT 10 (фирма Tektronix)—распространенного пакета, имеющегося на машине VAX и других компьютерах. На большинстве компьютеров PLOT 10 может также вызываться из Паскаля. Краткое описание графических подпрограмм нз пакета PLOT 10, задействованных в наших программах, приведено в табл. П.

ГЛАВА 10

PROGRAM integ

*          вычисляется интеграл от f(x) в пределах от х = а до х = b CALL initial(a,b, h, n)

CALL rectang!e(a, b, h, n, area)

CALL output(area)

STOP

END

SUBROUTINE initia!( a, b,h,n)

•           a - нижний предел интегрирования, b - верхний предел а = 0.0

Ь = 0.5*( 3.14159)

WRITE(6,*) 'число интервалов = '

READ( 5, *) п

h = (Ь - а)/п

RETURN

END

SUBROUTINE rectangle( a, b, h. n, area) x = a sum = 0

DO 100 i = 0,n - 1

sum = sum + f(x)

x = x + h 100 CONTINUE

area - sum*h

RETURN

END

SUBROUTINE output(area) WRITE(6,13) area 13       FORMAT(2x, 'площадь = \F12.7) RETURN END

FUNCTION f(x) f = cos(x) RETURN END

ГЛАВА 11

PROGRAM rwalk

•           моделирование Моите-Карло одномерного случайного блуждания DIMENSION prob(-64:64)

CALL start(p,N,ntrial, iseed) DO 100 itrial = l.ntria! CALL walk( ix, p, N, iseed)

*          сбор данных после N шагов CALL data(ix,xcum, x2cum, prob)

100 CONTINUE

CALL aver( N, ntrial, xcum, x2cum, prob)

STOP

END

SUBROUTINE start( р, N, ntrial, iseed)

WRITE(6,*) 'начальное случайное число (положительное целое) = ' READ( 5,*) iseed

число испытаний ntrial = 100

вероятность шага вправо р = 0.5

WRITE(6,*) 'максимальное число шагов = '

READ(5,*) N

RETURN

END

SUBROUTINE wa!k(ix,p,N, iseed)

N-шаговое блуждание

начальное положение для каждого испытания ix = 0

DO 100 istep = I.N

IF (ran( iseed). le.p) then

ix = ix + 1 ELSE

ix = ix - 1 ENDIF 100 CONTINUE RETURN END

SUBROUTINE data( ix, xcum, x2cum, prob)

DIMENSION prob(-64:64)

xcum = xcum + ix

x2cum = x2cum + ix*ix

prob(ix) = prob(ix) + 1

RETURN

END

SUBROUTINE aver( N, ntrial, xcum, x2cum, prob)*           средние значения для N-шагового блуждания

DIMENSION prob(-64:64)

znorm = 1.0/ntrial

xbar = xcum*znorm

x2bar = x2cum*znorm

DO 10 i = -N,N

prob(i) = prob( i)*znorm WRITE(6,13) i,prob(i) 13 FORMAT(2x,i6,f10.5) 10 CONTINUE

var = x2bar - xbar*xbar

sigma = sqrt(var)

WRITE(6,*) 'среднее смещение = ', xbar

WRITE(6,*) 'sigma = ', sigma

RETURN

END

ГЛАВА 12

PROGRAM percl

DIMENSION lat( 100,100), r( 100,100), icl( 100,100), np( 5000)DIMENSION mass(5000)CALL start(L, p, ntrial, iseed)DO 100 itrial = 1, ntrial35         CALL assign(L,r)

CALL occupy( L, p, lat, r)

CALL cluster( L, lat, icl, np)

проверка вертикального протекания CALL span(L, icl, np, l.ivert)

проверка горизонтального протекания CALL span(L, icl, np, 2, ihori)

IF ((ivert + ihori). eq.0) then

нет протекания ни в одном направлении ispan = 0

ELSEIF ((ivert*ihori).eq.0) then

протекание только в одном направлении ispan = ivert + ihori

IF (ivert.ne.O) spvert = spvert + 1 IF( ihori. ne. 0) sphori = sphori + 1 ELSE

*          соединение в обоих направлениях IF (ivert.eq. ihori) then

iipan = ivert spvert = spvert + 1 sphori = sphori + 1 ELSE

*          два соединяющих кластера —> это испытание пропускаем GO ТО 35

ENDIF ENDIF

вычисление массы каждого кластера CALL iize( L, id, np, mass, labmax)

вычисление среднего размера иесоединяющих кластеров CALL mean( labmax, mass, ispan, clsize, sum)

dtot = dtot + clsize

вычисление вероятности принадлежности соединяющему кластеру CALL Pinf(mass, ispan, sum, zPinf)

zPtot = zPtot + zPinf

вычисление радиуса циркуляции наибольшего несоед. кластера CALL radius( L, id, np, mass, labmax, ispan, radgyr)

radtot = radtot + radgyr 100 CONTINUE

CALL output( ntrial, spvert, sphori,dtot, zPtot, radtot)

STOP

END

SUBROUTINE start( L, p. ntrial, iseed)

WRITE(6,«) 'введите линейный размер решетки'

READ( 5, *) L

WRITE(6,*) 'введите вероятность занятия ячейки' READ{5,*) р

WRITE(6,*) 'enter seed, positive integer' READ(5,*) iseed

WRITE(6,«) 'введите число испытаний'

READ( 5, *) ntrial

RETURN

END

SUBROUTINE assign(L,r)

присваивание   каждой   ячейке   решетки   случайного числа DIMENSION R( 100,100)

DO 10 i = 1.L DO 10 j = 1.L

r{i,j) = ran(iseed) 10 CONTINUE RETURN END

SUBROUTINE occupy(L,p,lat,r)

занятие ячеек решетки DIMENSION Ы( 100,100), r{ 100,100) DO 10 i = 1,L

DO 10 j = I.L

IF (r(i,j).le.p) then

lat(i.j) = -1 ELSE

lat(i,j) = 0 ENDIF 10 CONTINUE RETURN END

SUBROUTINE cluster( L, lat, id, np)

маркировка ячеек

метки хранятся в массиве id, массив правильных меток пр DIMENSION lat( 100,100), icl( 100,100), np( 5000)

DO 10 k = 1,1000 np(k) = 0 10 CONTINUE

DO 15 i = I.L DO 15 j = 1,L icKi.i) = 0 15 CONTINUE nclust = 0

IF (lat(1,1).lt.0) CALL newcl(nclust,id, np, 1,1) DO 20 i = 2,L

IF (lat(i,1).lt.0) then ileft = i - 1

IF (lat(ileft,1).lt.0) then icl(i,1) = icl(ileft,1) ELSE

CALL newcl( nclust, icl, np,i,1) ENDIF ENDIF

20 CONTINUE

DO 30 j = 2,L

IF (lat(1,j).lt.0) then jdown = j - 1 IF (lat(1, jdown). It. 0) then icl(l.j) = icl(1, jdown) ELSE

CALL newcl( nclust, icl, np,1,j) ENDIF ENDIF

DO 40 i = 2,L

IF (lat(i.j).lt.O) then jdown = j - 1 ileft = i - 1

IF (icl(i,jdown) + icl(ileft,j).eq.O) then

CALL newcl( nclust, icl, np,i,j) ELSE

CALL neigh( lat, icl, np, i, j) ENDIF ENDIF

40 CONTINUE 30 CONTINUE

RETURN

END

SUBROUTINE newcl( nclust, icl, np, i, j)*   инициируется новый кластер, указанный в подпрограмме cluster

DIMENSION icl( 100,100), пр{ 5000) nclust = nclust + 1 icl(i,j) = nclust np{ nclust) = 0 RETURN END

SUBROUTINE neigh( lat, icl, np, i, j)

присваивание кпастеру меток соседей DIMENSION lat(100,100), icl( 100,100), np( 5000) jdown = j - 1

Heft = i - 1

IF (lat(i,jdown)*lat(ileft,j).gt.0) then

CALL choice{icl, np, i,j, ileft,jdown)

RETURN ENDIF

IF (icl(i, jdown). gt.0) then

icl(i,j) = icl(i,jdown)

RETURN ENDIF

icl(i.j) = icl( ileft, j)

RETURN

END

SUBROUTINE choice(icl, np, i, j, ileft, jdown)

присваивание кпастеру меток при наличии выбора соседей DIMENSION icl( 100,100), пр( 5000)

IF (icl( ileft,j).eq.icl(i,jdown)) then

icl(i,j) = icl(ileft,j) ELSE

*          соседи имеют разные метки nleft = icl( ileft, j)

ndown = icl(i, jdown) CALL proper(np, nleft) CALL proper( np, ndown) kmax = max0( nleft, ndown) kmin = min0( nleft, ndown) icl(i.j) = kmin

IF (kmax.ne.kmin)np(kmax) = kmin ENDIF RETURN END

SUBROUTINE properfnp, label)

определяется массив правильных меток пр DIMENSION пр{5000)

10        IF (пр{ label). eq.O) RETURN label = np( label) GO TO 10 END

SUBROUTINE span( L, id, np, idir, ispan)

выясняется существование соединяющего кластера проверкой

наличия ячеек с одинаковой правильной меткой на верхней и

нижней (idir = 1 - вертикальное протекание) ипи левой и

правой (idir * 1 - горизонтальное протекание) границах

решетки. На выходе: ispan = метке соединяющего кластера DIMENSION icl(100,100),np(5000)

DO 100 И = 1.L

IF (idir.eq. 1) then

iarg = i1

jarg = 1 ELSE

iarg = 1

jarg = И ENDIF

IF (icl( iarg, jarg). gt.O) then nl = icl( iarg, jarg) CALL proper(np, nl) DO 20 i2 = 1,L

IF (idir.eq.1) then iarg = i2 jarg = L ELSE

iarg = L jarg = i2 ENDIF

IF (icl(iarg,jarg).gt.O) then n2 = icl(iarg, jarg) CALL proper( np, n2) IF (n1.eq.n2) then ispan = nl RETURN ENDIF ENDIF

20 CONTINUE ENDIF

100 CONTINUE ispan = 0 RETURN END

SUBROUTINE size{ L, id, np, mass, labmax) DIMENSION icl( 100,100), np( 5000), mass (5000) DO 10 к = 1,1000 mass(k) = 0 10 CONTINUE labmax = 0 DO 100 i = 1,L DO 100 j = 1,L

IF (icl(i,j).gt.0) then label = icl(i,j) CALL proper( np, label) IF (label, gt. labmax) labmax = label mass( label) = mass( label) + 1 ENDIF 100 CONTINUE RETURN END

SUBROUTINE mean( labmax, mass, ispan, clsize, sum)

вычисление среднего размера несоединяющих кластеров DIMENSION mass(SOOO)

sum = О sum2 = О

DO 10 n = 1, labmax IF (n.ne. ispan) then

sum = sum + mass(n)

sum2 = sum2 + mass(n)**2 ENDIF

10 CONTINUE

clsize = sum2/sum

RETURN

END

SUBROUTINE Pinf( mass, ispan, sum, zPinf)

вычисление вероятности принадлежности соединяющему кластеру DIMENSION mass(SOOO)

IF (ispan. eq.O) then

*          соединяющего кластера нет zPinf =0.0

RETURN ENDIF

zPinf = mass(ispan)/(sum + mass( ispan))

RETURN

END

SUBROUTINE radius( L, id, np, mass, labmax, ispan, radgyr)

расчет радиуса циркуляции наибольшего несоединяющего кластера

сначала находится наибольших несоединяющий кластер DIMENSION icl( 100,100), пр{ 5000), mass( 5000)

maxcl = 0 xcm = 0.0 уст =0.0 r2 = 0.0

DO 50 к = 1,labmax

IF ((mass(k).ne. 0).and.(k.ne. ispan)) then IF (mass(k).gt.maxcl) then maxcl = mass(k) maxk = к ENDIF ENDIF

50 CONTINUE

вычисление r**2 и центра масс DO 100 i = I.L

DO 100 j = I.L

IF (icl(i.j).gt.O) then label = icl(i.j) CALL proper(np, label) IF (label, eq. maxk) then xcm = xcm + i ycm = ycm + j r2 = r2 + i*i + j*j ENDIF ENDIF 100 CONTINUE

вычисляем радиус циркуляции (вычитая центр масс) r2cm = (хст*хст + ycm*ycm)/(maxcl*maxcl) radgyr = sqrt( r2/maxcl - г2ст)

RETURN END

SUBROUTINE output( ntrial, spvert, sphori, cltot, zPtot, radtot)

spvert = spvert/ntrial

sphori = sphori/ntrial

cltot = cltot/ntrial

zPtot = zPtot/ntrial

radtot = radtot/ntrial

WRITE(6,*) 'доля вертикального протекания = '.spvert WRITE(6,*) 'доля горизонтального протекания = '.sphori WRITE(6,*) 'средний размер кластера (иесоедиияющего) = '.cltot WRITE(6,*) 'вероятность принадлежности соед. кластеру = '.zPtot WRITE(6,*) 'радиус циркуляции наибольшего кластера = '.radtot RETURN END

PROGRAM invasion

расчет оккупирующего перколяционного кластера, рисование

кластера, печать относительной занятости и Р(г) DIMENSION г{ 0:100,0:20), регх( 500), регу( 500) CALL setup{ Lx, Ly, iseed)

CALL assign( Lx, Ly, r, perx, pery, iseed)

CALL invade{ perx, pery, r, Lx, Ly)

CALL aver( Lx, Ly, r)

STOP

END

SUBROUTINE sehip(Lx,Ly, iseed)

WRITE(6,*) 'Размер решетки в направлении у = ' READ(5,*) Ly

WRITE(6,*) 'начальное случайное число = ' READ(5,*) iseed Lx = 2*Ly

задание оконных координат CALL initt(1200)

CALL dwindo( 0.0, float( Lx+1), 0.0, float( Ly+2))

RETURN

END

SUBROUTINE assign( Lx, Ly, r, perx, pery, iseed)

присваивание   случайных   чисел   каждой ячейке

занятие первого столбца

ячейки второго столбца образуют начальный лериметр DIMENSION г( 0:100,0:20), регх( 500), регу( 500)

DO 200 j = 1,Ly

занимается первый столбец

г{1.|) = 1 CALL box(l.j)

присваивание случайных чисел остатку j-й строки DO 100 i = 2,Lx

r(i,j) = ran( iseed) 100 CONTINUE

для ячеек периметра r(i,j) устанавливается больше 2 r(2,j) = 2* r(2,j)

sort perimeter sites

CALL sort(perx, pery, r,2, j, j) 200 CONTINUE RETURN END

SUBROUTINE sort(perx,pery,r, iO, jO, nper)

двоичная сортировка: список делится пополам; определяется, к

какой половине относится новое число; этот список тоже делится

пополам и определяется, в какой половине новое число. Процесс

продолжается до определения точного места нового числа DIMENSION г( 0:100,0:20), регх( 500), регу( 500)

если только одна ячейка периметра IF (nper.eq. 1) then

perx(nper) = iO pery(nper) = jO RETURN ENDIF

если новая ячейка меньше всех предыдущих ячеек периметра i = perx(nper-l)

j = pery(nper-1)

IF (r(i0,j0).lt.r(i,j)) then

perx(nper) = iO

pery(nper) = jO

RETURN ENDIF

начало сортировки Шелла, k2 - середина списка, k1,k3 - концы И = 1

k3 = nper - 1 к2 = (kl + кЗ)/2 100 CONTINUE

цикл до отыскания точного положения i = perx(k2)

j = pery(k2)

определение в какой половине списка находится новая ячейка

IF (r(iO,jO).gt.r{i,j)) then

кЗ = к2 ELSE

И = к2 ENDIF

новая середина к2 = (к1 + кЗ)/2

IF ((k1.eq.k2).or.(k2.eq.k3)) then

точное место найдено и равно кЗ

сдвигаем все ячейки выше кЗ на одну позицию вверх DO 10 кк = прег,кЗ+1,-1

perx(kk) = perx(kk-l) pery{kk) = регу(кк-1) 10 CONTINUE

регх(кЗ) = 10 регу(кЗ) = J0 RETURN ENDIF

GO TO 100

RETURN

END

SUBROUTINE sort2( perx, pery, r, iO, j0, nper)

стандартная сортировка методом вставок, элементы perx и pery с

наибольшим аргументом соответствуют наименьшему значению DIMENSION г( 0:100,0:20), регх( 500), регу( 500)

DO 200 к = 1,прег-1 i = регх(к) j = регу(к)

IF (r{i0,j0) .gt. r(i,j)) then

*          вставляем новую ячейку DO 100 kk = nper,k+1,-1

perx(kk) = perx(kk-l) pery(kk) = pery(kk-l) 100 CONTINUE perx(k) = iO pery(k) = JO RETURN ENDIF

200 CONTINUE

новая ячейка меньше всех предыдущих ячеек периметра perx(k) = iO

pery(k) = JO

RETURN

END

SUBROUTINE invade( perx, pery, r, Lx, Ly)

массивы nnx and nny указывают положение ближайших соседей

относительно каждой ячейки

DIMENSION г( 0:100,0:20), регх( 500), регу( 500) DIMENSION nnx(4),nny(4) DATA ппх/1,-1,0,0/ DATA ппу/0,0,1,-1/ nper = Ly

DO 1000 idurnrny = 1,5000 i = perx(nper) j = pery(nper) nper = nper - 1

•           занятые ячейки имеют значения между 1 и 2 r{i.j) = r<i,j) - 1

CALL box(i,j) DO 100 nn = 1,4

находим новые ячейки периметра nx = i + nnx(nn)

ny = j + nny(nn)

периодические краевые условия no у IF (ny .gt. Ly) then

ny = 1 ELSEIF (ny.lt.1) then

ny = Ly ENDIF

новая ячейка периметра

IF (r{nx,ny).lt.1) then

r{nx, ny) = r{nx,ny) + 2 nper = nper + 1

CALL sort( perx, pery, r,nx,ny, nper) ENDIF

100 CONTINUE

•           выход по достижении кластером правой границы IF (i.gt.Lx) RETURN

1000 CONTINUE RETURN END

SUBROUTINE aver{Lx,Ly,r)

нахождение доли занятых ячеек н распределения вероятности

занятия для средней половины решетки DIMENSION г[ 0:100,0:20) DIMENSION р(0:20),п${0:20)

Lmin = Lx/3 Lmax = 2*Lmin

число ячеек в средней половине п = (Lmax - Lmin + 1)*Ly

DO 200 i = Lmin, Lmax DO 100 j = 1,Ly

iarg = 20*(amod(r(i, j),1.0)) ns(iarg) = ns(iarg) + 1 IF ((r(i,j) .gt. 1).and.(r(i,j) .It. 2) ) then occupied = occupied + 1

•           плотность вероятности P(r) p(iarg) = p(iarg) + 1

ENDIF

100 CONTINUE 200 CONTINUE

для вывода результатов в файл по каналу 9 ждем нажатия клавиши fraction = occupied/n

CALL finitt(1,1)

WRITE(9,») 'доля оккупированных ячеек = ',fraction DO 300 i = 0,20 rr = i/20.0

IF (ns(i) .gt. 0) WRITE(9,*) rr, p(i)/ns(i) 300 CONTINUE

RETURN END

SUBROUTINE box(i,j)

вычерчивание прямоугольника, состоящего нз 11 полос хО = i - 0.5

у0 = j - 0.5 DO 10 k = 0,10

у = у0 + k*0.1

CALL movea{x0,Y)

•           проводим прямую относитепьно х0,у CALL drawr(1.0,0.0)

10 CONTINUE CALL tsend RETURN END

ГЛАВА 14

PROGRAM entropy

вычисляется энтропия методом подсчета совпадений Ma DIMENSION rnleft( 10), mright( 10), rnicro( 0:2000)

CALL start( nl, nr, mleft, mright, micro, nexch, iseed)

обмен частнц

CALL exch( nl, nr, nexch, mleft, mright, micro, iseed)

расчет частоты совпадений н энтропии CALL output{ nexch, micro)

STOP END

SUBROUTINE start( nl, nr, mleft, mright, micro, nexch, iseed)

ввод параметров н выбор начальной конфигурации частиц DIMENSION mleft( 10), mright( 10), micro( 0:2000) WRITE(6,*) 'полное число частиц = '

READ( 5, *) N

WRITE(6,*) 'число частиц слева = ' READ( 5, *) nl

WRITE(6,*) 'начальное случайное число (целое положительное) = ' READ( 5, •) iseed

число частиц с правой стороны nr = N - nl

micro(O) = 0 DO 100 il = l.nl

список номеров частиц с левой стороны mleft(il) = il

начальное микросостояние micro(O) = micro(0)*2 + 2

100 CONTINUE

DO 200 ir = 1,nr

список номеров частиц с правой стороны mright(ir) = ir + nl

200 CONTINUE

WRITE(6,*) 'число обменов = '

READ( 5, •) nexch

RETURN

END

SUBROUTINE exch(nl, nr,nexch, mleft, mright, micro, iseed)

обмен частицы слева с номером i left

с частицей справа с номером fright DIMENSION mleft( 10), mright( 10), micro( 0:2000) DO 100 iexch = 1, nexch

индексы массивов выбираем случайным образом i left = int(ran(iseed)*nl + 1)

iright = int(ran(iseed)*nr + 1) jleft = mleft(ileft) jright= mright( iright)

новый номер частицы в левом массиве mlef»(ileft) = jright

новый номер частицы в правом массиве mright( iright) = jleH

определяем новое мнкросостояние micro(iexch) = micro(iexch - 1) + 2**jrighi micro(iexch) = micro(iexch) - 2**jleft

100 CONTINUE RETURN END

SUBROUTINE output(nexch, micro)

вычисляются частота совпадений и энтропия

полное число сравнений DIMENSION micro( 0:2000) ncomp = nexch*(nexch - 1)/2

сравниваем микросостояния DO 200 iexch = l.nexch - 1

DO 100 jexch = iexch + 1,nexch

IF (micro( iexch). eq.micro( jexch)) ncoin = ncoin + 1 100 CONTINUE 200 CONTINUE

частота совпадений

rate = float(ncoin)/float( ncomp) IF (rate.gt.0) S = alog(1.0/rate) WRITE(6,*) 'оценка энтропии = ', S RETURN END

ГЛАВА 15

PROGRAM conduct

алгоритм демона для одномерной модели Изинга

подвод и отвод тепла предусмотрен в концевых узлах DIMENSION s( 1000), edemon( 1000), edcum( 1000), zmgcum( 1000) CALL start( N, nmcs, s, zJ, iseed, nheat)

DO 100 imcs = l.nmcs DO 50 i = 1.N-2

istep = (imcs - 1)*(N-2) + i

IF ( mod( istep, nheat). eq.0) CALL heat( N, zJ, s, edemon) CALL change( N, edemon, s, zJ, accept, iseed) 50 CONTINUE

CALL data( N, s, edemon, zmgcum, edcum) 100 CONTINUE

CALL averg( N, nmcs, zmgcum, edcum, accept)

STOP

END

SUBROUTINE start( N, nmcs, s, zJ, iseed, nheat)

ввод параметров с клавиатуры, запись на устройство 9 DIMENSION s(1000)

WRITE(6,») 'число спинов = ' READ(5,*) N

WRITE(9,*) 'число спинов = ',N

WRITE(6,*) 'число шагов Монте-Карло на спин = ' READ(5,*) nmcs

WRITE(9,*) 'число шагов Монте-Карло на спин = '.nmcs WRITE(6,») 'константа обменного взаимодействия = ' READ( 5, •) zJ

WRITE(9,*) 'константа обменного взаимодействия = ',zJ WRITE(6,*) 'начальное случайное число = ' READ(5,*) iseed

WRITE(9,*) 'начальное случайное число = ',iseed WRITE(6,*) 'время между изменением концевых спинов = ' READ(5,») nheat

WRITE(9,*) 'время между изменением концевых спинов = ',nheat

DO 100 i = 1,N

•           задание случайной начальной конфигурации IF (ran( iseed). gt. 0.5) then

s(i) = 1 ELSE

s(i) = -1 END IF

100 CONTINUE RETURN END

SUBROUTINE change( N, edemon, s, zJ, accept, iseed) DIMENSION s( 1000), edemon( 1000)

вычисляется случайный спин от 2-го до (N-l)-ro ispin = int(ran(iseed)*(N-2) + 2)

дииамнка опрокидывания спина, пробное изменение энергии = de de = 2*zJ*s(ispin)*(s(ispin-1) + s(ispin+1))

IF (de.le.edemon(ispin)) then

s( ispin) = -s( ispin)

accept = accept + 1

edemon( ispin) = edemon( ispin) - de ENDIF RETURN END

SUBROUTINE data( N, s, edemon, zmgcum, edcum)

накопление данных

DIMENSION s( 1000), edemon( 1000), edcum( 1000), zmgcumj 1000) DO 10 i = 1.N

edcum(i) = edcum(i) + edemon(i)

zmgcum(i) = zmgcum(i) + s(i) 10 CONTINUE

edemon(l) = 0 edemon(n) = 0 RETURN END

SUBROUTINE averg( N, nmcs, zmgcum, edcum, accept)

•           нормирование накопленных данных н печать результатов DIMENSION zmgcum( 1000), edcum( 1000)

znorm = 1.0/nmcs

accept = accept*znorm/(N-2)

WRITE(6,*) 'коэффициент принятия = '.accept

WRITE(9,*) 'коэффициент принятия = '.accept

*          величина edcum для спина 1 (ипи N) равна скорости поступления

•           (отвода) тепла

WRITE(6,*) 'скорость подвода тепла = ',edcum(1)*znorm WRITE(6,*) 'скорость отвода тепла = ',edcum(N)*znorm DO 10 i = 2.N-1

edave = edcum( i)*znorm

IF (edave. ne. 0.0) temp = 4.0/alog(1.0 + 4.0/edave) zmag = zmgcum(i)*znorm WRITE(6,13) i, edave, temp, zmag WRITE(9,13) i, edave, temp, zmag 13 FORMAT(1x,i4,3(2x,f13.5)) 10 CONTINUE RETURN END

SUBROUTINE heat(N,zJ,s, edemon)

*          производится попытка добавить (отвести) тепло у спина 1 (N) DIMENSION s( 1000), edemon( 1000)

IF ((s(1)*s(2).gt.0).and.(s(N)*s(N-1).lt.0)) then edemon(l) = edemon(1) + 2*zJ s(1) = -s(1)

edemon(N) = edemon(N) - 2*zJ

s(N) = -s(N) ENDIF RETURN END

ГЛАВА 16

PROGRAM Ising2

алгоритм Метрополиса для двумерной модели Изинга DIMENSION spin(32,32),w(-4:4)

CALL start( N, L, T, nmcs, spin, E, M, w, iseed) DO 100 imcs = 1,nmcs

CALL Melrop{ N, L, spin, E, M, w, ratio, iseed) CALL data( E, M, ecum, e2cum, mcum, m2cum) 100 CONTINUE

CALL output{N, nmcs,ecum, e2cum, mcum, m2cum, ratio)

STOP

END

SUBROUTINE start(N, L,T, nmcs, spin, E,M,w, iseed) DIMENSION spin(32,32),w(-4 :4) WRITE(6,*) 'линейный размер решетки = ' READ'S,*) L

число спинов N = L»L

WRITE(6,*) 'число шагов Монте-Карло на спин = ' READ(5,*) nmcs

WRITE(6,*) 'приведенная температура = ' READ( 5, *) Т

WRITE(6,*) 'начальное случайное число = ' READ{ 5, *) iseed

случайная начальная конфигурация DO 200 j = 1,L

DO 100 i = 1,L

IF (ran(iseed) .It. 0.5) then

spin(i,j) = 1 ELSE

spin(i,j) = -1 ENDIF

•           суммарная намагниченность M = M + spin(i.j)

100 CONTINUE 200 CONTINUE

*          вычисление начальной энергии Е DO 300 j = 1,L

IF (j .eq. L) then

iup = 1 ELSE

iup = j + 1 ENDIF

DO 400 i = 1,L

IF (i .eq. L) then

iright = 1 ELSE

iright = i + I ENDIF

sum = spin(i,iup) + spin( iright, j)

•           полная энергия

E = E - spin( i, j)*sum 400 CONTINUE 300 CONTINUE

•           вычисление вероятностей перехода

*          индекс массива w равен сумме спинов четырех соседей е4 = ехр(-4.0/Т)

е8 = е4*е4

•           аргумент функции ехр = изменению энергии w(4) = е8

w(-4) = е8 w(2) = е4 w(-2) = е4 RETURN END

SUBROUTINE Metrop( N, L, spin, E,M,w, ratio, iseed) DIMENSION spin(32,32),w(-4 :4) DO 100 ispin = 1,N

•           случайный разыгрыш координат спина i = int(L*ran( iseed) + 1)

j = int(L*ran( iseed) + 1)

находим значения соседних спинов

периодические краевые условия CALL periodic(i, j,L, spin, isum)

IF ($pin(i,j)*i$um . le. 0 ) then

CALL accept( i, j, M, E, isum, spin, ratio)

ELSEIF (ran(iseed) .It. w(isum) ) then CALL accept{ i, j, M, E, isum, spin, ratio)

ENDIF

100 CONTINUE RETURN END

SUBROUTINE periodic(i,j,L,spin,isum)

находится   сумма   значений   соседних спинов;

используются периодические краевые условия DIMENSION spin(32,32)

IF (i .eq. 1) then left = spin(L, j) ELSE

left = spin(i - 1,j) ENDIF

IF (i .eq. L) then

iright = spin(1,j) ELSE

iright = spin(i + 1,j) ENDIF

IF (j .eq. 1) then

idown = spin( i, L) ELSE

idown = spin(i,j - 1) ENDIF

IF (j .eq. L) then iup = spin(i,1) ELSE

iup = spin(i,j + 1) ENDIF

isum = left + iright + iup + idown

RETURN

END

SUBROUTINE ассер*{ i, j, М, Е, isum, spin, ratio)

DIMENSION spin(32,32)

spin(i.j) = -spin(i,j)

ratio = ratio + 1

M = M + 2*spin(i,j)

E = E - 2*spin{i,j)*isum

RETURN

END

SUBROUTINE data( E, M, ecum, e2cum, mcum, m2cum)

накапливание данных после каждого шага Монте-Карло на спин ecum = ecum + Е

e2cum = e2cum + Е*Е mcum = mcum + M m2cum = m2cum + M*M RETURN END

SUBROUTINE output( N, nmcs, ecum, e2cum, mcum, m2cum, ratio) znorm = 1/float(nmcs*N)

средние на спин ratio = ratio*znorm eave = ecum*znorm e2ave = e2cum*znorm save = mcum*znorm s2ave = m2cum*znorm

WRITE{6,*) 'коэффициент принятия = '.ratio WRITE{6,*) "средняя энергия на спин - '.eave WRITE{6,*) 'средний квадрат энергии на спин = ',e2ave WRITE{6,*) 'средняя намагниченность - ', save WRITE{6,*) 'средний квадрат намагниченности = ',s2ave RETURN END

ГЛАВА 17

PROGRAM wp

график движения квантовомеханического волнового пакета DIMENSION аг{ -100:100), ai( -100:100), prob( -100:100) DIMENSION prl(-100:100,-100:100),pim(-100:100,-100:100) CALL param( zkO, width, xO, xmin, h, dk, dt, tmax, m, N)

CALL packet( width, xO, m, dk, ar, ai)

CALL free{m,N,zkO,dk,h,ar,ai,prl,pim) CALL barrie( m, N, zkO, dk, h. ar, ai, prl, pim) ntime = tmax/dt

DO 10 itime = 0, ntime t = itime*dt

CALL move(t, m, N, zkO, dk, h, prl, pim, prob) CALL plot{ prob, xmi n, xmax, N, h) 10 CONTINUE

конец рисования, выставление курсора в позицию 1,1 CALL finitt(1,1)

STOP END

SUBROUTINE param( zkO, width, xO, xmin, h, dk, dt, tmax, m, N) WRITE{6,*) 'максимальное значение x = ' READ{ 5, *) xmax

WRITE{6,») 'минимальное значение x = '

READ{ 5, *) xmin

WRITE{6,*) 'размер шага = '

READ( 5, *) h

N = количество вычисляемых кординат N = (xmax - xmin)/h

WRITE(6,*) 'центральный волновой вектор = ' READ( 5, *) zkO

WRITE{6,*) 'число волновых векторов = ' READ{5,») m

WRITE(6,*) 'ширина волнового пакета в х-пространстве = ' READ{5,») width

WRITE(6,*) 'начальная координата пакета = ' READ{ 5, *) xO

dk = 2.0/{width*m) m = m/2 N = N/2

WRITE{6,*) 'шаг no времени = ' READ{ 5, *) dt

WRITE{6,*) 'полное время = ' READ{ 5, *) tmax

WRITE{6,*) 'максимальное значение квадрата psi = READ( 5, *) ymax CALL initt{1200)

CALL dwindo{xmin,xmax,-0.1,Ymax)

RETURN

END

SUBROUTINE packet( width, xO, m, dk, ar, ai)*    вычисляются коэффициенты Фурье

DIMENSION ar{-100:100), ai{-100:100) sum2 = 0.0 DO 100 ik = -m,m dzk = ik*dk

a = exp{-dzk*dzk*width*width/4.0) ar{ik) = a*co${-dzk*x0) ai(ik) = a*sin(-dzk*x0)

$um2 = $um2 + ar{ik)*ar(ik) + ai(ik)*ai(ik) 100 CONTINUE

znorm = dk*$qrt(1.0/($um2)) DO 200 ik = -m,m

ar(ik) = ar(ik)*znorm ai(ik) = ai{ik)*znorm 200 CONTINUE RETURN END

SUBROUTINE free( m, N, zkO, dk, h, ar, ai, prl, pirn) DIMENSION prl{ -100:100, -100:100), pim( -100:100, -100:100) DIMENSION ar{-100:100),ai(-100:100) DO 100 ik = -m,m zk = zkO + ik*dk DO 50 ix = -N,N arg = zk*ix*dx fr = cos(arg) fi = sin(arg)

prl(ik, ix) = fr*ar(ik) - (i*ai(ik) pim(ik, ix) = fr*ai(ik) + fi*ar{ik) 50 CONTINUE 100 CONTINUE RETURN END

SUBROUTINE move{ », m, N, zkO, dk, h, prl, pirn, prob)DIMENSION prl{-100:100,-100:100), pim{-100:100,-100:100)DIMENSION prob( -100:100)DO 100 ix = -N, Npsir =0.0psii = 0.0DO 10 ik = -m,mzk = zkO + ik*dkarg = -0.5*zk*zk*tfr = cos(arg)fi = sin(arg)*           вещественная и мнимая части psi

psir = piir + fr*prl{ ik, ix) - fi*pim(ik, ix) psii = psii + fi*prl( ik, ix) + fr*pim(ik, ix) 10 CONTINUE

prob(ix) = psir*psir + psii* psii 100 CONTINUE RETURN END

SUBROUTINE barrie( m, N, zkO, dk, h, ar, ai, prl, pirn) DIMENSION ar{-100:100),ai{-100:100)

DIMENSION prl{-100:100,-100:100), pim{-100:100,-100:100)

вычисление a(k)(exp(ik1x) + b exp(-iklx)) для x > 0

и a(k)*c*exp(ik2x) для x > 0 WRITE(6,*) 'высота барьера = ' READ{5,*) VO

DO 100 ik = -m, m zk1 = zkO + ik*dk E = 0.5*zk1*zk1 IF (E.gt.VO) then

zk2 = $qrt{2.0*(E - VO))

br = (zk1 - zk2)/{zk1 + zk2)

bi = 0.0

cr = 2.0*zk1/{zk1 + zk2) ci = 0.0 ELSE

zk2 = $qrt{2.0*{V0 - E)) denom = zk1*zk1 + zk2*zk2 br = (zk1*zk1 - zk2*zk2)/denom bi = 2.0*zk1*zk2/denom cr = 2.0*zk1*zk1/denom ci = -2.0*zk1*zk2/denom ENDIF

DO 40 ix = -Nx,0 arg = zk1*ix*h cs = cos(arg) si = sin(arg)

*          далее использовано sin(-kx) =-sin(kx) и cos(-kx) = cos(kx) fr = cs + br*cs + bi*si

fi = si + bi*cs - br*si prl(ik, ix) = fr*ar(ik) - fi*ai(ik) pim(ik, ix) = fr*ai(ik) + fi*ar(ik) 40 CONTINUE

DO 60 ix = 1,N arg = zk2*ix*h IF (E.gt.VO) then

cs = cos(arg)

si = sin(arg) ELSE

cs = exp(-arg)

si = 0.0 ENDIF

(r - cr*cs - ci*si fi = ci *cs + cr*si pri(ik, ix) = fr*ar(ik) - fi*ai(ik) pim(ik, ix) = fr*ai(ik) + (i*ar(ik) 60 CONTINUE 100 CONTINUE RETURN END

SUBROUTINE plot(prob,xmin,N,h)

стирание экрана CALL newpag

нанесение делений DO 10 i = -N,N

x = i*h

CALL movea(x,0.0) CALL drawa(x.0.01*h) 10 CONTINUE

проведение горизонтальной прямой, после цикла х = xmax CALL movea(xmin, 0.0)

CALL drawa(x,0.0)

построение графика квадрата волновой функции CALL movea{-N*h,prob(-N))

DO 20 i = -N, N x = i*h

CALL drawa(x,prob(i)) 20 CONTINUE RETURN END