home *** CD-ROM | disk | FTP | other *** search
/ Gold Fish 2 / goldfish_vol2_cd1.bin / files / misc / sci / accrete / src / enviro.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-06-12  |  25.0 KB  |  653 lines
  1. #include <stdio.h>
  2. #include <string.h>
  3. #include <math.h>
  4. #include <stddef.h>
  5. #include <stdlib.h>
  6. #ifdef AMIGA
  7. #ifdef HUGE_VAL
  8. #undef HUGE_VAL
  9. #endif /* HUGE_VAL */
  10. #endif /* AMIGA */
  11. #include <float.h>
  12. #include <time.h>
  13.  
  14. #include "const.h"
  15. #include "structs.h"
  16. #include "accrete.h"
  17.  
  18. double luminosity(mass_ratio)
  19. double mass_ratio; 
  20. {
  21.      double n; 
  22.      
  23.      if (mass_ratio < 1.0)
  24.       n = 1.75 * (mass_ratio - 0.1) + 3.325;
  25.      else 
  26.       n = 0.5 * (2.0 - mass_ratio) + 4.4;
  27.      return(pow(mass_ratio,n));
  28. }
  29.  
  30.  
  31. /*--------------------------------------------------------------------------*/
  32. /*   This function, given the orbital radius of a planet in AU, returns     */
  33. /*   the orbital 'zone' of the particle.                                    */
  34. /*--------------------------------------------------------------------------*/
  35.  
  36. int orbital_zone(orbital_radius)
  37. double orbital_radius; 
  38. {
  39.      if (orbital_radius < (4.0 * sqrt(stellar_luminosity_ratio)))
  40.       return(1);
  41.      else 
  42.      {
  43.       if ((orbital_radius >= (4.0 * sqrt(stellar_luminosity_ratio))) && (orbital_radius < (15.0 * sqrt(stellar_luminosity_ratio))))
  44.            return(2);
  45.       else 
  46.            return(3);
  47.      }
  48. }
  49.  
  50.  
  51. /*--------------------------------------------------------------------------*/
  52. /*   The mass is in units of solar masses, and the density is in units      */
  53. /*   of grams/cc.  The radius returned is in units of km.                   */
  54. /*--------------------------------------------------------------------------*/
  55.  
  56. double volume_radius(mass, density)
  57. double mass, density;
  58. {
  59.      double volume; 
  60.      
  61.      mass = mass * SOLAR_MASS_IN_GRAMS;
  62.      volume = mass / density;
  63.      return(pow((3.0 * volume) / (4.0 * PI),(1.0 / 3.0)) / CM_PER_KM);
  64. }
  65.  
  66. /*--------------------------------------------------------------------------*/
  67. /*    Returns the radius of the planet in kilometers.                       */
  68. /*   The mass passed in is in units of solar masses, the orbital radius     */
  69. /*   in A.U.                                                                */
  70. /*   This formula is listed as eq.9 in Fogg's article, although some typos  */
  71. /*   crop up in that eq.  See "The Internal Constitution of Planets", by    */
  72. /*   Dr. D. S. Kothari, Mon. Not. of the Royal Astronomical Society, vol 96 */
  73. /*   pp.833-843, 1936 for the derivation.  Specifically, this is Kothari's  */
  74. /*   eq.23, which appears on page 840.                                      */
  75. /*--------------------------------------------------------------------------*/
  76.  
  77. double kothari_radius(mass, orbital_radius, giant, zone)
  78. double mass, orbital_radius;
  79. int giant, zone;
  80. {
  81.      double temp, temp2, atomic_weight, atomic_num;
  82.      
  83.      if (zone == 1)
  84.      {
  85.       if (giant)
  86.       {
  87.            atomic_weight = 9.5;
  88.            atomic_num = 4.5;
  89.       }
  90.       else 
  91.       {
  92.            atomic_weight = 15.0;
  93.            atomic_num = 8.0;
  94.       }
  95.      }
  96.      else 
  97.       if (zone == 2)
  98.       {
  99.            if (giant)
  100.            {
  101.             atomic_weight = 2.47;
  102.             atomic_num = 2.0;
  103.            }
  104.            else 
  105.            {
  106.             atomic_weight = 10.0;
  107.             atomic_num = 5.0;
  108.            }
  109.       }
  110.       else 
  111.       {
  112.            if (giant)
  113.            {
  114.             atomic_weight = 7.0;
  115.             atomic_num = 4.0;
  116.            }
  117.            else 
  118.            {
  119.             atomic_weight = 10.0;
  120.             atomic_num = 5.0;
  121.            }
  122.       }
  123.      temp = atomic_weight * atomic_num;
  124.      temp = (2.0 * BETA_20 * pow(SOLAR_MASS_IN_GRAMS,(1.0 / 3.0))) / (A1_20 * pow(temp,(1.0 / 3.0)));
  125.      temp2 = A2_20 * pow(atomic_weight,(4.0 / 3.0)) * pow(SOLAR_MASS_IN_GRAMS,(2.0 / 3.0));
  126.      temp2 = temp2 * pow(mass,(2.0 / 3.0));
  127.      temp2 = temp2 / (A1_20 * pow(atomic_num, 2.0));
  128.      temp2 = 1.0 + temp2;
  129.      temp = temp / temp2;
  130.      temp = (temp * pow(mass,(1.0 / 3.0))) / CM_PER_KM;
  131.      return(temp);
  132. }
  133.  
  134.  
  135. /*--------------------------------------------------------------------------*/
  136. /*  The mass passed in is in units of solar masses, and the orbital radius  */
  137. /*  is in units of AU.  The density is returned in units of grams/cc.       */
  138. /*--------------------------------------------------------------------------*/
  139.  
  140. double empirical_density(mass, orbital_radius, gas_giant)
  141. double mass, orbital_radius;
  142. int gas_giant;
  143. {
  144.      double temp; 
  145.      
  146.      temp = pow(mass * EARTH_MASSES_PER_SOLAR_MASS,(1.0 / 8.0));
  147.      temp = temp * pow(r_ecosphere / orbital_radius,(1.0 / 4.0));
  148.      if (gas_giant)
  149.       return(temp * 1.2);
  150.      else 
  151.       return(temp * 5.5);
  152. }
  153.  
  154.  
  155. /*--------------------------------------------------------------------------*/
  156. /*  The mass passed in is in units of solar masses, and the equatorial      */
  157. /*  radius is in km.  The density is returned in units of grams/cc.         */
  158. /*--------------------------------------------------------------------------*/
  159.  
  160. double volume_density(mass, equatorial_radius)
  161. double mass, equatorial_radius;
  162. {
  163.      double volume; 
  164.      
  165.      mass = mass * SOLAR_MASS_IN_GRAMS;
  166.      equatorial_radius = equatorial_radius * CM_PER_KM;
  167.      volume = (4.0 * PI * pow(equatorial_radius,3.0)) / 3.0;
  168.      return(mass / volume);
  169. }
  170.  
  171.  
  172. /*--------------------------------------------------------------------------*/
  173. /*  The separation is in units of AU, and both masses are in units of solar */
  174. /*  masses.  The period returned is in terms of Earth days.                 */
  175. /*--------------------------------------------------------------------------*/
  176.  
  177. double period(separation, small_mass, large_mass)
  178. double separation, small_mass, large_mass;
  179. {
  180.      double period_in_years; 
  181.      
  182.      period_in_years = sqrt(pow(separation,3.0) / (small_mass + large_mass));
  183.      return(period_in_years * DAYS_IN_A_YEAR);
  184. }
  185.  
  186.  
  187. /*--------------------------------------------------------------------------*/
  188. /*   Fogg's information for this routine came from Dole "Habitable Planets  */
  189. /* for Man", Blaisdell Publishing Company, NY, 1964.  From this, he came    */
  190. /* up with his eq.12, which is the equation for the base_angular_velocity   */
  191. /* below.  Going a bit further, he found an equation for the change in      */
  192. /* angular velocity per time (dw/dt) from P. Goldreich and S. Soter's paper */
  193. /* "Q in the Solar System" in Icarus, vol 5, pp.375-389 (1966).  Comparing  */
  194. /* to the change in angular velocity for the Earth, we can come up with an  */
  195. /* approximation for our new planet (his eq.13) and take that into account. */
  196. /*--------------------------------------------------------------------------*/
  197.  
  198. double day_length(mass, radius, orbital_period, eccentricity, giant)
  199. double mass, radius, orbital_period, eccentricity;
  200. int giant;
  201. {
  202.      double base_angular_velocity, planetary_mass_in_grams, k2, temp,
  203.      equatorial_radius_in_cm, change_in_angular_velocity, spin_resonance_period;
  204.      
  205.      spin_resonance = FALSE;
  206.      if (giant)
  207.       k2 = 0.24;
  208.      else 
  209.       k2 = 0.33;
  210.      planetary_mass_in_grams = mass * SOLAR_MASS_IN_GRAMS;
  211.      equatorial_radius_in_cm = radius * CM_PER_KM;
  212.      base_angular_velocity = sqrt(2.0 * J * (planetary_mass_in_grams) / (k2 * pow(equatorial_radius_in_cm, 2.0)));
  213.      /*   This next term describes how much a planet's rotation is slowed by    */
  214.      /*  it's moons.  Find out what dw/dt is after figuring out Goldreich and   */
  215.      /*  Soter's Q'.                                                            */
  216.      change_in_angular_velocity = 0.0;
  217.      temp = base_angular_velocity + (change_in_angular_velocity * age);
  218.      /*   'temp' is now the angular velocity. Now we change from rad/sec to     */
  219.      /*  hours/rotation.                               */
  220.      temp = 1.0 / ((temp / radians_per_rotation) * SECONDS_PER_HOUR);
  221.      if (temp >= orbital_period)
  222.      {
  223.       spin_resonance_period = ((1.0 - eccentricity) / (1.0 + eccentricity)) * orbital_period;
  224.       if (eccentricity > 0.1)
  225.       {
  226.            temp = spin_resonance_period;
  227.            spin_resonance = TRUE;
  228.       }
  229.       else 
  230.            temp = orbital_period;
  231.      }
  232.      return(temp);
  233. }
  234.  
  235.  
  236. /*--------------------------------------------------------------------------*/
  237. /*   The orbital radius is expected in units of Astronomical Units (AU).    */
  238. /*   Inclination is returned in units of degrees.                           */
  239. /*--------------------------------------------------------------------------*/
  240.  
  241. int inclination(orbital_radius)
  242. double orbital_radius; 
  243. {
  244.      int temp; 
  245.      
  246.      temp = (int)(pow(orbital_radius,0.2) * about(EARTH_AXIAL_TILT,0.4));
  247.      return(temp % 360);
  248. }
  249.  
  250.  
  251. /*--------------------------------------------------------------------------*/
  252. /*   This function implements the escape velocity calculation.  Note that   */
  253. /*  it appears that Fogg's eq.15 is incorrect.                              */
  254. /*  The mass is in units of solar mass, the radius in kilometers, and the   */
  255. /*  velocity returned is in cm/sec.                                         */
  256. /*--------------------------------------------------------------------------*/
  257.  
  258. double escape_vel(mass, radius)
  259. double mass, radius;
  260. {
  261.      double mass_in_grams, radius_in_cm;
  262.      
  263.      mass_in_grams = mass * SOLAR_MASS_IN_GRAMS;
  264.      radius_in_cm = radius * CM_PER_KM;
  265.      return(sqrt(2.0 * GRAV_CONSTANT * mass_in_grams / radius_in_cm));
  266. }
  267.  
  268.  
  269. /*--------------------------------------------------------------------------*/
  270. /*  This is Fogg's eq.16.  The molecular weight (usually assumed to be N2)  */
  271. /*  is used as the basis of the Root Mean Square velocity of the molecule   */
  272. /*  or atom.  The velocity returned is in cm/sec.                           */
  273. /*--------------------------------------------------------------------------*/
  274.  
  275. double rms_vel(molecular_weight, orbital_radius)
  276. double molecular_weight, orbital_radius;
  277. {
  278.      double exospheric_temp; 
  279.      
  280.      exospheric_temp = EARTH_EXOSPHERE_TEMP / pow(orbital_radius, 2.0);
  281.      return(sqrt((3.0 * MOLAR_GAS_CONST * exospheric_temp) / molecular_weight) * CM_PER_METER);
  282. }
  283.  
  284.  
  285. /*--------------------------------------------------------------------------*/
  286. /*   This function returns the smallest molecular weight retained by the    */
  287. /*  body, which is useful for determining the atmosphere composition.       */
  288. /*  Orbital radius is in A.U.(ie: in units of the earth's orbital radius),  *)
  289.     (*  mass is in units of solar masses, and equatorial radius is in units of  */
  290. /*  kilometers.                                                             */
  291. /*--------------------------------------------------------------------------*/
  292.  
  293. double molecule_limit(orbital_radius, mass, equatorial_radius)
  294. double orbital_radius, mass, equatorial_radius;
  295. {
  296.      double numerator, denominator1, denominator2, escape_velocity, temp;
  297.      
  298.      escape_velocity = escape_vel(mass,equatorial_radius);
  299.      return((3.0 * pow(GAS_RETENTION_THRESHOLD * CM_PER_METER, 2.0) * MOLAR_GAS_CONST * EARTH_EXOSPHERE_TEMP) / pow(escape_velocity, 2.0));
  300. }
  301.  
  302.  
  303. /*--------------------------------------------------------------------------*/
  304. /*   This function calculates the surface acceleration of a planet.  The    */
  305. /*  mass is in units of solar masses, the radius in terms of km, and the    */
  306. /*  acceleration is returned in units of cm/sec2.                           */
  307. /*--------------------------------------------------------------------------*/
  308.  
  309. double acceleration(mass, radius)
  310. double mass, radius;
  311. {
  312.      return(GRAV_CONSTANT * (mass * SOLAR_MASS_IN_GRAMS) / pow(radius * CM_PER_KM, 2.0));
  313. }
  314.  
  315.  
  316. /*--------------------------------------------------------------------------*/
  317. /*   This function calculates the surface gravity of a planet.  The         */
  318. /*  acceleration is in units of cm/sec2, and the gravity is returned in     */
  319. /*  units of Earth gravities.                                               */
  320. /*--------------------------------------------------------------------------*/
  321.  
  322. double gravity(acceleration)
  323. double acceleration; 
  324. {
  325.      return(acceleration / EARTH_ACCELERATION);
  326. }
  327.  
  328.  
  329. /*--------------------------------------------------------------------------*/
  330. /*  Note that if the orbital radius of the planet is greater than or equal  */
  331. /*  to R_inner, 99% of it's volatiles are assumed to have been deposited in */
  332. /*  surface reservoirs (otherwise, it suffers from the greenhouse effect).  */
  333. /*--------------------------------------------------------------------------*/
  334.  
  335. int greenhouse(zone, orbital_radius, greenhouse_radius)
  336. int zone; 
  337. double orbital_radius, greenhouse_radius;
  338. {
  339.      if ((orbital_radius < greenhouse_radius) && (zone == 1))
  340.       return(TRUE);
  341.      else 
  342.       return(FALSE);
  343. }
  344.  
  345.  
  346. /*--------------------------------------------------------------------------*/
  347. /*  This implements Fogg's eq.17.  The 'inventory' returned is unitless.    */
  348. /*--------------------------------------------------------------------------*/
  349.  
  350. double vol_inventory(mass, escape_vel, rms_vel, stellar_mass, zone, greenhouse_effect)
  351. double mass, escape_vel, rms_vel, stellar_mass;
  352. int zone, greenhouse_effect;
  353. {
  354.      double velocity_ratio, proportion_const, temp1, temp2, mass_in_earth_units;
  355.      
  356.      velocity_ratio = escape_vel / rms_vel;
  357.      if (velocity_ratio >= GAS_RETENTION_THRESHOLD)
  358.      {
  359.       switch (zone) {
  360.            case 1:
  361.             proportion_const = 100000.0;
  362.             break;
  363.            case 2:
  364.             proportion_const = 75000.0;
  365.             break;
  366.            case 3:
  367.             proportion_const = 250.0;
  368.             break;
  369.            default:
  370.             printf("Error: orbital zone not initialized correctly!\n");
  371.             break;
  372.            }
  373.       mass_in_earth_units = mass * EARTH_MASSES_PER_SOLAR_MASS;
  374.       temp1 = (proportion_const * mass_in_earth_units) / stellar_mass;
  375.       temp2 = about(temp1,0.2);
  376.       if (greenhouse_effect)
  377.            return(temp2);
  378.       else 
  379.            return(temp2 / 100.0);
  380.      }
  381.      else 
  382.       return(0.0);
  383. }
  384.  
  385.  
  386. /*--------------------------------------------------------------------------*/
  387. /*  This implements Fogg's eq.18.  The pressure returned is in units of     */
  388. /*  millibars (mb).  The gravity is in units of Earth gravities, the radius */
  389. /*  in units of kilometers.                                                 */
  390. /*--------------------------------------------------------------------------*/
  391.  
  392. double pressure(volatile_gas_inventory, equatorial_radius, gravity)
  393. double volatile_gas_inventory, equatorial_radius, gravity;
  394. {
  395.      equatorial_radius = EARTH_RADIUS_IN_KM / equatorial_radius;
  396.      return(volatile_gas_inventory * gravity / pow(equatorial_radius, 2.0));
  397. }
  398.  
  399. /*--------------------------------------------------------------------------*/
  400. /*   This function returns the boiling point of water in an atmosphere of   */
  401. /*   pressure 'surface_pressure', given in millibars.  The boiling point is */
  402. /*   returned in units of Kelvin.  This is Fogg's eq.21.                    */
  403. /*--------------------------------------------------------------------------*/
  404.  
  405. double boiling_point(surface_pressure)
  406. double surface_pressure; 
  407. {
  408.      double surface_pressure_in_bars; 
  409.      
  410.      surface_pressure_in_bars = surface_pressure / MILLIBARS_PER_BAR;
  411.      return(1.0 / (log(surface_pressure_in_bars) / -5050.5 + 1.0 / 373.0));
  412. }
  413.  
  414.  
  415. /*--------------------------------------------------------------------------*/
  416. /*   This function is Fogg's eq.22.  Given the volatile gas inventory and   */
  417. /*   planetary radius of a planet (in Km), this function returns the        */
  418. /*   fraction of the planet covered with water.                             */
  419. /*   I have changed the function very slightly:  the fraction of Earth's    */
  420. /*   surface covered by water is 71%, not 75% as Fogg used.                 */
  421. /*--------------------------------------------------------------------------*/
  422.  
  423. double hydrosphere_fraction(volatile_gas_inventory, planetary_radius)
  424. double volatile_gas_inventory, planetary_radius;
  425. {
  426.      double temp; 
  427.      
  428.      temp = (0.71 * volatile_gas_inventory / 1000.0) * pow(EARTH_RADIUS_IN_KM / planetary_radius, 2.0);
  429.      if (temp >= 1.0)
  430.       return(1.0);
  431.      else 
  432.       return(temp);
  433. }
  434.  
  435.  
  436. /*--------------------------------------------------------------------------*/
  437. /*   Given the surface temperature of a planet (in Kelvin), this function   */
  438. /*   returns the fraction of cloud cover available.  This is Fogg's eq.23.  */
  439. /*   See Hart in "Icarus" (vol 33, pp23 - 39, 1978) for an explanation.     */
  440. /*   This equation is Hart's eq.3.                                          */
  441. /*   I have modified it slightly using constants and relationships from     */
  442. /*   Glass's book "Introduction to Planetary Geology", p.46.                */
  443. /*   The 'CLOUD_COVERAGE_FACTOR' is the amount of surface area on Earth     */
  444. /*   covered by one Kg. of cloud.                        */
  445. /*--------------------------------------------------------------------------*/
  446.  
  447. double cloud_fraction(surface_temp, smallest_MW_retained, equatorial_radius, hydrosphere_fraction)
  448. double surface_temp, smallest_MW_retained, equatorial_radius,
  449.      hydrosphere_fraction;
  450. {
  451.      double water_vapor_in_kg, fraction, surface_area, hydrosphere_mass;
  452.      
  453.      if (smallest_MW_retained > WATER_VAPOR)
  454.       return(0.0);
  455.      else 
  456.      {
  457.       surface_area = 4.0 * PI * pow(equatorial_radius, 2.0);
  458.       hydrosphere_mass = hydrosphere_fraction * surface_area * EARTH_WATER_MASS_PER_AREA;
  459.       water_vapor_in_kg = (0.00000001 * hydrosphere_mass) * exp(Q2_36 * (surface_temp - 288.0));
  460.       fraction = CLOUD_COVERAGE_FACTOR * water_vapor_in_kg / surface_area;
  461.       if (fraction >= 1.0)
  462.            return(1.0);
  463.       else 
  464.            return(fraction);
  465.      }
  466. }
  467.  
  468.  
  469. /*--------------------------------------------------------------------------*/
  470. /*   Given the surface temperature of a planet (in Kelvin), this function   */
  471. /*   returns the fraction of the planet's surface covered by ice.  This is  */
  472. /*   Fogg's eq.24.  See Hart[24] in Icarus vol.33, p.28 for an explanation. */
  473. /*   I have changed a constant from 70 to 90 in order to bring it more in   */
  474. /*   line with the fraction of the Earth's surface covered with ice, which  */
  475. /*   is approximatly .016 (=1.6%).                                          */
  476. /*--------------------------------------------------------------------------*/
  477.  
  478. double ice_fraction(hydrosphere_fraction, surface_temp)
  479. double hydrosphere_fraction, surface_temp;
  480. {
  481.      double temp; 
  482.      
  483.      if (surface_temp > 328.0) 
  484.       surface_temp = 328.0;
  485.      temp = pow(((328.0 - surface_temp) / 90.0),5.0);
  486.      if (temp > (1.5 * hydrosphere_fraction))
  487.       temp = (1.5 * hydrosphere_fraction);
  488.      if (temp >= 1.0)
  489.       return(1.0);
  490.      else 
  491.       return(temp);
  492. }
  493.  
  494.  
  495. /*--------------------------------------------------------------------------*/
  496. /*  This is Fogg's eq.19.  The ecosphere radius is given in AU, the orbital */
  497. /*  radius in AU, and the temperature returned is in Kelvin.            */
  498. /*--------------------------------------------------------------------------*/
  499.  
  500. double eff_temp(ecosphere_radius, orbital_radius, albedo)
  501. double ecosphere_radius, orbital_radius, albedo;
  502. {
  503.      return(sqrt(ecosphere_radius / orbital_radius) * pow((1.0 - albedo) / 0.7,0.25) * EARTH_EFFECTIVE_TEMP);
  504. }
  505.  
  506.  
  507. /*--------------------------------------------------------------------------*/
  508. /*  This is Fogg's eq.20, and is also Hart's eq.20 in his "Evolution of     */
  509. /*  Earth's Atmosphere" article.  The effective temperature given is in     */
  510. /*  units of Kelvin, as is the rise in temperature produced by the          */
  511. /*  greenhouse effect, which is returned.                                   */
  512. /*--------------------------------------------------------------------------*/
  513.  
  514. double green_rise(optical_depth, effective_temp, surface_pressure)
  515. double optical_depth, effective_temp, surface_pressure;
  516. {
  517.      double convection_factor; 
  518.      
  519.      convection_factor = EARTH_CONVECTION_FACTOR * pow((surface_pressure / EARTH_SURF_PRES_IN_MILLIBARS),0.25);
  520.      return(pow((1.0 + 0.75 * optical_depth),0.25) - 1.0) * effective_temp * convection_factor;
  521. }
  522.  
  523.  
  524. /*--------------------------------------------------------------------------*/
  525. /*   The surface temperature passed in is in units of Kelvin.               */
  526. /*   The cloud adjustment is the fraction of cloud cover obscuring each     */
  527. /*   of the three major components of albedo that lie below the clouds.     */
  528. /*--------------------------------------------------------------------------*/
  529.  
  530. double planet_albedo(water_fraction, cloud_fraction, ice_fraction, surface_pressure)
  531. double water_fraction, cloud_fraction, ice_fraction, surface_pressure;
  532. {
  533.      double rock_fraction, cloud_adjustment, components, cloud_contribution,
  534.      rock_contribution, water_contribution, ice_contribution;
  535.      
  536.      rock_fraction = 1.0 - water_fraction - ice_fraction;
  537.      components = 0.0;
  538.      if (water_fraction > 0.0)
  539.       components = components + 1.0;
  540.      if (ice_fraction > 0.0)
  541.       components = components + 1.0;
  542.      if (rock_fraction > 0.0)
  543.       components = components + 1.0;
  544.      cloud_adjustment = cloud_fraction / components;
  545.      if (rock_fraction >= cloud_adjustment)
  546.       rock_fraction = rock_fraction - cloud_adjustment;
  547.      else 
  548.       rock_fraction = 0.0;
  549.      if (water_fraction > cloud_adjustment)
  550.       water_fraction = water_fraction - cloud_adjustment;
  551.      else 
  552.       water_fraction = 0.0;
  553.      if (ice_fraction > cloud_adjustment)
  554.       ice_fraction = ice_fraction - cloud_adjustment;
  555.      else 
  556.       ice_fraction = 0.0;
  557.      cloud_contribution = cloud_fraction * about(CLOUD_ALBEDO,0.2);
  558.      if (surface_pressure == 0.0)
  559.       rock_contribution = rock_fraction * about(AIRLESS_ROCKY_ALBEDO,0.3);
  560.      else 
  561.       rock_contribution = rock_fraction * about(ROCKY_ALBEDO,0.1);
  562.      water_contribution = water_fraction * about(WATER_ALBEDO,0.2);
  563.      if (surface_pressure == 0.0)
  564.       ice_contribution = ice_fraction * about(AIRLESS_ICE_ALBEDO,0.4);
  565.      else 
  566.       ice_contribution = ice_fraction * about(ICE_ALBEDO,0.1);
  567.      return(cloud_contribution + rock_contribution + water_contribution + ice_contribution);
  568. }
  569.  
  570.  
  571. /*--------------------------------------------------------------------------*/
  572. /*   This function returns the dimensionless quantity of optical depth,     */
  573. /*   which is useful in determining the amount of greenhouse effect on a    */
  574. /*   planet.                                                                */
  575. /*--------------------------------------------------------------------------*/
  576.  
  577. double opacity(molecular_weight, surface_pressure)
  578. double molecular_weight, surface_pressure;
  579. {
  580.      double optical_depth; 
  581.      
  582.      optical_depth = 0.0;
  583.      if ((molecular_weight >= 0.0) && (molecular_weight < 10.0))
  584.       optical_depth = optical_depth + 3.0;
  585.      if ((molecular_weight >= 10.0) && (molecular_weight < 20.0))
  586.       optical_depth = optical_depth + 2.34;
  587.      if ((molecular_weight >= 20.0) && (molecular_weight < 30.0))
  588.       optical_depth = optical_depth + 1.0;
  589.      if ((molecular_weight >= 30.0) && (molecular_weight < 45.0))
  590.       optical_depth = optical_depth + 0.15;
  591.      if ((molecular_weight >= 45.0) && (molecular_weight < 100.0))
  592.       optical_depth = optical_depth + 0.05;
  593.      if (surface_pressure >= (70.0 * EARTH_SURF_PRES_IN_MILLIBARS))
  594.       optical_depth = optical_depth * 8.333;
  595.      else 
  596.       if (surface_pressure >= (50.0 * EARTH_SURF_PRES_IN_MILLIBARS))
  597.            optical_depth = optical_depth * 6.666;
  598.       else 
  599.            if (surface_pressure >= (30.0 * EARTH_SURF_PRES_IN_MILLIBARS))
  600.             optical_depth = optical_depth * 3.333;
  601.            else 
  602.             if (surface_pressure >= (10.0 * EARTH_SURF_PRES_IN_MILLIBARS))
  603.              optical_depth = optical_depth * 2.0;
  604.             else 
  605.              if (surface_pressure >= (5.0 * EARTH_SURF_PRES_IN_MILLIBARS))
  606.                   optical_depth = optical_depth * 1.5;
  607.      return(optical_depth);
  608. }
  609.  
  610.  
  611. /*--------------------------------------------------------------------------*/
  612. /*   The temperature calculated is in degrees Kelvin.                       */
  613. /*   Quantities already known which are used in these calculations:         */
  614. /*     planet->molecule_weight                        */
  615. /*     planet->surface_pressure                        */
  616. /*       R_ecosphere                                                        */
  617. /*     planet->a                                */
  618. /*     planet->volatile_gas_inventory                        */
  619. /*     planet->radius                                */
  620. /*     planet->boil_point                            */
  621. /*--------------------------------------------------------------------------*/
  622.  
  623. void iterate_surface_temp(planet)
  624. planet_pointer *planet; 
  625. {
  626.      double surface_temp, effective_temp, greenhouse_rise, previous_temp,
  627.      optical_depth, albedo, water, clouds, ice;
  628.      
  629.      optical_depth = opacity((*planet)->molecule_weight,(*planet)->surface_pressure);
  630.      effective_temp = eff_temp(r_ecosphere,(*planet)->a,EARTH_ALBEDO);
  631.      greenhouse_rise = green_rise(optical_depth,effective_temp,(*planet)->surface_pressure);
  632.      surface_temp = effective_temp + greenhouse_rise;
  633.      previous_temp = surface_temp - 5.0;        /* force the while loop the first time */
  634.      while ((fabs(surface_temp - previous_temp) > 1.0)) {
  635.            previous_temp = surface_temp;
  636.            water = hydrosphere_fraction((*planet)->volatile_gas_inventory,(*planet)->radius);
  637.            clouds = cloud_fraction(surface_temp,(*planet)->molecule_weight,(*planet)->radius,water);
  638.            ice = ice_fraction(water,surface_temp);
  639.            if ((surface_temp >= (*planet)->boil_point) || (surface_temp <= FREEZING_POINT_OF_WATER))
  640.             water = 0.0;
  641.            albedo = planet_albedo(water,clouds,ice,(*planet)->surface_pressure);
  642.            optical_depth = opacity((*planet)->molecule_weight,(*planet)->surface_pressure);
  643.            effective_temp = eff_temp(r_ecosphere,(*planet)->a,albedo);
  644.            greenhouse_rise = green_rise(optical_depth,effective_temp,(*planet)->surface_pressure);
  645.            surface_temp = effective_temp + greenhouse_rise;
  646.       }
  647.      (*planet)->hydrosphere = water;
  648.      (*planet)->cloud_cover = clouds;
  649.      (*planet)->ice_cover = ice;
  650.      (*planet)->albedo = albedo;
  651.      (*planet)->surface_temp = surface_temp;
  652. }
  653.