--
-- SMN model of a 1D spring
-- uses either FD Data or Euphoria atomic data
-- John Ringland 2005/08/27
--
include SMN.e
include test.e

puts(1,"Note: There are other simulation scenarios available by changing the source code comments\n")
puts(1,"\n{ time, spring potential, equilibrium length, velocity_1, position_1, velocity_2, position_2 }\n\n")

---------------------------------------------------------
-- create and test a 1D_oscilator, which is a spring with two 1D_particles
-- using FD data
---------------------------------------------------------
---------------------------------------------------------
global function create_spring(atom k, atom dt, atom mass1, atom mass2, sequence iSV)
    SMElement m0, m1
    SystemMatrix M
    StateVector V

    m0 = {{fdZero},{fdZero}}
    m1 = {{fdOne},{fdOne}}
    M = {{m0,{{defaultFDV(k)},{fdOne}},m0,{{defaultFDV(k)},{fdOne}}, m0, {{defaultFDV(-k)},{fdOne}}},
          {m0,m1,m0,m0,m0,m0},
          {{{defaultFDV(-dt/mass1)},{fdOne}},m0, m1,m0,m0,m0},
          {{{defaultFDV(power(dt,2)/(-2*mass1))},{fdOne}}, m0, {{defaultFDV(dt)},{fdOne}}, m1, m0, m0},
          {{{defaultFDV(dt/mass2)},{fdOne}}, m0, m0, m0, m1, m0},
          {{{defaultFDV(power(dt,2)/(2*mass2))},{fdOne}}, m0, m0, m0, {{defaultFDV(dt)},{fdOne}}, m1} }
    V = {}
    for i = 1 to length(iSV) do
        V = V & {{defaultFDV(iSV[i])}}
    end for
    
    return {M,V}
end function

--Test Spring Model using FD Data
smModel spring

-- slack spring
--spring = create_spring(.1,.1,2,2,{0,2,0,-1,-0.2,1})
--test_FLModel(spring, 100, 10, .1)

-- stiff spring
--spring = create_spring(.9,.1,2,2,{0,2,0,-1,-0.2,1})
--test_FLModel(spring, 100, 10, .1)


---------------------------------------------------------
-- create and test a 1D_oscilator, which is a spring with two 1D_particles
-- using euphoria atomic data
---------------------------------------------------------
---------------------------------------------------------
global function create_scalar_spring(atom k, atom dt, atom F, atom mass1, atom mass2, sequence iSV)
    SystemMatrix M
    StateVector V

    M = {{0,k,0,k,0,-k},
          {0,1,0,0,0,0},
          {-dt/mass1,0,1-F,0,0,0},
          {power(dt,2)/(-2*mass1),0,dt,1,0,0},
          {dt/mass2,0,0,0,1-F,0},
          {power(dt,2)/(2*mass2),0,0,0,dt,1} }
    V = {}
    for i = 1 to length(iSV) do
        V = V & {iSV[i]}
    end for
    
    return {M,V}
end function

-- both particles have identical initial velocities
-- so spring moves without oscilation
--spring = create_scalar_spring(.1,.1,0.015,2,2,{0,2,-0.1,-1,-0.1,1})
--test_FLModel(spring, 100, 100, .1)

-- spring moves with oscilation and settles into new equilibrium
--spring = create_scalar_spring(.99,.1,0.02,2,2,{0,2,0,-1,0,0.9})
--test_FLModel(spring, 1000, 10, .1)

-- the power seemed to amplify until I put in some friction and used a small time step
-- that is an interesting phenomenon to explore some time, i.e. map out where and why it explodes.
-- e.g. try: 
--spring = create_scalar_spring(.99,.2,0.015,2,2,{0,2,0,-1,0,0.9})
--test_FLModel(spring, 1000, 1, .1)

-- 
-- 
puts(1,"\nthe second mass is displaced slightly inwards, the spring moves with oscilation and settles into a new equilibrium\n")
spring = create_scalar_spring(.2,.2,0.015,2,2,{0,2,0,-1,0,0.9})
test_FLModel(spring, 100, 100, .2)