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

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include “const.h”
#include “structs.h”
#include “accrete.h”
#include “stargen.h”
#include “utils.h”

#include “Dumas.h”

/ Now for some variables global to the accretion process: */
int dust_left;
long double r_inner;
long double r_outer;
long double reduced_mass;
long double dust_density;
long double cloud_eccentricity;
dust_pointer dust_head = NULL;
planet_pointer planet_head = NULL;
gen_pointer hist_head = NULL;

void set_initial_conditions(long double inner_limit_of_dust,
long double outer_limit_of_dust)
{
gen_pointer hist;
hist = (gen_pointer)malloc(sizeof(generation));
hist->dusts = dust_head;
hist->planets = planet_head;
hist->next = hist_head;
hist_head = hist;

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;
}

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

long double nearest_planet(long double stell_mass_ratio)
{
return(0.3 * pow(stell_mass_ratio,(1.0 / 3.0)));
}

long double farthest_planet(long double stell_mass_ratio)
{
return(50.0 * pow(stell_mass_ratio,(1.0 / 3.0)));
}

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

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

int dust_available(long double inside_range, long 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(long double min, long double max, long double mass,
long double crit_mass, long double body_inner_bound,
long double body_outer_bound)
{
int gas;
dust_pointer node1;
dust_pointer node2;
dust_pointer 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;
}
}

long double collect_dust(long double last_mass, long double *new_dust,
long double *new_gas,
long double a, long double e,
long double crit_mass, dust_pointer dust_band)
{
long double mass_density;
long double temp1;
long double temp2;
long double temp;
long double temp_density;
long double bandwidth;
long double width;
long double volume;
long double gas_density = 0.0;
long double new_mass;
long double next_mass;
long double next_dust = 0;
long double next_gas = 0;


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));
gas_density = mass_density - temp_density;
}

if (((dust_band->outer_edge <= r_inner)
|| (dust_band->inner_edge >= r_outer)))
{
return(collect_dust(last_mass, new_dust, new_gas,
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;

new_mass = volume * mass_density;
*new_gas = volume * gas_density;
*new_dust = new_mass - *new_gas;

next_mass = collect_dust(last_mass, &next_dust, &next_gas,
a,e,crit_mass, dust_band->next_band);

new_gas = new_gas + next_gas;
new_dust = new_dust + next_dust;

return(new_mass + next_mass);
}
}
}


/--------------------------------------------------------------------------/
/ 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. /
/--------------------------------------------------------------------------/

long double critical_limit(long double orb_radius, long double eccentricity,
long double stell_luminosity_ratio)
{
long double temp;
long double perihelion_dist;

perihelion_dist = (orb_radius - orb_radius * eccentricity);
temp = perihelion_dist * sqrt(stell_luminosity_ratio);
return(B * pow(temp,-0.75));
}



void accrete_dust(long double *seed_mass, long double *new_dust, long double *new_gas,
long double a, long double e, long double crit_mass,
long double body_inner_bound, long double body_outer_bound)
{
long double new_mass = (*seed_mass);
long double temp_mass;

do
{
temp_mass = new_mass;
new_mass = collect_dust(new_mass, new_dust, new_gas,
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(long double a, long double e, long double mass, long double crit_mass,
long double dust_mass, long double gas_mass,
long double stell_luminosity_ratio,
long double body_inner_bound, long double body_outer_bound,
int do_moons)
{
planet_pointer the_planet;
planet_pointer next_planet;
planet_pointer prev_planet;
int finished;
long double temp;
long double diff;
long double dist1;
long double dist2;

finished = FALSE;
prev_planet = NULL;

// First we try to find an existing planet with an over-lapping orbit.

for (the_planet = planet_head;
the_planet != NULL;
the_planet = the_planet->next_planet)
{
diff = the_planet->a - a;

if ((diff > 0.0))
{
dist1 = (a * (1.0 + e) * (1.0 + reduced_mass)) - a;
/ x aphelion /
reduced_mass = pow((the_planet->mass / (1.0 + the_planet->mass)),(1.0 / 4.0));
dist2 = the_planet->a
- (the_planet->a * (1.0 - the_planet->e) * (1.0 - reduced_mass));
}
else
{
dist1 = a - (a * (1.0 - e) * (1.0 - reduced_mass));
/ x perihelion */
reduced_mass = pow((the_planet->mass / (1.0 + the_planet->mass)),(1.0 / 4.0));
dist2 = (the_planet->a * (1.0 + the_planet->e) * (1.0 + reduced_mass))
- the_planet->a;
}

if (((fabs(diff) <= fabs(dist1)) || (fabs(diff) <= fabs(dist2))))
{
long double new_dust = 0;
long double new_gas = 0;
long double new_a = (the_planet->mass + mass) /
((the_planet->mass / the_planet->a) + (mass / a));

temp = the_planet->mass * sqrt(the_planet->a) * sqrt(1.0 - pow(the_planet->e,2.0));
temp = temp + (mass * sqrt(a) * sqrt(sqrt(1.0 - pow(e,2.0))));
temp = temp / ((the_planet->mass + mass) * sqrt(new_a));
temp = 1.0 - pow(temp,2.0);
if (((temp < 0.0) || (temp >= 1.0)))
temp = 0.0;
e = sqrt(temp);

if (do_moons)
{
long double existing_mass = 0.0;

if (the_planet->first_moon != NULL)
{
planet_pointer m;

for (m = the_planet->first_moon;
m != NULL;
m = m->next_planet)
{
existing_mass += m->mass;
}
}

if (mass < crit_mass)
{
if ((mass * SUN_MASS_IN_EARTH_MASSES) < 2.5
&& (mass * SUN_MASS_IN_EARTH_MASSES) > .0001
&& existing_mass < (the_planet->mass * .05)
)
{
planet_pointer the_moon = (planets )malloc(sizeof(planets));

the_moon->type = tUnknown;
/ the_moon->a = a; /
/ the_moon->e = e; */
the_moon->mass = mass;
the_moon->dust_mass = dust_mass;
the_moon->gas_mass = gas_mass;
the_moon->atmosphere = NULL;
the_moon->next_planet = NULL;
the_moon->first_moon = NULL;
the_moon->gas_giant = FALSE;
the_moon->atmosphere = NULL;
the_moon->albedo = 0;
the_moon->gases = 0;
the_moon->surf_temp = 0;
the_moon->high_temp = 0;
the_moon->low_temp = 0;
the_moon->max_temp = 0;
the_moon->min_temp = 0;
the_moon->greenhs_rise = 0;
the_moon->minor_moons = 0;

if ((the_moon->dust_mass + the_moon->gas_mass)
> (the_planet->dust_mass + the_planet->gas_mass))
{
long double temp_dust = the_planet->dust_mass;
long double temp_gas = the_planet->gas_mass;
long double temp_mass = the_planet->mass;

the_planet->dust_mass = the_moon->dust_mass;
the_planet->gas_mass = the_moon->gas_mass;
the_planet->mass = the_moon->mass;

the_moon->dust_mass = temp_dust;
the_moon->gas_mass = temp_gas;
the_moon->mass = temp_mass;
}

if (the_planet->first_moon == NULL)
the_planet->first_moon = the_moon;
else
{
the_moon->next_planet = the_planet->first_moon;
the_planet->first_moon = the_moon;
}

finished = TRUE;

if (flag_verbose & 0x0100)
fprintf (stderr, "Moon Captured… "
“%5.3Lf AU (%.2LfEM) <- %.2LfEM\n”,
the_planet->a, the_planet->mass * SUN_MASS_IN_EARTH_MASSES,
mass * SUN_MASS_IN_EARTH_MASSES
);
}
else
{
if (flag_verbose & 0x0100)
fprintf (stderr, "Moon Escapes… “
“%5.3Lf AU (%.2LfEM)%s %.2LfEM%s\n”,
the_planet->a, the_planet->mass * SUN_MASS_IN_EARTH_MASSES,
existing_mass < (the_planet->mass * .05) ? “” : " (big moons)”,
mass * SUN_MASS_IN_EARTH_MASSES,
(mass * SUN_MASS_IN_EARTH_MASSES) >= 2.5 ? “, too big” :
(mass * SUN_MASS_IN_EARTH_MASSES) <= .0001 ? “, too small” : “”
);
}
}
}

if (!finished)
{
if (flag_verbose & 0x0100)
fprintf (stderr, "Collision between two planetesimals! "
“%4.2Lf AU (%.2LfEM) + %4.2Lf AU (%.2LfEM = %.2LfEMd + %.2LfEMg [%.3LfEM])-> %5.3Lf AU (%5.3Lf)\n”,
the_planet->a, the_planet->mass * SUN_MASS_IN_EARTH_MASSES,
a, mass * SUN_MASS_IN_EARTH_MASSES,
dust_mass * SUN_MASS_IN_EARTH_MASSES, gas_mass * SUN_MASS_IN_EARTH_MASSES,
crit_mass * SUN_MASS_IN_EARTH_MASSES,
new_a, e);

temp = the_planet->mass + mass;
accrete_dust(&temp, &new_dust, &new_gas,
new_a,e,stell_luminosity_ratio,
body_inner_bound,body_outer_bound);

the_planet->a = new_a;
the_planet->e = e;
the_planet->mass = temp;
the_planet->dust_mass += dust_mass + new_dust;
the_planet->gas_mass += gas_mass + new_gas;
if (temp >= crit_mass)
the_planet->gas_giant = TRUE;

while (the_planet->next_planet != NULL && the_planet->next_planet->a < new_a)
{
next_planet = the_planet->next_planet;

if (the_planet == planet_head)
planet_head = next_planet;
else
prev_planet->next_planet = next_planet;

the_planet->next_planet = next_planet->next_planet;
next_planet->next_planet = the_planet;
prev_planet = next_planet;
}
}

finished = TRUE;
break;
}
else
{
prev_planet = the_planet;
}
}

if (!(finished)) // Planetesimals didn’t collide. Make it a planet.
{
the_planet = (planets *)malloc(sizeof(planets));

the_planet->type = tUnknown;
the_planet->a = a;
the_planet->e = e;
the_planet->mass = mass;
the_planet->dust_mass = dust_mass;
the_planet->gas_mass = gas_mass;
the_planet->atmosphere = NULL;
the_planet->first_moon = NULL;
the_planet->atmosphere = NULL;
the_planet->albedo = 0;
the_planet->gases = 0;
the_planet->surf_temp = 0;
the_planet->high_temp = 0;
the_planet->low_temp = 0;
the_planet->max_temp = 0;
the_planet->min_temp = 0;
the_planet->greenhs_rise = 0;
the_planet->minor_moons = 0;

if ((mass >= crit_mass))
the_planet->gas_giant = TRUE;
else
the_planet->gas_giant = FALSE;

if ((planet_head == NULL))
{
planet_head = the_planet;
the_planet->next_planet = NULL;
}
else if ((a < planet_head->a))
{
the_planet->next_planet = planet_head;
planet_head = the_planet;
}
else if ((planet_head->next_planet == NULL))
{
planet_head->next_planet = the_planet;
the_planet->next_planet = NULL;
}
else
{
next_planet = planet_head;
while (((next_planet != NULL) && (next_planet->a < a)))
{
prev_planet = next_planet;
next_planet = next_planet->next_planet;
}
the_planet->next_planet = next_planet;
prev_planet->next_planet = the_planet;
}
}
}


planet_pointer dist_planetary_masses(long double stell_mass_ratio,
long double stell_luminosity_ratio,
long double inner_dust,
long double outer_dust,
long double outer_planet_limit,
long double dust_density_coeff,
planet_pointer seed_system,
int do_moons)
{
long double a;
long double e;
long double mass;
long double dust_mass;
long double gas_mass;
long double crit_mass;
long double planet_inner_bound;
long double planet_outer_bound;
planet_pointer seeds = seed_system;

set_initial_conditions(inner_dust,outer_dust);
planet_inner_bound = nearest_planet(stell_mass_ratio);

if (outer_planet_limit == 0)
planet_outer_bound = farthest_planet(stell_mass_ratio);
else
planet_outer_bound = outer_planet_limit;

while (dust_left)
{
if (seeds != NULL)
{
a = seeds->a;
e = seeds->e;
seeds = seeds->next_planet;
}
else
{
a = random_number(planet_inner_bound,planet_outer_bound);
e = random_eccentricity( );
}

mass = PROTOPLANET_MASS;
dust_mass = 0;
gas_mass = 0;

if (flag_verbose & 0x0200)
fprintf (stderr, “Checking %Lg AU.\n”,a);

if (dust_available(inner_effect_limit(a, e, mass),
outer_effect_limit(a, e, mass)))
{
if (flag_verbose & 0x0100)
fprintf (stderr, “Injecting protoplanet at %Lg AU.\n”, a);

dust_density = dust_density_coeff * sqrt(stell_mass_ratio)
* exp(-ALPHA * pow(a,(1.0 / N)));
crit_mass = critical_limit(a,e,stell_luminosity_ratio);
accrete_dust(&mass, &dust_mass, &gas_mass,
a,e,crit_mass,
planet_inner_bound,
planet_outer_bound);

dust_mass += PROTOPLANET_MASS;

if (mass > PROTOPLANET_MASS)
coalesce_planetesimals(a,e,mass,crit_mass,
dust_mass, gas_mass,
stell_luminosity_ratio,
planet_inner_bound,planet_outer_bound,
do_moons);
else if (flag_verbose & 0x0100)
fprintf (stderr, “… failed due to large neighbor.\n”);
}
else if (flag_verbose & 0x0200)
fprintf (stderr, “… failed.\n”);
}
return(planet_head);
}

void free_dust (dust_pointer head)
{
dust_pointer node;
dust_pointer next;

for(node = head;
node != NULL;
node = next)
{
next = node->next_band;
free (node);
}

}

void free_planet (planet_pointer head)
{
planet_pointer node;
planet_pointer next;

for(node = head;
node != NULL;
node = next)
{
next = node->next_planet;

free (node);
}
}

void free_generations()
{
gen_pointer node;
gen_pointer next;

for(node = hist_head;
node != NULL;
node = next)
{
next = node->next;

if (node->dusts)
free_dust (node->dusts);

if (node->planets)
free_planet (node->planets);

free (node);
}

if (dust_head != NULL)
free_dust (dust_head);

if (planet_head != NULL)
free_planet (planet_head);

dust_head = NULL;
planet_head = NULL;
hist_head = NULL;
}

void free_atmosphere(planet_pointer head)
{
planet_pointer node;

for (node = head;
node != NULL;
node = node->next_planet)
{
if (node->atmosphere != NULL)
{
free(node->atmosphere);

node->atmosphere = NULL;
}

if (node->first_moon != NULL)
{
free_atmosphere(node->first_moon);
}
}
}