//---------------------------------------------------------------------------

#include "stdafx.h"
#include "Solver.h"

//---------------------------------------------------------------------------


nbsolverobj solver;


class rndgeneratorobj
{
  __int32 y1rnd,y2rnd,y3rnd;
  __int32 rndvalue;

  public:

  void init(__int32 rndseed);

  mfloat random();
};


mfloat rndgeneratorobj::random()
{
  static const __int32 m1rnd=30269;
  static const __int32 m2rnd=30307;
  static const __int32 m3rnd=30323;
  static const __int32 k1rnd=171;
  static const __int32 k2rnd=172;
  static const __int32 k3rnd=170;

  mfloat x1,x2,x3;
  double buf;


//  rndvalue=(((rndvalue+1)*75)%65537)-1;
//  return (mfloat)rndvalue/65536.0;


  y1rnd=(k1rnd*y1rnd)%m1rnd;
  y2rnd=(k2rnd*y2rnd)%m2rnd;
  y3rnd=(k3rnd*y3rnd)%m3rnd;

  x1=(mfloat)y1rnd/(mfloat)m1rnd;
  x2=(mfloat)y2rnd/(mfloat)m2rnd;
  x3=(mfloat)y3rnd/(mfloat)m3rnd;

  return modf(x1+x2+x3,&buf);
}

void rndgeneratorobj::init(__int32 rndseed)
{
  int i;

  y1rnd=1023;
  y2rnd=11;
  y3rnd=3007;

  rndvalue=4096+31;

  for (i=1; i<rndseed; i++)
    random();
}

rndgeneratorobj rndgenerator;



inline mfloat r()
{
  return 1.0-2.0*rndgenerator.random();
}

inline mfloat r01()
{
  return rndgenerator.random();
}



inline void mad(vec &x, const vec &y, mfloat my)
{
  x.x+=y.x*my;
  x.y+=y.y*my;
  x.z+=y.z*my;

  return;

/* 
//  /G7 /O3 -Qip -QxW 

  __m128d r0 = _mm_load_pd((const double*)&y);
  __m128d r1 = _mm_load_sd(&y.z);
  __m128d r3 = _mm_load1_pd(&my); 

  r0 = _mm_mul_pd(r0,r3);
  r1 = _mm_mul_sd(r1,r3);

  __m128d x1 = _mm_load_pd((const double*)&x);
  __m128d x2 = _mm_load_sd(&x.z);

  r0 = _mm_add_pd(r0,x1);
  r1 = _mm_add_sd(r1,x2);

  _mm_store_pd((double*)&x,r0);
  _mm_store_sd(&x.z,r1);
*/

  //  return;

/* 
  __asm
  {
    mov eax,x
    mov ecx,y

    movapd xmm0,[ecx]

    movsd xmm3,my
    shufpd xmm3,xmm3,0

    movsd xmm1,[ecx+16]
    mulpd xmm0,xmm3
    mulsd xmm1,xmm3
    addpd xmm0,[eax]
    addsd xmm1,[eax+16]

    movapd [eax],xmm0
    movsd [eax+16],xmm1
  }
*/
}

inline void sub(vec &x, vec &y, vec &z)
{
  x.x=y.x-z.x;
  x.y=y.y-z.y;
  x.z=y.z-z.z;

/* 
  __asm
  {
    mov eax,x
    mov ecx,y
    mov edx,z

    movapd xmm0,[ecx]
    movsd xmm1,[ecx+16]
    subpd xmm0,[edx]
    subsd xmm1,[edx+16]
    
    movapd [eax],xmm0
    movsd [eax+16],xmm1 
  }
*/
}

inline void mul(vec &x, mfloat my)
{
  x.x*=my;
  x.y*=my;
  x.z*=my;
}


inline void coeffs(vec &x, mfloat &sqlen, mfloat &nnlen, mfloat &invlen)
{
  sqlen=SQ(x.x)+SQ(x.y)+SQ(x.z)+0.0000000001;
  nnlen=sqrt(sqlen);
  invlen=1.0/nnlen;
  sqlen=invlen*invlen;
}

inline mfloat dot(vec &x, vec &y)
{
  return x.x*y.x+x.y*y.y+x.z*y.z;
}

/*
2 order of accuracy
void nbsolverobj::step(mfloat step)
{
  int i;
  mfloat newcine1=0;
  mfloat newcine2=0;
  mfloat energy;
  mfloat dt=step;
  mfloat idt;
  mfloat k=200.0*(1.0-1.0/vfactor);
  idt=1.0/dt;
  vec imp;
  vec idvel;
  vec nvel;

  imp.x=0;
  imp.y=0;
  imp.z=0;
  imp.w=0;
 

  for (i=0; i<len; i++)
  {
    if (!pts[i].moveable) continue;

    nvel=pts[i].acc;
    idvel=pts[i].vel;
    mad(idvel,nvel,dt);
    mad(nvel,pts[i].vel,2.0*k+idt);
    mul(nvel,1.0/(idt-2.0*k));
    pts[i].vel=nvel;
    mad(pts[i].pos,pts[i].vel,dt);
 
    energy=0.5*pts[i].mas*dot(pts[i].vel,pts[i].vel);
    newcine1+=energy;
    energy=0.5*pts[i].mas*dot(idvel,idvel);
    newcine2+=energy;
    mad(imp,pts[i].vel,pts[i].mas);
  }
 
  energyleak=(newcine1-newcine2); // kinetic energy loss due to friction
  potenergy+=-newcine1+cinenergy; // potential energy increase
  cinenergy=newcine2;
  impulse=imp;
}
*/

void nbsolverobj::step(mfloat step)
{
    int i;
  mfloat newcine1=0;
  mfloat newcine2=0;
  mfloat energy;
  vec imp;

  imp.x=0;
  imp.y=0;
  imp.z=0;

    for (i=0; i<len; i++)
    {
    if (!pts[i].moveable) continue;

      mad(pts[i].vel,pts[i].acc,step);

    energy=0.5*pts[i].mas*dot(pts[i].vel,pts[i].vel);
    newcine1+=energy;
    newcine2+=energy*SQ(vfactor);

    mul(pts[i].vel,vfactor);
    mad(pts[i].pos,pts[i].vel,step);

    mad(imp,pts[i].vel,pts[i].mas);
    }

  energyleak+=(newcine1-newcine2); // kinetic energy loss due to friction
  potenergy+=-newcine1+cinenergy; // potential energy increase
  cinenergy=newcine2;
  impulse=imp;
}

void nbsolverobj::evaluate()
{
    int i,j;
    vec dir,cpos,cacc;
    mfloat sqlen,coeff,nnlen,invlen,crad,cmas,acc1;

    for (i=0; i<len; i++)
    {
      cacc.x=0;
      cacc.y=0;
      cacc.z=0;

      cpos.x=pts[i].pos.x;
      cpos.y=pts[i].pos.y;
      cpos.z=pts[i].pos.z;

    crad=pts[i].rad;
    cmas=pts[i].mas;

      for (j=0; j<i; j++)
      {
      sub(dir,pts[j].pos,cpos);

          coeffs(dir,sqlen,nnlen,invlen);
        coeff=sqlen*(nnlen-crad-pts[j].rad)*
           (quadfactor*sqlen+linlfactor*invlen);

      mul(dir,coeff);
      mad(cacc,dir,pts[j].mas);
      mad(pts[j].acc,dir,-cmas);
       }

      pts[i].acc.x=cacc.x;
      pts[i].acc.y=cacc.y;
      pts[i].acc.z=cacc.z;
    }
}


void nbsolverobj::drawobjects(sceneguidt sceneguid, float time)
{
  sphguidt sphguid;
  int i;

  initlight(gmodeltype,sceneguid,time);
//  clearscene(sceneguid);

  for (i=0; i<len; i++)
  {
    sphguid=sceneallocatesphere(sceneguid);

    setspherepos3f(sphguid,15000.0*pts[i].pos.x,
                           15000.0*pts[i].pos.y,
                           15000.0*pts[i].pos.z);
    setsphereradius(sphguid,15000.0*pts[i].prd);
    setspherecolorfp(sphguid,pts[i].clr);
    setspherevisible(sphguid,true);

    if (gmodeltype==6)
    {
      if (i<3)
      {
        enablesphereownbgillum(sphguid,true);
        setspherebgillum3f(sphguid,1,1,1);
      }
    }
    else
      if (gmodeltype==5)
        if (i==0)
        {
          enablesphereownbgillum(sphguid,true);
          setspherebgillum3f(sphguid,1,1,1);
        }
  }
}

mfloat nbsolverobj::precision()
{
    int i,j;
    vec dir,cpos,cacc;
    mfloat sqlen,coeff,nnlen,invlen,crad,acc1;
  mfloat merr=0;
  mfloat berr=0.0000000001;


    for (i=0; i<len; i++)
    {
      cacc.x=0;
      cacc.y=0;
      cacc.z=0;

      cpos.x=pts[i].pos.x;
      cpos.y=pts[i].pos.y;
      cpos.z=pts[i].pos.z;

    crad=pts[i].rad;


      for (j=0; j<len; j++)
      {
      if (i==j) continue;

      dir.x=pts[j].pos.x-cpos.x;
      dir.y=pts[j].pos.y-cpos.y;
      dir.z=pts[j].pos.z-cpos.z;

      sqlen=SQ(dir.x)+SQ(dir.y)+SQ(dir.z)+0.0000000001;

      nnlen=sqrt(sqlen);
      invlen=1.0/nnlen;
      coeff=invlen*invlen*(nnlen-crad-pts[j].rad)*
      (quadfactor*invlen*invlen+linlfactor*invlen);
      dir.x*=coeff;
      dir.y*=coeff;
      dir.z*=coeff;
      cacc.x+=dir.x*pts[j].mas;
      cacc.y+=dir.y*pts[j].mas;
      cacc.z+=dir.z*pts[j].mas;
      }

    merr+=SQ(pts[i].acc.x-cacc.x)+
          SQ(pts[i].acc.y-cacc.y)+
          SQ(pts[i].acc.z-cacc.z);

    berr+=SQ(pts[i].acc.x)+
             SQ(pts[i].acc.y)+
          SQ(pts[i].acc.z);
  }

  return sqrt(merr/berr);
}


void nbsolverobj::setmodel1()
{
  int i;

/////////////////////////////////////////////////////////////////////////////
// This is a model of interacting bodies. Loss by friction is accounted, 
// i.e. the bodies move in a viscous medium. The interaction is implemented 
// according to the inverse-square law, and at small distances attraction 
// is replaced by repulsion.
// In course of the interaction the bodies gather around the common center 
// of mass, small objects being closer to the center of mass, 
// - such disposition is dictated by the minimum of the potential energy.
// As you can see, the potential energy stabilizes.
//

  len=120;
  quadfactor=1.0f;
  linlfactor=0.0f;
  vfactor=0.995;

  for (i=0; i<len; i++)
  {
    float rad=r01()*r01()*r01()*3.0f+0.01;

    pts[i].clr[0]=0.8;
    pts[i].clr[1]=0.5;
    pts[i].clr[2]=0.2;

    pts[i].pos.x=r()*8.0f;
    pts[i].pos.y=r()*8.0f;
    pts[i].pos.z=r()*8.0f;

    pts[i].vel.x=r()*8.0f;;
    pts[i].vel.y=r()*8.0f;;
    pts[i].vel.z=r()*8.0f;;

    pts[i].rotaxis[0]=r();
    pts[i].rotaxis[1]=r();
    pts[i].rotaxis[2]=r();
    pts[i].moveable=true;

    pts[i].rad=rad*3.0;
    pts[i].prd=rad;
    pts[i].mas=500.0f*rad*rad*rad;
  }
}

void nbsolverobj::setmodel2()
{
  int i;

/////////////////////////////////////////////////////////////////////////////
// This is a model of interacting bodies. Loss by friction is accounted, 
// i.e. the bodies move in a viscous medium. The interaction is implemented 
// according to the inverse-radius law, and at small distances attraction 
// is replaced by repulsion.
// The configuration of the objects becomes spherical, - such disposition 
// corresponds to the minimum of the potential energy.
// 

  len=110;
  quadfactor=0.0f;
  linlfactor=8.0f;
  vfactor=0.995;

  for (i=0; i<len; i++)
  {
    pts[i].clr[0]=1;
    pts[i].clr[1]=1;
    pts[i].clr[2]=1;

      pts[i].pos.x=r()*1.0f;
      pts[i].pos.y=r()*1.0f;
      pts[i].pos.z=r()*1.0f;

    pts[i].vel.x=0;
      pts[i].vel.y=0;
      pts[i].vel.z=0;

    pts[i].rotaxis[0]=r();
    pts[i].rotaxis[1]=r();
    pts[i].rotaxis[2]=r();
    pts[i].moveable=true;

    float rad=0.6;
    pts[i].rad=rad*1.2;
    pts[i].prd=rad;
      pts[i].mas=100.0f*rad*rad*rad;
  }

  for (i=0; i<40; i++)
  {
    pts[i].clr[0]=0.2;
    pts[i].clr[1]=1;
    pts[i].clr[2]=1;

      pts[i].pos.x=r()*1.0f;
      pts[i].pos.y=r()*1.0f;
      pts[i].pos.z=r()*1.0f;

    pts[i].vel.x=0;
      pts[i].vel.y=0;
      pts[i].vel.z=0;

    pts[i].rotaxis[0]=r();
    pts[i].rotaxis[1]=r();
    pts[i].rotaxis[2]=r();
    pts[i].moveable=true;

    mfloat rad=0.25;
    pts[i].rad=rad*10.0;
    pts[i].prd=rad;
      pts[i].mas=2000.0f*rad*rad*rad;
  }
}

void nbsolverobj::setmodel3()
{

/////////////////////////////////////////////////////////////////////////////
// This is a model of interacting bodies. Loss by friction is absent. 
// The interaction is implemented according to the inverse-radius law, 
// and at small distances attraction is replaced by repulsion.
// The configuration of the objects doesn't stabilize with time.
//

  len=80;
  quadfactor=0.0f;
  linlfactor=1.0f;
  vfactor=1.0;

  for (int i=0; i<8; i++)
  {
    float c1=r()*4.0f,c2=r()*4.0f,c3=r()*4.0f;

    for (int j=0; j<10; j++)
    {
      mfloat rad=r01()*r01()+0.1;

      pts[i*10+j].clr[0]=1;
      pts[i*10+j].clr[1]=1;
      pts[i*10+j].clr[2]=1;

      pts[i*10+j].vel.x=0;
      pts[i*10+j].vel.y=0;
      pts[i*10+j].vel.z=0;

      pts[i*10+j].pos.x=c1+r()*0.5f;
      pts[i*10+j].pos.y=c2+r()*0.5f;
      pts[i*10+j].pos.z=c3+r()*0.5f;

      pts[i*10+j].rotaxis[0]=r();
      pts[i*10+j].rotaxis[1]=r();
      pts[i*10+j].rotaxis[2]=r();
      pts[i*10+j].moveable=true;

      rad=0.3;
      pts[i*10+j].rad=rad*2.0;
      pts[i*10+j].prd=rad;
      pts[i*10+j].mas=50.0f*rad*rad*rad;
    }
  }
}

void nbsolverobj::setmodel4()
{

/////////////////////////////////////////////////////////////////////////////
// This is a model of interacting bodies. Dissipative forces, i.e. loss 
// by friction, are absent. The interaction is implemented according 
// to the inverse-square law, and at small distances attraction is replaced 
// by repulsion.
// The configuration of the objects doesn't stabilize with time.
//

  len=100;
  quadfactor=1.0f;
  linlfactor=0.0f;
  vfactor=1.000;

  for (int i=0; i<len; i++)
  {
    mfloat rad=r01()*r01()+0.1;

    pts[i].clr[0]=1.0f-rad*0.2;
    pts[i].clr[1]=0.6f+rad*0.2;
    pts[i].clr[2]=0.3f+0.1*r();
    pts[i].moveable=true;

    pts[i].vel.x=0;
    pts[i].vel.y=0;
    pts[i].vel.z=0;

    pts[i].pos.x=r()*4.0f;
    pts[i].pos.y=r()*4.0f;
    pts[i].pos.z=r()*4.0f;

    pts[i].rotaxis[0]=r();
    pts[i].rotaxis[1]=r();
    pts[i].rotaxis[2]=r();

    rad=0.2;
    pts[i].rad=rad*2.0;
    pts[i].prd=rad;
    pts[i].mas=100.0f*rad*rad*rad;
  }

  for (i=0; i<4; i++)
  {
    for (int j=0; j<4; j++)
    {
      pts[i*4+j].clr[0]=0.3f;
      pts[i*4+j].clr[1]=1.0f;
      pts[i*4+j].clr[2]=0.6f;

       pts[i*4+j].pos.x=j-1.5;
       pts[i*4+j].pos.y=i-1.5;
       pts[i*4+j].pos.z=0;

      pts[i*4+j].rotaxis[0]=r();
      pts[i*4+j].rotaxis[1]=r();
      pts[i*4+j].rotaxis[2]=r();
      pts[i*4+j].moveable=true;

      mfloat rad=0.4;
      pts[i*4+j].rad=rad*3.0;
      pts[i*4+j].prd=rad;
       pts[i*4+j].mas=100.0f*rad*rad*rad;
    }
  }
}

void nbsolverobj::setmodel5()
{
  float col;
  i32 i;

/////////////////////////////////////////////////////////////////////////////
// This is a model of interacting bodies. Loss by friction is absent. 
// The interaction is implemented according to the inverse-square law, 
// and at small distances attraction is replaced by repulsion. 
// The configuration is a star with a belt of asteroids.
// The configuration of the objects doesn't stabilize with time. 
// Because of unstable configuration almost in any starting conditions 
// some asteroids fly away to infinity due to interaction. Such behavior 
// is typical of a system of many bodies which interact according to 
// the inverse-square law without friction (this is a physically accurate 
// law of gravitational interaction). But stability of a configuration 
// with one central body is quite high and increases as the mass of 
// the central body grows.
//

  len=100;
  quadfactor=1.0f;
  linlfactor=0.0f;
  vfactor=1.0;

  for (i=0; i<len; i++)
  {
    mfloat q=r()*3.14156;

    col=0.65+0.35*r01();
    pts[i].clr[0]=col;
    pts[i].clr[1]=col;
    pts[i].clr[2]=col;
    pts[i].moveable=true;

    pts[i].pos.x=0.4*r()+3.0f*sin(q);
    pts[i].pos.y=0.4*r()+3.0f*cos(q);
    pts[i].pos.z=0.4*r()+3.0f*cos(q)*sin(q);

    pts[i].vel.x=0.1*r()-8.0f*cos(q);
    pts[i].vel.y=0.1*r()+8.0f*sin(q);
    pts[i].vel.z=0.0;

    pts[i].rotaxis[0]=r();
    pts[i].rotaxis[1]=r();
    pts[i].rotaxis[2]=r();

    mfloat rad=0.3;
    pts[i].rad=rad*1.0;
    pts[i].prd=rad;
    pts[i].mas=50.0f*rad*rad*rad;
  }

  for (i=0; i<1; i++)
  {
    pts[i].clr[0]=1;
    pts[i].clr[1]=1;
    pts[i].clr[2]=0;

    pts[i].vel.x=0;
      pts[i].vel.y=0;
      pts[i].vel.z=0;

    pts[i].pos.x=0;
      pts[i].pos.y=0;
      pts[i].pos.z=0;

    float rad=1.0;
    pts[i].rad=rad*1.0;
    pts[i].prd=rad;
      pts[i].mas=300.0f*rad*rad*rad;
  }
}

void nbsolverobj::setmodel6()
{
  int i;

/////////////////////////////////////////////////////////////////////////////
// This is a model of interacting bodies. Loss by friction is absent. 
// The interaction is implemented according to the inverse-radius law, 
// and at small distances attraction is replaced by repulsion.
// The configuration is three stars with a common belt of asteroids. 
// The configuration of the objects doesn't stabilize with time. 
// The configuration is not stable, however no objects fly away to infinity. 
// Such behavior is typical and connected with the fact that 
// the corresponding potential increases without limit at infinity. 
// This model is an example of the system which is located in an infinitely 
// deep potential well.
//

  len=100;
  quadfactor=0.0f;
  linlfactor=1.0f;
  vfactor=1.000000;


  for (i=0; i<len; i++)
  {
    mfloat q=r()*3.14156;
    mfloat rad;

    pts[i].clr[0]=1;
    pts[i].clr[1]=1;
    pts[i].clr[2]=1;
    pts[i].moveable=true;

    pts[i].pos.x=0.4*r()+4.0f*sin(q);
    pts[i].pos.y=0.4*r()+4.0f*cos(q);
    pts[i].pos.z=0.4*r();

    pts[i].vel.x=0.1*r()-8.0f*cos(q);
    pts[i].vel.y=0.1*r()+8.0f*sin(q);
    pts[i].vel.z=0.0;

    pts[i].rotaxis[0]=r();
    pts[i].rotaxis[1]=r();
    pts[i].rotaxis[2]=r();

    rad=0.4;
    pts[i].rad=rad*2.0;
    pts[i].prd=rad;
    pts[i].mas=20.0f*rad*rad*rad;
  }

  for (i=0; i<3; i++)
  {
    pts[i].vel.x=0;
    pts[i].vel.y=0;
    pts[i].vel.z=0;

    pts[i].pos.x=r()*7.0;
    pts[i].pos.y=r()*7.0;
    pts[i].pos.z=r()*7.0;

    mfloat rad=1.0;
    pts[i].rad=rad*2.0;
    pts[i].prd=rad;
    pts[i].mas=30.0f*rad*rad*rad;
  }

  pts[0].clr[0]=1;
  pts[0].clr[1]=0;
  pts[0].clr[2]=0;

  pts[1].clr[0]=0;
  pts[1].clr[1]=1;
  pts[1].clr[2]=0;

  pts[2].clr[0]=0;
  pts[2].clr[1]=0;
  pts[2].clr[2]=1;
}

void nbsolverobj::setmodel7()
{
  int i;

/////////////////////////////////////////////////////////////////////////////
// This is a model of interacting bodies. Loss by friction is accounted, 
// i.e. the bodies move in a viscous medium. The interaction is implemented 
// according to the inverse-radius law, and at small distances attraction 
// is replaced by repulsion.
// The configuration of the objects becomes spherical, - such disposition 
// corresponds to the minimum of the potential energy. Such behavior can 
// be considered a demonstration of surface tension. And such model can 
// describe a drop of water.
//

  len=100;
  quadfactor=0.0f;
  linlfactor=1.0f;
  vfactor=0.99;

  for (i=0; i<len; i++)
  {
    pts[i].clr[0]=0.2f+0.1*r();
    pts[i].clr[1]=0.5f+0.1*r();
    pts[i].clr[2]=0.9f+0.1*r();

    pts[i].pos.x=r()*1.0f;
    pts[i].pos.y=r()*1.0f;
    pts[i].pos.z=r()*1.0f;

    pts[i].vel.x=0;
    pts[i].vel.y=0;
    pts[i].vel.z=0;

    pts[i].rotaxis[0]=r();
    pts[i].rotaxis[1]=r();
    pts[i].rotaxis[2]=r();
    pts[i].moveable=true;

    float rad=1.2;
    pts[i].rad=rad*1.2;
    pts[i].prd=rad;
    pts[i].mas=20.0f*rad*rad*rad;
  }

  for (i=0; i<10; i++)
  {
    pts[i].clr[0]=0.2f+0.1*r();
    pts[i].clr[1]=0.9f+0.1*r();
    pts[i].clr[2]=0.9f+0.1*r();

    pts[i].pos.x=r()*1.0f;
    pts[i].pos.y=r()*1.0f;
    pts[i].pos.z=r()*1.0f;

    pts[i].vel.x=0;
    pts[i].vel.y=0;
    pts[i].vel.z=0;

    pts[i].rotaxis[0]=r();
    pts[i].rotaxis[1]=r();
    pts[i].rotaxis[2]=r();
    pts[i].moveable=true;

    mfloat rad=1.4;
    pts[i].rad=rad*1.2;
    pts[i].prd=rad;
    pts[i].mas=20.0f*rad*rad*rad;
  }
}



void nbsolverobj::initlight1(sceneguidt sceneguid, float time)
{
  float x,y,z;
  float l;

  
  l=20000.0;
  
  lightguid1=sceneallocatelight(sceneguid);
  setlightcenter3f(lightguid1,-15*l,0,-15*l);
  setlightdistance(lightguid1,10000000);
  setlightcoef(lightguid1,1,0.0000001,0);
  setlightintense(lightguid1,0.85);
  setlightcolor3f(lightguid1,1,1,1);
  enablelight(lightguid1,true);


  lightguid2=sceneallocatelight(sceneguid);
  setlightcenter3f(lightguid2,0,15*l,0);
  setlightdistance(lightguid2,10000000);
  setlightcoef(lightguid2,1,0.0000001,0);
  setlightintense(lightguid2,0.85);
  setlightcolor3f(lightguid2,1,0.5,1);
  enablelight(lightguid2,true);


  lightguid3=sceneallocatelight(sceneguid);
  setlightcenter3f(lightguid3,0,-15*l,0);
  setlightdistance(lightguid3,10000000);
  setlightcoef(lightguid3,1,0.0000001,0);
  setlightintense(lightguid3,0.6);
  setlightcolor3f(lightguid3,1,0.5,1);
  enablelight(lightguid3,true);

  /*
  lightguid1=sceneallocatelight(sceneguid);
  setlightcenter3f(lightguid1,15000.0*4.0f,
                              15000.0*0.0f,
                              15000.0*0.0f);
  setlightdistance(lightguid1,1000000);
  setlightcoef(lightguid1,1,0.00001,0);
  setlightintense(lightguid1,1);
  setlightcolor3f(lightguid1,0.2,1,0.2);
  enablelight(lightguid1,true);


  lightguid2=sceneallocatelight(sceneguid);
  setlightcenter3f(lightguid2,-15000.0*4.0f,
                              15000.0*0.0f,
                              15000.0*0.0f);
  setlightdistance(lightguid2,1000000);
  setlightcoef(lightguid2,1,0.00005,0);
  setlightintense(lightguid2,1);
  setlightcolor3f(lightguid2,1,0.2,0.2);
  enablelight(lightguid2,true);
  */
}

void nbsolverobj::initlight2(sceneguidt sceneguid, float time)
{
  float x1,y1,z1;
  float x2,y2,z2;
  float l;

  l=15000.0;


  x1=l*25*cos(time*0.1);
  y1=0;
  z1=l*25*sin(time*0.1);

  lightguid1=sceneallocatelight(sceneguid);
  setlightcenter3f(lightguid1,x1,y1,z1);
  setlightdistance(lightguid1,1000000);
  setlightcoef(lightguid1,1,0.0000001,0);
  setlightintense(lightguid1,0.85);
  setlightcolor3f(lightguid1,0,1,0.2);
  enablelight(lightguid1,true);


  x2=l*25*cos(time*0.13+PI);
  y2=0;
  z2=l*25*sin(time*0.13+PI);

  lightguid2=sceneallocatelight(sceneguid);
  setlightcenter3f(lightguid2,x2,y2,z2);
  setlightdistance(lightguid2,1000000);
  setlightcoef(lightguid2,1,0.0000001,0);
  setlightintense(lightguid2,0.85);
  setlightcolor3f(lightguid2,0.2,0,1);
  enablelight(lightguid2,true);
}

void nbsolverobj::initlight3(sceneguidt sceneguid, float time)
{
  lightguid1=sceneallocatelight(sceneguid);
  setlightcenter3f(lightguid1,75000.0*4.0f,
                              15000.0*0.0f,
                              15000.0*0.0f);
  setlightdistance(lightguid1,1000000);
  setlightcoef(lightguid1,1,0.00001,0);
  setlightintense(lightguid1,3.5);
  setlightcolor3f(lightguid1,0,1,0.2);
  enablelight(lightguid1,true);


  lightguid2=sceneallocatelight(sceneguid);
  setlightcenter3f(lightguid2,-75000.0*4.0f,
                              15000.0*0.0f,
                              15000.0*0.0f);
  setlightdistance(lightguid2,1000000);
  setlightcoef(lightguid2,1,0.00001,0);
  setlightintense(lightguid2,3.5);
  setlightcolor3f(lightguid2,0,0,1);
  enablelight(lightguid2,true);


  lightguid3=sceneallocatelight(sceneguid);
  setlightcenter3f(lightguid3,75000.0*0.0f,
                              15000.0*0.0f,
                              -75000.0*4.0f);
  setlightdistance(lightguid3,1000000);
  setlightcoef(lightguid3,1,0.00001,0);
  setlightintense(lightguid3,1.0);
  setlightcolor3f(lightguid3,1,0,1);
  enablelight(lightguid3,true);
}

void nbsolverobj::initlight4(sceneguidt sceneguid, float time)
{
  lightguid1=sceneallocatelight(sceneguid);
  setlightcenter3f(lightguid1,15000.0*4.0f,
                              15000.0*0.0f,
                              15000.0*0.0f);
  setlightdistance(lightguid1,1000000);
  setlightcoef(lightguid1,1,0.00001,0);
  setlightintense(lightguid1,1);
  setlightcolor3f(lightguid1,0.2,1,0.2);
  enablelight(lightguid1,true);


  lightguid2=sceneallocatelight(sceneguid);
  setlightcenter3f(lightguid2,-15000.0*4.0f,
                              15000.0*0.0f,
                              15000.0*0.0f);
  setlightdistance(lightguid2,1000000);
  setlightcoef(lightguid2,1,0.00005,0);
  setlightintense(lightguid2,1);
  setlightcolor3f(lightguid2,1,0.2,0.2);
  enablelight(lightguid2,true);
}

void nbsolverobj::initlight5(sceneguidt sceneguid, float time)
{
//  setscenebgillum3b(sceneguid,60,60,60);
  
  lightguid1=sceneallocatelight(sceneguid);
  setlightcenter3f(lightguid1,15000.0*pts[0].pos.x,
                              15000.0*pts[0].pos.y,
                              15000.0*pts[0].pos.z);
  setlightdistance(lightguid1,1000000);
  setlightcoef(lightguid1,0.0001,0,0.000000001);
  setlightintense(lightguid1,3);
  setlightcolorfp(lightguid1,pts[0].clr);
  enablelight(lightguid1,true);


  lightguid2=sceneallocatelight(sceneguid);
  setlightcenter3f(lightguid2,-15000.0*10.0f,
                              15000.0*0.0f,
                              -15000.0*20.0f);
  setlightdistance(lightguid2,1000000);
  setlightcoef(lightguid2,1,0.00001,0);
  setlightintense(lightguid2,1);
  setlightcolor3f(lightguid2,0,0,1);
  enablelight(lightguid2,true);
}

void nbsolverobj::initlight6(sceneguidt sceneguid, float time)
{
  lightguid1=sceneallocatelight(sceneguid);
  setlightcenter3f(lightguid1,15000.0*pts[0].pos.x,
                              15000.0*pts[0].pos.y,
                              15000.0*pts[0].pos.z);
  setlightdistance(lightguid1,1000000);
  setlightcoef(lightguid1,0.0001,0.000035,0);
  setlightintense(lightguid1,2.25);
  setlightcolorfp(lightguid1,pts[0].clr);
  enablelight(lightguid1,true);


  lightguid2=sceneallocatelight(sceneguid);
  setlightcenter3f(lightguid2,15000.0*pts[1].pos.x,
                              15000.0*pts[1].pos.y,
                              15000.0*pts[1].pos.z);
  setlightdistance(lightguid2,1000000);
  setlightcoef(lightguid2,0.0001,0.000035,0);
  setlightintense(lightguid2,2.25);
  setlightcolorfp(lightguid2,pts[1].clr);
  enablelight(lightguid2,true);


  lightguid3=sceneallocatelight(sceneguid);
  setlightcenter3f(lightguid3, 15000.0*pts[2].pos.x,
                               15000.0*pts[2].pos.y,
                               15000.0*pts[2].pos.z);
  setlightdistance(lightguid3,1000000);
  setlightcoef(lightguid3,0.0001,0.000035,0);
  setlightintense(lightguid3,2.25);
  setlightcolorfp(lightguid3,pts[2].clr);
  enablelight(lightguid3,true);
}

void nbsolverobj::initlight7(sceneguidt sceneguid, float time)
{
  lightguid1=sceneallocatelight(sceneguid);
  setlightcenter3f(lightguid1,-15000.0*10.0f,
                              15000.0*0.0f,
                              -15000.0*20.0f);
  setlightdistance(lightguid1,1000000);
  setlightcoef(lightguid1,1,0.000001,0);
  setlightintense(lightguid1,1);
  setlightcolor3f(lightguid1,0,1,1);
  enablelight(lightguid1,true);
}

void nbsolverobj::initlight(int modeltype, sceneguidt sceneguid, float time)
{
//  gmodeltype=modeltype;
  switch (modeltype)
  {
    case 1:
    {
      initlight1(sceneguid,time);
      break;
    }
    case 2:
    {
      initlight2(sceneguid,time);
      break;
    }
    case 3:
    {
      initlight3(sceneguid,time);
      break;
    }
    case 4:
    {
      initlight4(sceneguid,time);
      break;
    }
    case 5:
    {
      initlight5(sceneguid,time);
      break;
    }
    case 6:
    {
      initlight6(sceneguid,time);
      break;
    }
    case 7:
    {
      initlight7(sceneguid,time);
      break;
    }
    default:
    {
      initlight1(sceneguid,time);
      break;
    }
  }
}

void nbsolverobj::init(int modeltype, int randseed)
{
  gmodeltype=modeltype;
  rndgenerator.init(randseed);

  potenergy=0;
  energyleak=0;

  switch (modeltype)
  {
    case 1:
    {
      setmodel1();
      break;
    }
    case 2:
    {
      setmodel2();
      break;
    }
    case 3:
    {
      setmodel3();
      break;
    }
    case 4:
    {
      setmodel4();
      break;
    }
    case 5:
    {
      setmodel5();
      break;
    }
    case 6:
    {
      setmodel6();
      break;
    }
    case 7:
    {
      setmodel7();
      break;
    }
    default:
    {
      setmodel1();
      break;
    }
  }
}
