home *** CD-ROM | disk | FTP | other *** search
/ Gold Fish 2 / goldfish_vol2_cd1.bin / files / misc / sci / accrete / src / accrete.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-02-01  |  13.3 KB  |  455 lines
  1.  
  2. /*----------------------------------------------------------------------*/
  3. /*                           BIBLIOGRAPHY                               */
  4. /*  Dole, Stephen H.  "Formation of Planetary Systems by Aggregation:   */
  5. /*      a Computer Simulation"  October 1969,  Rand Corporation Paper   */
  6. /*    P-4226.                                */
  7. /*----------------------------------------------------------------------*/
  8.  
  9. #include <stdio.h>
  10. #include <string.h>
  11. #include <math.h>
  12. #include <stddef.h>
  13. #include <stdlib.h>
  14. #ifdef AMIGA
  15. #ifdef HUGE_VAL
  16. #undef HUGE_VAL
  17. #endif /* HUGE_VAL */
  18. #endif /* AMIGA */
  19. #include <float.h>
  20. #include <time.h>
  21.  
  22. #include "const.h"
  23. #include "structs.h"
  24. #include "accrete.h"
  25.  
  26.  
  27. /* Variables global to the accretion process:      */
  28. int dust_left;
  29. double r_inner, r_outer, reduced_mass, dust_density, cloud_eccentricity;
  30. dust_pointer dust_head;
  31.  
  32.  
  33. void set_initial_conditions(inner_limit_of_dust, outer_limit_of_dust)
  34. double inner_limit_of_dust, outer_limit_of_dust;
  35. {
  36.      dust_head = (dust_bands *)malloc(sizeof(dust_bands));
  37.      planet_head = NULL;
  38.      dust_head->next_band = NULL;
  39.      dust_head->outer_edge = outer_limit_of_dust;
  40.      dust_head->inner_edge = inner_limit_of_dust;
  41.      dust_head->dust_present = TRUE;
  42.      dust_head->gas_present = TRUE;
  43.      dust_left = TRUE;
  44.      cloud_eccentricity = 0.2;
  45. }
  46.  
  47. double stellar_dust_limit(stellar_mass_ratio)
  48. double stellar_mass_ratio; 
  49. {
  50.      return(200.0 * pow(stellar_mass_ratio,(1.0 / 3.0)));
  51. }
  52.  
  53. double innermost_planet(stellar_mass_ratio)
  54. double stellar_mass_ratio; 
  55. {
  56.      return(0.3 * pow(stellar_mass_ratio,(1.0 / 3.0)));
  57. }
  58.  
  59. double outermost_planet(stellar_mass_ratio)
  60. double stellar_mass_ratio; 
  61. {
  62.      return(50.0 * pow(stellar_mass_ratio,(1.0 / 3.0)));
  63. }
  64.  
  65. double inner_effect_limit(a, e, mass)
  66. double a, e, mass;
  67. {
  68.      return (a * (1.0 - e) * (1.0 - mass) / (1.0 + cloud_eccentricity));
  69. }
  70.  
  71. double outer_effect_limit(a, e, mass)
  72. double a, e, mass;
  73. {
  74.      return (a * (1.0 + e) * (1.0 + reduced_mass) / (1.0 - cloud_eccentricity));
  75. }
  76.  
  77. int dust_available(inside_range, outside_range)
  78. double inside_range, outside_range;
  79. {
  80.      dust_pointer current_dust_band;
  81.      int dust_here;
  82.      
  83.      current_dust_band = dust_head;
  84.      while ((current_dust_band != NULL)
  85.         && (current_dust_band->outer_edge < inside_range))
  86.       current_dust_band = current_dust_band->next_band;
  87.      if (current_dust_band == NULL)
  88.       dust_here = FALSE;
  89.      else dust_here = current_dust_band->dust_present;
  90.      while ((current_dust_band != NULL)
  91.         && (current_dust_band->inner_edge < outside_range)) {
  92.            dust_here = dust_here || current_dust_band->dust_present;
  93.            current_dust_band = current_dust_band->next_band;
  94.       }
  95.      return(dust_here);
  96. }
  97.  
  98. void update_dust_lanes(min, max, mass, crit_mass,
  99.                body_inner_bound, body_outer_bound)
  100. double min, max, mass, crit_mass, body_inner_bound, body_outer_bound;
  101. {
  102.      int gas; 
  103.      dust_pointer node1, node2, node3;
  104.      
  105.      dust_left = FALSE;
  106.      if ((mass > crit_mass))
  107.       gas = FALSE;
  108.      else 
  109.       gas = TRUE;
  110.      node1 = dust_head;
  111.      while ((node1 != NULL))
  112.      {
  113.       if (((node1->inner_edge < min) && (node1->outer_edge > max)))
  114.       {
  115.            node2 = (dust_bands *)malloc(sizeof(dust_bands));
  116.            node2->inner_edge = min;
  117.            node2->outer_edge = max;
  118.            if ((node1->gas_present == TRUE))
  119.             node2->gas_present = gas;
  120.            else 
  121.             node2->gas_present = FALSE;
  122.            node2->dust_present = FALSE;
  123.            node3 = (dust_bands *)malloc(sizeof(dust_bands));
  124.            node3->inner_edge = max;
  125.            node3->outer_edge = node1->outer_edge;
  126.            node3->gas_present = node1->gas_present;
  127.            node3->dust_present = node1->dust_present;
  128.            node3->next_band = node1->next_band;
  129.            node1->next_band = node2;
  130.            node2->next_band = node3;
  131.            node1->outer_edge = min;
  132.            node1 = node3->next_band;
  133.       }
  134.       else 
  135.            if (((node1->inner_edge < max) && (node1->outer_edge > max)))
  136.            {
  137.             node2 = (dust_bands *)malloc(sizeof(dust_bands));
  138.             node2->next_band = node1->next_band;
  139.             node2->dust_present = node1->dust_present;
  140.             node2->gas_present = node1->gas_present;
  141.             node2->outer_edge = node1->outer_edge;
  142.             node2->inner_edge = max;
  143.             node1->next_band = node2;
  144.             node1->outer_edge = max;
  145.             if ((node1->gas_present == TRUE))
  146.              node1->gas_present = gas;
  147.             else 
  148.              node1->gas_present = FALSE;
  149.             node1->dust_present = FALSE;
  150.             node1 = node2->next_band;
  151.            }
  152.            else 
  153.             if (((node1->inner_edge < min) && (node1->outer_edge > min)))
  154.             {
  155.              node2 = (dust_bands *)malloc(sizeof(dust_bands));
  156.              node2->next_band = node1->next_band;
  157.              node2->dust_present = FALSE;
  158.              if ((node1->gas_present == TRUE))
  159.                   node2->gas_present = gas;
  160.              else 
  161.                   node2->gas_present = FALSE;
  162.              node2->outer_edge = node1->outer_edge;
  163.              node2->inner_edge = min;
  164.              node1->next_band = node2;
  165.              node1->outer_edge = min;
  166.              node1 = node2->next_band;
  167.             }
  168.             else 
  169.              if (((node1->inner_edge >= min) && (node1->outer_edge <= max)))
  170.              {
  171.                   if ((node1->gas_present == TRUE))
  172.                    node1->gas_present = gas;
  173.                   node1->dust_present = FALSE;
  174.                   node1 = node1->next_band;
  175.              }
  176.              else 
  177.                   if (((node1->outer_edge < min) || (node1->inner_edge > max)))
  178.                    node1 = node1->next_band;
  179.      }
  180.      node1 = dust_head;
  181.      while ((node1 != NULL))
  182.      {
  183.       if (((node1->dust_present)
  184.            && (((node1->outer_edge >= body_inner_bound)
  185.             && (node1->inner_edge <= body_outer_bound)))))
  186.            dust_left = TRUE;
  187.       node2 = node1->next_band;
  188.       if ((node2 != NULL))
  189.       {
  190.            if (((node1->dust_present == node2->dust_present)
  191.             && (node1->gas_present == node2->gas_present)))
  192.            {
  193.             node1->outer_edge = node2->outer_edge;
  194.             node1->next_band = node2->next_band;
  195.             free(node2);
  196.            }
  197.       }
  198.       node1 = node1->next_band;
  199.      }
  200. }
  201.  
  202. double collect_dust(last_mass, a, e, crit_mass, dust_band)
  203. double last_mass, a, e, crit_mass;
  204. dust_pointer dust_band;
  205. {
  206.      double mass_density, temp1, temp2, temp, temp_density, bandwidth, width, volume;
  207.      
  208.      temp = last_mass / (1.0 + last_mass);
  209.      reduced_mass = pow(temp,(1.0 / 4.0));
  210.      r_inner = inner_effect_limit(a, e, reduced_mass);
  211.      r_outer = outer_effect_limit(a, e, reduced_mass);
  212.      if ((r_inner < 0.0))
  213.       r_inner = 0.0;
  214.      if ((dust_band == NULL))
  215.       return(0.0);
  216.      else 
  217.      {
  218.       if ((dust_band->dust_present == FALSE))
  219.            temp_density = 0.0;
  220.       else 
  221.            temp_density = dust_density;
  222.       if (((last_mass < crit_mass) || (dust_band->gas_present == FALSE)))
  223.            mass_density = temp_density;
  224.       else 
  225.            mass_density = K * temp_density / (1.0 + sqrt(crit_mass / last_mass)
  226.                           * (K - 1.0));
  227.       if (((dust_band->outer_edge <= r_inner)
  228.            || (dust_band->inner_edge >= r_outer)))
  229.            return(collect_dust(last_mass,a,e,crit_mass, dust_band->next_band));
  230.       else
  231.       {
  232.            bandwidth = (r_outer - r_inner);
  233.            temp1 = r_outer - dust_band->outer_edge;
  234.            if (temp1 < 0.0)
  235.             temp1 = 0.0;
  236.            width = bandwidth - temp1;
  237.            temp2 = dust_band->inner_edge - r_inner;
  238.            if (temp2 < 0.0)
  239.             temp2 = 0.0;
  240.            width = width - temp2;
  241.            temp = 4.0 * PI * pow(a,2.0) * reduced_mass
  242.             * (1.0 - e * (temp1 - temp2) / bandwidth);
  243.            volume = temp * width;
  244.            return(volume * mass_density
  245.               + collect_dust(last_mass,a,e,crit_mass,
  246.                      dust_band->next_band));
  247.       }
  248.      }
  249. }
  250.  
  251.  
  252. /*--------------------------------------------------------------------------*/
  253. /*   Orbital radius is in AU, eccentricity is unitless, and the stellar     */
  254. /*  luminosity ratio is with respect to the sun.  The value returned is the */
  255. /*  mass at which the planet begins to accrete gas as well as dust, and is  */
  256. /*  in units of solar masses.                                               */
  257. /*--------------------------------------------------------------------------*/
  258.  
  259. double critical_limit(orbital_radius, eccentricity, stellar_luminosity_ratio)
  260. double orbital_radius, eccentricity, stellar_luminosity_ratio;
  261. {
  262.      double temp, perihelion_dist;
  263.      
  264.      perihelion_dist = (orbital_radius - orbital_radius * eccentricity);
  265.      temp = perihelion_dist * sqrt(stellar_luminosity_ratio);
  266.      return(B * pow(temp,-0.75));
  267. }
  268.  
  269.  
  270.  
  271. void accrete_dust(seed_mass, a, e, crit_mass,
  272.           body_inner_bound, body_outer_bound)
  273. double *seed_mass, a, e, crit_mass,
  274.      body_inner_bound, body_outer_bound;
  275. {
  276.      double perihelion_dist, new_mass, temp_mass;
  277.      
  278.      new_mass = (*seed_mass);
  279.      do
  280.      {
  281.       temp_mass = new_mass;
  282.       new_mass = collect_dust(new_mass,a,e,crit_mass,
  283.                   dust_head);
  284.      }
  285.      while (!(((new_mass - temp_mass) < (0.0001 * temp_mass))));
  286.      (*seed_mass) = (*seed_mass) + new_mass;
  287.      update_dust_lanes(r_inner,r_outer,(*seed_mass),crit_mass,body_inner_bound,body_outer_bound);
  288. }
  289.  
  290.  
  291.  
  292. void coalesce_planetesimals(a, e, mass, crit_mass,
  293.                 stellar_luminosity_ratio,
  294.                 body_inner_bound, body_outer_bound)
  295. double a, e, mass, crit_mass, stellar_luminosity_ratio,
  296.      body_inner_bound, body_outer_bound;
  297. {
  298.      planet_pointer node1, node2, node3;
  299.      int coalesced; 
  300.      double temp, dist1, dist2, a3;
  301.      
  302.      coalesced = FALSE;
  303.      node1 = planet_head;
  304.      while ((node1 != NULL))
  305.      {
  306.       node2 = node1;
  307.       temp = node1->a - a;
  308.       if ((temp > 0.0))
  309.       {
  310.            dist1 = (a * (1.0 + e) * (1.0 + reduced_mass)) - a;
  311.            /* x aphelion   */
  312.            reduced_mass = pow((node1->mass / (1.0 + node1->mass)),(1.0 / 4.0));
  313.            dist2 = node1->a
  314.             - (node1->a * (1.0 - node1->e) * (1.0 - reduced_mass));
  315.       }
  316.       else 
  317.       {
  318.            dist1 = a - (a * (1.0 - e) * (1.0 - reduced_mass));
  319.            /* x perihelion */
  320.            reduced_mass = pow(node1->mass / (1.0 + node1->mass),(1.0 / 4.0));
  321.            dist2 = (node1->a * (1.0 + node1->e) * (1.0 + reduced_mass))
  322.             - node1->a;
  323.       }
  324.       if (((fabs(temp) <= fabs(dist1)) || (fabs(temp) <= fabs(dist2))))
  325.       {
  326. #ifdef VERBOSE
  327.            printf("Collision between two planetesimals!\n");
  328. #endif
  329.            a3 = (node1->mass + mass) / ((node1->mass / node1->a) + (mass / a));
  330.            temp = node1->mass * sqrt(node1->a) * sqrt(1.0 - pow(node1->e,2.0));
  331.            temp = temp + (mass * sqrt(a) * sqrt(sqrt(1.0 - pow(e,2.0))));
  332.            temp = temp / ((node1->mass + mass) * sqrt(a3));
  333.            temp = 1.0 - pow(temp,2.0);
  334.            if (((temp < 0.0) || (temp >= 1.0)))
  335.             temp = 0.0;
  336.            e = sqrt(temp);
  337.            temp = node1->mass + mass;
  338.            accrete_dust(&(temp),a3,e,stellar_luminosity_ratio,
  339.                 body_inner_bound,body_outer_bound);
  340.            node1->a = a3;
  341.            node1->e = e;
  342.            node1->mass = temp;
  343.            node1 = NULL;
  344.            coalesced = TRUE;
  345.       }
  346.       else 
  347.            node1 = node1->next_planet;
  348.      }
  349.      if (!(coalesced))
  350.      {
  351.       node3 = (planets *)malloc(sizeof(planets));
  352.       node3->a = a;
  353.       node3->e = e;
  354.       if ((mass >= crit_mass))
  355.            node3->gas_giant = TRUE;
  356.       else 
  357.            node3->gas_giant = FALSE;
  358.       node3->mass = mass;
  359.       if ((planet_head == NULL))
  360.       {
  361.            planet_head = node3;
  362.            node3->next_planet = NULL;
  363.       }
  364.       else 
  365.       {
  366.            node1 = planet_head;
  367.            if ((a < node1->a))
  368.            {
  369.             node3->next_planet = node1;
  370.             planet_head = node3;
  371.            }
  372.            else 
  373.             if ((planet_head->next_planet == NULL))
  374.             {
  375.              planet_head->next_planet = node3;
  376.              node3->next_planet = NULL;
  377.             }
  378.             else 
  379.             {
  380.              while (((node1 != NULL) && (node1->a < a)))
  381.              {
  382.                   node2 = node1;
  383.                   node1 = node1->next_planet;
  384.              }
  385.              node3->next_planet = node1;
  386.              node2->next_planet = node3;
  387.             }
  388.       }
  389.      }
  390. }
  391.  
  392.  
  393. planet_pointer distribute_planetary_masses(stellar_mass_ratio,
  394.                        stellar_luminosity_ratio, inner_dust, outer_dust)
  395. double stellar_mass_ratio, stellar_luminosity_ratio, inner_dust, outer_dust;
  396. {
  397.      double a, e, mass, crit_mass,
  398.      planetesimal_inner_bound, planetesimal_outer_bound;
  399.      
  400.      set_initial_conditions(inner_dust,outer_dust);
  401.      planetesimal_inner_bound = innermost_planet(stellar_mass_ratio);
  402.      planetesimal_outer_bound = outermost_planet(stellar_mass_ratio);
  403.      while (dust_left)
  404.      {
  405.       a = random_number(planetesimal_inner_bound,planetesimal_outer_bound);
  406.       e = random_eccentricity( );
  407.       mass = PROTOPLANET_MASS;
  408. #ifdef VERBOSE
  409.       printf("Checking %f AU.\n", a);
  410. #else
  411. # ifndef SILENT
  412.       printf(".");
  413. # endif /* SILENT */
  414. #endif /* VERBOSE */
  415.  
  416.       if (dust_available(inner_effect_limit(a, e, mass),
  417.                  outer_effect_limit(a, e, mass))) {
  418. #ifdef VERBOSE
  419.             printf(".. Injecting protoplanet.\n");
  420. #endif
  421.             dust_density = DUST_DENSITY_COEFF * sqrt(stellar_mass_ratio)
  422.              * exp(-ALPHA * pow(a,(1.0 / N)));
  423.             crit_mass = critical_limit(a,e,stellar_luminosity_ratio);
  424.             accrete_dust(&(mass),a,e,crit_mass,
  425.                  planetesimal_inner_bound,
  426.                  planetesimal_outer_bound);
  427.             if ((mass != 0.0) && (mass != PROTOPLANET_MASS))
  428.              coalesce_planetesimals(a,e,mass,crit_mass,
  429.                         stellar_luminosity_ratio,
  430.                         planetesimal_inner_bound,planetesimal_outer_bound);
  431. #ifdef VERBOSE
  432.             else printf(".. failed due to large neighbor.\n");
  433. #endif
  434.            }
  435. #ifdef VERBOSE
  436.       else printf(".. failed.\n");
  437. #endif
  438.      }
  439.      return(planet_head);
  440. }
  441.  
  442.  
  443.  
  444. planet_pointer distribute_moon_masses(planetary_mass, stellar_luminosity_ratio,
  445.                       planet_eccentricity, inner_dust, outer_dust)
  446. double planetary_mass, stellar_luminosity_ratio, planet_eccentricity,
  447.      inner_dust, outer_dust;
  448. {
  449. /*
  450.      double a, e, mass, crit_mass,
  451.      planetesimal_inner_bound, planetesimal_outer_bound;
  452. */
  453.      return(NULL);
  454. }
  455.