class Complex {

                  public:

                                                                                        var var1,var2,var3, Re, Im;

                                                                                        Multiply(var1,var2);

                                                                                        Divide(var1,var2);

                                                                                        Real(var1);

                                                                                        Imaginary(var1);

                                                                                        Conjugate(var1);

                                                                                        RealToComplex(var1);

}

Complex::Multiply(v1,v2)

{

                  var3 = vector(v1.x*v2.x - v1.y*v2.y, v1.x*v2.y + v2.x*v1.y, 0);

                  return var3;

}

Complex::Divide(v1,v2)

{

                  var denominator = pow(v2.x,2) + pow(v2.y,2);

                  var3 = vector((v1.x*v2.x + v1.y*v2.y)/denominator, (-v1.x*v2.y + v2.x*v1.y)/denominator, 0);

                  return var3;

}

Complex::Real(v1)

{

                  Re = v1.x;

                  return Re;

}

Complex::Imaginary(v1)

{

                  Im = v1.y;

                  return Im;

}

Complex::Conjugate(v1)

{

                  v1.y = -v1.y;

                  return v1;

}

Complex::RealToComplex(v1)

{

                  if(typeof(v1)!=2)

                  {

                                   return null;

                  }

                  return vector(v1,0,0);

}

main(doc,op)

{

                  // Step 1

                  var obj = op->GetDown();

                  var frame = doc->GetTime()->GetFrame(doc->GetFps());

                  frame = frame/doc->GetFps();

                  var i,j,pos;

                  var m = op#ID_USERDATA:1; // spring mass

                  var c = op#ID_USERDATA:2; // damping coefficient

                  var k = op#ID_USERDATA:3; // spring coefficient

                  var x0 = op#ID_USERDATA:4; // initial displacement

                  var f1 = op#ID_USERDATA:5; // sine harmonic force

                  var f2 = op#ID_USERDATA:6; // cosine harmonic force

                  var w = op#ID_USERDATA:7; // harmonic frequency

                  var v0 = op#ID_USERDATA:8; // initial velocity

                  // For complex numbers

                  var cM = vector(m,0,0);

                  var cC = vector(c,0,0);

                  var cK = vector(k,0,0);

                  var cX0 = vector(x0,0,0);

                  var cF1 = vector(f1,0,0);

                  var cF2 = vector(f2,0,0);

                  var cW = vector(w,0,0);

                  var cV0 = vector(v0,0,0);

                  // Step 2

                  var a = new(Complex);

                  var var1 = vector(f1,f2,0);

                  var var2 = vector(k-m*pow(w,2),c*w,0);

                  var var3 = a->Divide(var1,var2);

                  var Re = a->Real(var3);

                  var Im = a->Imaginary(var3);

                  var cA = var3;

                  var tmp = pow(c,2) - 4.0*m*k;

                  var r = vector(0,0,0);

                  if(tmp >= 0)

                  {

                                   r = vector(sqrt(tmp), 0, 0);

                  }

                  else

                  {

                                   r = vector(0, sqrt(-tmp), 0);

                  }

                  var s1 = (-cC+r)/(2.0*m);

                  var s2 = (-cC-r)/(2.0*m);

                  // Step 3

                  var X0 = Re;

                  var V0 = a->Imaginary(a->Multiply(var3,vector(0,w,0)));

                  // Step 4

                  var denominator = a->Real(a->Multiply(s2-s1, a->Conjugate(s2-s1)));

                  var tmp1 = a->Multiply(s2, a->RealToComplex(x0-V0));

                  var tmp2 = a->Multiply(s1, a->RealToComplex(x0-V0));

                  var numerator1 = a->Multiply((tmp1 - a->RealToComplex(v0-V0)),a->Conjugate(s2-s1));

                  var numerator2 = a->Multiply((a->RealToComplex(v0-V0) - tmp2),a->Conjugate(s2-s1));

                  var p1 = numerator1/denominator;

                  var p2 = numerator2/denominator;

                  // Step 5

                  var X = a->Real(a->Multiply(cA, vector(cos(w*frame),sin(w*frame),0)));

                  var V = a->Real(a->Multiply(cA, vector(-w*sin(w*frame),w*cos(w*frame),0)));

                  // Step 6

                  var xx = X + a->Real(a->Real(a->Multiply(p1,a->Multiply(a->RealToComplex(exp(s1.x*frame)), vector(cos(s1.y*frame),sin(s1.y*frame),0)))) + a->Multiply(p2,a->Multiply(a->RealToComplex(exp(s2.x*frame)),vector(cos(s2.y*frame),sin(s2.y*frame),0))));

                  // Step 7: Set the new y-position of the selected point

                  var p_sel = obj->GetPointSelection();

                  var p_sel_cnt = p_sel->GetCount();

                  var pnt = obj->GetPoints();

                  var pnt_cnt = obj->GetPointCount();

                  if(pnt_cnt == 0 || p_sel_cnt == 0) return false;

                  var pos;

                  for(i=0 ; i<pnt_cnt ; i++)

                  {

                                   pos = pnt[i];

                                   if(p_sel->IsSelected(i))

                                   {

                                                     pos.y = xx;

                                   }

                                   else

                                   {

                                                     pos.y = 0;

                                   }

                                   pnt[i] = pos;

                  }

                  obj->SetPoints(pnt);

  obj->Message(MSG_UPDATE);

}