PROGRAM shot
      REAL r(2),v(2),a(2)
      CALL init
      DO WHILE(loopin(t,r,v,a) .EQ. 1)
        DO i = 1, 100
          CALL euler(t,r,v,a)
          CALL prtout(t,r,v,a)
          IF(r(2) .LE. 0.0) GOTO 100
        ENDDO
  100   CONTINUE
      ENDDO
      CALL GWQUIT(IR)
      END
C
      SUBROUTINE init
      PARAMETER(WX = 50.0, R2 = 1.414213, PI=3.141592)
      CHARACTER ttl*80, str*80
      COMMON /gparam/XPOS1,XPOS2,YPOS,FH,SMRK
      COMMON /param/dt,g,c,ncalc,height,v0,th
*
      CALL GWOPEN(ID,0)
      CALL GWVPORT(IR,0.0,0.0,R2,1.0)
      xmin = -WX*0.1
      xmax =  WX*1.1
      ymin = -WX/R2*0.1
      ymax =  WX/R2*1.1
      CALL GWINDOW(IR,xmin,ymin,xmax,ymax)
      CALL GWLINE(IR,xmin,0.0,xmax,0.0)
      CALL GWLINE(IR,0.0,ymin,0.0,ymax)
*
      xtw = (xmax-xmin)/100
      ytw = (ymax-ymin)/100
*
      CALL GWSETTXT(IR,3*xtw, 0.0, 8, -1, -1, ' ')
      x = 0
      DO WHILE(x .LT. xmax-11)
        x = x + 10
        CALL GWLINE(IR,x,-xtw,x,xtw)
        WRITE(str, '(I2)') nint(x)
        CALL GWPUTTXT(IR,x, -xtw, str)
      ENDDO
*
      CALL GWSETTXT(IR,-1.0, -1.0, 7, -1, -1, ' ')
      y = 0
      DO WHILE(y .LT. ymax-11)
        y = y + 10
        CALL GWLINE(IR,-ytw,y,ytw,y)
        WRITE(str, '(I2)') nint(y)
        CALL GWPUTTXT(IR,-2*ytw, y, str)
      ENDDO
      CALL GWSETTXT(IR,-1.0, -1.0, 3, -1, -1, ' ')
      CALL GWPUTTXT(IR,xmax-ytw, -xtw, '[m]')
      CALL GWPUTTXT(IR,-ytw, ymax-xtw, '[m]')
      CALL GWPUTTXT(IR,-2*ytw, -xtw, '0')
*
      g = 9.80665
      prper = 0.2
      dt = 0.001
      ncalc = prper/dt + 0.5
      WRITE(ttl,'(A,$)') 
     +   '初期条件： 高さ[m], 速さ[m/s], 角[deg] = '
      CALL GWINPUT(l,ttl, str)
      READ(str,*) height,v0,th
*
      CALL GWSETTXT(IR,-1.0, -1.0, 1, -1, -1, ' ')
      XPOS = 5
      YPOS = ymax-3
      SMRK = 0.7
      WRITE(str, '(A,F5.1,A)') '　高さ =',height,' [m]'
      CALL GWGETTXT(IR,DMY, FH, DMY, DMY, str)
      CALL GWPUTTXT(IR,XPOS, YPOS, str)
      YPOS = YPOS - FH
      WRITE(str, '(A,F5.1,A)') '初速度 =',v0,' [m/s]'
      CALL GWPUTTXT(IR,XPOS, YPOS, str)
      YPOS = YPOS - FH
      WRITE(str, '(A,F5.1,A)') '発射角 =',th,' [deg]'
      CALL GWPUTTXT(IR,XPOS, YPOS, str)
*
      CALL GWSETTXT(IR,-1.0, -1.0, 5, -1, -1, ' ')
      CALL GWSETSYM(IR,4*xtw, 0.0, 0, -1, -1, 'Symbol')
      CALL GWGETSYM(IR,FW, DMY, 103)               ! 103 = 'γ'
      XPOS = xmax-10
      YPOS = ymax-3
      XPOS1 = XPOS + FW/2
      XPOS2 = XPOS+2*FW
      CALL GWPUTSYM(IR,XPOS1, YPOS, 103)
      CALL GWPUTTXT(IR,XPOS2, YPOS, '[1/m]')
      YPOS = YPOS - FH/3
*
      th = th/180*PI
      END
C
      FUNCTION loopin(t,r,v,a)
      CHARACTER ttl*80, str*80
      REAL r(2),v(2),a(2)
      COMMON /gparam/XPOS1,XPOS2,YPOS,FH,S
      COMMON /param/dt,g,c,ncalc,height,v0,th
      DATA M,K/0,0/
      WRITE(ttl,'(A,$)') '抵抗係数 γ[1/m] (STOP if < 0) = '
      CALL GWINPUT(l,ttl, str)
      READ(str,*) c
      IF(c .LT. 0.0) THEN
        loopin = 0
        RETURN
      ENDIF
      t = 0
      v(1) = v0*cos(th)
      v(2) = v0*sin(th)
      r(1) = 0
      r(2) = height
      vv = sqrt(v(1)**2 + v(2)**2)
      a(1) = - c*vv*v(1)
      a(2) = - g -c*vv*v(2)
      M=MOD(M,6)+1
      K=MOD(K,18)+1
      CALL GWSETMRK(IR,M,S,K,-1,-1)
      WRITE(str,'(F5.3)') c
      YPOS = YPOS - FH
      CALL GWPUTMRK(IR,XPOS1,YPOS)
      CALL GWPUTTXT(IR,XPOS2, YPOS, str)
*
      WRITE(*,'(/2(1X,2A/))') 
     + '   t       x       y      v       a   ',
     + '    vx      vy      ax      ay  ',
     + '  (s)     (m)     (m)   (m/s)  (m/s^2)',
     + '  (m/s)   (m/s)  (m/s^2) (m/s^2)'
      CALL prtout(t,r,v,a)
      loopin = 1
      END
C
      SUBROUTINE prtout(t,r,v,a)
      REAL r(2),v(2),a(2)
      va = sqrt(v(1)**2+v(2)**2)
      aa = sqrt(a(1)**2+a(2)**2)
      WRITE(*,'(1X,F6.2,8F8.3)') t,r,va,aa,v,a
      CALL GWPUTMRK(IR,r(1),r(2))
      END
C
      SUBROUTINE euler(t,r,v,a)
      COMMON /param/dt,g,c,ncalc,height,v0,th
      REAL r(2),v(2),a(2)
      DO i = 1, ncalc
        r(1) = r(1) + v(1)*dt
        r(2) = r(2) + v(2)*dt
        vv = sqrt(v(1)**2 + v(2)**2)
        a(1) = - c*vv*v(1)
        a(2) = - g -c*vv*v(2)
        v(1) = v(1) + a(1)*dt
        v(2) = v(2) + a(2)*dt
        t = t + dt
        IF(r(2) .LE. 0.0) RETURN
      ENDDO
      END