c---------------------------------------------------------------------- c Integrates equations of motion for orbiting dumbbell. c---------------------------------------------------------------------- program dumb implicit none character*2 itoc character*64 catsqz integer iargc, indlnb, vsxynt integer vsrc real*8 x_vs(3), y_vs(3) integer stderr, stdin, stdout parameter ( stderr = 0, stdin = 5, stdout = 6 ) real*8 s2r8 real*8 r8arg integer s2i4 integer i4arg character*128 sarg real*8 r8_never parameter ( r8_never = -1.0d-60 ) integer i4_never parameter ( i4_never = -2 000 000 000 ) integer neq parameter ( neq = 6 ) integer maxout parameter ( maxout = 50 000 ) real*8 y0(neq) real*8 vtout(maxout), vyout(maxout,neq), & work(maxout), ketrans(maxout), & kerot(maxout), pegrav(maxout), etot(maxout) integer ntout, itout, ntout_succ integer ieq logical ltrace parameter ( ltrace = .true. ) integer maxdump parameter ( maxdump = 50 ) logical vstrace parameter ( vstrace = .true. ) logical lsodatrace parameter ( lsodatrace = .false. ) logical ivs_ok c---------------------------------------------------------------------- c LSODA Variables. c---------------------------------------------------------------------- external fcn, jac real*8 y(neq), yprime(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 Hard-code some of the parameters c---------------------------------------------------------------------- G = 1.0d0 MM = 1.0d0 c---------------------------------------------------------------------- c Parse command line arguments (initial values) followed c by m1/m2, d, tol c---------------------------------------------------------------------- if( iargc() .lt. neq ) then write(0,*) 'dumb: 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,*) 'dumb: Initial value ', ieq, ' unrecognizable ', & 'as a real number.' ivs_ok = .false. end if end do if( .not. ivs_ok ) go to 900 m1bym2 = r8arg(neq+1,0.5d0) mu = 1.0d0 / (1.0d0 + m1bym2) d = r8arg(neq+2,0.1d0) tol = r8arg(neq+3,default_tol) m1 = 1.0d0 m2 = m1bym2 d1 = 1.0d0 / (1.0d0 + m1bym2) * d d2 = d - d1 if( ltrace ) then write(0,2000) G, MM, m1bym2, d, d1, d2, tol 2000 format('Gravitational constant: ',F8.2/ & 'Central mass: ',F8.2/ & 'm1/m2: ',F12.5 / & 'Dumbbell separation: ',F12.5/ & 'd1: ',F12.5/ & 'd2: ',F12.5/ & 'LSODA tolerance: ',1P,E10.2) call dvdump(y0,neq,'dumb: initial values',0) 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,*) 'dumb: No output times read from standard in' write(0,*) 'dumb: Use (e.g.) dvmesh or dvgmesh to '// & 'generate output times.' write(0,*) write(0,*) 'dumb: Sample usages:' write(0,*) write(0,*) 'dumb: dvmesh 0.0 1.0 101 | dumb ... ' write(0,*) 'dumb: dvgmesh 0.0 1.0 101 10 | dumb ...' write(0,*) stop end if if( ltrace ) then if( ntout .le. maxdump ) then call dvdump(vtout,ntout,'dumb: output times',0) else call dvdmp1(vtout,work, & int(1.0d0 * ntout / maxdump + 0.5d0), & ntout,'dumb: selected output times',0) end if write(0,*) write(0,*) '***** 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----------------------------------------------------------------------- call dvcopy(y0,y,neq) c---------------------------------------------------------------------- c And call the external routine to initialize the auxiliary c variables. c----------------------------------------------------------------------- call fcn(neq,vtout(1),y,yprime) c---------------------------------------------------------------------- c Do the integration ... c---------------------------------------------------------------------- do itout = 1 , ntout do ieq = 1 , neq vyout(itout,ieq) = y(ieq) end do c---------------------------------------------------------------------- c Compute the total energy ... normalized for m2 = 1, c m1 = m2bym1. c---------------------------------------------------------------------- ketrans(itout) = 0.5d0 * (m1 + m2) * (y(2)**2 + y(4)**2) kerot(itout) = 0.5d0 * (m1 * m2) / (m1 + m2) * & (d * y(6))**2 pegrav(itout) = - G * MM * (m1 / sqrt(x1**2 + y1**2) + & m2 / sqrt(x2**2 + y2**2)) etot(itout) = ketrans(itout) + kerot(itout) + pegrav(itout) if( itout .eq. ntout ) go to 500 istate = 1 c---------------------------------------------------------------------- c Need these temporaries since lsoda overwrites tend ... c---------------------------------------------------------------------- tbgn = vtout(itout) tend = vtout(itout+1) 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) itout, ntout, vtout(itout), vtout(itout+1) 1000 format(' dumb: Step ',i4,' t = ',1pe10.3,' .. ',1pe10.3) write(0,*) 'dumb: lsoda reurns ', istate end if if( istate .lt. 0 ) then write(0,1500) istate, itout, ntout, vtout(itout), & vtout(itout+1) 1500 format(/' dumb: Error return ',i2,' from integrator LSODA.'/ & ' dumb: At output time ',i5,' of ',i5/ & ' dumb: Interval ',1pe11.3,' .. ',1pe11.3/) ntout_succ = itout - 1 go to 500 end if if( vstrace ) then x_vs(1) = x1 x_vs(2) = x2 y_vs(1) = y1 y_vs(2) = y2 vsrc = vsxynt('dumb',tend,x_vs,y_vs,2) end if end do ntout_succ = ntout 500 continue if( vstrace ) then vsrc = vsxynt('th',0.0d0,vtout,vyout(1,5),itout) vsrc = vsxynt('om',0.0d0,vtout,vyout(1,6),itout) vsrc = vsxynt('etot',0.0d0,vtout,etot,itout) vsrc = vsxynt('ketrans',0.0d0,vtout,ketrans,itout) vsrc = vsxynt('kerot',0.0d0,vtout,kerot,itout) vsrc = vsxynt('pegrav',0.0d0,vtout,pegrav,itout) end if write(*,*) vtout(itout), ( vyout(itout,ieq) , ieq = 1 , neq ) stop 900 continue write(0,*) 'usage: dumb '// & ' [ tol]' stop end