c=========================================================== c orbit: Integrates restricted gravitational 2-body c problem. See 'fcn.f' for more details on the equations c of motion. Initial values (specified on command line) c are c c x(0) v_x(0) y(0) v_y(0) c c Output to standard out is c c t x(t) v_x(t) y(t) v_y(t) c c for all requested output times 't' read from standard c input. c=========================================================== program orbit implicit none c----------------------------------------------------------- c Codename: Change this parameter string in accord with c the program name. Use a fixed-length character c string of the appropriate length. c----------------------------------------------------------- character*5 cdnm parameter ( cdnm = 'orbit' ) 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 times. c----------------------------------------------------------- integer maxout parameter ( maxout = 10 000 ) real*8 y0(neq) real*8 vtout(maxout), vyout(maxout,neq), & work(maxout) integer ntout, itout, ntout_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-6 ) c----------------------------------------------------------- c Common communication with routine 'fcn' in 'fcn.f' ... c----------------------------------------------------------- include 'fcn.inc' c---------------------------------------------------------------------- c Initialize parameters defined in common block ... c---------------------------------------------------------------------- G = 1.0d0 M = 1.0d0 c----------------------------------------------------------- c Parse command line arguments (initial values) ... c----------------------------------------------------------- if( iargc() .lt. neq ) then write(0,*) cdnm//': Insufficient initial values '// & ' specified on command line.' go to 900 end if ivs_ok = .true. do ieq = 1 , neq y0(ieq) = r8arg(ieq,r8_never) if( y0(ieq) .eq. r8_never ) then write(0,*) cdnm//': Initial value ', ieq, & ' unrecognizable as a real number.' ivs_ok = .false. end if end do if( .not. ivs_ok ) go to 900 c----------------------------------------------------------- c Get tolerance from command-line if specified, c otherwise use default ... c----------------------------------------------------------- tol = r8arg(neq+1,default_tol) c----------------------------------------------------------- c Echo command line arguments if "local tracing" c enabled ... c----------------------------------------------------------- if( ltrace ) then call dvdump(y0,neq,cdnm//': initial values',0) write(0,*) cdnm//': Tolerance ', tol end if c----------------------------------------------------------- c Get output times from standard input ... c----------------------------------------------------------- call dvfrom('-',vtout,ntout,maxout) if( ntout .le. 0 ) then write(0,*) write(0,*) cdnm//': No output times read from standard input' write(0,*) cdnm//': Use (e.g.) dvmesh or dvgmesh to '// & 'generate output times.' 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( ntout .le. maxdump ) then call dvdump(vtout,ntout,cdnm//': output times',0) else call dvdmp1(vtout,work, & int(1.0d0 * ntout / maxdump + 0.5d0), & ntout,cdnm//': selected output times',0) end if write(0,*) write(0,*) cdnm//': Initial time: ', vtout(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 itout = 2 , ntout istate = 1 c----------------------------------------------------------- c Need these temporaries since lsoda overwrites c tend ... c----------------------------------------------------------- tbgn = vtout(itout-1) tend = vtout(itout) 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, itout, ntout, vtout(itout), & vtout(itout+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, itout, ntout, & vtout(itout), vtout(itout+1) 1500 format(/' ',a,': Error return ',i2, & ' from integrator LSODA.'/ & ' At output time ',i5,' of ',i5/ & ' Interval ',1pe11.3,' .. ', & 1pe11.3/) ntout_succ = itout - 1 go to 500 end if do ieq = 1 , neq vyout(itout,ieq) = y(ieq) end do end do ntout_succ = ntout 500 continue if( neq .gt. 10 ) then write(0,*) '***** WARNING: Output lines will '// & 'be broken !!!! *****' end if do itout = 1 , ntout_succ write(*,2000) vtout(itout), & ( vyout(itout,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 times 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