class psystem{
  particle[] particles;
  int nump;

  psystem(){
  }

  psystem(int n, float r){
    float x, y, z,yl,zl;
    nump=n;
    particles = new particle[nump+1];
    for (int i=0;i<nump;i++){
      x=random(-r,r);
      yl=sqrt(r*r-x*x);
      y=random(-yl,yl);
      zl=sqrt(r*r-x*x-y*y);
      z=(random(1.0)<0.5)?-zl:zl;

      particles[i]=new particle(x,y,z);
    }
    particles[nump]=new particle(0.0,0.0,0.0);
    particles[nump].vel.assign(0.0,0.0,0.0);

  }

  void update(float dt){
    vector3D temp = new vector3D(0.0,0.0,0.0);
    for (int i=0;i<nump;i++){
      particles[i].force.assign(vector3D.mul(-0.0012,particles[i].pos));
      for(int j=i+1;j<nump;j++){
        temp.assign(vector3D.div(vector3D.unit(vector3D.sub(particles[i].pos,particles[j].pos)),0.001*vector3D.dist2(particles[i].pos,particles[j].pos)));
        particles[i].force.add(temp);
        particles[j].force.sub(temp);
      }
    }
    for (int i=0;i<nump;i++){
      particles[i].update(dt);
    }
  }

  void draw3D(vector3D offset, int relativeTo, float s, float ay, float ax){
    int ii=0;
    float dx,dy,dz,dr;
    int hu,st,br,al;
    vector3D relPos=new vector3D(0,0,0);
    vector3D relVel=new vector3D(0,0,0);
    pushMatrix();
    translate(-offset.x,-offset.y,-offset.z);
    rotateY(ay);
    rotateX(ax);
    for (int i=0;i<nump;i++){        
      if(rot){
        beginShape();
        for(int index=1;index<particles[i].bl;index++){
          ii=(index+particles[i].cb)%particles[i].bl;
          relPos.assign(vector3D.sub(particles[i].pbuffer[ii],particles[relativeTo].pbuffer[ii]));
          relVel.assign(vector3D.sub(particles[i].vbuffer[ii],particles[relativeTo].vbuffer[ii]));
          br=2*index;
          hu=-15+int(50*relVel.sqlength()/400.0);
          if (hu<0) hu+=256;
          st=250;
          al=constrain(100-int(100*relPos.length()/750.0),0,100);
          stroke(hu,st, br,al*(1.0-(float)index/512));
          dx=s*relPos.x;
          dy=s*relPos.y;
          dz=s*relPos.z;
          if (screenZ(dx,dy,dz)<500) vertex(dx,dy,dz);
        }
        endShape();
        
        noStroke();
        relPos.assign(vector3D.sub(particles[i].pbuffer[0],particles[relativeTo].pbuffer[0]));
        relVel.assign(vector3D.sub(particles[i].vbuffer[0],particles[relativeTo].vbuffer[0]));
        dx=s*relPos.x;
        dy=s*relPos.y;
        dz=s*relPos.z;
        br=20;
        hu=-15+int(50*relVel.sqlength()/400.0);
        if (hu<0) hu+=256;
        st=250;
        al=constrain(150-int(150*relPos.length()/750.0),0,255);
        fill(hu,st,br, al);
        if (screenZ(dx,dy,dz)<500){
          pushMatrix();
          translate(dx,dy,dz);
          rotateX(-ax);
          rotateY(-ay);
          ellipse(0,0,5.0,5.0);
          popMatrix();
        }
      }

      if(!rot){
        noStroke();
        beginShape();
        for(int index=1;index<particles[i].bl;index++){
          ii=(index+particles[i].cb)%particles[i].bl;
          relPos.assign(vector3D.sub(particles[i].pbuffer[ii],particles[relativeTo].pbuffer[ii]));
          relVel.assign(vector3D.sub(particles[i].vbuffer[ii],particles[relativeTo].vbuffer[ii]));
          br=2*index;
          hu=-15+int(50*relVel.sqlength()/400.0);
          if (hu<0) hu+=256;
          st=250;
          al=constrain(50-int(50*relPos.length()/750.0),0,50);
          fill(hu,st, br,al*(1.0-(float)index/particles[i].bl));
          dx=s*relPos.x;
          dy=s*relPos.y;
          dz=s*relPos.z;
          dr=0.1-index*0.1/particles[i].bl;
          if (screenZ(dx,dy,dz)<500)vertex((1.0-dr)*dx,(1.0-dr)*dy,(1.0-dr)*dz);
          if (screenZ(dx,dy,dz)<500)vertex((0.8+dr)*dx,(0.8+dr)*dy,(0.8+dr)*dz);
        }
        endShape();
        beginShape();
        for(int index=1;index<particles[i].bl;index++){
          ii=(index+particles[i].cb)%particles[i].bl;
          relPos.assign(vector3D.sub(particles[i].pbuffer[ii],particles[relativeTo].pbuffer[ii]));
          relVel.assign(vector3D.sub(particles[i].vbuffer[ii],particles[relativeTo].vbuffer[ii]));
          br=2*index;
          hu=-15+int(50*relVel.sqlength()/400.0);
          if (hu<0) hu+=256;
          st=250;
          al=constrain(100-int(100*relPos.length()/750.0),0,100);
          stroke(hu,st, br,al*(1.0-(float)index/512));
          dx=s*relPos.x;
          dy=s*relPos.y;
          dz=s*relPos.z;
          dr=0.1-index*0.1/particles[i].bl;
          if (screenZ(dx,dy,dz)<500)vertex((1.0-dr)*dx,(1.0-dr)*dy,(1.0-dr)*dz);
        }
        endShape();
        beginShape();
        for(int index=1;index<particles[i].bl;index++){
          ii=(index+particles[i].cb)%particles[i].bl;
          relPos.assign(vector3D.sub(particles[i].pbuffer[ii],particles[relativeTo].pbuffer[ii]));
          relVel.assign(vector3D.sub(particles[i].vbuffer[ii],particles[relativeTo].vbuffer[ii]));
          br=2*index;
          hu=-15+int(50*relVel.sqlength()/400.0);
          if (hu<0) hu+=256;
          st=250;
          al=constrain(100-int(100*relPos.length()/750.0),0,100);
          stroke(hu,st, br,al*(1.0-(float)index/512));
          dx=s*relPos.x;
          dy=s*relPos.y;
          dz=s*relPos.z;
          dr=0.1-index*0.1/particles[i].bl;
          if (screenZ(dx,dy,dz)<500)vertex((0.8+dr)*dx,(0.8+dr)*dy,(0.8+dr)*dz);
        }
        endShape();
        drawing=false;
      }    
    }
     popMatrix();
  }
}