      PROGRAM projetil
!Animates a projectile in vacuum using xmgrace
!Firstly, compile the FORTRAN program
!>gfortran filename.f90 -o filename.exe
!Then run the filename.exe program feeding xmgrace
!>./filename.exe | xmgrace -pipe -free
      IMPLICIT NONE
      INTEGER           :: i,j
      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
      npts    = 100
      nframes = 250
      pifrac  = 1.d0/4.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 vacuum / dt
      
      ymax = (v0y)**2.d0 / (2.d0*g)!highest point of the projectile in vacuum
      xmax = 2.d0*v0x*v0y / g!reach of the projectile in vacuum
      
      framestep = max(1,imax / nframes)
      
      xtick = int(xmax/5.d0) + 1
      ytick = int(ymax/5.d0) + 1
      
      write(*,*)"@focus off"!write the xmgrace command lines
      write(*,*)"@g0 on"
      write(*,*)"@with g0"
      write(*,*)'@xaxis label "x (m)"'
      write(*,*)'@yaxis label "y (m)"'
      write(*,*)"@ xaxis ticklabel on"
      write(*,*)"@ yaxis ticklabel on"
      write(*,*)"@ xaxis tick on"
      write(*,*)"@ yaxis tick on"
!      write(*,*)"@ xaxis tick major",xtick
      write(*,*)"@ xaxis tick major 3000"
      write(*,*)"@ xaxis tick minor ticks 1"
!      write(*,*)"@ yaxis tick major",ytick
      write(*,*)"@ yaxis tick major 1000"
      write(*,*)"@ yaxis tick minor ticks 1"
!      write(*,*)"@world xmin ",int(xmin)
      write(*,*)"@world xmax ",(1+int(xmax/1.d3))*1000
      write(*,*)"@world xmax 15000"
      write(*,*)"@world ymin ",int(ymin)
!      write(*,*)"@world ymax 4000"
      write(*,*)"@world ymax ",(1+int(ymax/1.d3))*1000
      
      write(*,*)"@title ",'"Fisica Computacional - Projeto 3"'
      write(*,*)"@title font 0"
      write(*,*)"@title size 1.200000"
      write(*,*)"@title color 1"
      write(*,2)"@subtitle ",'"m=',int(m),' Kg, Cd=',CD,', g=',g,' m/s\S2\N, \xr\0=',rho,' Kg/m\S3"'
2     format(a10,a3,I2,a8,f5.3,a4,f3.1,a17,f3.1,a9)
      write(*,*)"@subtitle font 0"
      write(*,*)"@subtitle size 1.000000"
      write(*,*)"@subtitle color 1"
      
      do i=1,imax,framestep
         write(*,*)"@with string"
         write(*,*)"@string on"
         write(*,*)"@string loctype view"
         write(*,*)"@string 0.17, 0.2"
         write(*,*)"@string color 1"
         write(*,*)"@string rot 0"
         write(*,*)"@string font 0"
         write(*,*)"@string just 0"
         write(*,*)"@string char size 1.000000"
         write(*,4)"@string def ",'"\xq\0\s0\N = ',pifrac,'\xp"'
4        format(a12,a13,f4.2,a4)
         
         write(*,*)"@with string"
         write(*,*)"@string on"
         write(*,*)"@string loctype view"
         write(*,*)"@string 0.18, 0.8"
         write(*,*)"@string color 1"
         write(*,*)"@string rot 0"
         write(*,*)"@string font 0"
         write(*,*)"@string just 0"
         write(*,*)"@string char size 1.000000"
         write(*,1)"@string def",' "t =',real(i-1)*dt,' s"'
1         format(a11,a5,f5.1,a3)
         
         write(*,*)"@S0 on"!data set for the track of free projectile
         write(*,*)"@type xy"
         write(*,*)"@s0 legend",'"vacuo"'
         write(*,*)"@S0 line type 1"
         write(*,*)"@S0 linestyle 1"
         write(*,*)"@S0 linewidth 1.0"
         write(*,*)"@S0 line color 1"
         write(*,*)"@S0 symbol 0"
         write(*,*)"@S0 symbol color 1 "
         write(*,*)"@S0 symbol size .25"
         write(*,*)"@S1 on"!data set for the free projectile
         write(*,*)"@type xy"
         write(*,*)"@S1 line type 0    "
         write(*,*)"@S1 symbol 1       "
         write(*,*)"@S1 symbol color 1 "
         write(*,*)"@S1 symbol fill pattern 1"
         write(*,*)"@S1 symbol fill color 1"
         write(*,*)"@S1 symbol size .5"
         
         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(*,*)"@g0.S0 point ",xvac,',',yvac
         enddo
         write(*,*) "@g0.S1 point ",xvac,',',yvac
         
         write(*,*)"&"
         write(*,*)"@redraw"
         write(*,*)"@sleep 0.02"!wait this number of seconds
         if(i.le.imax-framestep)then
            write(*,*)"@kill S0"
            write(*,*)"@kill S1"
            write(*,*)'@clear string'
         endif
      end do
      STOP
      END PROGRAM projetil
