program anyshape_falling
      real x(100),y(100),x_rot(100),y_rot(100),xold(100),yold(100)
      real k,m
      character*30 fname

      write(*,*)'File name of your rigid-body data? (30 char max)'
      read(*,*) fname

      open(unit=10,file=fname)
      read(10,*) n
      do i=1,n
        read(10,*) x(i), y(i)
      end do

      goto 300
      x(1)= 1.0
      y(1)= 0.0
      x(2)= 1.0
      y(2)= 0.3
      x(3)= 2.0
      y(3)= 0.0
      x(4)=-2.0
      y(4)= 0.0
      x(5)=-2.5
      y(5)= 0.5
      x(6)=-1.5
      y(6)= 0.5
      x(7)=-1.0
      y(7)= 0.0
 300  continue

      write(*,*)'Vx=?'
      read(*,*)vx
      write(*,*)'Vy=?'
      read(*,*)vy
      write(*,*)'k=?'
      read(*,*)k
      write(*,*)'m=?'
      read(*,*)m
      xcm=0.0
      ycm=1.0
      dt=0.00001
      i_dummy=pgopen('/xwin')
      call pgpap(3.0,0.75)
      call pgenv(0.0,10.0,0.0,10.0,1,0)
      do while (ycm>=0)
        v=sqrt(vx**2+vy**2)
        fd=(-1)*k*v
        ax=fd*(vx/v)/m
        ay=-9.8+fd*(vy/v)/m
        vx=vx+ax*dt
        vy=vy+ay*dt

        theta=asin(vy/v)

      do i=1,n
        x_rot(i) = x(i)*cos(theta) + (-1)*y(i)*sin(theta)
        y_rot(i) = x(i)*sin(theta) +      y(i)*cos(theta)
      end do

        xcm=xcm+vx*dt
        ycm=ycm+vy*dt
        do i=1,n
          x_rot(i)=x_rot(i)+xcm
          y_rot(i)=y_rot(i)+ycm
        end do
        call pgbbuf
          call pgsci(0)
          call pgline(n,xold,yold)
          call pgsci(1)
          call pgline(n,x_rot,y_rot)
        call pgebuf
        xcmold=xcm
        ycmold=ycm
        do i=1,n
          xold(i)=x_rot(i)
          yold(i)=y_rot(i)
        end do
      end do
      write(*,*)'Distance=',xcm,';Hight=',ycm
      call pgclos
      end