// Subroutine for particle-based simulations
// Au: Cagatay Basdogan, basdogan@mit.edu
// note that there are references to Open Inventor constructs in this routine

void ObjectStruct::calculate_force_applied_to_vertex(SbVec3f probe_tip)
{
        //called every time thru the loop - moves all vertices in the directions of forces applied on them

         int j, p,vertex_no, line_no;
         int vertex_num;
         SbVec3f  temp_vec;
         SbVec3f force_on_v0,force_on_v1,force_on_v2,force_on_v3,force_on_v;
         double mag_vec,diff_length;
         PolyhedronStruct *bone = lores;

         temp_vec.setValue(0.0,0.0,0.0);
 
        for (p=0; p<bone->GetNumberOfVertex(); p++)
         {
              if (bone->vertex[p].frozen = = no)
              {
                    force_on_v.setValue(0.0,0.0,0.0);
                    force_on_v0.setValue(0.0,0.0,0.0);
                    force_on_v1.setValue(0.0,0.0,0.0);
                    force_on_v2.setValue(0.0,0.0,0.0);
                    force_on_v3.setValue(0.0,0.0,0.0);
 

                    // FORCES THAT BRING THE VERTICES BACK (use for stability)
                    force_on_v0 = -1.0*bDrag*(bone->vertex[p].current_vel);
 
                    // FORCES BETWEEN THE CURRENT POS. OF VERTEX AND ITS NEIGHBORS
                    for (j=0; j<bone->vertex[p].GetNumberOfVertex(); j++)
                   {
                        // get the vertex and line indices of neighbors
                        vertex_no = bone->vertex[p].neighbor_vertex_index[j];
                        temp_vec = bone->vertex[vertex_no].new_coord - bone->vertex[p].new_coord;
                        mag_vec = temp_vec.normalize();
                        line_no = bone->vertex[p].neighbor_line_index[j];
                        diff_length = mag_vec - bone->line[line_no].length;

                        // spring force
                        force_on_v1 = force_on_v1 + (kNeighbor*diff_length*temp_vec);
                        // damping force
                        temp_vec = (bone->vertex[temp].current_vel) - (bone->vertex[p].current_vel);
                        force_on_v2 = force_on_v2 + (bNeighbor*temp_vec);
                   }

                  // FORCES DUE TO GRAVITY
                  force_on_v3.setValue(0.0,-(MASS_OF_PARTICLE * 9.8), 0.0);
 
                // THE TOTAL FORCE APPLIED ON THE VERTEX
                force_on_v = force_on_v0 + force_on_v1 + force_on_v2 + force_on_v3;
                bone->vertex[p].force_on_vertex = force_on_v;
              }

         }

        for (p=0; p<bone->GetNumberOfVertex(); p++)
        {
            //UPDATE ACC,VEL,POS (state of the particle)
            bone->vertex[p].coord = bone->vertex[p].new_coord;
            bone->vertex[p].previous_vel = bone->vertex[p].current_vel;
            bone->vertex[p].previous_acc = bone->vertex[p].current_acc;

           //EULER INTEGRATION
           bone->vertex[p].current_acc = (bone->vertex[p].force_on_vertex) / MASS_OF_PARTICLE;
           bone->vertex[p].current_vel = bone->vertex[p].previous_vel + DeltaTIME * bone->vertex[p].current_acc;
           bone->vertex[p].new_coord = bone->vertex[p].coord + DeltaTIME * bone->vertex[p].current_vel;
         }
};