home *** CD-ROM | disk | FTP | other *** search
/ Aminet 10 / aminetcdnumber101996.iso / Aminet / misc / sci / StarCollapse.lha / grav.c < prev    next >
C/C++ Source or Header  |  1995-08-20  |  15KB  |  631 lines
  1. //----------------------------------------------------------------------------
  2. // 3-D Gravity Simulator                                    © 1995 L. Vanhelsuwé
  3. // ---------------------                                    --------------------
  4. //
  5. // Model a "cloud" of masses interacting in 3 dimensions under the influence
  6. // of gravitational attraction.
  7. //
  8. // History:
  9. // --------
  10. // 07-AUG-95: started this file from standard skeleton NEW.c file
  11. // 08-AUG-95: seem to have the bare bones done: physics OK.
  12. // 10-AUG-95: now generating data file with 3-D animation for 3D viewer
  13. // 12-AUG-95: modified algorithm to do vector sum of forces. Code now CRAWLS.
  14. // 20-AUG-95: rewrote bottleneck (calc_pull) in 68882 FPU assembler.
  15. //
  16. // BUGS:
  17. // -----
  18. //----------------------------------------------------------------------------
  19.  
  20.     // some standard ANSI includes
  21.  
  22. #include    <stdio.h>
  23. #include    <string.h>
  24. #include    <stdlib.h>
  25.  
  26.     // some Amiga specific includes
  27.  
  28. #include    <exec/types.h>
  29. #include    <exec/memory.h>
  30.  
  31. #include    <dos/dos.h>
  32. #include    <dos/datetime.h>
  33.  
  34. #include    <intuition/screens.h>
  35. #include    <intuition/intuitionbase.h>
  36.  
  37. #include    <graphics/gfx.h>
  38.  
  39.     // some Amiga includes containing system call prototypes
  40.  
  41. #include    <clib/alib_protos.h>            // ANSI function prototypes
  42. #include    <clib/exec_protos.h>
  43. #include    <clib/dos_protos.h>
  44. #include    <clib/intuition_protos.h>
  45. #include    <clib/graphics_protos.h>
  46.  
  47. //-----------------------------------------------------------------------------
  48. // Program Constants, Limits
  49. // -------------------------
  50.  
  51. #define    GRAVITY_CONST        (100)
  52. #define    MASS_MASS            (1)
  53.  
  54. #define    STOP_BUTTON            !(*((char*)0xBFE001) & 0x80)
  55.  
  56. #define    DUMP_INFO            if (options.audit) printf
  57.  
  58. #define    CENTER_X            (350)
  59. #define    CENTER_Y            (250)
  60.  
  61. #define    CLUSTER_RADIUS        (250)
  62.  
  63. //-----------------------------------------------------------------------------
  64. // User Types
  65. // ----------
  66.  
  67. typedef ULONG    Display;
  68.  
  69. typedef    double    num;
  70.  
  71. typedef struct Mass {
  72.     num    mass;
  73.     num    x,y,z;
  74.     num    Vx,Vy,Vz;
  75.  
  76. } Mass;
  77.  
  78. typedef struct AnimFileHeader {
  79.     short    startframe;
  80.     short    numframes;
  81.     short    framesize;
  82.  
  83. } AnimFileHeader;
  84.  
  85. struct RDAres {                    // array for ReadArgs()
  86.  
  87.     ULONG frames;                // number of animation frames wanted
  88.     ULONG masses;                // number of masses taking part in simulation
  89.     ULONG delay;                // optional delay to slow down animation
  90.     ULONG audit;                // dump intermediate calculation results to stdout
  91.     ULONG nodisplay;            // do not display simulation (for faster simulation)
  92.  
  93. } options;
  94.  
  95. //-----------------------------------------------------------------------------
  96. // Function prototypes
  97. // -------------------
  98.  
  99. void        main                    (void);
  100.  
  101. void        init_platform_specific    (void);
  102. void        close_platform_specific    (void);
  103.  
  104. Display        open_display            (void);
  105. void        close_display            (Display d);
  106. void        erase_display            (Display d);
  107. void        write_pixel                (int x, int y);
  108.  
  109. BOOL        init_simulation            (void);
  110.  
  111. void        simulate_gravity        (void);
  112. void        display_cloud            (void);
  113. void        dump_simulation_frame    (void);
  114.  
  115. void            calc_pull            (int m, num*, num*, num*, num*);
  116. void        asm_calc_pull            (int m, num*, num*, num*, num*);
  117.  
  118. void        init_defaults            (struct variable *varptr);
  119. void        print_frame_number        (int fr);
  120. void        write_fileheader        (void);
  121.  
  122. double        random                    (int range);
  123.  
  124. //-----------------------------------------------------------------------------
  125. // Global statics
  126. // --------------
  127.  
  128. struct    IntuitionBase    *IntuitionBase;        // for access to ActiveWindow
  129. struct    GfxBase            *GfxBase;
  130.  
  131. struct    Screen            *sc;
  132. struct    RastPort        *rp;
  133.  
  134. Mass     *Cloud;                                // Mass Cloud[ CloudSize ];
  135.  
  136. #define    NumFrames    options.frames
  137. #define    AUDIT        options.audit
  138. #define NODISPLAY    options.nodisplay
  139.  
  140. Display    window;
  141.  
  142. int        fh;                                // output animation file handle
  143. short    *buffer;                        //    buffer[3*MAX_MASSES +1]; full data frame buffer
  144.  
  145. int        frame;
  146.  
  147. ULONG    CloudSize;
  148.  
  149. char    GRAV_CLI_TEMPLATE[] = "FRAMES/K,MASSES/K,DELAY/K,AUDIT/S,NODISPLAY/S";
  150.  
  151. //-----------------------------------------------------------------------------
  152. //
  153. //-----------------------------------------------------------------------------
  154.  
  155. void main (void) {
  156.  
  157. BOOL            ok;
  158.  
  159.     init_platform_specific();
  160.  
  161.     CloudSize    = options.masses;
  162.  
  163.     ok = init_simulation();
  164.     if (! ok) {
  165.         close_platform_specific();
  166.         exit(10);
  167.     }
  168.  
  169.     if (! NODISPLAY) {
  170.         window = open_display();
  171.     }
  172.  
  173.     write_fileheader();
  174.  
  175.     frame = NumFrames;
  176.  
  177.     while ( frame-- != 0) {
  178.  
  179.         if (STOP_BUTTON) break;
  180.  
  181.         simulate_gravity();
  182.  
  183.         dump_simulation_frame();
  184.  
  185.         if ( ! NODISPLAY) display_cloud();
  186.  
  187.         Delay(options.delay);
  188.     }
  189.  
  190.     free(Cloud);
  191.     free(buffer);
  192.  
  193.     if (! NODISPLAY) {
  194.         close_display( window );
  195.     }
  196.  
  197.     close_platform_specific();
  198. }
  199. //-----------------------------------------------------------------------------
  200. //
  201. //-----------------------------------------------------------------------------
  202. void        write_fileheader        (void) {
  203.  
  204. AnimFileHeader    an;
  205.  
  206.     an.startframe    = 0;                        // start with frame #0
  207.     an.numframes    = NumFrames;                // max frame #
  208.     an.framesize    = 2 + 6 * CloudSize;        // size of a data frame
  209.     Write(fh, &an, sizeof(an));                    // write animation file header
  210. }
  211.  
  212. //-----------------------------------------------------------------------------
  213. // Create starting model of cloud.
  214. //-----------------------------------------------------------------------------
  215.  
  216. BOOL        init_simulation            (void) {
  217.  
  218. int        m;
  219.  
  220. num        x,y,z,R;
  221.  
  222.     Cloud = malloc( CloudSize * sizeof(Mass));
  223.  
  224.     if (! Cloud) {
  225.         printf("Failed to allocate Cloud data structure.\n");
  226.         return FALSE;
  227.     }
  228.  
  229.     buffer = malloc( sizeof(short) * 3 * (CloudSize+1) );
  230.  
  231.     if (! buffer) {
  232.         printf("Failed to allocate write buffer.\n");
  233.         return FALSE;
  234.     }
  235.  
  236.         // Generate ball-shaped cluster of randomly distributed masses
  237.  
  238.     for (m=0; m< CloudSize; m++) {
  239.  
  240.         R = CLUSTER_RADIUS + 1;            // ensure first while cond fails
  241.  
  242.         while ( R > CLUSTER_RADIUS ) {
  243.  
  244.             x = random(CLUSTER_RADIUS);
  245.             y = random(CLUSTER_RADIUS);
  246.             z = random(CLUSTER_RADIUS);
  247.  
  248. //            printf("\nTrying %f,%f,%f \n", x,y,z);
  249.  
  250.             R = sqrt(x*x + y*y + z*z);
  251.  
  252. //            printf("Radius = %f\n", R);
  253.         }
  254.  
  255.         Cloud[m].mass    = 1.0 + (rand()%10) * MASS_MASS;
  256.  
  257.         Cloud[m].x        = x;
  258.         Cloud[m].y        = y;
  259.         Cloud[m].z        = z;
  260.  
  261.         Cloud[m].Vx        = 0;
  262.         Cloud[m].Vy        = 0;
  263.         Cloud[m].Vz        = 0;
  264.  
  265.         DUMP_INFO ("Body %03d (x,y,z) = (%f,%f,%f)\n", m, x,y,z);
  266.     }
  267.  
  268. //    printf("Initial Mass cloud generated.\n");
  269.  
  270.     return TRUE;
  271. }
  272.  
  273. //-----------------------------------------------------------------------------
  274. //-----------------------------------------------------------------------------
  275. double random (int range) {
  276.  
  277. int throw;
  278.  
  279.     throw = rand() % (2*range) - range;
  280.     return (double) throw;
  281. }
  282. //-----------------------------------------------------------------------------
  283. //
  284. //-----------------------------------------------------------------------------
  285.  
  286. void        simulate_gravity        (void) {
  287.  
  288. int        m;
  289.  
  290. num        F,Fx,Fy,Fz;
  291. num        V,Vx,Vy,Vz;
  292.  
  293.     for (m=0; m < CloudSize; m++) {
  294.  
  295.     //        calc_pull (m, &F, &Fx, &Fy, &Fz);
  296.         asm_calc_pull (m, &F, &Fx, &Fy, &Fz);
  297.  
  298.         V = F / Cloud[m].mass;
  299.  
  300.         Vx = Vy = Vz = 0;
  301.  
  302.         if (F != 0) {
  303.             Vx = V * Fx / F;
  304.             Vy = V * Fy / F;
  305.             Vz = V * Fz / F;
  306.         }
  307.  
  308.         Cloud[m].Vx += Vx;
  309.         Cloud[m].Vy += Vy;
  310.         Cloud[m].Vz += Vz;
  311.  
  312.         Cloud[m].x += Cloud[m].Vx;
  313.         Cloud[m].y += Cloud[m].Vy;
  314.         Cloud[m].z += Cloud[m].Vz;
  315.     }
  316. }
  317.  
  318. //-----------------------------------------------------------------------------
  319. //
  320. // Calculate total force acting on 1 body by summing all VECTOR forces
  321. // stemming from individual interactions with all other bodies.
  322. //
  323. //-----------------------------------------------------------------------------
  324.  
  325. void    calc_pull (int star, num *F_ptr, num *Fx_ptr, num *Fy_ptr, num *Fz_ptr) {
  326.  
  327. int m;
  328. num    F,Fx,Fy,Fz;
  329. num Fcouple;                    // force magnitude between 2 masses
  330. num    Fxcouple,Fycouple,Fzcouple;    // force vector
  331. num    dx,dy,dz;
  332. num    R;
  333.  
  334.     DUMP_INFO ("\nCalculating pull on body %03d\n", star);
  335.  
  336.         // reset force accumulators
  337.  
  338.     Fx = Fy = Fz = 0;
  339.  
  340.         // forall bodies (except body under consideration)
  341.  
  342.     for (m=0; m < CloudSize; m++) {
  343.  
  344.         if (m == star) continue;            // skip oneself...
  345.  
  346.             // calc distance between bodies
  347.  
  348.         dx = Cloud[m].x - Cloud[star].x;
  349.         dy = Cloud[m].y - Cloud[star].y;
  350.         dz = Cloud[m].z - Cloud[star].z;
  351.     
  352.         R = dx*dx + dy*dy + dz*dz;
  353.  
  354.         DUMP_INFO ("Squared distance %f (R=%f)\n", R, sqrt((double)R));
  355.     
  356.             // If body collides with another, no force results
  357.     
  358.         if (R == 0) continue;
  359.     
  360.         Fcouple = GRAVITY_CONST * Cloud[star].mass * Cloud[m].mass / R;
  361.  
  362.         DUMP_INFO ("Resulting intermediate force = %f\n", Fcouple);
  363.  
  364.         R = sqrt(R);
  365.     
  366.         if (R == 0) continue;
  367.     
  368.         Fxcouple = Fcouple * dx/R;
  369.         Fycouple = Fcouple * dy/R;
  370.         Fzcouple = Fcouple * dz/R;
  371.  
  372.         DUMP_INFO ("Intermed Fx = %03f Fy = %03f Fz = %03f\n", Fxcouple, Fycouple, Fzcouple);
  373.  
  374.         Fx += Fxcouple;
  375.         Fy += Fycouple;
  376.         Fz += Fzcouple;
  377.     }
  378.  
  379.     F = sqrt( Fx*Fx + Fy*Fy + Fz*Fz );
  380.  
  381.     DUMP_INFO ("SUM      Fx = %03f Fy = %03f Fz = %03f\n", Fx, Fy, Fz);
  382.     DUMP_INFO ("Resulting Force = %f\n", F);
  383.  
  384.     *F_ptr    = F;
  385.     *Fx_ptr    = Fx;
  386.     *Fy_ptr    = Fy;
  387.     *Fz_ptr    = Fz;
  388. }
  389.  
  390. //-----------------------------------------------------------------------------
  391. //
  392. //
  393. //-----------------------------------------------------------------------------
  394.  
  395. void        dump_simulation_frame    (void) {
  396.  
  397. int        m;
  398. int        x,y,z;
  399. short    *ptr;
  400.  
  401.     ptr = buffer;
  402.  
  403.     *ptr++ = CloudSize;                    // store number of vertices in frame
  404.  
  405.     for (m=0; m < CloudSize; m++) {
  406.  
  407.         x = Cloud[m].x + 0.5;            // round correctly !!
  408.         y = Cloud[m].y + 0.5;
  409.         z = Cloud[m].z + 0.5;
  410.  
  411.         *ptr++ = (short) x;
  412.         *ptr++ = (short) y;
  413.         *ptr++ = (short) z;
  414.     }
  415.  
  416.     Write(fh, buffer, 2 + sizeof(short) * 3 * CloudSize);
  417. }
  418. //-----------------------------------------------------------------------------
  419. //
  420. //
  421. //-----------------------------------------------------------------------------
  422.  
  423. void        display_cloud            (void) {
  424.  
  425. int        m;
  426. int        x,y,z;
  427.  
  428.     erase_display( window );
  429.  
  430.     for (m=0; m < CloudSize; m++) {
  431.  
  432.         x = Cloud[m].x + 0.5;            // round correctly !!
  433.         y = Cloud[m].y + 0.5;
  434.         z = Cloud[m].z + 0.5;
  435.  
  436.         write_pixel( CENTER_X+x , CENTER_Y+y );
  437.     }
  438.  
  439.     print_frame_number(frame);
  440. }
  441.  
  442. //-----------------------------------------------------------------------------
  443. //
  444. //-----------------------------------------------------------------------------
  445. //-----------------------------------------------------------------------------
  446. //
  447. //-----------------------------------------------------------------------------
  448. //-----------------------------------------------------------------------------
  449. //
  450. //-----------------------------------------------------------------------------
  451. //-----------------------------------------------------------------------------
  452. //
  453. //-----------------------------------------------------------------------------
  454. //-----------------------------------------------------------------------------
  455. //
  456. //-----------------------------------------------------------------------------
  457.  
  458. //-----------------------------------------------------------------------------
  459. //
  460. //-----------------------------------------------------------------------------
  461.  
  462. // The Variable structure defines program globals which can be accessed or
  463. // changed via the command line and/or the AREXX interface.
  464.  
  465. struct variable {
  466.  
  467.     char    *name;            // name of a global parameter
  468.     ULONG    *address;        // its address
  469.  
  470.     ULONG    init_val;        // initial default value
  471.     ULONG    min_val;        // range minimum
  472.     ULONG    max_val;        // range maximum
  473. };
  474.  
  475. #define    PARAMETER        (1)
  476.  
  477. struct variable varlist[] = {
  478.  
  479. //  FLAGS/ATTR     NAME            ADDRESS            DEFAULT    MINIMUM    MAXIMUM
  480. //    ----------    -----            -------            -------    -------    -------
  481.  
  482.     {            "CLOUDSIZE"        ,&options.masses,50        ,2        ,1000    },
  483.     {            "NUM_FRAMES"    ,&NumFrames        ,50        ,1        ,1000    },
  484.     {            "DELAY"            ,&options.delay    ,0        ,0        ,200    },
  485.  
  486.     {            NULL            ,NULL            ,NULL    ,0        ,0        }
  487. };
  488.  
  489. //-----------------------------------------------------------------------------
  490.  
  491. void init_platform_specific        (void) {
  492.  
  493. struct RDArgs *rdargs;
  494. struct variable *var;
  495.  
  496. #ifndef DEBUG
  497.     if (! (rdargs = ReadArgs(GRAV_CLI_TEMPLATE, (ULONG*) &options, NULL))) {
  498.         printf("ERROR: ARGUMENTS INCORRECT.\n");
  499.         printf("Type Grav ? for a full syntax template.\n\n");
  500.  
  501.         exit(100);
  502.  
  503.     } // Free RDA only AFTER we've parsed ptrs in options pointing into RDA **!!
  504.  
  505. #endif
  506.  
  507. // Any non-initialized simulation globals get initialized now.
  508.  
  509.     init_defaults(&varlist[0]);            // others are intialized from CMD line.
  510.  
  511. // **!! NOW we can free RDA (AFTER we've parsed ptrs)
  512.  
  513.     FreeArgs(rdargs);
  514.  
  515.     var = &varlist[0];
  516.     while(var->name) {
  517.  
  518.         if (*var->address > var->max_val) {
  519.             printf("%s has to be in the range [%d..%d] ! (%d > %d)\n",
  520.                 var->name, var->min_val, var->max_val, *var->address, var->max_val);
  521.             exit(10);
  522.         } else
  523.  
  524.         if (*var->address < var->min_val) {
  525.             printf("%s has to be in the range [%d..%d] ! (%d < %d)\n",
  526.                 var->name, var->min_val, var->max_val, *var->address, var->min_val);
  527.             exit(10);
  528.         }
  529.  
  530.         var++;                    // goto next var in list
  531.     }
  532.  
  533.     IntuitionBase    =
  534.         (struct IntuitionBase *)        OpenLibrary("intuition.library",0);
  535.     GfxBase            =
  536.         (struct GfxBase *)                OpenLibrary("graphics.library",0);
  537.  
  538.     fh = Output();            // get access to stdout (no need to Close() )
  539. }
  540. //-----------------------------------------------------------------------------
  541. //
  542. //-----------------------------------------------------------------------------
  543.  
  544. void close_platform_specific    (void) {
  545.  
  546.     CloseLibrary( (struct Library *) IntuitionBase);
  547.     CloseLibrary( (struct Library *) GfxBase);
  548. }
  549.  
  550. //-----------------------------------------------------------------------------
  551. //
  552. //-----------------------------------------------------------------------------
  553.  
  554. Display    open_display            (void) {
  555.  
  556.     sc = OpenScreenTags(NULL,
  557.                         SA_LikeWorkbench, TRUE,
  558.                         TAG_END);
  559.  
  560.     rp = &sc->RastPort;
  561.  
  562.     SetAPen(rp, 1);
  563.  
  564.     return (Display) sc;
  565. }
  566.  
  567.  
  568. //-----------------------------------------------------------------------------
  569. //
  570. //-----------------------------------------------------------------------------
  571.  
  572. void close_display                (Display d) {
  573.  
  574.     CloseScreen( (struct Screen *) d);
  575.  
  576. }
  577. //-----------------------------------------------------------------------------
  578. //
  579. //-----------------------------------------------------------------------------
  580.  
  581. void write_pixel                (int x, int y) {
  582.  
  583.     if (x<0 || x>600 || y<0 || y>500) return;
  584.  
  585.     WritePixel(rp, x,y);
  586. }
  587.  
  588. //-----------------------------------------------------------------------------
  589. //
  590. //-----------------------------------------------------------------------------
  591. void        erase_display            (Display d) {
  592.  
  593.     SetRast    (rp, 0);
  594. }
  595.  
  596. //----------------------------------------------------------------
  597. // Set all non-initialized Simulation Global Parameters to a default value
  598. // IF the user hasn't given a value himself.
  599. //----------------------------------------------------------------
  600. void init_defaults    (struct variable *varptr) {
  601.  
  602. ULONG *var;
  603.  
  604.     while (varptr->name) {
  605.         var = varptr->address;
  606.  
  607.         if (*var) {                            // if RDA initialized a ptr to an arg,
  608.             *var = atoi((char*) *var);        // convert ASCII arg to number
  609.         } else {
  610.             *var = varptr->init_val;        // else use built-in default value.
  611.         }
  612.  
  613.         varptr++;                // goto next variable
  614.     }
  615. }
  616. //-----------------------------------------------------------------------------
  617. //
  618. //-----------------------------------------------------------------------------
  619. void        print_frame_number    (int fr) {
  620.  
  621. char    fr_str[11];
  622.  
  623.     sprintf(fr_str, "FRAME %04d", fr);
  624.  
  625.     Move(rp, 40,40);
  626.     Text(rp, fr_str, 10);
  627. }
  628. //-----------------------------------------------------------------------------
  629. //
  630. //-----------------------------------------------------------------------------
  631.