program pinball implicit none integer pgopen,i_again,n_nail,ii,num_of_nails_left parameter (n_nail = 150) real x,vx,ax, y,vy,ay, vx_ini,vy_ini, x_range,y_range,delta_t real mass,k,radius_2 real nail_x(n_nail),nail_y(n_nail),old_pos_x,old_pos_y,radius parameter (radius = 0.5) radius_2 = radius - 0.02 call input_n_init(vx_ini,vy_ini,delta_t,mass,k) if (pgopen('/xwin').le.0) stop c call pgpap(4.0, 2.0) call pgpap(2.0, 2.0) call pgsfs(2) 100 write(*,*) 'How far (x-dir) and how high (y-dir) is the plot ?' read (*,*) x_range, y_range call pgeras call pgenv(0.0,x_range,0.0,y_range,1,0) call nail_line_up (x_range,y_range,n_nail,nail_x,nail_y, & num_of_nails_left,radius) x = 0.0+radius y = 0.9*y_range vx = vx_ini vy = vy_ini old_pos_x = x old_pos_y = y do while (y.gt.0.0) call fly_drift (x,vx,ax,y,vy,ay,delta_t,mass,k) call pgbbuf call pgsci (0) call pgcirc (old_pos_x,old_pos_y,radius_2) call pgsci (1) call pgcirc (x,y,radius_2) call pgpt (n_nail-num_of_nails_left,nail_x,nail_y,-1) call pgebuf old_pos_x = x old_pos_y = y call test_collide(x,vx,y,vy,nail_x,nail_y,n_nail,radius,ii, & x_range,y_range) if (ii.ne.0) call bounch (x,vx,y,vy,nail_x(ii),nail_y(ii)) end do 102 write(*,*) 'Again ? (1=Yes;0=No)' read(*,*) i_again if (i_again.eq.1) goto 100 if (i_again.eq.0) stop goto 102 call pgclos end subroutine input_n_init(vx,vy,delta_t,mass,k) implicit none real vx,vy,v,angle,ke,delta_t,mass,k integer i_gun_powder write(*,*) 'How maqny packs of gun powder (integer) ?' read (*,*) i_gun_powder write (*,*) 'What is the angle of shooting ?' read (*,*) angle write(*,*) 'What is the small time-interval delta_t (< 0.001) ?' read (*,*) delta_t write(*,*) 'What is the mass of the particle ?' read (*,*) mass write(*,*) 'What is the coefficient of air resistence ?' read (*,*) k angle = angle * 3.14159268/180 ke = 1.0*i_gun_powder v = sqrt(2.0*ke/mass) vx = v*cos(angle) vy = v*sin(angle) end subroutine fly_drift(x,vx,ax,y,vy,ay,delta_t,mass,k) implicit none real x,vx,ax,y,vy,ay,delta_t,theda,g real v_old,v_old_x,v_new_x,v_old_y,v_new_y,mass,v_mid_x real a_mid_x,a_mid_y,a_drift_x,a_drift_y,a_drift,k,v_mid_y g = -9.8 v_old_x = vx v_old_y = vy v_old = sqrt(v_old_x**2 + v_old_y**2) a_drift = -1*k*v_old**2/mass c theda = atan(v_old_y/v_old_x) c theda = asin(v_old_y/v_old) c a_drift_x = a_drift*cos(theda) c a_drift_y = a_drift*sin(theda) a_drift_x = a_drift*(v_old_x/v_old) a_drift_y = a_drift*(v_old_y/v_old) a_mid_x = a_drift_x a_mid_y = a_drift_y + g v_mid_x = v_old_x + a_mid_x*0.5*delta_t v_mid_y = v_old_y + a_mid_y*0.5*delta_t v_new_x = v_old_x + a_mid_x*delta_t x = x + v_mid_x*delta_t v_new_y = v_old_y + a_mid_y*delta_t y = y + v_mid_y*delta_t vx = v_new_x vy = v_new_y end subroutine test_collide(x,vx,y,vy,nail_x,nail_y,n_nail,radius,ii &,x_range,y_range) implicit none integer n_nail,i,ii real nail_x(n_nail),nail_y(n_nail),x,vx,y,vy,radius real d(n_nail),x_range,y_range c assuming hitting to one of the nails only ii = 0 do i=1,n_nail d(i) = sqrt( (x-nail_x(i))**2 + (y-nail_y(i))**2 ) if (d(i).le.radius) ii = i end do if ((x.le.radius).or.(x.ge.x_range-radius)) vx = -1*vx if ((y.le.radius).or.(y.ge.y_range-radius)) vy = -1*vy end subroutine bounch (x,vx,y,vy,nail_x,nail_y) real x,vx,y,vy,nail_x,nail_y real normal_x,normal_y,v_para_x,v_para_y,v_perp_x,v_perp_y,temp real v_dot_n c N_vec = nail_vec - ball_vec, then normalise N_vec to become n_vec normal_x = nail_x - x normal_y = nail_y - y temp = sqrt(normal_x**2 + normal_y**2) normal_x = normal_x / temp normal_y = normal_y / temp c inner product of V_ball and n_vec v_dot_n = vx*normal_x + vy*normal_y c V_perpendicular = Vdotn n v_perp_x = v_dot_n*normal_x v_perp_y = v_dot_n*normal_y c V_parallele = V_ball - V_peprp v_para_x = vx - v_perp_x v_para_y = vy - v_perp_y c bouch back on n diection v_perp_x = -1*v_perp_x v_perp_y = -1*v_perp_y c New vx anc vy after bouch back vx = v_para_x + v_perp_x vy = v_para_y + v_perp_y end subroutine nail_line_up (x_range,y_range,n_nail,nail_x,nail_y, & num_of_nails_left,radius) implicit none integer n_nail real x_range, y_range, nail_x(n_nail), nail_y(n_nail),shift integer n_col, n_row, row_size,col_size,num_of_nails_left,i integer ix,iy real dx,dy,h_space,v_space,radius 200 write(*,*) "What is the horizontal and vertical nail spacing with &respect to ball diameter ?" read (*,*) h_space, v_space c if (h_space.le.1.0)) stop ('Nails too close !') c if (h_space.le.1.0) then c write(*,*) 'Nails too close !' c goto 200 c endif h_space = h_space * 2 * radius v_space = v_space * 2 * radius row_size = int(x_range/h_space) - 1 col_size = int(0.75*y_range/v_space) - 1 c col_size = int(y_range/v_space) - 1 c if ((row_size.eq.0).or.(col_size.eq.0)) stop ('table too small') if ((row_size.eq.0).or.(col_size.eq.0)) then write(*,*) 'table too small' stop endif dx = x_range / row_size dy = 0.75*y_range / col_size c dy = y_range / col_size num_of_nails_left=n_nail c shift = 1.0 shift = 2.0*radius do iy = 1, col_size do ix = 1, row_size i = (iy-1)*row_size + ix nail_x(i) = shift + (ix-1)*dx nail_y(i) = iy* dy c call pgpt (1,nail_x,nail_y,-1) num_of_nails_left = num_of_nails_left - 1 if (num_of_nails_left.eq.0) exit end do if (num_of_nails_left.eq.0) exit c shift = 1.0 - shift shift = 2.0*radius - shift end do call pgpt (n_nail-num_of_nails_left,nail_x,nail_y,-1) end