research/transmission_spectroscopy/universesim/src/main.c

   1 /*
   2  * Universe Simulator
   3  */
   4
   5 #include <stdlib.h>
   6 #include <stdio.h>
   7 #include <tParticle.h>
   8 #include <signal.h>
   9
  10 #define PROTON_MASS     1836
  11 #define NEUTRON_MASS    1839
  12 #define ELECTRON_MASS   1
  13
  14 // === CONSTANTS ===
  15 #define TPL_PROB_GRAN   100000
  16 #define TPL_MASS_SHIFT  0
  17 #define TPL_CHARGE_SHIFT        0
  18 #define TPL_LOCATION_SHIFT      (-16)
  19 //#define       NUM_INIT_PARTICLES      32
  20 //#define       NUM_INIT_PARTICLES      64
  21 //#define       NUM_INIT_PARTICLES      256
  22 #define NUM_INIT_PARTICLES      1024
  23 //#define       NUM_INIT_PARTICLES      0x10000
  24
  25 // === IMPORTS ===
  26 extern void     Visualiser_Initialise(void);
  27 extern void     Update_DoTimestep(void);
  28 extern void     Visualiser_Update(void);
  29
  30 // === PROTOTYPES ===
  31 void    SaveUniverse(int signal);
  32
  33 // === GLOBALS ===
  34 const struct sPartTpl {
  35         size_t  Mass;
  36         size_t  Charge;
  37          int    Probability;    // per `TPL_PROB_GRAN`
  38 }       gaTemplates[] = {
  39         {PROTON_MASS,    1, TPL_PROB_GRAN/4},   // Proton-Like
  40         {NEUTRON_MASS,   0, TPL_PROB_GRAN/4},   // Neutron-Like
  41         {ELECTRON_MASS, -1, TPL_PROB_GRAN/2},   // Electron-Like
  42 };
  43 const int       giNumTemplates = (sizeof(gaTemplates)/sizeof(gaTemplates[0]));
  44 tParticle       gaParticles[NUM_INIT_PARTICLES];
  45 tParticle       *gpParticles = gaParticles;
  46  int    giNumParticles = NUM_INIT_PARTICLES;
  47  int    giCurrentTimestamp = 0;
  48 volatile int    giUpdatePos;
  49 FILE    *gLogFile;
  50
  51 // === CODE ===
  52 int main(int argc, char *argv[])
  53 {
  54          int    i;
  55         FILE    *fp;
  56          int    tmp, len;
  57
  58         signal(SIGINT, SaveUniverse);
  59
  60         printf("Mathematical Universe Simulator\n");
  61         printf(" By John Hodge (thePowersGang)\n");
  62         printf("\n");
  63
  64         if( argc > 1 && (fp = fopen(argv[1], "rb")) )
  65         {
  66                 tmp = 0;
  67                 len = fread(&tmp, 2, 1, fp);
  68                 if(tmp != sizeof(tLargeInt)) {
  69                         fprintf(stderr, "Unable to restore from '%s', largeint size does not match (%i != %i)\n",
  70                                 argv[1], tmp, sizeof(tLargeInt));
  71                         exit(1);
  72                 }
  73                 len = fread(&tmp, 2, 1, fp);
  74                 if(tmp != sizeof(tVector)/sizeof(tLargeInt)) {
  75                         fprintf(stderr, "Unable to restore from '%s', vector order does not match (%i != %i)\n",
  76                                 argv[1], tmp, sizeof(tVector)/sizeof(tLargeInt));
  77                         exit(1);
  78                 }
  79
  80                 len = fread(&tmp, 4, 1, fp);
  81                 if( tmp > NUM_INIT_PARTICLES ) {
  82                         fprintf(stderr, "Unable to restore from '%s', there is not enough static space allocated (%i > %i)\n",
  83                                 argv[1], tmp, NUM_INIT_PARTICLES);
  84                         exit(1);
  85                 }
  86                 giNumParticles = tmp;
  87                 printf("giNumParticles = %i\n", tmp);
  88
  89                 // read timestamp
  90                 tmp = giNumParticles*(sizeof(tLargeInt)*2+sizeof(tVector));
  91                 fseek(fp, 0, SEEK_END);
  92                 if( (ftell(fp) - 8) % tmp ) {
  93                         fprintf(stderr,
  94                                 "Warning: File has trailing bytes (0x%x), may have been an unclean close (Len 0x%x)\n",
  95                                 (ftell(fp) - 8) % tmp, ftell(fp)
  96                                 );
  97                         return 0;
  98                 }
  99                 giCurrentTimestamp = (ftell(fp) - 8) / tmp;
 100                 printf("giCurrentTimestamp = %i\n", giCurrentTimestamp);
 101                 fseek(fp, 8 + (giCurrentTimestamp-1) * tmp, SEEK_SET);
 102                 for( i = 0; i < giNumParticles; i++ )
 103                 {
 104                         len = fread(&gaParticles[i].Location, sizeof(tVector), 1, fp);
 105                         len = fread(&gaParticles[i].Mass, sizeof(tLargeInt), 1, fp);
 106                         len = fread(&gaParticles[i].Charge, sizeof(tLargeInt), 1, fp);
 107                 }
 108
 109                 fclose(fp);
 110
 111                 printf("Restored from '%s' (%i particles)\n", argv[1], giNumParticles);
 112                 printf("WARNING: Restoring a universe state can lead to errors\n");
 113         }
 114         else
 115         {
 116                  int    j;
 117                 size_t  val;
 118                  int    aCounts[giNumTemplates];
 119                 for( j = 0; j < giNumTemplates; j++ )   aCounts[j] = 0;
 120
 121                 printf("Initialising with %i atomic particles\n", giNumParticles);
 122                 printf(" Basing off %i templates\n", giNumTemplates);
 123
 124                 if( argc > 1 ) {
 125                         size_t  seed;
 126                         seed = atoi(argv[1]);
 127                         printf(" Seeding random number generator with %i\n", seed);
 128                         srand( seed );
 129                 }
 130
 131                 printf("Generating particles...\n");
 132                 for( i = 0; i < giNumParticles; i ++ )
 133                 {
 134                         val = rand() % TPL_PROB_GRAN;
 135                         for( j = 0; j < giNumTemplates; j++ )
 136                         {
 137                                 if( val < gaTemplates[j].Probability )  break;
 138                                 val -= gaTemplates[j].Probability;
 139                         }
 140                         if( j == giNumTemplates ) {
 141                                 fprintf(stderr, "ERROR: Templates do not add up to %i\n", TPL_PROB_GRAN);
 142                                 return 1;
 143                         }
 144                         aCounts[j] ++;
 145                         #if TPL_LOCATION_SHIFT < 0
 146                         if( i == 0 )    printf("Setting coords for #0\n");
 147                         LargeInt_SetNative(&gaParticles[i].Location.C[0], rand()>>(-TPL_LOCATION_SHIFT), 0);
 148                         LargeInt_SetNative(&gaParticles[i].Location.C[1], rand()>>(-TPL_LOCATION_SHIFT), 0);
 149                         //LargeInt_SetNative(&gaParticles[i].Location.C[2], rand()>>(-TPL_LOCATION_SHIFT), 0);
 150                         #else
 151                         LargeInt_SetNative(&gaParticles[i].Location.C[0], rand(), TPL_LOCATION_SHIFT);
 152                         LargeInt_SetNative(&gaParticles[i].Location.C[1], rand(), TPL_LOCATION_SHIFT);
 153                         LargeInt_SetNative(&gaParticles[i].Location.C[2], rand(), TPL_LOCATION_SHIFT);
 154                         #endif
 155
 156                         LargeInt_SetNative(&gaParticles[i].Mass, gaTemplates[j].Mass, TPL_MASS_SHIFT);
 157                         LargeInt_SetNative(&gaParticles[i].Charge, gaTemplates[j].Charge, TPL_CHARGE_SHIFT);
 158                         printf("%2i%%\r", i*100/giNumParticles);
 159                 }
 160                 printf("   \n");
 161
 162                 printf("Final counts of particles:\n");
 163                 for( j = 0; j < giNumTemplates; j++ )
 164                 {
 165                         printf("%i: %i\n", j, aCounts[j]);
 166                 }
 167         }
 168
 169         printf("Starting simulation.\n");
 170
 171         Visualiser_Initialise();
 172         Update_Initialise();
 173
 174
 175         gLogFile = fopen("universe.blg", "ab");
 176         if( giCurrentTimestamp ) {
 177         //      fseek(gLogFile, giCurrentTimestamp*giNumParticles*(sizeof(tLargeInt)*2+sizeof(tVector)), SEEK_SET);
 178         }
 179         else {
 180                 tmp = sizeof(tLargeInt);        fwrite(&tmp, 2, 1, gLogFile);
 181                 tmp = sizeof(tVector)/sizeof(tLargeInt);        fwrite(&tmp, 2, 1, gLogFile);
 182                 tmp = giNumParticles;   fwrite(&tmp, 4, 1, gLogFile);
 183         }
 184
 185         giUpdatePos = giNumParticles;
 186         for(;;giCurrentTimestamp++)
 187         {
 188                 printf("Iteration %i\n", giCurrentTimestamp);
 189                 Update_DoTimestep();
 190                 Visualiser_Update();
 191                 for(i=0,giUpdatePos=0;i<giNumParticles;i++,giUpdatePos++)
 192                 {
 193                         fwrite(&gaParticles[i].Location, sizeof(tVector), 1, gLogFile);
 194                         fwrite(&gaParticles[i].Mass, sizeof(tLargeInt), 1, gLogFile);
 195                         fwrite(&gaParticles[i].Charge, sizeof(tLargeInt), 1, gLogFile);
 196                 }
 197         }
 198
 199         return 0;
 200 }
 201
 202 void SaveUniverse(int signal)
 203 {
 204         // Detect an interrupt during update
 205         if(giUpdatePos < giNumParticles)        return;
 206         fclose(gLogFile);
 207         exit(0);
 208 }