16
|
1
|
|
2 // written in the D programming language
|
|
3
|
|
4 module samples.MagnetsElectric;
|
|
5
|
|
6 import chipmunkd.chipmunk;
|
|
7
|
|
8 import samples.ChipmunkDemo;
|
|
9
|
|
10 import std.math;
|
|
11 import core.stdc.string:strcmp;
|
|
12 import core.stdc.stdio:sprintf;
|
|
13
|
|
14 enum WIDTH = 600;
|
|
15 enum HEIGHT = 400;
|
|
16
|
|
17 enum SINGMAX = 10; // Maximum number of singularities per body
|
|
18 enum NMAG = 10; // Number of magnets
|
|
19 enum NCHG = 10; // Number of charged bodies
|
|
20 enum NMIX = 10; // Number of charged magnets
|
|
21
|
|
22 enum COU_MKS = 8.987551787e9; // Some physical constants
|
|
23 enum MAG_MKS = 1e-7;
|
|
24
|
|
25 // Prototypes
|
|
26 alias void function(DataforForce* data)SingForceFunc;
|
|
27
|
|
28 // Structures
|
|
29 // Singularities
|
|
30 struct ActorSingularity{
|
|
31 // Number of singularities
|
|
32 int Nsing;
|
|
33 // Value of the singularities
|
|
34 cpFloat value[SINGMAX];
|
|
35 // Type of the singularities
|
|
36 char type[SINGMAX][100];
|
|
37 // Global position of the singularities
|
|
38 cpVect Gpos[SINGMAX];
|
|
39 // Local position of the singularities
|
|
40 cpVect position[SINGMAX];
|
|
41 // Angle of the singularities measured in the body axes
|
|
42 cpFloat angle[SINGMAX];
|
|
43 // Angle of the singularities measured from x
|
|
44 cpFloat Gangle[SINGMAX];
|
|
45 // Force function
|
|
46 SingForceFunc force_func[SINGMAX];
|
|
47 // Force function
|
|
48 SingForceFunc torque_func[SINGMAX];
|
|
49 }
|
|
50 alias ActorSingularity Sing;
|
|
51
|
|
52 // Data for the force functions
|
|
53 struct DataforForce{
|
|
54 //Everything in global coordinates
|
|
55 // Position of the source
|
|
56 cpVect p0;
|
|
57 // Observed position
|
|
58 cpVect p;
|
|
59 // Relative position source-observed
|
|
60 cpVect relp;
|
|
61 // distance, disntace^2, ditance ^3
|
|
62 cpFloat r[3];
|
|
63 // angle of the source
|
|
64 cpFloat ang0;
|
|
65 // angle of the observed singularity
|
|
66 cpFloat ang;
|
|
67 // Foce value
|
|
68 cpVect F;
|
|
69 // Torque value
|
|
70 cpFloat T;
|
|
71 }
|
|
72 alias DataforForce ForceData;
|
|
73
|
|
74 // Global Varibales
|
|
75 static cpSpace *space;
|
|
76
|
|
77
|
|
78 // **** Forces ****** //
|
|
79 // Calculate the forces between two bodies. all this functions requieres
|
|
80 // a pointer to an structure with the necessary fields.
|
|
81
|
|
82 // forces between charges
|
|
83 static void
|
|
84 CoulombForce(ForceData* data){
|
|
85 data.F=cpvmult(cpvnormalize(data.relp),COU_MKS/data.r[1]);
|
|
86 }
|
|
87
|
|
88 // forces between magnets
|
|
89 static void
|
|
90 MagDipoleForce(ForceData* data){
|
|
91 static cpFloat phi,alpha,beta,Fr,Fphi;
|
|
92
|
|
93 // Angle of the relative position vector
|
|
94 phi=cpvtoangle(data.relp);
|
|
95 alpha=data.ang0;
|
|
96 beta=data.ang;
|
|
97
|
|
98 alpha =phi - alpha;
|
|
99 beta = phi - beta;
|
|
100
|
|
101
|
|
102 // Components in polar coordinates
|
|
103 Fr=(2.0e0*cos(alpha)*cos(beta) - sin(alpha)*sin(beta));
|
|
104 Fphi=sin(alpha+beta);
|
|
105 // printf("%g %g %g %g %g\n",phi,alpha,beta,Fphi);
|
|
106
|
|
107 // Cartesian coordinates
|
|
108 data.F=cpv(Fr*cos(phi)-Fphi*sin(phi),Fr*sin(phi)+Fphi*cos(phi));
|
|
109 data.F=cpvmult(data.F,-3.e0*MAG_MKS/(data.r[1]*data.r[1]));
|
|
110 }
|
|
111
|
|
112 static void
|
|
113 MagDipoleTorque(ForceData* data){
|
|
114 static cpFloat phi,alpha,beta;
|
|
115
|
|
116 phi=cpvtoangle(data.relp);
|
|
117 alpha=data.ang0;
|
|
118 beta=data.ang;
|
|
119 alpha =phi - alpha;
|
|
120 beta = phi - beta;
|
|
121
|
|
122 // Torque. Though we could use a component of F to save some space,
|
|
123 // we use another variables for the sake of clarity.
|
|
124
|
|
125 data.T=(MAG_MKS/data.r[2])*(3.0e0*cos(alpha)*sin(beta) + sin(alpha-beta));
|
|
126 }
|
|
127 // ******* //
|
|
128
|
|
129 // This function fills the data structure for the force functions
|
|
130 // The structure Sing has the information about the singularity (charge or magnet)
|
|
131 static void
|
|
132 FillForceData(Sing* source,int inds, Sing* obs,int indo, ForceData* data)
|
|
133 {
|
|
134 // Global Position and orientation of the source singularity
|
|
135 data.p0=source.Gpos[inds];
|
|
136 data.ang0=source.Gangle[inds];
|
|
137
|
|
138 // Global Position and orientation of the observed singularity
|
|
139 data.p=obs.Gpos[indo];
|
|
140 data.ang=obs.Gangle[indo];
|
|
141
|
|
142 // Derived magnitudes
|
|
143 data.relp=cpvsub(data.p,data.p0); //Relative position
|
|
144 data.r[0]=cpvlength(data.relp); // Distance
|
|
145 data.r[1]=cpvlengthsq(data.relp); // Square Distance
|
|
146 data.r[2]=data.r[0]*data.r[1]; // Cubic distance
|
|
147
|
|
148 source.force_func[inds](data); // The value of the force
|
|
149 data.F= cpvmult(data.F,source.value[inds]*obs.value[indo]);
|
|
150 }
|
|
151
|
|
152 // Calculation of the interaction
|
|
153 static void
|
|
154 LRangeForceApply(cpBody *a, cpBody *b){
|
|
155
|
|
156 Sing* aux = cast(Sing*)a.data;
|
|
157 Sing* aux2 = cast(Sing*)b.data;
|
|
158 cpVect delta;
|
|
159 // General data needed to calculate interaction
|
|
160 static ForceData fdata;
|
|
161 fdata.F=cpvzero;
|
|
162
|
|
163 // Calculate the forces between the charges of different bodies
|
|
164 for (int i=0; i<aux.Nsing; i++)
|
|
165 {
|
|
166 for (int j=0; j<aux2.Nsing; j++)
|
|
167 {
|
|
168 if(!strcmp(aux.type[i].ptr,aux2.type[j].ptr))
|
|
169 {
|
|
170 //printf("%s %s\n",aux.type[i],aux2.type[j]);
|
|
171 FillForceData (aux2,j,aux,i,&fdata);
|
|
172
|
|
173 //Force applied to body A
|
|
174 delta=cpvsub(aux.Gpos[i], a.p);
|
|
175 cpBodyApplyForce(a,fdata.F, delta);
|
|
176
|
|
177 if(aux.torque_func[i] != null)
|
|
178 {
|
|
179 //Torque on A
|
|
180 aux.torque_func[i](&fdata);
|
|
181 a.t += aux.value[i]*aux2.value[j]*fdata.T;
|
|
182
|
|
183 }
|
|
184 }
|
|
185 }
|
|
186 }
|
|
187 }
|
|
188
|
|
189 // function for the integration of the positions
|
|
190 // The following functions are variations to the starndrd integration in Chipmunk
|
|
191 // you can go ack to the standard ones by doing the appropiate changes.
|
|
192 static void
|
|
193 ChargedBodyUpdatePositionVerlet(cpBody *_body, cpFloat dt)
|
|
194 {
|
|
195 // Long range interaction
|
|
196 cpArray *bodies = space.bodies;
|
|
197 static cpBody* B;
|
|
198 Sing* aux=cast(Sing*)_body.data;
|
|
199 Sing* aux2;
|
|
200
|
|
201 // General data needed to calculate interaction
|
|
202 static ForceData fdata;
|
|
203 fdata.F=cpvzero;
|
|
204
|
|
205 for(int i=0; i< bodies.num; i++)
|
|
206 {
|
|
207 B=cast(cpBody*)bodies.arr[i];
|
|
208 aux2=cast(Sing*)B.data;
|
|
209
|
|
210 if(B != _body)
|
|
211 {
|
|
212 // Calculate the forces between the singularities of different bodies
|
|
213 LRangeForceApply(_body, B);
|
|
214 }
|
|
215 }
|
|
216
|
|
217 cpVect dp = cpvmult(cpvadd(_body.v, _body.v_bias), dt);
|
|
218 dp = cpvadd(dp,cpvmult(cpvmult(_body.f, _body.m_inv), 0.5e0*dt*dt));
|
|
219 _body.p = cpvadd(_body.p, dp);
|
|
220
|
|
221 cpBodySetAngle(_body, cast(float)(_body.a + (_body.w + _body.w_bias)*dt
|
|
222 + 0.5*_body.t*_body.i_inv*dt*dt));
|
|
223
|
|
224 // Update position of the singularities
|
|
225 aux = cast(Sing*)_body.data;
|
|
226 for (int i=0; i<aux.Nsing; i++)
|
|
227 {
|
|
228 aux.Gpos[i]=cpvadd(_body.p,cpvrotate(cpv(aux.position[i].x,
|
|
229 aux.position[i].y), _body.rot));
|
|
230 aux.Gangle[i]= aux.angle[i] + _body.a;
|
|
231 }
|
|
232
|
|
233
|
|
234 _body.v_bias = cpvzero;
|
|
235 _body.w_bias = 0.0f;
|
|
236 }
|
|
237
|
|
238 // function for the integration of the velocities
|
|
239 static void
|
|
240 ChargedBodyUpdateVelocityVerlet(cpBody *_body, cpVect gravity, cpFloat damping, cpFloat dt)
|
|
241 {
|
|
242 _body.v = cpvadd(_body.v, cpvmult(cpvadd(gravity, cpvmult(_body.f, _body.m_inv)), 0.5e0*dt));
|
|
243 _body.w = _body.w + _body.t*_body.i_inv*0.5e0*dt;
|
|
244
|
|
245 _body.f = cpvzero;
|
|
246 _body.t = 0.0e0;
|
|
247
|
|
248 // Long range interaction
|
|
249 cpArray *bodies = space.bodies;
|
|
250 static cpBody* B;
|
|
251
|
|
252 // General data needed to calculate interaction
|
|
253 static ForceData fdata;
|
|
254 fdata.F=cpvzero;
|
|
255
|
|
256 for(int i=0; i< bodies.num; i++)
|
|
257 {
|
|
258 B=cast(cpBody*)bodies.arr[i];
|
|
259
|
|
260 if(B != _body)
|
|
261 {
|
|
262 // Calculate the forces between the singularities of different bodies
|
|
263 LRangeForceApply(_body, B);
|
|
264 }
|
|
265 }
|
|
266 _body.v = cpvadd(cpvmult(_body.v,damping), cpvmult(cpvadd(gravity, cpvmult(_body.f, _body.m_inv)), 0.5e0*dt));
|
|
267 _body.w = _body.w*damping + _body.t*_body.i_inv*0.5e0*dt;
|
|
268 }
|
|
269
|
|
270 static void
|
|
271 update(int ticks)
|
|
272 {
|
|
273 enum int steps = 10;
|
|
274 enum cpFloat dt = 1.0/60.0/cast(cpFloat)steps;
|
|
275
|
|
276 cpArray *bodies = space.bodies;
|
|
277
|
|
278 for(int i=0; i< bodies.num; i++)
|
|
279 cpBodyResetForces(cast(cpBody*)bodies.arr[i]);
|
|
280
|
|
281 for(int i=0; i<steps; i++){
|
|
282 cpSpaceStep(space, dt);
|
|
283 }
|
|
284
|
|
285 }
|
|
286
|
|
287 static void
|
|
288 make_mag(cpVect p, cpFloat ang, cpFloat mag)
|
|
289 {
|
|
290 int nverts=6;
|
|
291 cpVect verts[] = [
|
|
292 cpv(-10,-10),
|
|
293 cpv(-10, 10),
|
|
294 cpv( 10, 10),
|
|
295 cpv( 15, 5),
|
|
296 cpv( 15, -5),
|
|
297 cpv( 10,-10)
|
|
298 ];
|
|
299
|
|
300 cpBody *_body = cpBodyNew(1.0, cpMomentForPoly(1.0f, nverts, verts.ptr, cpvzero));
|
|
301 _body.p = p;
|
|
302 _body.v = cpvzero;
|
|
303 cpBodySetAngle(_body, ang);
|
|
304 _body.w = 0.0e0;
|
|
305
|
|
306 // Load the singularities
|
|
307 Sing *magnet=cast(Sing*)cpmalloc(Sing.sizeof);
|
|
308 magnet.Nsing=1;
|
|
309 magnet.value[0]=mag;
|
|
310 sprintf(magnet.type[0].ptr,"magdipole");
|
|
311
|
|
312 // The position and angle could be different form the one of the body
|
|
313 magnet.position[0]=cpvzero;
|
|
314 magnet.Gpos[0]=cpvadd(p,magnet.position[0]);
|
|
315 magnet.angle[0]=0.0f;
|
|
316 magnet.Gangle[0]=ang;
|
|
317
|
|
318 magnet.force_func[0]=&MagDipoleForce;
|
|
319 magnet.torque_func[0]=&MagDipoleTorque;
|
|
320
|
|
321 _body.data=magnet;
|
|
322
|
|
323 _body.position_func=&ChargedBodyUpdatePositionVerlet;
|
|
324 _body.velocity_func=&ChargedBodyUpdateVelocityVerlet;
|
|
325 cpSpaceAddBody(space, _body);
|
|
326
|
|
327 cpShape *shape = cpPolyShapeNew(_body, nverts, verts.ptr, cpvzero);
|
|
328 shape.e = 0.0; shape.u = 0.7;
|
|
329 cpSpaceAddShape(space, shape);
|
|
330 }
|
|
331
|
|
332 static void
|
|
333 make_charged(cpVect p, cpFloat chg)
|
|
334 {
|
|
335 int nverts=4;
|
|
336 cpVect verts[] = [
|
|
337 cpv(-10,-10),
|
|
338 cpv(-10, 10),
|
|
339 cpv( 10, 10),
|
|
340 cpv( 10,-10)
|
|
341 ];
|
|
342
|
|
343 cpBody *_body = cpBodyNew(1.0, cpMomentForPoly(1.0, nverts, verts.ptr, cpvzero));
|
|
344 _body.p = p;
|
|
345 _body.v = cpvzero;
|
|
346 cpBodySetAngle(_body, 0.0f);
|
|
347 _body.w = 0.0e0;
|
|
348
|
|
349 // Load the singularities
|
|
350 Sing *charge=cast(Sing*)cpmalloc(Sing.sizeof);;
|
|
351 charge.Nsing=1;
|
|
352 charge.value[0]=chg;
|
|
353 sprintf(charge.type[0].ptr,"electrical\0");
|
|
354
|
|
355 // The position and angle could be different form the one of the body
|
|
356 charge.position[0]=cpvzero;
|
|
357 charge.Gpos[0]=cpvadd(p,charge.position[0]);
|
|
358 charge.Gangle[0]=0;
|
25
|
359 charge.angle[0]=0.0f;
|
16
|
360
|
|
361 charge.force_func[0]=&CoulombForce;
|
|
362 charge.torque_func[0]=null;
|
|
363
|
|
364 _body.data=charge;
|
|
365
|
|
366 _body.position_func=&ChargedBodyUpdatePositionVerlet;
|
|
367 _body.velocity_func=&ChargedBodyUpdateVelocityVerlet;
|
|
368 cpSpaceAddBody(space, _body);
|
|
369
|
|
370 cpShape *shape = cpPolyShapeNew(_body, nverts, verts.ptr, cpvzero);
|
|
371 shape.e = 0.0; shape.u = 0.7;
|
|
372 cpSpaceAddShape(space, shape);
|
|
373 }
|
|
374 void
|
|
375 make_mix(cpVect p, cpFloat ang, cpFloat mag,cpFloat chg)
|
|
376 {
|
|
377 int nverts=5;
|
|
378 cpVect verts[] = [
|
|
379 cpv(-10,-10),
|
|
380 cpv(-10, 10),
|
|
381 cpv( 10, 10),
|
|
382 cpv( 20, 0),
|
|
383 cpv( 10,-10)
|
|
384 ];
|
|
385
|
|
386 cpBody *_body = cpBodyNew(1.0, cpMomentForPoly(1.0, nverts, verts.ptr, cpvzero));
|
|
387 _body.p = p;
|
|
388 _body.v = cpvzero;
|
|
389 cpBodySetAngle(_body, ang);
|
|
390 _body.w = 0.0e0;
|
|
391
|
|
392 // Load the singularities
|
|
393 Sing *mix=cast(Sing*)cpmalloc(Sing.sizeof);;
|
|
394 mix.Nsing=2;
|
|
395 mix.value[0]=mag;
|
|
396 mix.value[1]=chg;
|
|
397 sprintf(mix.type[0].ptr,"magdipole\0");
|
|
398 sprintf(mix.type[1].ptr,"electrical\0");
|
|
399
|
|
400 // The position and angle could be different form the one of the body
|
|
401 mix.position[0]=cpvzero;
|
|
402 mix.Gpos[0]=cpvadd(p,mix.position[0]);
|
|
403 mix.position[1]=cpvzero;
|
|
404 mix.Gpos[1]=cpvadd(p,mix.position[1]);
|
|
405 mix.Gangle[0]=ang;
|
|
406 mix.Gangle[1]=ang;
|
25
|
407 mix.angle[0]=0.0f;
|
|
408 mix.angle[1]=0.0f;
|
16
|
409
|
|
410 mix.force_func[0]=&MagDipoleForce;
|
|
411 mix.force_func[1]=&CoulombForce;
|
|
412 mix.torque_func[0]=&MagDipoleTorque;
|
|
413 mix.torque_func[1]=null;
|
|
414
|
|
415 _body.data=mix;
|
|
416
|
|
417 _body.position_func=&ChargedBodyUpdatePositionVerlet;
|
|
418 _body.velocity_func=&ChargedBodyUpdateVelocityVerlet;
|
|
419 cpSpaceAddBody(space, _body);
|
|
420
|
|
421 cpShape *shape = cpPolyShapeNew(_body, nverts, verts.ptr, cpvzero);
|
|
422 shape.e = 0.0; shape.u = 0.7;
|
|
423 cpSpaceAddShape(space, shape);
|
|
424 }
|
|
425
|
|
426
|
|
427 static cpSpace*
|
|
428 init()
|
|
429 {
|
|
430 cpResetShapeIdCounter();
|
|
431 space = cpSpaceNew();
|
|
432 space.iterations = 5;
|
|
433 space.gravity = cpvzero; //cpv(0,-100);
|
|
434
|
|
435 cpSpaceResizeActiveHash(space, 30.0, 2999);
|
|
436
|
|
437 // Screen border
|
|
438 /* shape = cpSegmentShapeNew(staticBody, cpv(-320,-240), cpv(-320,240), 0.0f);
|
|
439 shape.e = 1.0; shape.u = 1.0;
|
|
440 cpSpaceAddShape(space, shape);
|
|
441
|
|
442 shape = cpSegmentShapeNew(staticBody, cpv(320,-240), cpv(320,240), 0.0f);
|
|
443 shape.e = 1.0; shape.u = 1.0;
|
|
444 cpSpaceAddShape(space, shape);
|
|
445
|
|
446 shape = cpSegmentShapeNew(staticBody, cpv(-320,-240), cpv(320,-240), 0.0f);
|
|
447 shape.e = 1.0; shape.u = 1.0;
|
|
448 cpSpaceAddShape(space, shape);
|
|
449
|
|
450 // Reference line
|
|
451 // Does not collide with other objects, we just want to draw it.
|
|
452 shape = cpSegmentShapeNew(staticBody, cpv(-320,0), cpv(320,0), 0.0f);
|
|
453 shape.collision_type = 1;
|
|
454 cpSpaceAddShape(space, shape);
|
|
455 // Add a collision pair function to filter collisions
|
|
456 cpSpaceAddCollisionPairFunc(space, 0, 1, null, null);
|
|
457 */
|
|
458
|
|
459 //srand(cast(uint) time(null));
|
|
460 cpVect p;
|
|
461 cpFloat ang;
|
|
462
|
|
463 // Create magnets
|
|
464 for(int i=0; i<NMAG; i++)
|
|
465 {
|
|
466 p.x=(2.0e0*frand() - 1.0e0)*WIDTH/2.0f;
|
|
467 p.y=(2.0e0*frand() - 1.0e0)*HEIGHT/2.0f;
|
|
468 ang=(2.0e0*frand() - 1.0e0)*3.1415;
|
|
469 make_mag(p, ang,1.0e7);
|
|
470 }
|
|
471
|
|
472 // Create charged objects
|
|
473 for(int i=0; i<NCHG; i++)
|
|
474 {
|
|
475 p.x=(2.0e0*frand() - 1.0e0)*WIDTH/2.0f;
|
|
476 p.y=(2.0e0*frand() - 1.0e0)*HEIGHT/2.0f;
|
|
477 ang=(2.0e0*frand() - 1.0e0)*3.1415;
|
|
478 make_charged(p,1.0e-3*pow(-1.0,i%2));
|
|
479 }
|
|
480
|
|
481 // Create charged magnets objects
|
|
482 for(int i=0; i<NMIX; i++)
|
|
483 {
|
|
484 p.x=(2.0e0*frand() - 1.0e0)*WIDTH/2.0f;
|
|
485 p.y=(2.0e0*frand() - 1.0e0)*HEIGHT/2.0f;
|
|
486 ang=(2.0e0*frand() - 1.0e0)*3.1415;
|
|
487 make_mix(p, ang,1.0e7*pow(-1.0,i%2), 1.0e-3*pow(-1.0,i%2));
|
|
488 }
|
|
489
|
|
490 return space;
|
|
491 }
|
|
492
|
|
493 static void
|
|
494 destroy()
|
|
495 {
|
|
496 cpSpaceFreeChildren(space);
|
|
497 cpSpaceFree(space);
|
|
498 }
|
|
499
|
|
500 chipmunkDemo MagnetsElectric = {
|
|
501 "Magnets and Electric Charges (By: Juan Pablo Carbajal)",
|
|
502 null,
|
|
503 &init,
|
|
504 &update,
|
|
505 &destroy,
|
|
506 };
|