program sp211 ! 空気抵抗あり(速度の2乗に比例) implicit none ! x方向 y方向を独立に計算 integer i, imax real*8 pi, g, m, k2, tmax, dt parameter( pi=3.141592653 ) parameter( g=9.80217, m=7, k2=0.01, tmax=6, dt=0.01 ) real*8 t0, X0, Y0, V0, theta0 parameter( t0=0, X0=0, Y0=40, V0=30, theta0=45 ) real*8 t, X,Xmid, Y,Ymid, V,Vmid, Vx,Vxmid, Vy,Vymid, & Ax,Axmid, Ay,Aymid imax = tmax / dt t = t0 ! sec X = X0 ! meter Y = Y0 ! meter V = V0 ! m/s Vx = V0*cos(theta0*pi/180) ! radian Vy = V0*sin(theta0*pi/180) ! radian open(1,file='sp211.dat') write(1,'(a1,3x,a6,2(8x,a4),2(9x,a7),5x,a6)') & '#','T(sec)','X(m)','Y(m)','Vx(m/s)','Vy(m/s)','V(m/s)' write(1,'(1p6e11.3)') t, X, Y, Vx, Vy, V do i=1,imax Ax = - k2/m * Vx * Vx ! <-- k2/m*Vx*Vx Ay = -g - k2/m * abs(Vy) * Vy ! <-- k2/m*Vx*Vy Vxmid = Vx + Ax * dt/2 Vymid = Vy + Ay * dt/2 Vmid = sqrt( Vxmid**2 + Vymid**2 ) Xmid = X + Vx * dt/2 ymid = Y + Vy * dt/2 Axmid = - k2/m * Vxmid * Vxmid ! <-- k2/m*Vxmid*Vxmid Aymid = -g - k2/m * abs(Vymid) * Vymid ! <-- k2/m*Vymid*Vymid Vx = Vx + Axmid*dt Vy = Vy + Aymid*dt V = sqrt( Vx**2 + Vy**2 ) X = X + Vxmid*dt Y = Y + Vymid*dt t = t + dt write(1,'(1p6e11.3)') t, X, Y, Vx, Vy, V enddo close(1) stop end