Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Fortran Simple Harmonic Oscillator Problem

  1. Dec 8, 2011 #1
    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!!! :cry:

    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


    t= 0
    dt = 0.1

    dx = 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
    ! 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

    ! 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
  2. jcsd
  3. Dec 9, 2011 #2


    Staff: Mentor

    Is there some reason that you think they shouldn't overlap? I don't know what values you're putting in for the masses, and with some values for these, plus the values of the spring and damping constants and viscous drag, it might be that two or more of the particles collide from time to time.
  4. Dec 9, 2011 #3
    I just figured that conceptually the spring would compress to a point and never allow the particles to move past one another, mathematically I see how its possible but doesn't that just mean I'm modelling the physical situation inaccurately?
  5. Dec 12, 2011 #4


    Staff: Mentor

    I don't recall seeing anything that sets a minimum length for the springs. Your program doesn't understand the physical characteristics of springs such as minimum length when it is compressed, or maximum length at which the spring deforms or breaks when it is stretched. If you want your program to model that behavior, you need to have code that does this.
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook