- #1
Dylicious
- 5
- 0
Hello fellow computer physics nerds,
I'm trying to write a program to plot the positions of the three particles connected by two springs (one dimensional) in Fortran 90. I have a main program block and a module that calls a PGPLOT.
My problem is that the positions of the second and third particle overlap at some points as my code stands and I don't know why!
Here is my code please help meeeeeeeeeeeee
program springparticles
use ploty
implicit none
! Variable Declaration
real :: position_p1, position_p2, position_p3, velocity_p1, velocity_p2, velocity_p3, mass_p1, mass_p2, mass_p3, force_p1, force_p2, force_p3, dx, r
real, parameter :: spring_constant=0.1 , damping_constant=0.01 , viscous_drag=0.1
integer :: n
real, dimension(1000)::ap, bp, cp, tp
real :: t,dt
write (*,*)'Please enter the mass of the first particle: '
read *, mass_p1
write (*,*)'Second Particle: '
read *, mass_p2
write (*,*)'Third Particle: '
read *, mass_p3
position_p1=1.0
position_p2=2.0
position_p3=3.0
t= 0
dt = 0.1
dx = 0
force_p1=0
force_p2=0
force_p3=0
velocity_p1=0
velocity_p2=0
velocity_p3=0
print *,'Particle One |||||| Particle Two |||||| Particle Three'
do n = 1,1000
t = t +dt
dx = dx + 0.01
! force_p1 = force_p1-viscous_drag*velocity_p1
! force_p2 = force_p2-viscous_drag*velocity_p2
r= 0.5
force_p1 = spring_constant*(position_p2-position_p1-r) - damping_constant*(velocity_p1-velocity_p2)
force_p2 = - spring_constant*(position_p2-position_p1-r) - damping_constant*(velocity_p1-velocity_p2) + spring_constant*(position_p3-position_p2-r) - damping_constant*(velocity_p2-velocity_p3)
force_p3 = -spring_constant*(position_p3-position_p2-r) - damping_constant*(velocity_p2-velocity_p3)
velocity_p1 = velocity_p1+force_p1/mass_p1*dt
velocity_p2 = velocity_p2+force_p2/mass_p1*dt
velocity_p3 = velocity_p3+force_p3/mass_p3*dt
position_p1 = position_p1+velocity_p1*dt
position_p2 = position_p2+velocity_p2*dt
position_p3 = position_p3+velocity_p3*dt
! velocity_p1 = velocity_p1+dx/dt
! velocity_p2 = velocity_p2+force_p2/mass_p1p2*dt
print *,position_p1,'|| ',position_p2,'|| ',position_p3
tp(n) = t
ap(n)= position_p1
bp(n)= position_p2
cp(n)= position_p3
end do
call plotfunction(tp, ap, bp, cp)
end program springparticles
module ploty
implicit none
contains
!-----------------------------------------------------------------------
! This module uses pgplot to plot the positions of the particles
!-----------------------------------------------------------------------
subroutine plotfunction(t, a, b, c)
implicit none
integer, parameter :: d = 100
real :: x(100), y(100)
real, intent(in) :: t(:), a(:), b(:), c(:)
integer :: pgopen
! Open a plot window
IF (PGOPEN('/XWINDOW') .LE. 0) STOP
! Set-up plot axes
call PGENV(minval(t),maxval(t),min(minval(a),minval(b)),max(maxval(c),maxval(b)),0,0)
call PGLAB('Time', 'Distance from Equilibrium', 'Positions of Particles')
! Change plot colour to colour 6 ()
! Compute the function at the points
! x(d) = t(:)
! y(d) = a(:)
! write(6,*) t(::10)
! write(6,*) a(::10)
! Plot the curve
call PGSCI(6)
call PGLINE(size(a),t,a)
call PGSCI(4)
call PGLINE(size(b),t,b)
call PGSCI(7)
call PGLINE(size(c),t,c)
! Pause and then close plot window
call PGCLOS
end subroutine plotfunction
!-----------------------------------------------------------------------
end module ploty
I'm trying to write a program to plot the positions of the three particles connected by two springs (one dimensional) in Fortran 90. I have a main program block and a module that calls a PGPLOT.
My problem is that the positions of the second and third particle overlap at some points as my code stands and I don't know why!
Here is my code please help meeeeeeeeeeeee
program springparticles
use ploty
implicit none
! Variable Declaration
real :: position_p1, position_p2, position_p3, velocity_p1, velocity_p2, velocity_p3, mass_p1, mass_p2, mass_p3, force_p1, force_p2, force_p3, dx, r
real, parameter :: spring_constant=0.1 , damping_constant=0.01 , viscous_drag=0.1
integer :: n
real, dimension(1000)::ap, bp, cp, tp
real :: t,dt
write (*,*)'Please enter the mass of the first particle: '
read *, mass_p1
write (*,*)'Second Particle: '
read *, mass_p2
write (*,*)'Third Particle: '
read *, mass_p3
position_p1=1.0
position_p2=2.0
position_p3=3.0
t= 0
dt = 0.1
dx = 0
force_p1=0
force_p2=0
force_p3=0
velocity_p1=0
velocity_p2=0
velocity_p3=0
print *,'Particle One |||||| Particle Two |||||| Particle Three'
do n = 1,1000
t = t +dt
dx = dx + 0.01
! force_p1 = force_p1-viscous_drag*velocity_p1
! force_p2 = force_p2-viscous_drag*velocity_p2
r= 0.5
force_p1 = spring_constant*(position_p2-position_p1-r) - damping_constant*(velocity_p1-velocity_p2)
force_p2 = - spring_constant*(position_p2-position_p1-r) - damping_constant*(velocity_p1-velocity_p2) + spring_constant*(position_p3-position_p2-r) - damping_constant*(velocity_p2-velocity_p3)
force_p3 = -spring_constant*(position_p3-position_p2-r) - damping_constant*(velocity_p2-velocity_p3)
velocity_p1 = velocity_p1+force_p1/mass_p1*dt
velocity_p2 = velocity_p2+force_p2/mass_p1*dt
velocity_p3 = velocity_p3+force_p3/mass_p3*dt
position_p1 = position_p1+velocity_p1*dt
position_p2 = position_p2+velocity_p2*dt
position_p3 = position_p3+velocity_p3*dt
! velocity_p1 = velocity_p1+dx/dt
! velocity_p2 = velocity_p2+force_p2/mass_p1p2*dt
print *,position_p1,'|| ',position_p2,'|| ',position_p3
tp(n) = t
ap(n)= position_p1
bp(n)= position_p2
cp(n)= position_p3
end do
call plotfunction(tp, ap, bp, cp)
end program springparticles
module ploty
implicit none
contains
!-----------------------------------------------------------------------
! This module uses pgplot to plot the positions of the particles
!-----------------------------------------------------------------------
subroutine plotfunction(t, a, b, c)
implicit none
integer, parameter :: d = 100
real :: x(100), y(100)
real, intent(in) :: t(:), a(:), b(:), c(:)
integer :: pgopen
! Open a plot window
IF (PGOPEN('/XWINDOW') .LE. 0) STOP
! Set-up plot axes
call PGENV(minval(t),maxval(t),min(minval(a),minval(b)),max(maxval(c),maxval(b)),0,0)
call PGLAB('Time', 'Distance from Equilibrium', 'Positions of Particles')
! Change plot colour to colour 6 ()
! Compute the function at the points
! x(d) = t(:)
! y(d) = a(:)
! write(6,*) t(::10)
! write(6,*) a(::10)
! Plot the curve
call PGSCI(6)
call PGLINE(size(a),t,a)
call PGSCI(4)
call PGLINE(size(b),t,b)
call PGSCI(7)
call PGLINE(size(c),t,c)
! Pause and then close plot window
call PGCLOS
end subroutine plotfunction
!-----------------------------------------------------------------------
end module ploty