      PROGRAM projetil
! Generates many .png files on the projectile trajectory
! using the xmgrace and loading a previously constructed .par file.
! Then produces an animated gif and delete all files.
      IMPLICIT NONE
      INTEGER           :: i,j,icont,c,d,u
      INTEGER           :: nframes!no. of frames in the simulation
      INTEGER           :: npts,nptstep!no. of points in a data set
      INTEGER           :: framestep!step for discarding useless frames
      INTEGER           :: imax,xtick,ytick
      REAL*8,PARAMETER  :: pi = 4.d0*datan(1.d0)
      REAL*8            :: pifrac,theta0!initial angle of the projectile velocity
      REAL*8            :: A!projectile cross section
      REAL*8            :: m!projectile mass
      REAL*8            :: dt!time discretization
      REAL*8            :: rho!air density
      REAL*8            :: CD!drag coefficient
      REAL*8            :: g!gravity
      REAL*8            :: v0,v0x,v0y!initial speed
      REAL*8            :: xmin,ymin,xmax,ymax
      REAL*8            :: t!time
      REAL*8            :: xvac,yvac!position of the projectile in the vacuum
      CHARACTER(65)     :: comando
      CHARACTER(20)     :: comando2
      npts    = 100
      nframes = 100
      pifrac  = 1.d0/5.d0
      theta0  = pi*pifrac
      A       = pi*(0.1491/2.d0)**2.d0
      m       = 42.d0
      dt      = 1.d-2
      rho     = 1.3d0
      CD      = 0.295d0
      g       = 9.8d0
      v0      = 377.d0
      v0x     = v0*dcos(theta0)
      v0y     = v0*dsin(theta0)
      xmin    = 0.d0
      ymin    = 0.d0

      imax = int( 2.d0*v0y / (g*dt) )!no. of iterations = flight time in the vacuum / dt
      
      ymax = (v0y)**2.d0 / (2.d0*g)!highest point of the projectile in the vacuum
      xmax = 2.d0*v0x*v0y / g!reach of the projectile in the vacuum
      
      framestep = max(1,imax / nframes)
      icont = 0
      do i=1,imax,framestep
         open(20,file='a.d',status='replace')
         open(21,file='b.d',status='replace')
         
         icont   = icont + 1
         nptstep = max(1,i/npts)
         do j=1,i,nptstep
            t    = dt*real(j-1,8)
            xvac = v0x*t
            yvac = v0y*t - 0.5d0*g*t*t
            write(20,*)xvac,yvac
         enddo
         write(21,*)xvac,yvac
         close(20)
         close(21)
         c = icont/100
         d = (icont - 100*c)/10
         u = icont - 100*c - 10*d
         comando='gracebat a.d b.d -param par.par -printfile txxx.png -hdevice PNG'
         comando2=' -timestamp '
         comando(45:45) = achar(48+c)
         comando(46:46) = achar(48+d)
         comando(47:47) = achar(48+u)
         CALL execute_command_line(comando)
!         CALL execute_command_line(comando//comando2)
      end do
      CALL execute_command_line("rm a.d")
      CALL execute_command_line("rm b.d")
!      comando = 'convert -delay 12 -loop 0 *.png linux.gif'!infinite loops
      comando ='convert -delay 6 -loop 1 *.png projetil.gif'!1 loop only
      CALL execute_command_line(comando)
      CALL execute_command_line("rm *.png")
      STOP
      END PROGRAM projetil
