--- /dev/null
+/*
+ */
+#include <common.h>
+#include <stdio.h>
+#include <video.h> // For the thread functions
+
+// === CONSTANTS ===
+#define G_CONSTANT 200 // Gravitational Constant
+#define E_CONSTANT 400 // Electric Constant
+#define N_CONSTANT 50 // Nuclear Constant
+#define WN_CONSTANT 100 // Weak Nuclear Constant
+
+#define MAX_THREADS 32
+
+#define SHOW_FORCES 1
+
+#define DEBUG 0
+#define DEBUG_SCAN (DEBUG&&)
+#define DEBUG_CHECKS (DEBUG&&0)
+#define DEBUG_INNER (DEBUG&&0)
+#define DEBUG_STATE (DEBUG&&0)
+
+// === PROTOTYPES ===
+void Update_DoTimestep(void);
+void Update_ParticleForces(tParticle *Part);
+void Update_ParticleVelocity(tParticle *Part);
+void Update_ParticlePosition(tParticle *Part);
+
+// === GLOBALS ===
+ int giNumberOfJobs = 16;
+volatile int gaThreadStatuses[MAX_THREADS]; // -1 = Die, 0 = Completed, 1 = Working
+
+// === CODE ===
+int Update_ProcessThread(void *Offset)
+{
+ int ofs = (intptr_t)Offset; // Also Thread ID
+ int i;
+
+ while( gaThreadStatuses[ofs] != -1 )
+ {
+ while( gaThreadStatuses[ofs] == 0 )
+ Thread_Yield();
+ if( gaThreadStatuses[ofs] == -1 )
+ break;
+ for( i = ofs; i < giNumParticles; i += giNumberOfJobs )
+ {
+ Update_ParticleForces( &gpParticles[i] );
+ }
+ for( i = ofs; i < giNumParticles; i += giNumberOfJobs )
+ {
+ Update_ParticleVelocity( &gpParticles[i] );
+ }
+ gaThreadStatuses[ofs] = 0;
+ }
+ return 0;
+}
+
+/**
+ * Start update threads
+ */
+void Update_Initialise(void)
+{
+ int i;
+ for( i = 1; i < giNumberOfJobs; i++ ) {
+ Thread_Create(Update_ProcessThread, (void*)(intptr_t)i);
+ }
+}
+
+void Update_Cleanup(void)
+{
+ int i;
+ for( i = 1; i < giNumberOfJobs; i++ ) {
+ gaThreadStatuses[i] = -1;
+ }
+}
+
+void Update_DoTimestep(void)
+{
+ int i;
+
+ // Start jobs again
+ for( i = 1; i < giNumberOfJobs; i++ ) {
+ gaThreadStatuses[i] = 1; // You have work to do
+ }
+
+ for( i = 0; i < giNumParticles; i += giNumberOfJobs )
+ {
+ #if DEBUG | SHOW_FORCES
+ printf(" Forces for #%i \r", i);
+ fflush(stdout);
+ #endif
+ Update_ParticleForces( &gpParticles[i] );
+ }
+ #if SHOW_FORCES
+ printf(" \r");
+ #endif
+
+ #if DEBUG
+ printf(" Forces Done \n");
+ printf(" Velocities\n");
+ #endif
+ for( i = 0; i < giNumParticles; i += giNumberOfJobs )
+ {
+ #if DEBUG
+ printf(" Accel & Velocity for %i %c", i, DEBUG_SCAN?'\r':'\n');
+ #endif
+ Update_ParticleVelocity( &gpParticles[i] );
+ }
+
+ // Wait for jobs to terminate
+ for( i = 1; i < giNumberOfJobs; i++ ) {
+ while( gaThreadStatuses[i] == 1 )
+ Thread_Yield();
+ }
+
+ #if DEBUG
+ printf(" Positions \n");
+ #endif
+ for( i = 0; i < giNumParticles; i++ )
+ {
+ Update_ParticlePosition( &gpParticles[i] );
+ }
+}
+
+void Update_ParticleForces(tParticle *Part)
+{
+ int i;
+ tLargeInt dist, dist2, dist3, dist4;
+ tLargeInt force, tmpQuad;
+ tVector line;
+ tLargeInt *partMass = &Part->Mass;
+
+ Vector_Zero(&Part->Forces);
+
+ for( i = 0; i < giNumParticles; i++ )
+ {
+ if( &gpParticles[i] == Part ) continue;
+ #if DEBUG_CHECKS
+ printf(" Checking against %4i ", i);
+ # if DEBUG_INNER
+ printf("\n");
+ # else
+ printf("\r");
+ # endif
+ #endif
+
+ Vector_Set( &line, &gpParticles[i].Location );
+ Vector_Sub( &line, &Part->Location );
+ Vector_Magnitude( &dist, &line );
+ LargeInt_PowNative( &dist2, &dist, 2 );
+ LargeInt_PowNative( &dist3, &dist, 3 );
+ LargeInt_PowNative( &dist4, &dist, 4 );
+
+ #if DEBUG_INNER
+ printf(" Calculating Gravity\n");
+ #endif
+ // Gravitational Forces
+ // -(mM)*G/(dist^2)
+ LargeInt_Mul( &tmpQuad, NULL, partMass, &gpParticles[i].Mass );
+ //printf(" tmpQuad = %s\n", LargeInt_DumpEx(&tmpQuad, 0) );
+ LargeInt_MulNative( &tmpQuad, NULL, &tmpQuad, G_CONSTANT );
+ LargeInt_Div( &force, NULL, &tmpQuad, &dist2 );
+ LargeInt_Neg( &force, &force );
+
+ #if DEBUG_INNER
+ printf(" Calculating Electric\n");
+ #endif
+ // Electric Forces
+ // - Opposites Attract, Like Repels, so no need to negate
+ // (e1e2)*E/(dist^2)
+ //printf(" Charges %s and %s\n",
+ // LargeInt_DumpEx(&Part->Charge, 0),
+ // LargeInt_DumpEx(&gpParticles[i].Charge, 1)
+ // );
+ LargeInt_Mul( &tmpQuad, NULL, &Part->Charge, &gpParticles[i].Charge );
+ //printf(" tmpQuad = %s\n", LargeInt_DumpEx(&tmpQuad, 0) );
+ LargeInt_MulNative( &tmpQuad, NULL, &tmpQuad, E_CONSTANT );
+ LargeInt_Div( &tmpQuad, NULL, &tmpQuad, &dist2 );
+ LargeInt_Add( &force, &force, &tmpQuad );
+
+ #if DEBUG_INNER
+ printf(" Calculating Nuclear\n");
+ #endif
+ // Nuclear Forces
+ // > Strong Nuclear (Attractive, 1/d^3)
+ // -(mM)*N/(dist^3)
+ LargeInt_Mul( &tmpQuad, NULL, &Part->Mass, &gpParticles[i].Mass );
+ LargeInt_MulNative( &tmpQuad, NULL, &tmpQuad, N_CONSTANT );
+ LargeInt_Div( &tmpQuad, NULL, &tmpQuad, &dist3 );
+ LargeInt_Neg( &tmpQuad, &tmpQuad );
+ LargeInt_Add( &force, &force, &tmpQuad );
+ // > Weak Nuclear (Repulsive, d^4)
+ // (mM)*N/(dist^4)
+ LargeInt_Mul( &tmpQuad, NULL, &Part->Mass, &gpParticles[i].Mass );
+ LargeInt_MulNative( &tmpQuad, NULL, &tmpQuad, WN_CONSTANT );
+ LargeInt_Div( &tmpQuad, NULL, &tmpQuad, &dist4 );
+ LargeInt_Add( &force, &force, &tmpQuad );
+
+ #if DEBUG_INNER
+ printf(" Calculating final force\n");
+ #endif
+ #if DEBUG_INNER >= 2
+ printf(" Mul (%s,%s) *= %s\n",
+ LargeInt_DumpEx(&line.C[0], 0),
+ LargeInt_DumpEx(&line.C[1], 1),
+ LargeInt_DumpEx(&force, 2)
+ );
+ #endif
+ Vector_MulScalar( &line, &force );
+
+ Vector_DivScalar( &line, &dist );
+ Vector_Add(&Part->Forces, &line);
+ }
+
+ #if 0
+ printf(" Part->Forces = (%s,%s)\n",
+ LargeInt_DumpEx(&Part->Forces.C[0], 0),
+ LargeInt_DumpEx(&Part->Forces.C[1], 1)
+ );
+ #endif
+
+ #if DEBUG_CHECKS
+ # if DEBUG_INNER
+ printf("\n");
+ # endif
+ #endif
+}
+
+void Update_ParticleVelocity(tParticle *Part)
+{
+ tVector accel;
+ #if DEBUG_STATE
+ printf(" Set\n");
+ #endif
+ Vector_Set( &accel, &Part->Forces );
+ #if DEBUG_STATE
+ printf(" Div (%s,%s) /= %s\n",
+ LargeInt_DumpEx(&accel.C[0], 0),
+ LargeInt_DumpEx(&accel.C[1], 1),
+ LargeInt_DumpEx(&Part->Mass, 2)
+ );
+ #endif
+ Vector_DivScalar( &accel, &Part->Mass );
+ #if DEBUG_STATE
+ printf(" = (%s,%s)\n",
+ LargeInt_DumpEx(&accel.C[0], 0),
+ LargeInt_DumpEx(&accel.C[1], 1)
+ );
+ printf(" Add\n");
+ #endif
+ Vector_Add( &Part->Velocity, &accel );
+}
+
+void Update_ParticlePosition(tParticle *Part)
+{
+ Vector_Add( &Part->Location, &Part->Velocity );
+}
git.ucc.asn.au / matches / honours.git / blobdiff
UCC git Repository :: git.ucc.asn.au