SHARE
TWEET

A simple (and crappy) C++ program.

a guest Sep 9th, 2019 161 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. /*
  2. A simple program to calculate the solid angle of all stars in the observable universe.
  3.  
  4. It then compares this to the solid angle of the sky (4 pi steradians) to get the odds
  5. that any position on Earth is directly under a star.
  6.  
  7. Created in association with an answer on astronomy.stackexchange:
  8. https://astronomy.stackexchange.com/a/33300/10678
  9.  
  10. Written in Visual Studio 2017 on Windows 8.1. Compiled in 64-bit mode.
  11.  
  12. I didn't try to optimize it, error-proof it, or make it look pretty. The universe calculation took about 4 minutes in debug mode on an i5 6600k. The galaxy calculation just took a few seconds.
  13.  
  14. *** REQUIRES THIRD-PARTY LIBRARY ***
  15. https://www.ttmath.org/
  16. Make sure to assemble and include the .asm -> .obj file.
  17. Extra instructions can be found at:
  18. https://stackoverflow.com/a/29732728/5313933
  19.  
  20. You'll also need to add the directory with "ttmath.h" to your includes in whatever IDE you're using.
  21. */
  22.  
  23. #include <stdio.h>
  24. #include <math.h>
  25. #include <ttmath.h>
  26. using namespace ttmath;
  27.  
  28. // Comment this line to run whole universe calculations. Uncomment to run galaxy calculations.
  29. //#define MILKY_WAY
  30.  
  31. int main()
  32. {
  33.     const long double PI = 3.1415926535897932384626433832795l;
  34.     const int bigPrecision = 10;
  35.  
  36.     // Use this code for the whole universe.
  37. #ifndef MILKY_WAY
  38.     // Setup the character values of the bignum values so we can be precise with the number of digits. For a*10^b, set i = b+1.
  39.     int i = 12; // 2.213*10^11
  40.     char* cLeadingS = new char[i + 1];
  41.     for (int j = 0; j < i; j++)
  42.     {
  43.         cLeadingS[j] = '0';
  44.     }
  45.     cLeadingS[0] = '2'; // Set specific digits here.
  46.     cLeadingS[1] = '2';
  47.     cLeadingS[2] = '1';
  48.     cLeadingS[3] = '3';
  49.     cLeadingS[i] = 0; // Make sure the string is null-terminated.
  50.  
  51.     i = 45; // 10^44
  52.     char* cRadLeftS = new char[i + 1]();
  53.     for (int j = 0; j < i; j++)
  54.     {
  55.         cRadLeftS[j] = '0';
  56.     }
  57.     cRadLeftS[0] = '1';
  58.     cRadLeftS[i] = 0;
  59.  
  60.     i = 18; // 4.9*10^17
  61.     char* cRadRightS = new char[i + 1]();
  62.     for (int j = 0; j < i; j++)
  63.     {
  64.         cRadRightS[j] = '0';
  65.     }
  66.     cRadRightS[0] = '4';
  67.     cRadRightS[1] = '9';
  68.     cRadRightS[i] = 0;
  69.  
  70.     i = 23; // 10^22
  71.     char* cDenomS = new char[i + 1]();
  72.     for (int j = 0; j < i; j++)
  73.     {
  74.         cDenomS[j] = '0';
  75.     }
  76.     cDenomS[0] = '1';
  77.     cDenomS[i] = 0;
  78. #endif // !MILKY_WAY
  79.  
  80.     // Use this code for the Milky Way galaxy.
  81. #ifdef MILKY_WAY
  82.     // Setup the character values of the bignum values so we can be precise with the number of digits. For a*10^b, set i = b+1.
  83.     int i = 7; // 1.005*10^6
  84.     char* cLeadingS = new char[i + 1];
  85.     for (int j = 0; j < i; j++)
  86.     {
  87.         cLeadingS[j] = '0';
  88.     }
  89.     cLeadingS[0] = '1'; // Set specific digits here.
  90.     cLeadingS[3] = '5';
  91.     cLeadingS[i] = 0; // Make sure the string is null-terminated.
  92.  
  93.     i = 35; // 10^34
  94.     char* cRadLeftS = new char[i + 1]();
  95.     for (int j = 0; j < i; j++)
  96.     {
  97.         cRadLeftS[j] = '0';
  98.     }
  99.     cRadLeftS[0] = '1';
  100.     cRadLeftS[i] = 0;
  101.  
  102.     i = 18; // 4.9*10^17
  103.     char* cRadRightS = new char[i + 1]();
  104.     for (int j = 0; j < i; j++)
  105.     {
  106.         cRadRightS[j] = '0';
  107.     }
  108.     cRadRightS[0] = '4';
  109.     cRadRightS[1] = '9';
  110.     cRadRightS[i] = 0;
  111.  
  112.     i = 18; // 10^17
  113.     char* cDenomS = new char[i + 1]();
  114.     for (int j = 0; j < i; j++)
  115.     {
  116.         cDenomS[j] = '0';
  117.     }
  118.     cDenomS[0] = '1';
  119.     cDenomS[i] = 0;
  120. #endif // MILKY_WAY
  121.  
  122.     // Setup the variables, using the character strings to instantiate non-zero values.
  123.     ttmath::Big<1, bigPrecision> omegaTotal = 0;
  124.     ttmath::Big<1, bigPrecision> omegaCurrent = 1;
  125.     ttmath::Big<1, bigPrecision> omegaPercent = 0;
  126.     ttmath::Big<1, bigPrecision> cLeading = cLeadingS;
  127.     ttmath::Big<1, bigPrecision> cRadLeft = cRadLeftS;
  128.     ttmath::Big<1, bigPrecision> cRadRight = cRadRightS;
  129.     ttmath::Big<1, bigPrecision> cDenom = cDenomS;
  130.     ttmath::Big<1, bigPrecision> One = 1;
  131.     ttmath::Big<1, bigPrecision> FourPi = (double)(4 * PI);
  132.     ttmath::Big<1, bigPrecision> Thousand = 1000;
  133.     ttmath::Big<1, bigPrecision> kMax;
  134.     ttmath::Big<1, bigPrecision> tempVal1;
  135.     ttmath::Big<1, bigPrecision> tempVal2;
  136.     ttmath::UInt<bigPrecision> oneInVal;
  137.  
  138.     // Setup the shell count depending on which calculation we're doing.
  139. #ifndef MILKY_WAY
  140.     kMax = 44000;
  141. #endif // !MILKY_WAY
  142. #ifdef MILKY_WAY
  143.     kMax = 155;
  144. #endif // MILKY_WAY
  145.  
  146.     /*
  147.     We're calculating omegaTotal = sum from k=1 to 44000 of OmegaCurrent, where
  148.  
  149.     omegaCurrent = (1 - omegaTotal/(4*pi)) * cLeading * k^2 *
  150.         (1 - [sqrt(k^2 * cRadLeft - cRadRight] / [k * cDenom])
  151.  
  152.     With the result being units of radians.
  153.     */
  154.     for (ttmath::Big<1, bigPrecision> k = 1; k <= kMax; k++)
  155.     {
  156.         tempVal1 = k * k * cRadLeft;
  157.         tempVal1 -= cRadRight;
  158.         tempVal1 = ttmath::Sqrt(tempVal1);
  159.         tempVal1 /= k;
  160.         tempVal1 /= cDenom;
  161.         tempVal1 = One - tempVal1;
  162.         tempVal1 *= k * k;
  163.         tempVal1 *= cLeading;
  164.  
  165.         tempVal2 = One / FourPi;
  166.         tempVal2 *= omegaTotal;
  167.         tempVal2 = One - tempVal2;
  168.  
  169.         tempVal1 *= tempVal2;
  170.  
  171.         omegaCurrent = tempVal1;
  172.         omegaTotal += omegaCurrent;
  173.  
  174.         // Tell the command window every 1000 iterations so we don't get bored.
  175.         if (k >= Thousand)
  176.         {
  177.             Thousand += 1000;
  178.             std::cout << "Shell number " << k << " of 44000. omegaTotal is currently \n" << omegaTotal << "\n";
  179.         }
  180.     }
  181.  
  182.     tempVal1 = (double)(100 / (4 * PI));
  183.     tempVal2 = (double)(1 / (4 * PI));
  184.     omegaPercent = omegaTotal * tempVal1;
  185.     tempVal1 = One / (omegaTotal * tempVal2);
  186.  
  187.     // Output message should change depending on which math we're doing.
  188. #ifndef MILKY_WAY
  189.     std::cout << "The total solid angle subtended by all the stars in the observable universe is:\n";
  190. #endif // !MILKY_WAY
  191. #ifdef MILKY_WAY
  192.     std::cout << "The total solid angle subtended by all the stars in the Milky Way is:\n";
  193. #endif // MILKY_WAY
  194.  
  195.     std::cout << omegaTotal << " steradians, which equates to " << omegaPercent << "%.\n";
  196.     std::cout << "The odds of you standing under a star right now are one in " << tempVal1 << ".\n";
  197.  
  198.     return 0;
  199. }
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top