accrete.c

/*----------------------------------------------------------------------*/
/*                           BIBLIOGRAPHY                               */
/*  Dole, Stephen H.  "Formation of Planetary Systems by Aggregation:   */
/*      a Computer Simulation"  October 1969,  Rand Corporation Paper   */
/*	P-4226.								*/
/*----------------------------------------------------------------------*/

#include	<stdio.h>
#include	<math.h>
#include	"const.h"
#include	"structs.h"

/* externals from C library not elsewhere declared: */
extern char *malloc();
extern void free();

extern double random_number();
extern double random_eccentricity();
extern int flag_verbose, flag_lisp;

/* A few variables global to the entire program:		*/
planet_pointer planet_head;

/* Now for some variables global to the accretion process:      */
int dust_left;
double r_inner, r_outer, reduced_mass, dust_density, cloud_eccentricity;
dust_pointer dust_head;


void set_initial_conditions(inner_limit_of_dust, outer_limit_of_dust)
double inner_limit_of_dust, outer_limit_of_dust;
{
     dust_head = (dust *)malloc(sizeof(dust));
     planet_head = NULL;
     dust_head->next_band = NULL;
     dust_head->outer_edge = outer_limit_of_dust;
     dust_head->inner_edge = inner_limit_of_dust;
     dust_head->dust_present = TRUE;
     dust_head->gas_present = TRUE;
     dust_left = TRUE;
     cloud_eccentricity = 0.2;
}

double stellar_dust_limit(stellar_mass_ratio)
double stellar_mass_ratio; 
{
       return(200.0 * pow(stellar_mass_ratio,(1.0 / 3.0)));
}

double innermost_planet(stellar_mass_ratio)
double stellar_mass_ratio; 
{
     return(0.3 * pow(stellar_mass_ratio,(1.0 / 3.0)));
}

double outermost_planet(stellar_mass_ratio)
double stellar_mass_ratio; 
{
     return(50.0 * pow(stellar_mass_ratio,(1.0 / 3.0)));
}

double inner_effect_limit(a, e, mass)
double a, e, mass;
{
     return (a * (1.0 - e) * (1.0 - mass) / (1.0 + cloud_eccentricity));
}

/* In the original accrete.c, the return statement uses reduced_mass, not mass */
double outer_effect_limit(a, e, mass)
double a, e, mass;
{
     return (a * (1.0 + e) * (1.0 + mass) / (1.0 - cloud_eccentricity));
}

int dust_available(inside_range, outside_range)
double inside_range, outside_range;
{
     dust_pointer current_dust_band;
     int dust_here;
     
     current_dust_band = dust_head;
     while ((current_dust_band != NULL)
	    && (current_dust_band->outer_edge < inside_range))
	  current_dust_band = current_dust_band->next_band;
     if (current_dust_band == NULL)
	  dust_here = FALSE;
     else dust_here = current_dust_band->dust_present;
     while ((current_dust_band != NULL)
	    && (current_dust_band->inner_edge < outside_range)) {
	       dust_here = dust_here || current_dust_band->dust_present;
	       current_dust_band = current_dust_band->next_band;
	  }
     return(dust_here);
}

void update_dust_lanes(min, max, mass, crit_mass,
		       body_inner_bound, body_outer_bound)
double min, max, mass, crit_mass, body_inner_bound, body_outer_bound;
{
     int gas; 
     dust_pointer node1, node2, node3;
     
     dust_left = FALSE;
     if ((mass > crit_mass))
	  gas = FALSE;
     else 
	  gas = TRUE;
     node1 = dust_head;
     while ((node1 != NULL))
     {
	  if (((node1->inner_edge < min) && (node1->outer_edge > max)))
	  {
	       node2 = (dust *)malloc(sizeof(dust));
	       node2->inner_edge = min;
	       node2->outer_edge = max;
	       if ((node1->gas_present == TRUE))
		    node2->gas_present = gas;
	       else 
		    node2->gas_present = FALSE;
	       node2->dust_present = FALSE;
	       node3 = (dust *)malloc(sizeof(dust));
	       node3->inner_edge = max;
	       node3->outer_edge = node1->outer_edge;
	       node3->gas_present = node1->gas_present;
	       node3->dust_present = node1->dust_present;
	       node3->next_band = node1->next_band;
	       node1->next_band = node2;
	       node2->next_band = node3;
	       node1->outer_edge = min;
	       node1 = node3->next_band;
	  }
	  else 
	       if (((node1->inner_edge < max) && (node1->outer_edge > max)))
	       {
		    node2 = (dust *)malloc(sizeof(dust));
		    node2->next_band = node1->next_band;
		    node2->dust_present = node1->dust_present;
		    node2->gas_present = node1->gas_present;
		    node2->outer_edge = node1->outer_edge;
		    node2->inner_edge = max;
		    node1->next_band = node2;
		    node1->outer_edge = max;
		    if ((node1->gas_present == TRUE))
			 node1->gas_present = gas;
		    else 
			 node1->gas_present = FALSE;
		    node1->dust_present = FALSE;
		    node1 = node2->next_band;
	       }
	       else 
		    if (((node1->inner_edge < min) && (node1->outer_edge > min)))
		    {
			 node2 = (dust *)malloc(sizeof(dust));
			 node2->next_band = node1->next_band;
			 node2->dust_present = FALSE;
			 if ((node1->gas_present == TRUE))
			      node2->gas_present = gas;
			 else 
			      node2->gas_present = FALSE;
			 node2->outer_edge = node1->outer_edge;
			 node2->inner_edge = min;
			 node1->next_band = node2;
			 node1->outer_edge = min;
			 node1 = node2->next_band;
		    }
		    else 
			 if (((node1->inner_edge >= min) && (node1->outer_edge <= max)))
			 {
			      if ((node1->gas_present == TRUE))
				   node1->gas_present = gas;
			      node1->dust_present = FALSE;
			      node1 = node1->next_band;
			 }
			 else 
			      if (((node1->outer_edge < min) || (node1->inner_edge > max)))
				   node1 = node1->next_band;
     }
     node1 = dust_head;
     while ((node1 != NULL))
     {
	  if (((node1->dust_present)
	       && (((node1->outer_edge >= body_inner_bound)
		    && (node1->inner_edge <= body_outer_bound)))))
	       dust_left = TRUE;
	  node2 = node1->next_band;
	  if ((node2 != NULL))
	  {
	       if (((node1->dust_present == node2->dust_present)
		    && (node1->gas_present == node2->gas_present)))
	       {
		    node1->outer_edge = node2->outer_edge;
		    node1->next_band = node2->next_band;
		    free(node2);
	       }
	  }
	  node1 = node1->next_band;
     }
}

double collect_dust(last_mass, a, e, crit_mass, dust_band)
double last_mass, a, e, crit_mass;
dust_pointer dust_band;
{
     double mass_density, temp1, temp2, temp, temp_density, bandwidth, width, volume;
     
     temp = last_mass / (1.0 + last_mass);
     reduced_mass = pow(temp,(1.0 / 4.0));
     r_inner = inner_effect_limit(a, e, reduced_mass);
     r_outer = outer_effect_limit(a, e, reduced_mass);
     if ((r_inner < 0.0))
	  r_inner = 0.0;
     if ((dust_band == NULL))
	  return(0.0);
     else 
     {
	  if ((dust_band->dust_present == FALSE))
	       temp_density = 0.0;
	  else 
	       temp_density = dust_density;
	  if (((last_mass < crit_mass) || (dust_band->gas_present == FALSE)))
	       mass_density = temp_density;
	  else 
	       mass_density = K * temp_density / (1.0 + sqrt(crit_mass / last_mass)
						  * (K - 1.0));
	  if (((dust_band->outer_edge <= r_inner)
	       || (dust_band->inner_edge >= r_outer)))
	       return(collect_dust(last_mass,a,e,crit_mass, dust_band->next_band));
	  else
	  {
	       bandwidth = (r_outer - r_inner);
	       temp1 = r_outer - dust_band->outer_edge;
	       if (temp1 < 0.0)
		    temp1 = 0.0;
	       width = bandwidth - temp1;
	       temp2 = dust_band->inner_edge - r_inner;
	       if (temp2 < 0.0)
		    temp2 = 0.0;
	       width = width - temp2;
	       temp = 4.0 * PI * pow(a,2.0) * reduced_mass
		    * (1.0 - e * (temp1 - temp2) / bandwidth);
	       volume = temp * width;
	       return(volume * mass_density
		      + collect_dust(last_mass,a,e,crit_mass,
				     dust_band->next_band));
	  }
     }
}


/*--------------------------------------------------------------------------*/
/*   Orbital radius is in AU, eccentricity is unitless, and the stellar     */
/*  luminosity ratio is with respect to the sun.  The value returned is the */
/*  mass at which the planet begins to accrete gas as well as dust, and is  */
/*  in units of solar masses.                                               */
/*--------------------------------------------------------------------------*/

double critical_limit(orbital_radius, eccentricity, stellar_luminosity_ratio)
double orbital_radius, eccentricity, stellar_luminosity_ratio;
{
     double temp, perihelion_dist;
     
     perihelion_dist = (orbital_radius - orbital_radius * eccentricity);
     temp = perihelion_dist * sqrt(stellar_luminosity_ratio);
     return(B * pow(temp,-0.75));
}



void accrete_dust(seed_mass, a, e, crit_mass,
		  body_inner_bound, body_outer_bound)
double *seed_mass, a, e, crit_mass,
     body_inner_bound, body_outer_bound;
{
     double new_mass, temp_mass;
     
     new_mass = (*seed_mass);
     do
     {
	  temp_mass = new_mass;
	  new_mass = collect_dust(new_mass,a,e,crit_mass,
				  dust_head);
     }
     while (!(((new_mass - temp_mass) < (0.0001 * temp_mass))));
     (*seed_mass) = (*seed_mass) + new_mass;
     update_dust_lanes(r_inner,r_outer,(*seed_mass),crit_mass,body_inner_bound,body_outer_bound);
}



void coalesce_planetesimals(a, e, mass, crit_mass,
			    stellar_luminosity_ratio,
			    body_inner_bound, body_outer_bound)
double a, e, mass, crit_mass, stellar_luminosity_ratio,
     body_inner_bound, body_outer_bound;
{
     planet_pointer node1, node2, node3;
     int coalesced; 
     double temp, dist1, dist2, a3;
     
     coalesced = FALSE;
     node1 = planet_head;
     while ((node1 != NULL))
     {
	  node2 = node1;
	  temp = node1->a - a;
	  if ((temp > 0.0))
	  {
	       dist1 = (a * (1.0 + e) * (1.0 + reduced_mass)) - a;
	       /* x aphelion   */
	       reduced_mass = pow((node1->mass / (1.0 + node1->mass)),(1.0 / 4.0));
	       dist2 = node1->a
		    - (node1->a * (1.0 - node1->e) * (1.0 - reduced_mass));
	  }
	  else 
	  {
	       dist1 = a - (a * (1.0 - e) * (1.0 - reduced_mass));
	       /* x perihelion */
	       reduced_mass = pow(node1->mass / (1.0 + node1->mass),(1.0 / 4.0));
	       dist2 = (node1->a * (1.0 + node1->e) * (1.0 + reduced_mass))
		    - node1->a;
	  }
	  if (((fabs(temp) <= fabs(dist1)) || (fabs(temp) <= fabs(dist2))))
	  {
	       if (flag_verbose)
		       if (flag_lisp)
			       printf(";Collision between two planetesimals!\n");
		       else
			       printf("Collision between two planetesimals!\n");
	       a3 = (node1->mass + mass) / ((node1->mass / node1->a) + (mass / a));
	       temp = node1->mass * sqrt(node1->a) * sqrt(1.0 - pow(node1->e,2.0));
	       temp = temp + (mass * sqrt(a) * sqrt(sqrt(1.0 - pow(e,2.0))));
	       temp = temp / ((node1->mass + mass) * sqrt(a3));
	       temp = 1.0 - pow(temp,2.0);
	       if (((temp < 0.0) || (temp >= 1.0)))
		    temp = 0.0;
	       e = sqrt(temp);
	       temp = node1->mass + mass;
	       accrete_dust(&(temp),a3,e,stellar_luminosity_ratio,
			    body_inner_bound,body_outer_bound);
	       node1->a = a3;
	       node1->e = e;
	       node1->mass = temp;
	       node1 = NULL;
	       coalesced = TRUE;
	  }
	  else 
	       node1 = node1->next_planet;
     }
     if (!(coalesced))
     {
	  node3 = (planets *)malloc(sizeof(planets));
	  node3->a = a;
	  node3->e = e;
	  if ((mass >= crit_mass))
	       node3->gas_giant = TRUE;
	  else 
	       node3->gas_giant = FALSE;
	  node3->mass = mass;
	  if ((planet_head == NULL))
	  {
	       planet_head = node3;
	       node3->next_planet = NULL;
	  }
	  else 
	  {
	       node1 = planet_head;
	       if ((a < node1->a))
	       {
		    node3->next_planet = node1;
		    planet_head = node3;
	       }
	       else 
		    if ((planet_head->next_planet == NULL))
		    {
			 planet_head->next_planet = node3;
			 node3->next_planet = NULL;
		    }
		    else 
		    {
			 while (((node1 != NULL) && (node1->a < a)))
			 {
			      node2 = node1;
			      node1 = node1->next_planet;
			 }
			 node3->next_planet = node1;
			 node2->next_planet = node3;
		    }
	  }
     }
}


planet_pointer distribute_planetary_masses(stellar_mass_ratio,
          stellar_luminosity_ratio, inner_dust, outer_dust)
double stellar_mass_ratio, stellar_luminosity_ratio, inner_dust, outer_dust;
{
     double a, e, mass, crit_mass,
     planetesimal_inner_bound, planetesimal_outer_bound;
     
     set_initial_conditions(inner_dust,outer_dust);
     planetesimal_inner_bound = innermost_planet(stellar_mass_ratio);
     planetesimal_outer_bound = outermost_planet(stellar_mass_ratio);
     while (dust_left)
     {
          a = random_number(planetesimal_inner_bound,planetesimal_outer_bound);
          e = random_eccentricity( );
          mass = PROTOPLANET_MASS;
          if (flag_verbose)
               if (flag_lisp) printf(";Checking %lg AU.\n",a);
               else printf("Checking %lg AU.\n",a);
          if (dust_available(inner_effect_limit(a, e, mass),
                             outer_effect_limit(a, e, mass)))
          {
               if (flag_verbose)
                    if (flag_lisp) printf(";.. Injecting protoplanet.\n");
                    else printf(".. Injecting protoplanet.\n");
               dust_density = DUST_DENSITY_COEFF * sqrt(stellar_mass_ratio)
                         * exp(-ALPHA * pow(a,(1.0 / N)));
               crit_mass = critical_limit(a,e,stellar_luminosity_ratio);
               accrete_dust(&(mass),a,e,crit_mass, planetesimal_inner_bound,
                         planetesimal_outer_bound);
               if ((mass != 0.0) && (mass != PROTOPLANET_MASS))
                    coalesce_planetesimals(a,e,mass,crit_mass, stellar_luminosity_ratio,
                              planetesimal_inner_bound,planetesimal_outer_bound);
               else if (flag_verbose)
                    if (flag_lisp) printf(";.. failed due to large neighbor.\n");
                    else printf(".. failed due to large neighbor.\n");
          } else if (flag_verbose)
               if (flag_lisp) printf(";.. failed, no dust available.\n");
               else printf(".. failed, no dust available.\n");
     }
     return(planet_head);
}



#ifdef	MOON
/*
 * distribute_moon_masses does the same thing to a 
 * planetary system as distribute_planetary_masses
 * does to the whole solar system.
 */
planet_pointer distribute_moon_masses(planetary_mass, stellar_luminosity_ratio,
                                      planet_eccentricity, inner_dust, outer_dust)
double planetary_mass, stellar_luminosity_ratio, planet_eccentricity,
     inner_dust, outer_dust;
{
     double a, e, mass, crit_mass, reduced_mass, cloud_eccentricity,
        planetesimal_inner_bound, planetesimal_outer_bound;
     
     set_initial_conditions(inner_dust, outer_dust);

     while (dust_left) {

          /* it appears that set_initial_conditions sets this */
          cloud_eccentricity = 0.2;
          /* should these next 2 lines be before while dust_left? */
          planetesimal_inner_bound = innermost_planet(planetary_mass);
          planetesimal_outer_bound = outermost_planet(planetary_mass);
          /* what is protomoon_mass? */
          mass = PROTOMOON_MASS;
          a = random_number(planetesimal_inner_bound,planetesimal_outer_bound);
          e = random_eccentricity;
          crit_mass = critical_limit(a, planet_eccentricity,
               stellar_luminosity_ratio);
          accrete_dust (mass, reduced_mass, a, e, crit_mass, stellar_mass_ratio,
               cloud_eccentricity, planetesimal_inner_bound,
               planetesimal_outer_bound);
          /* original code was "mass <> 0.0" and "mass <> protoplanet_mass" */
          if ( mass != 0.0 && mass != PROTOMOON_MASS)
               coalesce_planetesimals (reduced_mass, a, e, mass, crit_mass,
                    planetary_mass, stellar_luminosity_ratio, 
                    cloud_eccentricity, planetesimal_inner_bound, 
                    planetesimal_outer_bound);
     }

     return(planet_head);
}

double planet_dust_limit(planetary_mass)
double planetary_mass;
{
       /* this is what stellar_dust_limit returns */
       /* return(200.0 * pow(stellar_mass_ratio,(1.0 / 3.0))); */
       /* Since I cannot find the original formula, I guess: */
       return(0.1 * pow(planetary_mass,(1.0 / 3.0)));
}
#endif	/* MOON */