/*----------------------------------------------------------------------*/
/* BIBLIOGRAPHY */
/* Dole, Stephen H. "Formation of Planetary Systems by Aggregation: */
/* a Computer Simulation" October 1969, Rand Corporation Paper */
/* P-4226. */
/*----------------------------------------------------------------------*/
/* A few variables global to the entire program: */
#include "structs.h"
/* 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 (double inner_limit_of_dust, double outer_limit_of_dust)
{
dust_head = (dust_bands *)malloc(sizeof(dust_bands));
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 (double stellar_mass_ratio)
{
return(200.0 * pow(stellar_mass_ratio,(1.0 / 3.0)));
}
double
innermost_planet (double stellar_mass_ratio)
{
return(0.3 * pow(stellar_mass_ratio,(1.0 / 3.0)));
}
double
outermost_planet (double stellar_mass_ratio)
{
return(50.0 * pow(stellar_mass_ratio,(1.0 / 3.0)));
}
double
inner_effect_limit (double a, double e, double mass)
{
return (a * (1.0 - e) * (1.0 - mass) / (1.0 + cloud_eccentricity));
}
double
outer_effect_limit (double a, double e, double mass)
{
return (a * (1.0 + e) * (1.0 + reduced_mass) / (1.0 - cloud_eccentricity));
}
int
dust_available (double inside_range, double 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 (double min, double max, double mass, double crit_mass, double body_inner_bound, double 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_bands *)malloc(sizeof(dust_bands));
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_bands *)malloc(sizeof(dust_bands));
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_bands *)malloc(sizeof(dust_bands));
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_bands *)malloc(sizeof(dust_bands));
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 (double orbital_radius, double eccentricity, double 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 (double *seed_mass, double a, double e, double crit_mass, double body_inner_bound, double 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 (double a, double e, double mass, double crit_mass, double stellar_luminosity_ratio, double body_inner_bound, double body_outer_bound)
{
planet_pointer node1, node2, node3;
int coalesced;
double temp, dist1, dist2, a3;
coalesced = FALSE;
node1 = planet_head;
node2 = NULL;
node3 = NULL;
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))))
{
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 (double stellar_mass_ratio, double stellar_luminosity_ratio, double inner_dust, double 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;
#ifdef VERBOSE
printf("Checking %f AU.\n",a);
#endif
if (dust_available(inner_effect_limit(a, e, mass),
outer_effect_limit(a, e, mass))) {
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 printf(".. failed due to large neighbor.\n");
}
#ifdef VERBOSE
else printf(".. failed.\n");
#endif
}
return(planet_head);
}
syntax highlighted by Code2HTML, v. 0.9.1