void generate_planets(primary*		sun,
					  int 			random_tilt, 
					  char			flag_char,
					  int			sys_no,
					  char*			system_name,
					  int			do_gases,
					  int			do_moons)
{
	planet_pointer	planet;
	int				planet_no = 0;

	for (planet = innermost_planet, planet_no = 1; 
		 planet != NULL; 
		 planet = planet->next_planet, planet_no++)
	{
		planet->first_moon		= NULL;
		planet->atmosphere		= NULL;
		planet->gases			= 0;
		planet->high_temp		= 0;
		planet->low_temp		= 0;
		planet->max_temp		= 0;
		planet->min_temp		= 0;
		planet->greenhs_rise	= 0;
		planet->planet_no		= planet_no;
		planet->primary_body			= sun;
		planet->resonant_period = FALSE;

#ifdef MOONS
		if (do_moons)
		{
			long double m1 = planet->mass;
			long double m2 = sun->mass;
			long double mu = m2 / (m1 + m2);
			long double e = planet->e;
			long double a = planet->a;
			long double outer_limit = (0.464 + (-0.380 * mu) + (-0.631 * e) +
									   (0.586 * mu * e) + (0.150 * pow2(e)) +
									   (-0.198 * mu * pow2(e))) * a;
//			long double farthest = farthest_planet(planet->mass);
			long double outer_dust = stellar_dust_limit(planet->mass);

//			if (outer_limit > farthest)
//				outer_limit = farthest;
			
		    planet->first_moon = dist_moon_masses (planet->mass,
												   sun->luminosity,
												   planet->a, 
												   planet->e,
												   0.0, 
												   outer_dust,
												   1.0,
												   outer_limit);
			if (planet->first_moon)
			{
				int 			n;
				planet_pointer	ptr;
				for (n=0, ptr=planet->first_moon;
					 ptr != NULL;
					 ptr=ptr->next_planet)
				{
					if (ptr->mass * SUN_MASS_IN_EARTH_MASSES > .0001)
					{
						n++;
						fprintf (stderr, "%d: %d %6.3Lg %6.3Lg\n",	
								planet_no,
								n,
								ptr->mass * SUN_MASS_IN_EARTH_MASSES,
								ptr->mass / planet->mass);
					}
				}
			}
		}
#endif

		planet->orbit_zone 		= orb_zone(sun->luminosity, planet->a);

		planet->orb_period		= period(planet->a,planet->mass,sun->mass);
		if (random_tilt)
			planet->axial_tilt 	= inclination(planet->a);
		planet->exospheric_temp = EARTH_EXOSPHERE_TEMP / pow2(planet->a / sun->r_ecosphere);
		planet->rms_velocity 	= rms_vel(MOL_NITROGEN,planet->exospheric_temp);
		planet->core_radius 	= kothari_radius(planet->dust_mass,FALSE,planet->orbit_zone);

		// Calculate the radius as a gas giant, to verify it will retain gas.
		// Then if mass > Earth, it's at least 5% gas and retains He, it's
		// some flavor of gas giant.
		
		planet->density 		= empirical_density(planet->mass,planet->a, sun->r_ecosphere, TRUE);
		planet->radius 			= volume_radius(planet->mass,planet->density);
		
		planet->surf_accel   	= acceleration(planet->mass,planet->radius);
		planet->surf_grav 	 	= gravity(planet->surf_accel);

		planet->molec_weight	= min_molec_weight(planet);

		if (((planet->mass * SUN_MASS_IN_EARTH_MASSES) > 1.0)
		  && ((planet->gas_mass / planet->mass)        > 0.05)
		  && (min_molec_weight(planet)				  <= 4.0))
		{
			if ((planet->gas_mass / planet->mass) < 0.20)
				planet->type = tSubSubGasGiant;
			else if ((planet->mass * SUN_MASS_IN_EARTH_MASSES) < 20.0)
				planet->type = tSubGasGiant;
			else
				planet->type = tGasGiant;
		}
		else // If not, it's rocky.
		{
			planet->radius 		= kothari_radius(planet->mass,FALSE,planet->orbit_zone);
			planet->density 	= volume_density(planet->mass,planet->radius);

			planet->surf_accel  = acceleration(planet->mass,planet->radius);
			planet->surf_grav 	= gravity(planet->surf_accel);
			
			if ((planet->gas_mass / planet->mass)        > 0.000001)
			{
				long double h2_mass = planet->gas_mass * 0.85;
				long double he_mass = (planet->gas_mass - h2_mass) * 0.999;
				
				long double h2_loss = 0.0;
				long double he_loss = 0.0;
				
							 
				long double h2_life = gas_life(MOL_HYDROGEN, planet);
				long double he_life = gas_life(HELIUM, planet);
								
				if (h2_life < sun->age)
				{
					h2_loss = ((1.0 - (1.0 / exp(sun->age / h2_life))) * h2_mass);
					
					planet->gas_mass -= h2_loss;
					planet->mass     -= h2_loss;
					
					planet->surf_accel  = acceleration(planet->mass,planet->radius);
					planet->surf_grav 	= gravity(planet->surf_accel);
				}
				
				if (he_life < sun->age)
				{
					he_loss = ((1.0 - (1.0 / exp(sun->age / he_life))) * he_mass);
					
					planet->gas_mass -= he_loss;
					planet->mass     -= he_loss;
					
					planet->surf_accel  = acceleration(planet->mass,planet->radius);
					planet->surf_grav 	= gravity(planet->surf_accel);
				}
				
				if (((h2_loss + he_loss) > .000001) && (flag_verbose & 0x0080))
					fprintf (stderr, "%s (-s%ld -%c%d) %d\tLosing gas: H2: %5.3Lf EM, He: %5.3Lf EM\n",
							 system_name,
							 flag_seed, flag_char, sys_no,
							 planet_no,
							 h2_loss * SUN_MASS_IN_EARTH_MASSES, he_loss * SUN_MASS_IN_EARTH_MASSES);
			}
		}

		planet->day = day_length(planet);	/* Modifies planet->resonant_period */
		planet->esc_velocity 	= escape_vel(planet->mass,planet->radius);

		if ((planet->type == tGasGiant)
		 || (planet->type == tSubGasGiant) 
		 || (planet->type == tSubSubGasGiant))
		{
			planet->greenhouse_effect 		= FALSE;
			planet->volatile_gas_inventory 	= INCREDIBLY_LARGE_NUMBER;
			planet->surf_pressure 			= INCREDIBLY_LARGE_NUMBER;
			
			planet->boil_point 				= INCREDIBLY_LARGE_NUMBER;
			
			planet->surf_temp				= INCREDIBLY_LARGE_NUMBER;
			planet->greenhs_rise 			= 0;
			planet->albedo 					= about(GAS_GIANT_ALBEDO,0.1);
			planet->hydrosphere 			= 1.0;
			planet->cloud_cover	 			= 1.0;
			planet->ice_cover	 			= 0.0;
			planet->surf_grav			 	= gravity(planet->surf_accel);
			planet->molec_weight			= min_molec_weight(planet);
			planet->surf_grav 				= INCREDIBLY_LARGE_NUMBER;
 			planet->estimated_temp			= est_temp(sun->r_ecosphere, planet->a,  planet->albedo);
			planet->estimated_terr_temp		= est_temp(sun->r_ecosphere, planet->a,  EARTH_ALBEDO);
 			
			{
				long double temp = planet->estimated_terr_temp;

				if ((temp >= FREEZING_POINT_OF_WATER)
				 && (temp <= EARTH_AVERAGE_KELVIN + 10.)
				 && (sun->age > 2.0E9))
				{
					habitable_jovians++;
				
					if (flag_verbose & 0x8000)
					{
						fprintf (stderr, "%s (-s%ld -%c%d) %d\t%s (%4.2LfEM %5.3Lf By) with earth-like temperature (%.1Lf C, %.1Lf F, %+.1Lf C Earth).\n",
								 system_name,
								 flag_seed, flag_char, sys_no,
								 planet_no,
								 planet->type == tGasGiant ? "Jovian" :
								 planet->type == tSubGasGiant ? "Sub-Jovian" :
								 planet->type == tSubSubGasGiant ? "Gas Dwarf" :
								 "Big",
								 planet->mass * SUN_MASS_IN_EARTH_MASSES,
								 sun->age /1.0E9,
								 temp - FREEZING_POINT_OF_WATER,
								 32 + ((temp - FREEZING_POINT_OF_WATER) * 1.8),
								 temp - EARTH_AVERAGE_KELVIN);
					}
				}
			}
		}
		else
		{
			planet->estimated_temp			= est_temp(sun->r_ecosphere, planet->a,  EARTH_ALBEDO);
			planet->estimated_terr_temp		= est_temp(sun->r_ecosphere, planet->a,  EARTH_ALBEDO);

			planet->surf_grav 				= gravity(planet->surf_accel);
			planet->molec_weight			= min_molec_weight(planet);
	
			planet->greenhouse_effect 		= grnhouse(sun->r_ecosphere, planet->a);
			planet->volatile_gas_inventory 	= vol_inventory(planet->mass,
														    planet->esc_velocity,
														    planet->rms_velocity,
														    sun->mass,
														    planet->orbit_zone,
														    planet->greenhouse_effect,
														    (planet->gas_mass
														     / planet->mass) > 0.000001);
			planet->surf_pressure 			= pressure(planet->volatile_gas_inventory,
													   planet->radius,
													   planet->surf_grav);

			if ((planet->surf_pressure == 0.0))
				planet->boil_point 			= 0.0;
			else 
				planet->boil_point 			= boiling_point(planet->surf_pressure);
			
			iterate_surface_temp(planet);		/*	Sets:
												 *		planet->surf_temp
												 *		planet->greenhs_rise
												 *		planet->albedo
												 *		planet->hydrosphere
												 *		planet->cloud_cover
												 *		planet->ice_cover
												 */

			if (do_gases &&
				(planet->max_temp >= FREEZING_POINT_OF_WATER) &&
				(planet->min_temp <= planet->boil_point))
				calculate_gases(sun, planet, flag_char, sys_no, planet_no, system_name);
			
			/*
			 *	Next we assign a type to the planet.
			 */
			 
			if (planet->surf_pressure < 1.0)
			{
				if ((planet->mass * SUN_MASS_IN_EARTH_MASSES) < ASTEROID_MASS_LIMIT)
					planet->type 			= tAsteroids;
				else
					planet->type 			= tRock;
			}
			else if ((planet->surf_pressure > 6000.0) &&
					 (planet->molec_weight <= 2.0))	// Retains Hydrogen
			{
				planet->type = tSubSubGasGiant;
				planet->gases = 0;
				free(planet->atmosphere);
				planet->atmosphere = NULL;
			}
			else 
			{										// Atmospheres:
				if (((int)planet->day == (int)(planet->orb_period * 24.0) || 
										 (planet->resonant_period)))
					planet->type = t1Face;
				else if (planet->hydrosphere >= 0.95)
					planet->type = tWater;				// >95% water
				else if (planet->ice_cover >= 0.95)
					planet->type = tIce;				// >95% ice
				else if (planet->hydrosphere > 0.05)
					planet->type = tTerrestrial;		// Terrestrial
														// else <5% water
				else if (planet->max_temp > planet->boil_point)
					planet->type = tVenusian;			// Hot = Venusian
				else if ((planet->gas_mass / planet->mass) > 0.0001)
				{										// Accreted gas
					planet->type = tIce;				// But no Greenhouse
					planet->ice_cover = 1.0;			// or liquid water
				}										// Make it an Ice World
				else if (planet->surf_pressure <= 250.0)// Thin air = Martian
					planet->type = tMartian;
				else if (planet->surf_temp < FREEZING_POINT_OF_WATER)
					planet->type = tIce;
				else
				{
					planet->type = tUnknown;
					
					if (flag_verbose & 0x0001)
						fprintf (stderr, "%12s\tp=%4.2Lf\tm=%4.2Lf\tg=%4.2Lf\tt=%+.1Lf\t-s%ld -%c%d\t%s %d\t Unknown %s\n", 
										type_string (planet->type),
										planet->surf_pressure,
										planet->mass * SUN_MASS_IN_EARTH_MASSES,
										planet->surf_grav,
										planet->surf_temp  - EARTH_AVERAGE_KELVIN,
										flag_seed, 
										flag_char,
										sys_no,
										system_name,
										planet_no,
										((int)planet->day == (int)(planet->orb_period * 24.0) || 
										 (planet->resonant_period)) ? "(1-Face)" : ""
								 );
				}
			}
		}
		
		/*
		 *	Now we're ready to test for habitable planets,
		 *	so we can count and log them and such
		 */
{
	int tIndex = 0;

	switch (planet->type)
	{
		case tUnknown:			tIndex = 0;		break;
		case tRock:				tIndex = 1;		break;
		case tVenusian:			tIndex = 2;		break;
		case tTerrestrial:		tIndex = 3;		break;
		case tSubSubGasGiant:	tIndex = 4;		break;
		case tSubGasGiant:		tIndex = 5;		break;
		case tGasGiant:			tIndex = 6;		break;
		case tMartian:			tIndex = 7;		break;
		case tWater:			tIndex = 8;		break;
		case tIce:				tIndex = 9;		break;
		case tAsteroids: 		tIndex = 10;	break;
		case t1Face:			tIndex = 11;	break;
	}
	
	if (type_counts[tIndex] == 0)
		++type_count;
	
	++type_counts[tIndex];
	
}
		{
			unsigned int breathe = breathability (planet);
			
			if ((breathe == BREATHABLE) &&
				(!planet->resonant_period) &&		// Option needed?
				((int)planet->day != (int)(planet->orb_period * 24.0)))
			{
				int	list_it	= FALSE;
				
				habitable++;
				
				if (min_breathable_temp > planet->surf_temp)
				{
					min_breathable_temp = planet->surf_temp;

					if (flag_verbose & 0x0002)
						list_it = TRUE;
				}

				if (max_breathable_temp < planet->surf_temp)
				{
					max_breathable_temp = planet->surf_temp;

					if (flag_verbose & 0x0002)
						list_it = TRUE;
				}
			
				if (min_breathable_g > planet->surf_grav)
				{
					min_breathable_g = planet->surf_grav;

					if (flag_verbose & 0x0002)
						list_it = TRUE;
				}
				
				if (max_breathable_g < planet->surf_grav)
				{
					max_breathable_g = planet->surf_grav;

					if (flag_verbose & 0x0002)
						list_it = TRUE;
				}
				
				if (planet->type == tTerrestrial)
				{
					if (min_breathable_terrestrial_g > planet->surf_grav)
					{
						min_breathable_terrestrial_g = planet->surf_grav;

						if (flag_verbose & 0x0002)
							list_it = TRUE;
					}

					if (max_breathable_terrestrial_g < planet->surf_grav)
					{
						max_breathable_terrestrial_g = planet->surf_grav;

						if (flag_verbose & 0x0002)
							list_it = TRUE;
					}
				}

				if (min_breathable_p > planet->surf_pressure)
				{
					min_breathable_p = planet->surf_pressure;

					if (flag_verbose & 0x0002)
						list_it = TRUE;
				}
				
				if (max_breathable_p < planet->surf_pressure)
				{
					max_breathable_p = planet->surf_pressure;

					if (flag_verbose & 0x0002)
						list_it = TRUE;
				}
									
				if (flag_verbose & 0x0004)
					list_it = TRUE;
				
				if (list_it)
					fprintf (stderr, "%12s\tp=%4.2Lf\tm=%4.2Lf\tg=%4.2Lf\tt=%+.1Lf\t-s%ld -%c%d\t%s %d\n", 
							type_string (planet->type),
							planet->surf_pressure,
							planet->mass * SUN_MASS_IN_EARTH_MASSES,
							planet->surf_grav,
							planet->surf_temp  - EARTH_AVERAGE_KELVIN,
							flag_seed, 
							flag_char,
							sys_no,
							system_name,
							planet_no);
			}
		}


		if ((flag_verbose & 0x0800)
		 && (planet->dust_mass * SUN_MASS_IN_EARTH_MASSES >= 0.0006)
		 && (planet->gas_mass * SUN_MASS_IN_EARTH_MASSES >= 0.0006)
		 && (planet->type != tGasGiant) 
		 && (planet->type != tSubGasGiant)
		   )
		{
			int core_size = (int)((50. * planet->core_radius) / planet->radius);
			
			if (core_size <= 49)
			{
				fprintf (stderr, "%12s\tp=%4.2Lf\tr=%4.2Lf\tm=%4.2Lf\t-s%ld -%c%d\t%d\t%s %d\n", 
						type_string (planet->type),
						planet->core_radius,
						planet->radius,
						planet->mass * SUN_MASS_IN_EARTH_MASSES,
						flag_seed, 
						flag_char,
						sys_no,
						50-core_size,
						system_name,
						planet_no);
			}
		}
		
		{
			long double  rel_temp   = (planet->surf_temp -  FREEZING_POINT_OF_WATER) - 
						 			   EARTH_AVERAGE_CELSIUS;
			long double	 seas       = (planet->hydrosphere * 100.0);
			long double	 clouds     = (planet->cloud_cover * 100.0);
			long double	 pressure   = (planet->surf_pressure / 
									   EARTH_SURF_PRES_IN_MILLIBARS);
			long double	 ice        = (planet->ice_cover * 100.0);
			long double	 gravity    = planet->surf_grav;
			unsigned int breathe    = breathability (planet);
		
			if ((gravity 	>= .8) &&
				(gravity 	<= 1.2) &&
				(rel_temp 	>= -2.0) &&
				(rel_temp 	<= 3.0) &&
				(ice 		<= 10.) &&
				(pressure   >= 0.5) &&
				(pressure   <= 2.0) &&
				(clouds		>= 40.) &&
				(clouds		<= 80.) &&
				(seas 		>= 50.) &&
				(seas 		<= 80.) &&
				(planet->type != tWater) &&
				(breathe    == BREATHABLE))
			{
				earthlike++;
				
				if (flag_verbose & 0x0008)
					fprintf (stderr, "%12s\tp=%4.2Lf\tm=%4.2Lf\tg=%4.2Lf\tt=%+.1Lf\t%d -s%ld -%c%d\t%s %d\tEarth-like\n", 
									type_string (planet->type),
									planet->surf_pressure,
									planet->mass * SUN_MASS_IN_EARTH_MASSES,
									planet->surf_grav,
									planet->surf_temp  - EARTH_AVERAGE_KELVIN,
									habitable,
									flag_seed, 
									flag_char,
									sys_no,
									system_name,
									planet_no);
			} else if ((flag_verbose & 0x0008) &&
					 (breathe   == BREATHABLE) &&
					 (gravity	 > 1.3) &&
					 (habitable	 > 1) &&
					 ((rel_temp  < -2.0) ||
					  (ice		 > 10.))
					)
			{
				fprintf (stderr, "%12s\tp=%4.2Lf\tm=%4.2Lf\tg=%4.2Lf\tt=%+.1Lf\t%d -s%ld -%c%d\t%s %d\tSphinx-like\n", 
						type_string (planet->type),
						planet->surf_pressure,
						planet->mass * SUN_MASS_IN_EARTH_MASSES,
						planet->surf_grav,
						planet->surf_temp  - EARTH_AVERAGE_KELVIN,
						habitable,
						flag_seed, 
						flag_char,
						sys_no,
						system_name,
						planet_no);
			}
		}
	}
}