5 #include <video.h> // For the thread functions
8 #define G_CONSTANT 200 // Gravitational Constant
9 #define E_CONSTANT 400 // Electric Constant
10 #define N_CONSTANT 50 // Nuclear Constant
11 #define WN_CONSTANT 100 // Weak Nuclear Constant
13 #define MAX_THREADS 32
18 #define DEBUG_SCAN (DEBUG&&)
19 #define DEBUG_CHECKS (DEBUG&&0)
20 #define DEBUG_INNER (DEBUG&&0)
21 #define DEBUG_STATE (DEBUG&&0)
24 void Update_DoTimestep(void);
25 void Update_ParticleForces(tParticle *Part);
26 void Update_ParticleVelocity(tParticle *Part);
27 void Update_ParticlePosition(tParticle *Part);
30 int giNumberOfJobs = 16;
31 volatile int gaThreadStatuses[MAX_THREADS]; // -1 = Die, 0 = Completed, 1 = Working
34 int Update_ProcessThread(void *Offset)
36 int ofs = (intptr_t)Offset; // Also Thread ID
39 while( gaThreadStatuses[ofs] != -1 )
41 while( gaThreadStatuses[ofs] == 0 )
43 if( gaThreadStatuses[ofs] == -1 )
45 for( i = ofs; i < giNumParticles; i += giNumberOfJobs )
47 Update_ParticleForces( &gpParticles[i] );
49 for( i = ofs; i < giNumParticles; i += giNumberOfJobs )
51 Update_ParticleVelocity( &gpParticles[i] );
53 gaThreadStatuses[ofs] = 0;
59 * Start update threads
61 void Update_Initialise(void)
64 for( i = 1; i < giNumberOfJobs; i++ ) {
65 Thread_Create(Update_ProcessThread, (void*)(intptr_t)i);
69 void Update_Cleanup(void)
72 for( i = 1; i < giNumberOfJobs; i++ ) {
73 gaThreadStatuses[i] = -1;
77 void Update_DoTimestep(void)
82 for( i = 1; i < giNumberOfJobs; i++ ) {
83 gaThreadStatuses[i] = 1; // You have work to do
86 for( i = 0; i < giNumParticles; i += giNumberOfJobs )
88 #if DEBUG | SHOW_FORCES
89 printf(" Forces for #%i \r", i);
92 Update_ParticleForces( &gpParticles[i] );
99 printf(" Forces Done \n");
100 printf(" Velocities\n");
102 for( i = 0; i < giNumParticles; i += giNumberOfJobs )
105 printf(" Accel & Velocity for %i %c", i, DEBUG_SCAN?'\r':'\n');
107 Update_ParticleVelocity( &gpParticles[i] );
110 // Wait for jobs to terminate
111 for( i = 1; i < giNumberOfJobs; i++ ) {
112 while( gaThreadStatuses[i] == 1 )
117 printf(" Positions \n");
119 for( i = 0; i < giNumParticles; i++ )
121 Update_ParticlePosition( &gpParticles[i] );
125 void Update_ParticleForces(tParticle *Part)
128 tLargeInt dist, dist2, dist3, dist4;
129 tLargeInt force, tmpQuad;
131 tLargeInt *partMass = &Part->Mass;
133 Vector_Zero(&Part->Forces);
135 for( i = 0; i < giNumParticles; i++ )
137 if( &gpParticles[i] == Part ) continue;
139 printf(" Checking against %4i ", i);
147 Vector_Set( &line, &gpParticles[i].Location );
148 Vector_Sub( &line, &Part->Location );
149 Vector_Magnitude( &dist, &line );
150 LargeInt_PowNative( &dist2, &dist, 2 );
151 LargeInt_PowNative( &dist3, &dist, 3 );
152 LargeInt_PowNative( &dist4, &dist, 4 );
155 printf(" Calculating Gravity\n");
157 // Gravitational Forces
159 LargeInt_Mul( &tmpQuad, NULL, partMass, &gpParticles[i].Mass );
160 //printf(" tmpQuad = %s\n", LargeInt_DumpEx(&tmpQuad, 0) );
161 LargeInt_MulNative( &tmpQuad, NULL, &tmpQuad, G_CONSTANT );
162 LargeInt_Div( &force, NULL, &tmpQuad, &dist2 );
163 LargeInt_Neg( &force, &force );
166 printf(" Calculating Electric\n");
169 // - Opposites Attract, Like Repels, so no need to negate
171 //printf(" Charges %s and %s\n",
172 // LargeInt_DumpEx(&Part->Charge, 0),
173 // LargeInt_DumpEx(&gpParticles[i].Charge, 1)
175 LargeInt_Mul( &tmpQuad, NULL, &Part->Charge, &gpParticles[i].Charge );
176 //printf(" tmpQuad = %s\n", LargeInt_DumpEx(&tmpQuad, 0) );
177 LargeInt_MulNative( &tmpQuad, NULL, &tmpQuad, E_CONSTANT );
178 LargeInt_Div( &tmpQuad, NULL, &tmpQuad, &dist2 );
179 LargeInt_Add( &force, &force, &tmpQuad );
182 printf(" Calculating Nuclear\n");
185 // > Strong Nuclear (Attractive, 1/d^3)
187 LargeInt_Mul( &tmpQuad, NULL, &Part->Mass, &gpParticles[i].Mass );
188 LargeInt_MulNative( &tmpQuad, NULL, &tmpQuad, N_CONSTANT );
189 LargeInt_Div( &tmpQuad, NULL, &tmpQuad, &dist3 );
190 LargeInt_Neg( &tmpQuad, &tmpQuad );
191 LargeInt_Add( &force, &force, &tmpQuad );
192 // > Weak Nuclear (Repulsive, d^4)
194 LargeInt_Mul( &tmpQuad, NULL, &Part->Mass, &gpParticles[i].Mass );
195 LargeInt_MulNative( &tmpQuad, NULL, &tmpQuad, WN_CONSTANT );
196 LargeInt_Div( &tmpQuad, NULL, &tmpQuad, &dist4 );
197 LargeInt_Add( &force, &force, &tmpQuad );
200 printf(" Calculating final force\n");
203 printf(" Mul (%s,%s) *= %s\n",
204 LargeInt_DumpEx(&line.C[0], 0),
205 LargeInt_DumpEx(&line.C[1], 1),
206 LargeInt_DumpEx(&force, 2)
209 Vector_MulScalar( &line, &force );
211 Vector_DivScalar( &line, &dist );
212 Vector_Add(&Part->Forces, &line);
216 printf(" Part->Forces = (%s,%s)\n",
217 LargeInt_DumpEx(&Part->Forces.C[0], 0),
218 LargeInt_DumpEx(&Part->Forces.C[1], 1)
229 void Update_ParticleVelocity(tParticle *Part)
235 Vector_Set( &accel, &Part->Forces );
237 printf(" Div (%s,%s) /= %s\n",
238 LargeInt_DumpEx(&accel.C[0], 0),
239 LargeInt_DumpEx(&accel.C[1], 1),
240 LargeInt_DumpEx(&Part->Mass, 2)
243 Vector_DivScalar( &accel, &Part->Mass );
245 printf(" = (%s,%s)\n",
246 LargeInt_DumpEx(&accel.C[0], 0),
247 LargeInt_DumpEx(&accel.C[1], 1)
251 Vector_Add( &Part->Velocity, &accel );
254 void Update_ParticlePosition(tParticle *Part)
256 Vector_Add( &Part->Location, &Part->Velocity );