C (C) Copr. 05/2012 Alexander Moroz
C Wave-scattering.com
C http://www.wave-scattering.com/
C wavescattering@yahoo.com
C
C Last update: May 2012
C
C This is rabif, a computer code to generate data for
C the function F_0(x) introduced in my article
C
C [1] A. Moroz, Comment on ``Integrability of the Rabi model"
C http://arxiv.org/abs/1205.3139
C =====================================================================
C
C This code is made available on the following conditions:
C
C - You can use and modify it, as long as acknowledgment is given
C to the author.
C
C - You can not sell it or use it in any way to get money out of it
C (save for the research grants you might get thanks to it)
C
C - I'm making this code available to the best of my knowledge.
C It has been tested as throughly as possible.
C It work for Windows and Linux platforms. However, it's
C your responsibility to check for its accuracy in the environment
C and size/shape ranges you use.
C
C DISCLAIMER:
C
C WHILE THE COMPUTER PROGRAMS HAVE BEEN TESTED FOR A VARIETY OF CASES,
C IT IS NOT INCONCEIVABLE THAT THEY MAY CONTAIN UNDETECTED ERRORS. ALSO,
C INPUT PARAMETERS CAN BE USED WHICH ARE OUTSIDE THE ENVELOPE OF
C VALUES FOR WHICH RESULTS ARE COMPUTED ACCURATELY. FOR THIS REASON,
C THE AUTHOR DISCLAIM ALL LIABILITY FOR ANY DAMAGES THAT MAY RESULT
C FROM THE USE OF THE PROGRAMS.
C
C I would highly appreciate informing me of any problems encountered
C with this code. Improvements, comments, additions, all welcome at
C
C wavescattering@yahoo.com
C
C =====================================================================
program rabif
C--------/---------/---------/---------/---------/---------/---------/--
C Generates data for the function F_0(x) introduced in my article
C
C [1] A. Moroz, Comment on ``Integrability of the Rabi model"
C http://arxiv.org/abs/1205.3139
C
C by means of the "third method" by Gautschi described on p. 43,
C which makes use of Euler theorem.
C
C Once the respective roots have been bracketed, determines solutions
C of Rabi model by a solving a quantization
C condition (A.16) of Schweber by using routine zbrent (not supplied)
C from Numerical Recipes. If you are in possession of zbrent, you
C can comment out the code lines beginning with "ct" below.
C
C Tested for xg,dt,om= 0.7, 0.4, 1.0 corresponding to Fig 1 of Braak
C
C The relevant output data can be downloaded from
C http://www.wave-scattering.com/rabi.html
C--------/---------/---------/---------/---------/---------/---------/--
implicit none
integer NOUT
real*8 tol
c Parameters:
C ::: number of the output unit
*
PARAMETER (NOUT=35)
*
*
PARAMETER (tol=1.d-7)
*
integer*4 nmax,j
real*8 dt,om,xx,xr,xg,xdel,x0,f0,zbrent
complex*16 zrn,zz,za
external za,f0,zbrent
*
common/prm/ xg,dt,om
*
* Initial length of recurrence:
nmax=200
OPEN(UNIT=NOUT,FILE='f0f.dat')
rewind(NOUT)
write(nout,*)'#F0-function'
5 write(6,*) 'Read in xg,dt,om'
read(5,*) xg,dt,om
ca xg=0.7d0
ca dt=0.4d0
ca om=1.d0
write(nout,*)'#g=',xg
write(nout,*)'#Delta=', dt
write(nout,*)'#omega=', om
write(nout,*)
xdel=0.01d0
xx=-1.d0
C ====
C Here you can determine the root of f0 with a tolerance tol once
C the respective roots have been bracketed, i.e.
C you know that interval (x1,x2) with x1*x2<0 contains a single root.
C You have use either zbrent from your own Numerical Recipes or
C you can replace zbrent with your own root searching routine.
ct xr=zbrent(f0,1.8d0,1.9d0,tol)
ct write(6,*) 'Root=', xr
ct write(20,*) 'Root=', xr
ct go to 10
C ====
C Here you can generate the data to plot f0 as a function of
C the energy parameter xx=epsilon=x/omega, where "epsilon" was used
C by Schweber and "x" by Braak
do j=0,600
xx=xx+xdel
x0=f0(xx)
write(6,*) 'x', xx
write(6,*) 'f0', x0
write(nout,*) xx, x0
enddo
close(nout)
10 continue
c pause
end
real*8 function f0(x)
C--------/---------/---------/---------/---------/---------/---------/--
C Returns a quantization condition (A.16) of Schweber by means
C of the "third method" by Gautschi described on p. 43,
C which makes use of Euler theorem.
C
C Tested for xg,dt,om= 0.7, 0.4, 1.0 corresponding to Fig 1 of Braak
C--------/---------/---------/---------/---------/---------/---------/--
implicit none
integer*4 nmax
real*8 dt,om,x,xg,xx,xg1,dt1,om1
complex*16 zrn,zz,za
external za
*
common/prm/ xg,dt,om
*
common/ttt/ xx,xg1,dt1,om1
*
xx=x
xg1=xg
dt1=dt
om1=om
* Initial length of recurrence:
nmax=200
call rnbcf(nmax,0,zrn)
zz=za(0)
f0=dble(zrn+zz)
return
end
subroutine rnbcf(nmax,n,zrn)
C--------/---------/---------/---------/---------/---------/---------/--
C >>> nmax,n
C <<< zrn
C nmax can be rewritten on the exit via an intermediary nn
C
C rnbcf= rn by continued fraction via Euler theorem:
C For a predetermined n it computes the ratio
C r(n)=y(n+1)/y(n)
C of a minimal solution y(n) of three-term recurrence
C
C y(n+1)+a(n)*y(n)+b(n)*y(n-1)=0,
C
C where b(n).neq.0, by means of a continued fraction
C according to the third method by Gautschi. The method uses
C that (i) r(n) can be related to a continued fraction and
C that (ii) a continued fraction can be related by Euler theorem
C to infinite series.
C
C Calculation is based on recurrence (4.5) described on p. 43 of
C Gautschi.
C r(n) is determined within tolerance specified by TOL
C--------/---------/---------/---------/---------/---------/---------/--
implicit none
real*8 tol
C Relative error allowed - if the convergence
C within TOL is not reached, program issues warning
PARAMETER (TOL=1.d-16)
integer*4 nmax,n,nn,j
complex*16 cone,zrn,za,zb,zw,zv,zu,zw1,zv1,zu1,zaj,zajp1,
& zbj,zbjp1
external za,zb
c DATA PI/3.141592653589793d0/
DATA cone/(1.d0,0.d0)/
*
C Recurrence initialization (first line of 4.5 of Gautschi):
zu=cone
zaj = za(n+1)
zbj = zb(n+1)
zv =-zbj/zaj
zw = zv
C Recurrence (2nd to 4th lines of 4.5 of Gautschi):
nn=0 !used to fine tune nmax in order
!to reach the convergence
10 do j=1,nmax
zajp1 = za(n+j+1+nn)
zbjp1 = zb(n+j+1+nn)
zu1=cone/(cone-zbjp1*zu/(zaj*zajp1))
zv1=zv*(zu1-cone)
zw1=zw+zv1
if (abs(zv1/zw1).lt.tol) then
zrn=zw1
go to 20
end if
zu =zu1
zv =zv1
zw =zw1
zaj =zajp1
enddo
*
write(6,*) 'Recurrence not terminated '
write(6,*) 'NMAX increased to reach convergence'
nn=nmax
go to 10
20 write(6,*) 'Recurrence terminated successfully'
return
c stop
end
* Recurrence coefficients of of three-term recurrence
*
* y(n+1)+a(n)*y(n)+b(n)*y(n-1)=0
*
* are defined as external complex functions za and zb
*
complex*16 function za(j)
integer*4 j
real*8 xx,xg,dt,om
complex*16 ci,cone,zg,zx
common/ttt/ xx,xg,dt,om
DATA ci/(0.d0,1.d0)/,cone/(1.d0,0.d0)/
zg=2.d0*xg
zx=xx-j*om
za=zg/om +(dt**2/zx -zx)/zg !f_n(x) in Braak's (5)
za=-za/dble(j+1)
return
end
*
*
complex*16 function zb(j)
integer*4 j
complex*16 cone
DATA cone/(1.d0,0.d0)/
zb=cone/dble(j+1) !Braak's (4)
return
end