program HW02_03 implicit none * Program to solve a double pendulum problem from initial conditions * Arguments for the program are passed through runsting: * % HW02_03 [OutOPt] [error tolerance] * where and are lengths (m) * and are masses (kg) * are initial angles (deg) (Program assumes * a stationary start to pedulums) * and are duration of integation * and the output interval (sec) * [OutOpt] Combination of letters; option are * A - Animation, * E - Energy, * R - Runge-Kutta integration * S - simple integration * T - Slightly better integration.(but not Runge-Kutta) * At least RST must be specified. (RAE gives Runge-Kutta, animation * and energy). * is fractional error tolerance (rads) * include 'DPend.h' * MAIN Program variables real*8 t ! Current time (seconds) ., step ! Current step size (seconds) ., new_step ! Updated step size that will yield desired accuracy real*8 ths(2), oms(2) ! Theta and rate at start of integration step real*8 thc(2), omc(2) ! Theta and rate at end ! of integration step when done as one step real*8 tha(2), oma(2) ! Theta and rate at end of first half/step integration real*8 thb(2), omb(2) ! Theta and rate at end ! of second half/step integration (thb should = thc ! within tolerance. real*8 minstep ! Minimum step size needed. real*8 next_out, prev_step ! Time of next output step and previous step ! before stepping onto output time integer*4 sc ! Counter for step so that we don't get stuck integer*4 i,j ! Loop counters integer*4 rots(2) ! Turns of angle real*8 rtheta(2) ! Residual after removing rotations real*8 KE, PE ! Kinetic and potential energy logical step_not_OK ! Logical set true while step size is still being ! evaluated ., no_step_increase ! Set true to stop step from inceasing. ., save_step ! Save the values from this step when time falls ! on value ., energy_out ! Set true to output energy **** The command line arguments passed by user call read_commline **** OK Start the integration energy_out = .false. if( index(out_options,'E').gt.0 ) energy_out = .true. step = out_dt/2 ! Start at 1/2 output interval minstep = step * Initializa time t = 0.d0 next_out = t + out_dt * Initialize starting Theta and rate (We are using f90 array * features to assign the values in the arrays ths = thetaI oms = OmegaI num_save = 1 timeS(num_save) = t thetaS(:,num_save) = ths OmegaS(:,num_save) = oms next_out = t + out_dt * Report initial conditions call report_ICS(t, ths, oms, 'Initial Conditions') if( index(out_options,'A').gt.0 ) call animate(t, step, ths) do while ( t.lt. dur ) * Nake sure that the next step will nor overshoot the next * output epoch (this can happen because are are varying the * step size during the run.) no_step_increase = .false. save_step = .false. prev_step = step if ( t + step .gt.dur ) then step = dur - t no_step_increase = .true. save_step = .true. else ! See if we will miss an output time if( t+step.gt.next_out ) then step = next_out - t no_step_increase = .true. save_step = .true. elseif ( t+step.eq.next_out ) then save_step = .true. no_step_increase = .true. endif 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.20 ) sc = sc + 1 * Step with curret step size call int_step( thc, omc, ths, oms, step, t ) * Now take two steps with half size call int_step( tha, oma, ths, oms, step/2, t ) call int_step( thb, omb, tha, oma, 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, thc, thb, step, new_step, sc) * Make sure we have not increased step at wrong time if( .not.save_step ) then 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 else step_not_OK = .false. end if end do * save smallest step size if( step.lt. minstep .and. .not.save_step ) minstep = step **** OK step hase been correctly taken, update state and time t = t + step ths = thb oms = omb if( save_step .or. t.ge. dur ) then num_save = num_save + 1 timeS(num_save) = t thetaS(:,num_save) = ths OmegaS(:,num_save) = oms next_out = num_save*out_dt ! t + out_dt step = prev_step end if ***** Use more man's graphics with vt100 sequences if( index(out_options,'A').gt.0 ) call animate(t, step, ths) end do write(*,210) num_save 210 format('There are ',i4,' saved values ',/, . '. # Time Beam 1 Angle Rev Remain ', . 'Rate Beam 2 Angle Rev Remain Rate',/, . '. sec deg deg ', . 'deg/s deg deg deg/s') do i = 1, num_save rots = thetaS(:,i)/(2*pi) rtheta = thetaS(:,i)-rots*2*pi call comp_energy(thetaS(1,i),OmegaS(1,i),KE,PE) write(*, 230) i, timeS(i), (thetaS(j,i)*180/pi, . rots(j), rtheta(j)*180/pi, OmegaS(j,i)*180/pi,j=1,2) 230 format(' ',I4,1x,F10.3,2(F15.4,1x,I4,F11.5,1x,F11.5,2x)) if( energy_out ) write(*,240) i, timeS(i), PE, KE, PE+KE 240 format('ENERGY ',I4,1x,F10.3,' PE, KE, PE+KE ',3(1x,F12.5)) end do **** Output the final values call report_ICS(t, ths, oms,'Final conditions') write(*,320) minstep 320 format('Smallest step size needed ',F9.6,' seconds') end *----------------------------------------------------------------------- subroutine read_commline implicit none include 'DPend.h' * Routine to read the runsting for hw02_03 * % HW02_03 * where and are lengths (m) * and are masses (kg) * are initial angles (deg) (Program assumes * a stationary start to pedulums) * and are duration of integation * and the output interval (sec) * is fractional error tolerance (rads) character*80 argl ! String run for each argument in command line. * LOCAL VARIABLES integer*4 ierr ! IOSTAT errors ., lenline ! Length of string **** Get the input file name (len_trim could be used in f90 instead of lenline) call getarg(1,argl) if( lenline(argl).eq.0 ) then ! No arguments; tell user command line help call proper_commline( 0 ) else * See if we can read length 1 read(argl,*,iostat=ierr) lens(1) if( ierr.ne.0 ) call proper_commline( 1 ) endif **** Get length for beam 2 call getarg(2,argl) read(argl,*,iostat=ierr) lens(2) if( ierr.ne.0 ) call proper_commline( 2 ) **** Get mass for element 1 call getarg(3,argl) read(argl,*,iostat=ierr) mass(1) if( ierr.ne.0 ) call proper_commline( 3 ) **** Get mass for element 2 call getarg(4,argl) read(argl,*,iostat=ierr) mass(2) if( ierr.ne.0 ) call proper_commline( 4 ) **** Get Angle for beam 1. Input in deg, convert to rad call getarg(5,argl) read(argl,*,iostat=ierr) thetaI(1) if( ierr.ne.0 ) call proper_commline( 5 ) thetaI(1) = thetaI(1)*pi/180.d0 OmegaI(1) = 0.0d0 ! Initialize now (later versions may allow input) **** Get Angle for beam 1. Input in deg, convert to rad call getarg(6,argl) read(argl,*,iostat=ierr) thetaI(2) if( ierr.ne.0 ) call proper_commline( 6 ) thetaI(2) = thetaI(2)*pi/180.d0 OmegaI(2) = 0.0d0 ! Initialize now (later versions may allow input) **** Get duration of integration call getarg(7,argl) read(argl,*,iostat=ierr) dur if( ierr.ne.0 ) call proper_commline( 7 ) **** Get output interval call getarg(8,argl) read(argl,*,iostat=ierr) out_dt if( ierr.ne.0 ) call proper_commline( 8 ) * See if we can save all the output if( int(dur/out_dt).le.0 .or. int(dur/out_dt)+1.gt.max_save ) then write(*,810) dur, out_dt, int(dur/out_dt), max_save 810 format('HW02_03: FATAL Error. For duration ',F10.2, . ' sec at output interval ',f8.3,' sec, ',i8, . ' values need to be saved which is greater than ',i4, . ' allowed.',/, . 'Modify max_save in DPend.h to allow more saved values') stop 'HW02_03: To many saved values' end if **** See if output options passed call getarg(9,out_options) **** Get the error tolerance. Default to 1.d-4 call getarg(10,argl) if( lenline(argl).eq.0 ) then ! No arguments; set default tol = 1.d-4 else * See if we can read length 1 read(argl,*,iostat=ierr) tol if( ierr.ne.0 ) call proper_commline( 10 ) endif ***** Thats all return end *----------------------------------------------------------------------- subroutine proper_commline( missed_arg ) implicit none * Routine to output the correct command line for the double * pendulum problem. narg is the argument number that is * missing or incorrect integer*4 missed_arg ! Argument number missing write(*,120) 120 format('HW02_03: Double pendulum. Correct command line is:',/, . '% HW02_03 ', . ' [OutOpt] [error tolerance]',/, . 'where and are lengths (m)',/, . ' and are masses (kg)',/, . ' are initial angles (deg) ', . '(Program assumes',/, . ' a stationary start to pedulums)',/, . ' and are duration', . ' of integation',/, . ' and the output interval (sec)',/, . ' [OutOpt] Combination of letters; option are',/, . ' A - Animation,',/, . ' E - Energy,',/, . ' R - Runge-Kutta integration',/, . ' S - simple integration',/, . ' T - Slightly better integration.(but not' . ' Runge-Kutta)',/, . ' At least RST must be specified. (RAE gives', . ' Runge-Kutta, animation',/, . ' and energy).',/, . ' [Error tolerance] is fractional error', . ' tolerance (rads, Default 1e-4)',/) if( missed_arg.gt. 0 ) then write(*,140) missed_arg 140 format(/,'Argument number ',i2,' not correct or missing') else write(*,160) 160 format(/'No arguments given') end if stop 'HW02_03: Program stopped. Error in command line' * Thats all return end *----------------------------------------------------------------------- subroutine report_ICS(t, thr, omr, type) implicit none * Routine to report status of bodies at time t ond of type type include 'DPend.h' * PASSED VARIABLES real*8 t ! Current time days real*8 thr(2) ! Positions of bodies (km) ., omr(2) ! Velocities (km/s) character*(*) type ! Type of output.: String passed by user * LOCAL VARIABLES integer*4 i integer*4 rots(2) ! Number of rotations real*8 dthr(2) ! Angle after removing rotation write(*,120) type, t 120 format('+ 12.010 HW 02 Q 03: ',a,' at time ',F14.5,' seconds') rots = thr/(2*pi) dthr = thr - rots*2*pi write(*,140) (i, rots(i), dthr(i)*180/pi, omr(i)*180/pi,i = 1,2) 140 format('Beam ',i1,' Turns ',i4,' Angle ',F8.3, . ' deg, Rate ',F11.5,' deg/sec') 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, th) implicit none * Routine to animate the changes in position include 'DPend.h' * PASSED Variables real*8 t ! Current time ., step ! Current step size ., th(2) ! Positions of bodies (also use mass to ! get size real*8 maxx, maxy, minx,miny ! Min and max coordintes in initial state real*8 scalex, scaley ! Scales to convert x and y to screen coordinates ., offx, offy ! Offsets for origin in X and Y ., x(2), y(2) ! Coordinates of masses integer*4 SWX, SHY ! Screen width and height ., IX, IY ! Coorinates in screen coordinatea ., i ! Loop counter character*1 sym(2) ! Symbols for each element 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, 2 sym(i) = '*' end do sym(2) = '+' * Now get the scale SWX = 85 ! 85 characeter accross (5/3 ratio) SHY = 45 ! 45 lines down maxx = lens(1)+lens(2) minx = -maxx maxy = maxx miny = minx * Set Scale in X and Y Scalex = SWX/(maxx-minx) * 0.95d0 Scaley = SHY/(maxy-miny) * 0.95d0 offx = minx * 1.05d0 offy = miny * 1.05d0 write(mv45,100) char(27) 100 format(a1,'[47A') else * Move the cursor back 52-lines write(*,'(a)') mv45 endif ! print *,'Scale ',scalex, scaley, offx, offy write(*,120) t, step 120 format('TIME ',F10.5,' secs, Step ',F10.5, ' secs') * OK form up the output lines do i = 1,45 screen(i) = ' ' enddo x(1) = lens(1)*sin(th(1)) y(1) = -lens(1)*cos(th(1)) x(2) = x(1) + lens(2)*sin(th(2)) y(2) = y(1) - lens(2)*cos(th(2)) do i = 1,2 IX = 42 + nint(x(i)*scalex) IY = 22 - nint(y(i)*scaley) ! print *,'LOC',i,x(i),y(i), 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 screen(22)(42:42) = 'V' * write lines do i = 1,45 write(*,'(a)') screen(i) enddo **** Thats all return end *----------------------------------------------------------------------- subroutine int_step( the, ome, ths, oms, ss, time ) implicit none * Routine to integrate 1 time step of duration step in days include 'DPend.h' * Passed Variables real*8 the(2) ! Thetas at end of step (rad) ., ome(2) ! d(Theta)/dt at end of step (rad/s) ., ths(2) ! Thetas at start of step (rad) ., oms(2) ! d(Theta)/dt at emd of step (rad/s) ., ss ! Step size in seconds ., time ! Current time (seconds) **** Local variables real*8 T(2), W(2) ! Thetas and first derivatives for the two ! pendulums (rads rad/s) ., A(2) ! Second derivative of T (thetas) (rad/s^2) ., k1(2), k2(2), k3(2), k4(2) ! Intermediate integration ! points in Runge-Kutta integration. **** Get the accelerations of all bodies at start. Note: Here * we are using the f90 feature for array multiplication. * Use Runga-Kutta if S to T not specified for simplier integration. if( index(out_options,'S').eq.0 .and. . index(out_options,'T').eq.0) then call get_accel(time, ths, oms, A) k1 = SS*A * Step system to mid-point T = ths + SS*oms/2.d0 + SS*k1/8.d0 W = oms + k1/2.d0 * Get acceleration at this new point call get_accel(time+ss/2, T, W, A) k2 = SS*A T = ths + SS*oms/2.d0 + SS*k1/8.d0 W = oms + k2/2.d0 * Nexy value at mid point call get_accel(time+ss/2, T, W, A) k3 = SS*A T = ths + SS*oms + SS*k3/2.d0 W = oms + k3 * Now do final end point call get_accel(time+ss, T, W, A) k4 = SS*A the = ths + SS*(oms + (k1+k2+k3)/6.d0) ome = oms + (k1+2*k2+2*k3+k4)/6.d0 elseif ( index(out_options,'S').ne.0 ) then * Use a simple integration call get_accel(time, ths, oms, A) the = ths + oms*ss + A*ss**2/2 ome = oms + A*ss * Try a slightly better method than simple Euler else * Use the Euler eastimates to refine the acceleration * and use average value call get_accel(time, ths, oms, A) k1 = SS*A T = ths + oms*ss + k1*SS/2 W = oms + k1 call get_accel(time+ss, T, W, A) k2 = SS*A * This last step simply uses the average of the two acceleration * values. the = ths + SS*(oms + (k1+k2)/4.d0) ome = oms + (k1+k2)/2.d0 endif **** That all return end *----------------------------------------------------------------------- subroutine get_accel(time, T, W, A) implicit none * Routine to the second derivatives of angles Theta (T) given * Theta (T), d(Theta)/dt (W). The second derivatives are * returned in variable A include 'DPend.h' * PASSED VARIABLES real*8 time ! Time in seconds (not used) ., T(2), W(2) ! Thetas and first derivatives for the two ! pendulums (rads rad/s) ., A(2) ! Second derivative of T (rad/s^2) * LOCAL VARIABLES * none * Evalutate with expressions fo myphysicslab.com A(1) = (-g*(2*mass(1)+mass(2))*sin(T(1)) - . g*mass(2)*sin(T(1)-2*T(2)) - . 2*sin(T(1)-T(2))*mass(2)* . (W(2)**2*lens(2)+W(1)**2*lens(1)*cos(T(1)-T(2)))) / . (lens(1)*(2*mass(1)+mass(2)-mass(2)*cos(2*T(1)-2*T(2)))) A(2) = (2*sin(T(1)-T(2))* . (W(1)**2*lens(1)*(mass(1)+mass(2)) + . g*(mass(1)+mass(2))*cos(T(1)) + . W(2)**2*lens(2)*mass(2)*cos(T(1)-T(2)))) / . (lens(2)*(2*mass(1)+mass(2)-mass(2)*cos(2*T(1)-2*T(2)))) **** Thats all return end *----------------------------------------------------------------------- subroutine eval_step( t, thc, thb, step, new_step, sc) implicit none * Routine to see if step size meetds tolerances include 'DPend.h' * PASSED VARIABLES real*8 t ! current time in days ., thc(2) ! theta with 1-step ., thb(2) ! theta 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 dTh(2) ! Difference in angle estimates ., err ! Error in angles ., tolsc ! Tolerance scaled to keep error bounded integer*4 i, j ! Loop counters **** Get the difference and see how low dTh = thc - thb **** Project error to end of integration assuming that the error * level remains the same and will grow as sqrt(2). Take the * abs((dur-t)/step) because step can be negative if we overshoot * an output time. if( step.le.0 ) print *,'Time ',t,' Negative step ',step, 's' err = sqrt((dTh(1)**2+dTh(2)**2)*abs((dur-t)/step)) tolsc = tol/100 **** See if we have meet tolerance if( sc.gt.15 ) then write(*,120) t, sc, step, err, tol 120 format('UpateStep Time ',F12.6,1x, i3,' Size ',F9.6,' Err ', . E13.3,' Tol ',E13.3) endif new_step = step if( err.lt.tolsc/100 .and. step.lt. 1.d0 .and. sc.lt.10 ) then new_step = step*2 endif if( err.gt.tolsc .and. step.gt. 1.d0/2.d0**20 ) then new_step = step/2 endif **** Thats all return end *----------------------------------------------------------------------- subroutine comp_energy(TS,OS,KE,PE) implicit none * Routine to compute Kinetic and Potential energy include 'DPend.h' * PASSED VARIABLES real*8 TS(2), OS(2) ! Theta and rates of change real*8 KE, PE ! Kinetic and Potential Energy PE = -(mass(1)+mass(2))*g*lens(1)*cos(TS(1)) - . mass(2)*g*lens(2)*cos(TS(2)) KE = (mass(1)* lens(1)**2*OS(1)**2 + . mass(2)*(lens(1)**2*OS(1)**2 + lens(2)**2*OS(2)**2 + . 2*lens(1)*lens(2)*OS(1)*OS(2)*cos(TS(1)-TS(2))))/2 return end