cimprunge.f c
```Index
c    file imprunge.f
c
c    Using Runge Kutta Method
c
c    For Time-Dependent Schroedinger Equation
c    Impact Parameter Dependent
c
c    Written for Stockholm Course    March 1995
c
c    Based on standard routines used by
c    Atomic Physics Group in Bergen
c
c    Collected and tested by L. Kocbach
c
c    (based in part on standard functions and some routines
c         by e.g. J.P. Hansen, IBM Scientific Package (1972)
c         etc.....  )
c
c commons:
complex vc(32,32)
real  akof(32,32),vv(32,32)
c
c     akof(*,*) and vv(*,*) define the matrix elements
c     see vcalc()
c
complex y(32)
dimension vcc(32,32)
dimension vdiag(32)
dimension trafo(32,32)
common /blockb/ trafo
common /blocka/ vcc
common /block4/ vc, akof, vv,velo ,vdiag, rr
common /block5/ itest,iout,nout
common /block6/ bimpac
c
open(11,file='values',status='old')
open(15,file='b-prob')
c
write(*,*) nout , '       nout '
write(*,*) itest, '       itest '
write(*,*) velo,  '       velo '
write(*,*) tstep, '       tstep   dt   '
write(*,*) nstat, '       no of states, '
c
c Example case:  exp(-a(i,j)*t*t)  Read from 'values'
c
do 13 i = 1,nstat
c
do 3 i = 1,nstat
do 4 j = i,nstat
write(*,*) vv(i,j),akof(i,j) ,'  V a ',i,j
vv(j,i) = vv(i,j)
akof(j,i) = akof(i,j)
4        continue
3       continue
c
c      Here we are adding the LOOP over
c      Impact parameters
c
c
do 20000 kb=1,nb
c
bimpac=0.05+deltab*float(kb-1)
c
c      Starting Time-loop
c
do 6 i = 1,nstat
y(i) = cmplx(0.0,0.0)
6       continue
c
y(1) = cmplx(1.0,0.0)
c
tmin=-10.0/velo
t=tmin
tmax=-tmin
c
c      Adjusting the values to avoid BIG FILES
c
c      Time-loop
c
1000  continue
ttt=t
call runge(y,nstat,ttt,tstep)
t=t+tstep
if (t.lt.tmax) goto 1000
call output(nstat,y,bimpac)
c
20000 continue
c              Here ends the impact parameter loop
END
c
subroutine runge(y,mm,t,dt)
c     -----------------------------
c
c performs Runge Kutta 4-order for LINEAR systems
c for TDSE  :  complex amplitudes necesssary
c
c input y(t) , output y(t+dt)
c matrix elements calculated in vcalc
c
integer mm
real t,dt
complex y(32),k1(32),k2(32),k3(32),k4(32)
complex res(32)
c commons:
complex vc(32,32)
real  akof(32,32),vv(32,32)
dimension vdiag(32)
dimension trafo(32,32)
common /blockb/ trafo
common /block4/ vc, akof, vv,velo ,vdiag,rr
common /block6/ bimpac
c
call prdkt(k1,y,dt,mm)
do 10 i = 1,mm
res(i) = y(i)+ k1(i)*0.5
10   continue

call vcalc(mm,t+0.5*dt,dt/2.0)

call prdkt(k2,res,dt,mm)
do 11 i = 1,mm
res(i) = y(i)+ k2(i)*0.5
11   continue

call prdkt(k3,res,dt,mm)
do 12 i = 1,mm
res(i) = y(i)+ k3(i)
12   continue

t = t + dt
call vcalc(mm,t,dt/2.0)
call prdkt(k4,res,dt,mm)
c
c iteration step:
c
do 20 i = 1,mm
y(i) = y(i) + k1(i)/6 + k2(i)/3 + k3(i)/3 + k4(i)/6
20   continue
end
c
subroutine prdkt(ki,y,dt,mm)
c     -----------------------------
c
c returns ki = coup.matrix * y - vector * dt
c
real dt
complex ki(*),y(*)
integer mm
c commons:
complex vc(32,32)
real  akof(32,32),vv(32,32)
dimension vdiag(32)
dimension trafo(32,32)
common /blockb/ trafo
common /block4/ vc, akof, vv,velo ,vdiag, rr
common /block6/ bimpac
do 10 i = 1,mm
ki(i) = cmplx(0.,0.)
do 15 j = 1,mm
ki(i) = ki(i) + vc(i,j)*y(j)
15     continue
ki(i) = ki(i)*dt
10   continue
end
c
subroutine vcalc(nstat,t,dtx)
c     -----------------------------
c
c evaluates coupling matrix and multiplies by i
c
real dtx,drx,t,tt,velo
c commons:
complex vc(32,32)
real  akof(32,32),vv(32,32)
dimension vcc(32,32)
dimension vdiag(32)
dimension trafo(32,32)
common /blockb/ trafo
common /blocka/ vcc
common /block4/ vc, akof, vv,velo ,vdiag,rr
common /block6/ bimpac
c
c Build the matrix
c
c Example case:  exp(-a(i,j)*t*t)
c
tt=t+dtx
c
c
c       originally: rr= velo*tt, now modified
c
rr=sqrt(bimpac*bimpac+velo*tt*velo*tt)
c
do 3 i = 1,nstat
do 4 j = i,nstat
c
vcc(i,j) = vv(i,j)*exp(-akof(i,j)*rr*rr)
c
if (abs(vcc(i,j)) .lt. 0.000001)
.                           vcc(i,j) = 0.000001
vc(i,j) = cmplx(1.0,0.0) * vcc(i,j)
vc(j,i) = vc(i,j)
vcc(j,i) = vcc(i,j)
4        continue
vcc(i,i) = vcc(i,i) + vdiag(i)
vc(i,i) = vc(i,i) + cmplx(vdiag(i),0.0)
3       continue
c
c Multiply matrix by i
c
do 8 i = 1,nstat
do 9 j = 1,nstat
vc(i,j) = cmplx(0.0,-1.0) * vc(i,j)
9        continue
8       continue
return
end
c
subroutine output(nn,yy,bimp)
c     -----------------------------
c
complex yy(32)
c commons:
complex vc(32,32)
real  akof(32,32),vv(32,32), ach(32)
dimension vdiag(32)
dimension trafo(32,32)
common /blockb/ trafo
common /block4/ vc, akof, vv,velo ,vdiag, rr
common /block5/ itest,iout,nout
c
sum=0.0
do 333 k=1,nn
bach=cabs(yy(k))
ach(k) = bach*bach
sum=sum+ach(k)
333      continue
c
write(15,802) bimp,(ach(l),l=1,nn), rr
c
802        format(f10.3,'   ',5f15.6)
c
return
end
c
c
```