// by Greg M. Johnson http://members.xoom.com/gregjohn/animation.html. // Copyright 2001. /* I'm posting to p.b.s.-f. a pov file which models the effects of gravity on particles in our own solar system. The image below is from a scene file where the Sun, nine planets, and 10 dinosaur-killer-size asteroids are orbiting for a period of 50 earth years. The planets are the light blue/green streaks, the asteroids are red & yellow. Jupiter is the outermost orbit in this field of view. I was motivated to start working on this kind of model when someone said Jupiter was a blessing to earthlings in that its gravitational field protected us fom much asteroid bombardment. (The attached file doesn't answer this question yet.) The code shares its roots with my boid system, which is explained at : http://www.geocities.com/pterandon/boids.html The strengths of the model (and pov file) include: ----------------------------------------------------------------- 1. Actual use of gravity-based equations to determine instantaneous velocities of particles, rather than simple rotation transforms of povray objects. 2. Gravitational interactions between all particles are calculated, not just merely wrt/ the Sun, although it should be possible to see how this isn't too difficult to comment out, thus speeding up the calcs. 3. It is versatile in what you can to with the information: a) display a trace of the movement of particles (as shown in this JPG), or b) write to a file the positions and velocities of each particle at any increment desired (Note: some tweaking on the code on your part is required for writing data to file, but basic #fopen commands, etc. are there in the file. See my boids algo for examples of read & write). 4. The orbits for planets from earth and increasing radius are stable enough for most considerations. 5. The file contains information on the radius of the planets & Sun scaled to the scale of this system (1 pov unit= 1 Gigameter). Weaknesses of the model include: --------------------------------------------- 1. It took me 24 hours to compute this scene spanning 50 years and 20 particles on a 450 MHz Pentium II. Restating, if you hit the "RUN" button on this file right now, it will take 24 hours to get an image on my PII machine! 2. Planets drift. This slow speed is required for stable orbits in the inner planets. You can note in this photo that Mercury (Medium Blue) is starting to drift out into the orbit of Venus (Light Blue). Venus is also starting to drift into the orbit of Earth, which BTW has a nice stable orbit (Light Torquoise). This is a fundamental problem with all of the orbits: THEY ALL DRIFT, the inner ones more severely. Weeks of experimentation plus my own understanding confirm that this is an unavoidable problem with this kind of calculation, and IS NOT a problem with my gravitational constant, placing planets at incorrect velocities or positions. The algorithm essentially computes the forces acting on every particle at increments of 38 MINUTES (i.e., one earth orbit or year = 13600 ticks of the count variable "tt"). To get a perfect system (perfect elliptical orbit), one must take the limit of the system as it approaches zero, or in other words, compute the forces every ZERO SECONDS. Drawing particles that orbit the sun on perfect ellipses prevents one from considering inter-particle gravitational interactions, the original purpose of this POV file. To the extent of the problems listed herein, the file is your own "fixer-upper"! 3. Since it will rarely be the case that you will be interested in the position of particles every 38 minutes, care must be taken to set up the code for display or writing of variables to occur only once every 100, 1000, etc., counts as desired. */ #version unofficial MegaPov 0.5; #declare Grav= 0.0000113; //0.0000105*(10^(.15*clock)); // 0.00113 //#declare fudge=0.0001*(10^(3.5*clock)); #declare fudge=0; #include "colors.inc" /* text { ttf "timrom.ttf" str(Grav,8,7) .1, 0 pigment { Red } finish {ambient 0.8} //rotate<0,180,0> rotate <90,0,0> scale 50 translate <150,20,220> } */ #macro Reorient(Axis1,Axis2) #local vX1=vnormalize(Axis1); //john vansickle #local vX2=vnormalize(Axis2); #local vY=vnormalize(vcross(vX1,vX2)); #local vZ1=vnormalize(vcross(vX1,vY)); #local vZ2=vnormalize(vcross(vX2,vY)); matrix < vX1.x, vY.x,vZ1.x, vX1.y,vY.y,vZ1.y, vX1.z,vY.z, vZ1.z, 0,0,0 > matrix < vX2.x,vX2.y,vX2.z, vY.x,vY.y, vY.z, vZ2.x,vZ2.y,vZ2.z, 0,0,0 > #end #declare num=20 ; //number of actors #declare framenum=13600*50; //number of frames 13600 is one year. #declare RRR=seed(10012) ; #declare actorp= array[num+10] //the current position of each actor #declare actorv= array[num+10] //the current velocity of each actor #declare actoravoid= array[num+10] //the net instantaneous force from gravity on each actor #declare mass= array[num+10] //the mass of each actor #declare colorp= array[num+10] //the color of each actor //#declare radiuss=array[num+1] //the radius of each actor (optional, for some kinds of display only). //#fopen MyFile2 "c:\pov31g\include\boid13c06.inc" write //#write(MyFile2,"#declare Position=array[", framenum-1, "][", num, "][6] {\n") //initialize #declare n=0 ; #while (n<10) #declare colorp[n]=; #declare n=n+1 ; #end /* //radius of each actor, in Gm. Does not enter into gravity calcs. // From The New York Times 1998 Almanac, John W. Wright, ed. #declare radiuss[0]= 6.9599e-1; #declare radiuss[1]= 4880-6; // Gm #declare radiuss[2]= 12104e-6; #declare radiuss[3]= 12756-6; #declare radiuss[4]= 6794e-6; #declare radiuss[5]= 142984e-6; #declare radiuss[6]= 120536e-6; #declare radiuss[7]= 51118e-6; #declare radiuss[8]= 49532e-6; #declare radiuss[9]= 2340e-6; */ //masses divied by 10E24 kg // From The New York Times 1998 Almanac, John W. Wright, ed. #declare mass[0]=2000000; //The Sun #declare mass[1]=.33; //Mercury #declare mass[2]=4.87; //Venus #declare mass[3]=5.97; //Earth #declare mass[4]=.64; //Mars #declare mass[5]=1898; //Jupiter #declare mass[6]=568; //Saturn #declare mass[7]=86; //Uranus #declare mass[8]=102; //Neptune #declare mass[9]=.015; //Pluto //actor positions are textbook values for radius of orbit, in Gm. // From The New York Times 1998 Almanac, John W. Wright, ed. #declare adjust=1;//+fudge; #declare actorp[0]=<0,0,0>*adjust; #declare actorp[1]=<57.98,0,0>*adjust; #declare actorp[2]=<108.2,0,0>*adjust; #declare actorp[3]=<149.6,0,0>*adjust; #declare actorp[4]=<227.98,0,0>*adjust; #declare actorp[5]=<778.3,0,0>*adjust; #declare actorp[6]=<1427.0,0,0>*adjust; #declare actorp[7]=<2869.6,0,0>*adjust; #declare actorp[8]=<4496.6,0,0>*adjust; #declare actorp[9]=<5913.5,0,0>*adjust; //velocity of each actor in m/sec times a factor. // From The New York Times 1998 Almanac, John W. Wright, ed. #declare sf=0.1;//2*0.91;//fudge; #declare actorv[0]=z*0; #declare actorv[1]=sf*z*29.76; // mi/sec #declare actorv[2]=sf*z*21.76; #declare actorv[3]=sf*z*18.51; #declare actorv[4]=sf*z*14.99; #declare actorv[5]=sf*z*8.12; #declare actorv[6]=sf*z*6.00; #declare actorv[7]=sf*z*3.41; #declare actorv[8]=sf*z*2.95; #declare actorv[9]=sf*z*2.95; //catastrophic asteroid size //#declare mass[11]=1.05e-18; //#declare radiuss[11]=500e-6; //This creates random positions, velocities, and masses if you want some asteroids or n>10. #if (num>10) #declare n=10 ; #while (n; #declare n=n+1 ; #end #end //THis optional section rotates particles to give random starting position // ie., not all aligned in one straight line. #declare n=1 ; #while (n; #declare n=n+1 ; #end //Calculate effect of gravity on each of the actors. #declare i=0 ; #while (i look_at <1,0,0> location <.675,2500,-8> look_at <1,0,0> //see pluto //location <.675,5000,-8> look_at <1,0,0> //location <.675,600,-8> look_at <1,0,0> //see first four planets angle 65 }