readability.f90

!-----------
! Some code readability/performance improvements
!   _ indices don't have to be called with different names if they are not nested (i,j,k,l,... are most common names);
!   _ treating the boundary points should not be done inside of the loop, because the same operation will be repeated without being
!   necessary; it also improves the ease of read, since it is more logical;
!   _ defining a 'program' allows to have a better defined structure, and to use constants instead of passing them;
!   _ always use 'implicit none'! Is saves the situation sometimes;
!   _ if you are not defining a variable at runtime, better to declare it a parameter (=constant).



program main

  implicit none

  real, parameter :: k = 100.0
  real, parameter :: a = 10.0
  real, parameter :: h = 0.001
  real, parameter :: mass = 1.0

  real::x(1000),p(1000),k1(1000),k2(1000),k3(1000),k4(1000),l1,l2,l3,l4,t,tf
  integer::i,j,m

  write(*,*) "serial no of middlemost particle,time of observation"
  read(*,*) m,tf


  !define initial momentum
  do i=1,2*m+1
    p(i)=0
  end do
  !define initial position(gaussian tilt in middle zone)
  do i=1,m-10
    x(i)=0 
  end do  
  do i=m+10,2*m+1
    x(i)=0 
  end do  
  do i=m-9,m+9
    x(i)=exp(-(0.1*(i-m))**2)
  end do  


  t=0

  do while(t.le.tf)  !stop at time of observation

    do i=2,2*m

      !(k1=first runge kutta coefficient for ith particle moving under force field of (i-1) and (i+1)th particle,k2=second and so on)
      !started with j=2, not 1, as, end particles are fixed, same for 2m+1th particle. k(1),k(2m+1) zero for same reason 

      k1(1)=0
      k1(2*m+1)=0 
      do j=2,2*m
        k1(j)=f(k,a,x(j-1),x(j),x(j+1)) 
      end do

      k2(1)=0
      k2(2*m+1)=0
      do j=2,2*m
        k2(j)=f(k,a,x(j-1)+1*h*k1(j-1)/2,x(j)+h*k1(j)/2,x(j+1)+1*h*k1(j+1)/2)    
      end do

      k3(1)=0
      k3(2*m+1)=0
      do j=2,2*m
        k3(j)=f(k,a,x(j-1)+1*h*k2(j-1)/2,x(j)+h*k2(j)/2,x(j+1)+1*h*k2(j+1)/2)    
      end do

      k4(1)=0
      k4(2*m+1)=0
      do j=2,2*m
        k4(j)=f(k,a,x(j-1)+1*h*k3(j-1)/2,x(j)+h*k3(j)/2,x(j+1)+1*h*k3(j+1)/2)
      end do

      l1=g(p(i))
      l2=g((p(i)+h*l1/2))
      l3=g(p(i)+h*l2/2)
      l4=g(p(i)+h*l3)

      x(i)=x(i)+h*(l1+2*l2+2*l3+l4)/(6*mass)        !final runge kutta result
      p(i)=p(i)+h*(k1(i)+2*k2(i)+2*k3(i)+k4(i))/6
      t=t+h
    end do !for all particles(i) at time t
  end do !for time=tf=time of observation



  !write to file

  open(1,file="f.dat")

  do i=1,2*m+1
    write(1,*) i,x(i)
  end do

  close(1)

contains
  !f=dp/dt=Force,force on nth particle=f(k,a,x[n-1],x[n],x[n+1]),
  !1d lattice(2m+1) particles, all particles connected with each other in harmonic potential, with slight x^3 nonlinearity.

  function f(k,a,q,y,z) result(s)
    real::k,a,q,y,z,s
    s=-k*((y-q)+(y-z))-a*((y-q)**3+(y-z)**3)
  end function f

  !g=dx/dt=p/mass(mass=1),p=momentum, particles at lattice ends cant move

  function g(q) result(s)
    real::s,q
    s=q
  end function g


end program main