c======================================================================= c history: twsplt.f c twcmap: This program generates a map between successive c collisions with the infinite wall in the driven c triangle well system. This system has a Hamiltonian c H = p**2/(2*m) - V0*x - V*x*cos(w*t+theta) x < 600 c Because of the presence of the infinite wall the program c do2cje from the NAG library is used. This program integrates c the equations of motion until the particle hits the wall. c Values of J (action) and phi (phase) are written to standard out. c The momentum is then reversed and the equations of motion c are again integrated until another collision, etc. c======================================================================= program twcmap implicit none c----------------------------------------------------------------------- c n: number of odes c iwork: size of work array used by d02ejf c----------------------------------------------------------------------- integer n, iwork parameter (n = 2, iwork = 28 + 21*n) c----------------------------------------------------------------------- c tbgn: initial time of integration (overwritten by do2cje) c tend: final time c tol: tolerance for d02cje c----------------------------------------------------------------------- real tbgn, tend, tdt, dt parameter (dt = 1.0e-4) real tol, default_tol parameter (default_tol = 1.0e-8) c----------------------------------------------------------------------- c i: index used to increment output times c ifail: error return from do2cje c----------------------------------------------------------------------- integer ifail c----------------------------------------------------------------------- c work: work array used by d02cje c y: solution vector c g: function that monitors collisions with the wall c J: canonical action variable c phi: phase of driving field c----------------------------------------------------------------------- real work(iwork), y(n) real g, j, phi c----------------------------------------------------------------------- c Jmax: maximum value of J used as initial condition c Jmin: minimum value of J " c nphi: number of initial values of phi c nJ: " J c niter: number of strobe times (approx) c----------------------------------------------------------------------- real jmax, jmin integer nphi, nj, niter integer iphi, ij, iter c----------------------------------------------------------------------- c argument parsing variables c----------------------------------------------------------------------- real r8arg, r8_never parameter (r8_never = -1.0e-60) integer iargc c----------------------------------------------------------------------- c fcn: subroutine that calculates deriviatives c pederv: " " " Jacobian c out: subroutine that outputs solution at specified times c----------------------------------------------------------------------- external g, d02cje, fcn, cjwd02, cjxd02 c----------------------------------------------------------------------- c Common communication with subroutines c----------------------------------------------------------------------- include 'fcn.inc' real pi c----------------------------------------------------------------------- c parse command line arguments c----------------------------------------------------------------------- if (iargc() .lt. 2) go to 900 w = r8arg(1, r8_never) v = r8arg(2, r8_never) if (w .eq. r8_never .or. v .eq. r8_never) go to 900 tol = r8arg(3, default_tol) if (tol .le. 0.0e0) tol = default_tol c----------------------------------------------------------------------- c Set values of parameters c----------------------------------------------------------------------- v0 = 4.0e-1 m = 1.0e0 niter = 100 pi = 4.0e0*atan(1.0e0) per = 2.0e0*pi/w tend = 1.0e3*niter*per l = 6.0e2 c----------------------------------------------------------------------- c set iteration parameters c----------------------------------------------------------------------- jmax = 80 jmin = 50 nphi = 2 nj = 30 c----------------------------------------------------------------------- c iterate over initial conditions c----------------------------------------------------------------------- do iphi = 1, nphi do ij = 1, nj c----------------------------------------------------------------------- c initialize values of theta and J c----------------------------------------------------------------------- theta = 0.0e0+real(iphi-1)*pi/real(nphi-1) j = jmin + real(ij-1)*(jmax-jmin)/real(nj-1) c----------------------------------------------------------------------- c Initialize solution at t = 0 (y(1) = x, y(2) = p) c----------------------------------------------------------------------- tbgn = 0.0e0 y(1) = l y(2) = (3.0e0*pi*v0*m*j)**(1.0e0/3.0e0) c----------------------------------------------------------------------- c begin loop over collisions with wall c----------------------------------------------------------------------- do iter = 1, niter c----------------------------------------------------------------------- c Call d02ejf to integrate y until collsion with wall c----------------------------------------------------------------------- ifail = 0 call d02cje(tbgn,tend,n,y,fcn,tol,'D',cjxd02,g, & work,ifail) c----------------------------------------------------------------------- c Error message if d02cje fails c----------------------------------------------------------------------- if (ifail .gt. 0) then write(0,*) 'd02cje returns ifail = ',ifail stop end if c----------------------------------------------------------------------- c Reverse momentum and return to d02cje unless at the end c----------------------------------------------------------------------- j = y(2)*y(2)*y(2)/(3.0e0*pi*v0*m) phi = mod(w*tbgn + theta, 2.0e0*pi) write(*,*) phi, j y(2) = -y(2) ifail = 0 tdt = tbgn + dt call d02cje(tbgn,tdt,n,y,fcn,tol,'D',cjxd02,cjwd02, & work,ifail) end do end do end do stop c----------------------------------------------------------------------- c usage message c----------------------------------------------------------------------- 900 continue write(0,*) 'usage: twsplt []' write(0,*) write(0,*) 'w: frequency of driving field' write(0,*) 'V: strength of driving field' write(0,*) 'tol: optional tolerance for d02cje' stop end