c=========================================================== c History: sode.f c c Driver routine which integrates ODEs defining c static relativistic boson star. c c See class notes and Colpi et al, Phys Rev Lett, c vol 57, 2485--2488 for more details c=========================================================== program bstars implicit none character*6 cdnm parameter ( cdnm = 'bstars' ) integer iargc, indlnb, i4arg real*8 r8arg real*8 r8_never parameter ( r8_never = -1.0d-60 ) c----------------------------------------------------------- c Order of system. c----------------------------------------------------------- integer neq parameter ( neq = 4 ) c----------------------------------------------------------- c Storage for solution at requested output radii. c----------------------------------------------------------- integer maxout parameter ( maxout = 10 000 ) real*8 y0(neq) real*8 vxout(maxout), vyout(maxout,neq), & work(maxout) integer nxout, ixout, nxout_succ integer ieq logical ltrace parameter ( ltrace = .true. ) integer maxdump parameter ( maxdump = 50 ) logical lsodatrace parameter ( lsodatrace = .false. ) logical ivs_ok c----------------------------------------------------------- c LSODA Variables. c----------------------------------------------------------- external fcn, jac real*8 y(neq) real*8 tbgn, tend integer itol real*8 rtol, atol integer itask, istate, iopt integer lrw parameter ( lrw = 22 + neq * 16 ) real*8 rwork(lrw) integer liw parameter ( liw = 20 + neq ) integer iwork(liw) integer jt real*8 tol real*8 default_tol parameter ( default_tol = 1.0d-8 ) c----------------------------------------------------------- c Common communication with routine 'fcn' in 'fcn.f' ... c----------------------------------------------------------- include 'fcn.inc' c----------------------------------------------------------- c Parse command line arguments. Deviation from c 'sode'. Usage is c c bstars [] c----------------------------------------------------------- if( iargc() .lt. 3 ) go to 900 Lambda = r8arg(1,r8_never) y0(3) = r8arg(2,r8_never) Omega = r8arg(3,r8_never) tol = r8arg(4,default_tol) if( Lambda .eq. r8_never .or. y0(3) .eq. r8_never & .or. Omega .eq. r8_never ) go to 900 Omegasq = Omega ** 2 c----------------------------------------------------------- c Initialize those "boundary" conditions which are c fixed by regularity, or which are arbitrary. c----------------------------------------------------------- y0(1) = 0.0d0 y0(2) = 1.0d0 y0(4) = 0.0d0 c----------------------------------------------------------- c Echo command line arguments if "local tracing" c enabled ... c----------------------------------------------------------- if( ltrace ) then write(0,*) 'Lambda: ', Lambda write(0,*) 'sigma(0): ', y0(3) write(0,*) 'Omega: ', Omega end if c----------------------------------------------------------- c Get output radii from standard input ... c----------------------------------------------------------- call dvfrom('-',vxout,nxout,maxout) if( nxout .le. 0 ) then write(0,*) write(0,*) cdnm//': No output radii read from standard input' write(0,*) cdnm//': Use (e.g.) dvmesh or dvgmesh to '// & 'generate output radii.' write(0,*) write(0,*) cdnm//': Sample usages:' write(0,*) write(0,*) cdnm//': dvmesh 0.0 1.0 101' write(0,*) cdnm//': dvgmesh 0.0 1.0 101 10' write(0,*) stop end if if( ltrace ) then if( nxout .le. maxdump ) then call dvdump(vxout,nxout,cdnm//': output radii',0) else call dvdmp1(vxout,work, & int(1.0d0 * nxout / maxdump + 0.5d0), & nxout,cdnm//': selected output radii',0) end if write(0,*) write(0,*) cdnm//': Initial time: ', vxout(1) write(0,*) end if c----------------------------------------------------------- c Set LSODA parameters ... c----------------------------------------------------------- itol = 1 rtol = tol atol = tol itask = 1 iopt = 0 jt = 2 c----------------------------------------------------------- c Initialize the solution ... c----------------------------------------------------------- do ieq = 1 , neq y(ieq) = y0(ieq) vyout(1,ieq) = y0(ieq) end do c----------------------------------------------------------- c Do the integration ... c----------------------------------------------------------- do ixout = 2 , nxout istate = 1 c----------------------------------------------------------- c Need these temporaries since lsoda overwrites c tend ... c----------------------------------------------------------- tbgn = vxout(ixout-1) tend = vxout(ixout) call lsoda(fcn,neq,y,tbgn,tend, & itol,rtol,atol,itask, & istate,iopt,rwork,lrw,iwork,liw,jac,jt) if( lsodatrace ) then write(0,1000) cdnm, ixout, nxout, vxout(ixout), & vxout(ixout+1) 1000 format(' ',a,': Step ',i4,' t = ',1pe10.3, & ' .. ',1pe10.3) write(0,*) cdnm//': lsoda reurns ', istate end if if( istate .lt. 0 ) then write(0,1500) cdnm, istate, ixout, nxout, & vxout(ixout-1), vxout(ixout) 1500 format(/' ',a,': Error return ',i2, & ' from integrator LSODA.'/ & ' At output time ',i5,' of ',i5/ & ' Interval ',1pe11.3,' .. ', & 1pe11.3/) nxout_succ = ixout - 1 go to 500 end if do ieq = 1 , neq vyout(ixout,ieq) = y(ieq) end do end do nxout_succ = nxout 500 continue do ixout = 1 , nxout_succ write(*,2000) vxout(ixout), & ( vyout(ixout,ieq) , ieq = 1 , neq ) c----------------------------------------------------------- c NOTE: This format will not dump all values on a c single line if neq > 10 !! c----------------------------------------------------------- 2000 format(1P,10E25.16) end do stop 900 continue write(0,*) 'usage: '//cdnm// & ' []' write(0,*) ' ' write(0,*) ' Output radii read from standard input' stop end c=========================================================== c Some double precision vector utility routines. c c Originally from ~matt/vutil/dveclib.f c=========================================================== c----------------------------------------------------------- c Dumps vector labelled with LABEL on UNIT. c----------------------------------------------------------- subroutine dvdump(v,n,label,unit) implicit none real*8 v(*) character*(*) label integer i, n, st, unit if( n .lt. 1 ) go to 130 write(unit,100) label 100 format(/' <<< ',A,' >>>'/) st = 1 110 continue write(unit,120) ( v(i) , i = st , min(st+3,n)) 120 FORMAT(' ',4(1PE19.10)) st = st + 4 if( st .le. n ) go to 110 130 continue return end c----------------------------------------------------------- c Extension of dvdump which dumps every 'inc'th element. c----------------------------------------------------------- subroutine dvdmp1(v,w,inc,n,label,unit) implicit none real*8 v(*), w(*) character*(*) label integer inc, n, unit call dvinj(v,w,inc,n) call dvdump(w,1+(n-1)/inc,label,unit) return end c--------------------------------------------------------- c Injects every inc'th element of v1 into v2 c--------------------------------------------------------- subroutine dvinj(v1,v2,inc,n) implicit none real*8 v1(*), v2(*) integer i, inc, j, n j = 1 do i = 1 , n , inc v2(j) = v1(i) j = j + 1 end do return end