k=1000.0
m=1.0
c=1.0
u0=0.0
v0=1.0
f={sin(2.0*sqrt(k/m)*it)}; f={0.0}

g=1.0/2.0 // gamma
b=1.0/6.0 // beta

om=sqrt(k/m); xsi=c/(2*m*om)
dtcr=1/om*( xsi*(g-0.5)+sqrt(0.5*g-b+(xsi*xsi*(g-0.5)*(g-0.5))) )/(0.5*g-b)

dt=0.1*dtcr
N=1000

u=new double[N+1]; u[0]=u0
v=new double[N+1]; v[0]=v0
a=new double[N+1]; a[0]=(f(0)-c*v[0]-k*u[0])/m 

//beta-Newmark coefficients
w0=1/(b*dt*dt);	w1=g/(b*dt);			w2=1/(b*dt);	w3=1.0/(2*b)-1.0
w4=g/b-1.0;		w5=(g/b-2.0)*dt/2.0;	w6=dt*(1.0-g);	w7=dt*g

// K_effective
kef=k+w0*m+w1*c

// iterative computations
(1..N).each{
	// effective load
	fef=f(it*dt)+m*(w0*u[it-1]+w2*v[it-1]+w3*a[it-1])+c*(w1*u[it-1]+w4*v[it-1]+w5*a[it-1])

	// solve for current displacement
	u[it]=fef/kef

	// current acceleration and velocity
	a[it]=w0*(u[it]-u[it-1])-w2*v[it-1]-w3*a[it-1]
	v[it]=v[it-1]+w6*a[it-1]+w7*a[it]
}

respPlot= new PlotFrame()
respPlot.addFunction(new plotfunction(dt,u))
//respPlot.addFunction(new plotfunction(u,v))
respPlot.setMarker(false)
respPlot.show()

/*
import gendomain.*
theDomain=theUniverse.GENDomain()

aNode=new Node(1,0.0)
aNode.setNumDOFs_Steps(2,u.length)
theDomain.putNode(aNode)
(0..<u.size()).each{
	udof = [u[it], v[it]] as double[]
	theDomain.updateNodes(udof, it) 
}


theGP.stop()

theUniverse.getGENDomain(1).ColorNodesDEF=Color.blue
theUniverse.getGENDomain(1).fillNodes=false

if(theUniverse.getGENDomain(2)!=null){
	theUniverse.getGENDomain(2).ColorNodesDEF=Color.red
	theUniverse.getGENDomain(2).fillNodes=false
}

theGP.set_size(-0.05,0.05,-1.2,1.2)
theGP.trace=true
theGP.nodeScale=100
theGP.plotDeform(1, 1.0, 0,u.length,1)

*/

