program HW02_03 implicit none * Program to solve the N-body gravitational problem using Runge-Kutta * numerical integration. * Arguments for the program are passed through runsting: * % HW02_03 * include 'NBody.h' * MAIN Program variables real*8 t ! Current time (days) ., step ! Current step size (days) ., new_step ! Updated step size that will yield desired accuracy real*8 poss(2,max_body), vels(2,max_body) ! Position and velocity at start ! of integration step real*8 posc(2,max_body), velc(2,max_body) ! Position and velocity at end ! of integration step when done as one step real*8 posa(2,max_body), vela(2,max_body) ! Position and velocity at end ! of first half/step integration real*8 posb(2,max_body), velb(2,max_body) ! Position and velocity at end ! of second half/step integration (posb should = posc ! within tolerance. real*8 minstep ! Minum step size needed. integer*4 sc ! Counter for step so that we don't get stuck logical step_not_OK ! Logical set true while step size is still being ! evaluated ., no_step_increase ! Set true to stop step from inceasing. **** The runstring arguments passed by user call read_runstring( infile, dur, tol) **** Get the information on the initial conditions call read_nbody **** OK Start the integration step = 1.0d0 minstep = step * Initializa time t = 0.d0 closest = 1.d20 * Initialize starting position and velocity poss(:,1:num_body) = posi(:,1:num_body) vels(:,1:num_body) = veli(:,1:num_body) * Report initial conditions call report_ICS(t, poss, vels, 'Initial Conditions') call animate(t, step, poss) do while ( t.lt. dur ) * Nake sure that the next step will nor overshoot the end * of the time (this can happen because are are varying the * step size during the run. no_step_increase = .false. if ( t + step .gt.dur ) then step = dur - t no_step_increase = .true. endif * Iterate on step size to get the value needed to maintain the * requested tolerance step_not_OK = .true. sc = 0 do while ( step_not_OK .and. sc.lt.15 ) sc = sc + 1 * Step with curret step size call int_step( posc, velc, poss, vels, step, t ) * Now take two steps with half size call int_step( posa, vela, poss, vels, step/2, t ) call int_step( posb, velb, posa, vela, step/2, t+step/2 ) * Compare the two results and see if the meet tolerances or if * results too good can we increase step call eval_step( t, posc, posb, step, new_step, sc) * Make sure we have not increased step at wrong time if( no_step_increase .and. new_step.gt.step ) then new_step = step endif * If the new_step is the same as step then we are done. if( new_step.eq.step ) then step_not_OK = .false. else step = new_step endif end do * save smallest step size if( step.lt. minstep ) minstep = step **** OK step hase been correctly taken, update state and time t = t + step poss(:,1:num_body) = posb(:,1:num_body) vels(:,1:num_body) = velb(:,1:num_body) ***** Use more man's graphics with vt100 sequences call animate(t, step, poss) end do **** Output the final values call report_ICS(t, poss, vels,'Final conditions') write(*,320) minstep, closest, close_body 320 format('Smallest step size needed ',F9.6,' days',/, . 'Closest approach distance was ',ES14.5, . ' km Bodies ',I2,' and ',i2) end subroutine read_runstring( infile, dur, tol ) implicit none * Routine to read the runsting for hw02_03 * % hw02_03 real*8 dur ! Duration of integration days ., tol ! Relative tolerence of integration character*(*) infile ! Name of input file * LOCAL VARIABLES integer*4 ierr ! IOSTAT errors ., lenline ! Length of string character*64 runstring **** Get the input file name call getarg(1,infile) if( lenline(infile).eq.0 ) then write(*,120) 120 format('HW02_03:N-body problem. Correct runstring:',/, . '% HW02_03 ',/, . 'Default duration is 550 days, ', . 'default tolerance 1.e-6') stop 'HW02_03: No file name passed' endif **** See if other arguments are passed call getarg(2,runstring) if( lenline(runstring).gt.0 ) then read(runstring,*,iostat=ierr) dur if( ierr.ne.0 ) then write(*,220) ierr, runstring(1:lenline(runstring)) 220 format('IOSTAT error ',i5,' reading duration from ',a) stop 'HW02_3: Error reading duraton' endif else dur = 515.d0 endif * Check tolerance call getarg(3,runstring) if( lenline(runstring).gt.0 ) then read(runstring,*,iostat=ierr) tol if( ierr.ne.0 ) then write(*,240) ierr, runstring(1:lenline(runstring)) 240 format('IOSTAT error ',i5,' reading tolerance from ',a) stop 'HW02_3: Error reading tolerance' endif else tol = 1.d-6 endif ***** Thats all return end subroutine read_nbody implicit none * Routine to read the input data file. * File format lines are * Mass (kg) Xpos (km) Ypos (km) Xvel (km/s) Yvel (km/s) * include 'NBody.h' integer*4 n ! Counter for body number ., ierr ! IOSTAT error ., jerr ! IOSTAT error reading data line ., lenline ! Fuction to return length of line character*128 line ***** Open the file and readlines open(20,file=infile,iostat=ierr,status='old') if( ierr.ne. 0 ) then call report_error(ierr, . 'Opening N-body file ' // infile(1:lenline(infile))) stop 'HW02_03: Error opening N-body file' endif ***** Read file n = 0 do while ( ierr.eq.0 .and. n.lt.max_body ) read(20,'(a)',iostat=ierr) line if( ierr.eq.0 .and. line(1:1).eq.' ' ) then * read the line n = n + 1 read(line,*,iostat=jerr) mass(n), posi(:,n), veli(:,n) if( jerr.ne.0 ) then call report_error(jerr, . 'Reading N-body line ' // line(1:lenline(line))) stop 'HW02_03: Error reading N-body line' endif endif enddo if( n.eq. max_body ) then write(*,220) max_body 220 format('HW02_03: Maximum number of bodies ',i5,' reached') endif num_body = n **** Thats all return end subroutine report_ICS(t, posr, velr, type) implicit none * Routine to report status of bodies at time t ond of type type include 'NBody.h' * PASSED VARIABLES real*8 t ! Current time days real*8 posr(2,num_body) ! Positions of bodies (km) ., velr(2,num_body) ! Velocities (km/s) character*(*) type ! Type of output.: String passed by user * LOCAL VARIABLES integer*4 i write(*,120) type, t 120 format('+ 12.010 HW 02 Q 03: ',a,' At time ',F14.5,' days',/, . 'Body Mass (kg) PosX (km) PosY (km) ', . 'VelX (km/s) VelY (km/s)') do i = 1, num_body write(*,140) i, mass(i), posr(:,i),velr(:,i) 140 format(I4,2x,ES18.6,2x,2(ES18.7,1x),2x,2(ES18.7,1x)) end do return end subroutine report_error(ierr,mess) implicit none * Routine to report IOSTAT errors. Initially developed * for 12.010 HW2 Problem 2. * PASSED VARIABLES integer*4 ierr character*(*) mess * LOCAL VARIABLES * None ***** See if the IOSTAT error was non-zero if( ierr.ne.0 ) then write(*,120) ierr, mess 120 format(/,'IOSTAT Error ',i4,' occurred ',a) end if return end integer*4 function lenline(in) implicit none * Returns length of used portion of string. Originally * developed for HW2 Problem 2. f90 has an intrisic called * len_trim() that performs the same function. * PASSED VARIABLES character*(*) in ! The string passed that we need to find ! the last character of * LOCAL VARIABLES integer*4 len_in ! Declared length of in string integer*4 i ! Counter used to work backwards through string ***** Get the declared length of string len_in = LEN(in) i = len_in do while ( i.gt.0 .and. in(i:i).eq.' ') i = i - 1 end do lenline = i return end subroutine animate(t, step, pos) implicit none * Routine to animate the changes in position include 'NBody.h' * PASSED Variables real*8 t ! Current time ., step ! Current step size ., pos(2,num_body) ! Positions of bodies (also use mass to ! get size real*8 maxx, maxy, minx,miny ! Min and max coordintes in initial state ., maxa, mina ! Max and min of both x and y real*8 scalex, scaley ! Scales to convert x and y to screen coordinates ., offx, offy ! Offsets for origin in X and Y integer*4 SWX, SHY ! Screen width and height ., IX, IY ! Coorinates in screen coordinatea ., i ! Loop counter character*1 sym(max_body) ! Symbols for each body character*85 screen(45) ! 45-lines of 85-characters to output character*6 mv45 ! Escape sequence to move up 51 lines save sym, scalex, scaley, offx, offy, mv45 **** Initialize when t = 0; if( t .eq. 0 ) then * First call; so don't move cursor back up and compute size we need * Get o or . based on mass do i = 1, num_body if( mass(i).gt. 1.e30 ) then sym(i) = 'O' elseif( mass(i).gt. 1.e28 ) then sym(i) = 'o' else sym(i) = '.' endif end do * Now get the scale SWX = 85 ! 85 characeter accross (5/3 ratio) SHY = 45 ! 45 lines down maxx = pos(1,1) minx = pos(1,1) maxy = pos(2,1) miny = pos(2,1) do i = 2,num_body maxx = max(maxx, pos(1,i)) minx = min(minx, pos(1,i)) maxy = max(maxy, pos(2,i)) miny = min(miny, pos(2,i)) enddo maxa = max(maxx, maxy) mina = min(minx, miny) * Set Scale in X and Y Scalex = SWX/(maxa-mina) * 0.95d0 Scaley = SHY/(maxa-mina) * 0.95d0 offx = minx* 1.05d0 offy = (miny - 1.5e8 ) * 1.05d0 write(mv45,100) char(27) 100 format(a1,'[47A') else * Move the cursor back 52-lines write(*,'(a)') mv45 endif c print *,'Scale ',scalex, scaley, offx, offy write(*,120) t, step 120 format('TIME ',F13.6,' days, Step ',F10.6, ' days') * OK form up the output lines do i = 1,45 screen(i) = ' ' enddo do i = 1,num_body IX = nint((pos(1,i)-offx)*scalex) IY = 45-nint((pos(2,i)-offy)*scaley) c print *,'Pos ', IX,IY, ' ',sym(i) if( IX.ge.1 .and. IX.le.85 .and. . IY.ge.1 .and. IY.le.45 ) then if( screen(IY)(IX:IX).eq.' ' ) then screen(IY)(IX:IX) = sym(i) endif endif enddo * write lines do i = 1,45 write(*,'(a)') screen(i) enddo **** Thats all return end subroutine int_step( pose, vele, poss, vels, step, t ) implicit none * Routine to integrate 1 time step of duration step in days include 'NBody.h' * Passed Variables real*8 pose(2,num_body) ! Position at end of step ., vele(2,num_body) ! velocity at end of step ., poss(2,num_body) ! Position at start of step ., vels(2,num_body) ! Velocity at emd of step ., step ! Step size in days ., t ! Current time **** Local variables real*8 SS ! Step size in seconds ., P(2,num_body), V(2,num_body), A(2,num_body) ! Positon, ! velcity and acceleration at intermediate points ., k1(2,num_body), k2(2,num_body), k3(2, num_body) ., k4(2,num_body) ! Intermediate integration points in ! Runge-Kutta integration. **** Step in seconds SS = step*86400.d0 **** Get the accelerations of all bodies at start call get_accel(t, poss, vels, A) k1 = SS*A * Step system to mid-point P = poss + SS*vels/2.d0 + SS*k1/8.d0 V = vels + k1/2.d0 * Get acceleration at this new point call get_accel(t+step/2, P, V, A) k2 = SS*A P = poss + SS*vels/2.d0 + SS*k1/8.d0 V = vels + k2/2.d0 * Nexy value at mid point call get_accel(t+step/2, P, V, A) k3 = SS*A P = poss + SS*Vels + SS*k3/2.d0 V = vels + k3 * Now do final end point call get_accel(t+step, P, V, A) k4 = SS*A pose = poss + SS*(vels + (k1+k2+k3)/6.d0) vele = vels + (k1+2*k2+2*k3+k4)/6.d0 **** That all return end subroutine get_accel(t, P, V, A) implicit none * Routine to compute accelerations include 'NBody.h' * PASSED VARIABLES real*8 t ! Time in days (not used) ., P(2,num_body), V(2,num_body) ! Position and Velocity (not used) ! to compute accelerations ., A(2,num_body) ! Accerations at each body * LOCAL VARIABLES integer*4 i, j ! Counters over the bodies. real*8 dP(2,num_body) ! Position difference for each bodu ., R(num_body) ! Distances bewteen current body and others **** Loop over bodies do i = 1, num_body dP(1,:) = P(1,:) - P(1,i) dP(2,:) = P(2,:) - P(2,i) R = sqrt(dP(1,:)**2+dP(2,:)**2) *** Now get vector sum of accelerations A(:,i) = 0.d0 do j = 1,num_body if ( i.ne.j ) then A(:,i) = A(:,i) + G_univ*mass(j)*dP(:,j)/R(j)**3 endif enddo enddo **** Thats all return end subroutine eval_step( t, posc, posb, step, new_step, sc) implicit none * Routine to see if step size meetds tolerances include 'NBody.h' * PASSED VARIABLES real*8 t ! current time in days ., posc(2,num_body) ! Position with 1-step ., posb(2,num_body) ! Position with 2-step ., step, new_step ! Current and new estimate of step size ! new is factor of two different or the same integer*4 sc ! Step size evalution counter. If two any iterations then ! we stop (so debug out as well) * LOCAL VARIABLE real*8 dPe(2,num_body) ! Difference in position estimates ., R ! distance betweeb bodies ., Rmin ! Smallest difference ., err ! Position error in km ., minerr, maxerr ! Min and max relative error integer*4 i, j ! Loop counters **** Get the difference and see how low dPe = posc - posb * See if meets tolerance maxerr = 0 minerr = 1.d20 do i = 1,num_body err = sqrt(dPe(1,i)**2+dPe(2,i)**2) Rmin = 1.d20 do j = 1,num_body if( j.ne.i ) then R = sqrt((posb(1,j)-posb(1,i))**2 + . (posb(2,j)-posb(2,i))**2) if( R .lt.Rmin ) Rmin = R if( Rmin.lt.closest ) then closest = Rmin close_body(1) = i close_body(2) = j endif endif enddo * Now see if this smallest or largest error if( err/Rmin.lt.minerr) minerr = err/Rmin if( err/Rmin.gt.maxerr) maxerr = err/Rmin enddo **** See if we have meet tolerance if( sc.gt.11 ) then write(*,120) t, sc, step, maxerr, minerr, tol 120 format('UpateStep Time ',F12.6,1x, i3,' Size ',F9.6,' Errs ', . 2E13.3,' Tol ',E13.2) endif new_step = step if( maxerr.lt.tol/40 .and. step.lt. 1.d0 ) then new_step = step*2 endif if( maxerr.gt.tol .and. step.gt. 1.d0/2.d0**16 ) then new_step = step/2 endif **** Thats all return end