makemass.C Source File

00001 
00002        //=======================================================//    _\|/_
00003       //  __  _____           ___                    ___       //      /|\ ~
00004      //  /      |      ^     |   \  |         ^     |   \     //          _\|/_
00005     //   \__    |     / \    |___/  |        / \    |___/    //            /|\ ~
00006    //       \   |    /___\   |  \   |       /___\   |   \   // _\|/_
00007   //     ___/   |   /     \  |   \  |____  /     \  |___/  //   /|\ ~
00008  //                                                       //            _\|/_
00009 //=======================================================//              /|\ ~
00010 
00030 //++ Note: The conversion factor for scaling between dynamical and stellar masss
00031 //++        is properly set in the output snapshot.
00032 
00033 //              Steve McMillan, July 1996
00034 //              Simon Portegies Zwart, Tokyo, December 1997
00035 
00036 // need to repair bugs, see PJT comments
00037 
00038 #include "node.h"
00039 
00040 #define  MAXIMUM_ZAMS_MASS 100
00041 #define  SEED_STRING_LENGTH  256
00042 static char  tmp_string[SEED_STRING_LENGTH];
00043 
00044 #ifndef TOOLBOX
00045 
00046 
00047 //see: Eggleton, Fitchet & Tout 1989, ApJ 347, 998
00048 local real mf_Miller_Scalo(real m_lower, real m_upper) {
00049 
00050     real m, rnd;
00051     do {
00052         rnd = randinter(0,1);
00053         m = 0.19*rnd
00054             / (pow(1-rnd, 0.75) + 0.032*pow(1-rnd, 0.25));
00055     }
00056     while(m_lower>m || m>m_upper);
00057     return m;
00058 }
00059 
00060 //Tapered power-law (from Guido de Marchi private communication Oct. 2001)
00061 //dN/dLogm= m^(-1.35)*(1-exp(-m/0.25)^2.4)
00062 local real tapered_power_law(real m, 
00063                              const real x1 = -1.35, 
00064                              const real x2 = 2.4, 
00065                              const real mbreak = 0.25) {
00066 
00067   real dNdm = pow(m, x1)*(1-pow(exp(-m/mbreak), x2));
00068   return dNdm;
00069 }
00070 
00071 local real mf_GdeMarchi(real m_lower, real m_upper) {
00072 
00073 //dN/dLogm= m^(-1.35)*(1-exp(-m/0.15)^2.4)
00074 // for now we do it the poor way by Monte Carlo
00075 
00076   real x1 = -2.35;
00077   real x2 = 2.40;
00078   real mbreak = 0.15;
00079     real dNdm_min = tapered_power_law(m_lower, x1, x2, mbreak);
00080     real dNdm_max = tapered_power_law(m_upper, x1, x2, mbreak);
00081 //    PRC(dNdm_min);PRL(dNdm_max);
00082     real m, rnd;
00083     real dNdm_try, dNdm;
00084     do {
00085         rnd = randinter(0,1);
00086         m = m_lower + rnd*(m_upper-m_lower);
00087         dNdm_try = randinter(dNdm_min,dNdm_max);
00088         dNdm = tapered_power_law(m, x1, x2, mbreak);
00089     }
00090     while(dNdm_try>dNdm);
00091     return m;
00092 }
00093 
00094 // see: de la Fuente Marcos, Aarseth, Kiseleva, Eggleton
00095 // 1997, in Docobo, Elipe, McAlister (eds.), Visual Double
00096 // Stars: Formation, dynamics and Evolutionary Tracks, KAP: ASSL
00097 // Series vol. 223,  165
00098 local real mf_Scalo(real m_lower, real m_upper) {
00099     real m, rnd;
00100     do {
00101         rnd = randinter(0,1);
00102         m = 0.3*rnd
00103             / pow(1-rnd, 0.55);
00104     }
00105     while(m_lower>m || m>m_upper);
00106     return m;
00107 }
00108 
00109 local real Kroupa_Tout_Gilmore(real m_lower, real m_upper) {
00110 
00111     real m, rnd;
00112     do {
00113         rnd = randinter(0,1);
00114         m = 0.08 + (0.19*pow(rnd, 1.55) + 0.05*pow(rnd, 0.6))
00115             /        pow(1-rnd, 0.58);
00116     }
00117     while(m_lower>m || m>m_upper);
00118     return m;
00119 }
00120 
00121 real get_random_stellar_mass(real m_lower, real m_upper,
00122                              mass_function mf, real exponent) {
00123 
00124     real m;
00125     switch(mf) {
00126        case Equal_Mass:
00127           if (m_lower==0) {
00128              cerr << "get_random_stellar_mass:"<<endl;
00129              cerr << "unambiguous choise of Equal_Mass."<<endl;
00130              cerr << "Use -m option to set fixed stellar mass."<<endl;
00131              exit(1);
00132           }
00133           m =  m_lower;
00134                break;
00135        case mf_Power_Law:
00136           m =  general_power_law(m_lower, m_upper, exponent);
00137                break;
00138        case Miller_Scalo:
00139            m = mf_Miller_Scalo(m_lower, m_upper);
00140                break;
00141        case Scalo:
00142           m = mf_Scalo(m_lower, m_upper);
00143               break;
00144        case Kroupa: // Kroupa, Tout & Gilmore 1993, MNRAS 262, 545
00145           m = Kroupa_Tout_Gilmore(m_lower, m_upper);
00146               break;
00147        case Tapered_Power_Law:
00148           m = mf_GdeMarchi(m_lower, m_upper);
00149               break;
00150        default:
00151           cerr << "WARNING: \n"
00152                << "        real get_random_stellar_mass:\n"
00153                << "        Mass function parameters not properly defined.\n";
00154           exit(1);
00155     }
00156 
00157     return m;
00158 }
00159 
00160 char* type_string(mass_function mf) {
00161 
00162     local char  mf_name[SEED_STRING_LENGTH];    // permanent
00163     switch(mf) {
00164 
00165        case Equal_Mass:
00166             sprintf(mf_name, "Equal_Mass");
00167             break;      
00168        case mf_Power_Law:
00169             sprintf(mf_name, "Power_Law");
00170             break;      
00171        case Miller_Scalo:
00172             sprintf(mf_name, "Miller_Scalo");
00173             break;      
00174        case Scalo:
00175             sprintf(mf_name, "Scalo");
00176             break;      
00177        case Kroupa:
00178             sprintf(mf_name, "Kroupa");
00179             break;
00180        case Tapered_Power_Law:
00181             sprintf(mf_name, "Tapered_Power_Law");
00182             break;
00183        default:
00184             sprintf(mf_name, "Unknown");
00185             break;
00186     }
00187     return mf_name;
00188 }
00189 
00190 mass_function extract_mass_function_type_string(char* type_string) {
00191 
00192      mass_function type = Unknown_MF;
00193 
00194      if (!strcmp(type_string, "Equal_Mass"))
00195         type = Equal_Mass;
00196      else if (!strcmp(type_string, "Power_Law"))
00197         type = mf_Power_Law;
00198      else if (!strcmp(type_string, "Miller_Scalo"))
00199         type = Miller_Scalo;
00200      else if (!strcmp(type_string, "Scalo"))
00201         type = Scalo;
00202      else if (!strcmp(type_string, "Kroupa"))
00203         type = Kroupa;
00204      else if (!strcmp(type_string, "Tapered_Power_Law"))
00205         type = Tapered_Power_Law;
00206      else if (!strcmp(type_string, "Unknown"))
00207         type = Unknown_MF;
00208      else {
00209          cerr << "No proper mass function indicated in mkmass."<<endl;
00210          exit(-1);
00211      }
00212 
00213      return type;
00214    }
00215 
00216 real general_power_law(real lowerl, real upperl, real exponent) {
00217 
00218     real x, norm;
00219 
00220     // Normalizing factor.
00221     if (exponent == -1)
00222         norm = log(upperl/lowerl);
00223     else
00224         norm = pow(upperl/lowerl, 1+exponent) - 1;
00225 
00226     if (exponent == -1)
00227         x = lowerl*exp(norm*randinter(0,1));
00228     else
00229         x = lowerl*pow(norm*randinter(0,1) + 1, 1/(1+exponent));
00230 
00231     return x;
00232 }
00233 
00234 #else
00235 
00236 typedef  struct
00237 {
00238     node* str;
00239     real  mass;
00240 } nm_pair, *nm_pair_ptr;
00241 
00242 //-----------------------------------------------------------------------------
00243 //  compare_mass  --  compare the masses of two particles
00244 //-----------------------------------------------------------------------------
00245 
00246 local int compare_mass(const void * pi, const void * pj)
00247 {
00248     if (((nm_pair_ptr) pi)->mass < ((nm_pair_ptr) pj)->mass)
00249         return(1);
00250     else if (((nm_pair_ptr)pi)->mass > ((nm_pair_ptr)pj)->mass)
00251         return(-1);
00252     else
00253         return(0);
00254 }
00255 
00256 local void mkmass(node* b, mass_function mf,
00257                   real m_lower, real m_upper,
00258                   real exponent, real total_mass, bool renumber_stars) {
00259 
00260     real m, m_sum = 0;
00261     int n=0;
00262     for_all_daughters(node, b, bi) {
00263         m = get_random_stellar_mass(m_lower, m_upper, mf, exponent);
00264         n++;
00265         bi->set_mass(m);
00266         m_sum += bi->get_mass();
00267     }
00268 
00269     b->set_mass(m_sum);
00270 
00271     // Renumber the stars in order of mass.
00272     // Highest mass gets smallest number (strange choise, but).
00273     if (renumber_stars) {
00274       int istart = 1;
00275       bool M_flag = true;
00276       renumber(b, istart, M_flag);
00277     }
00278 
00279 #if 0
00280       nm_pair_ptr nm_table = new nm_pair[n];
00281       if (nm_table == NULL) {
00282         cerr << "mkmass: "
00283              << "not enough memory left for nm_table\n";
00284         return;
00285       }
00286 
00287       int i=0;
00288       for_all_daughters(node, b, bi) {
00289         nm_table[i].str = bi;
00290         nm_table[i].mass = bi->get_mass();
00291         i++;
00292       }
00293 
00294       qsort((void *)nm_table, (size_t)n, sizeof(nm_pair), compare_mass);
00295 
00296       for (i=0; i<n; i++) {
00297         nm_table[i].str->set_index(i+1);  // Ok, lets number from 1 to n
00298       }
00299 
00300       delete []nm_table;
00301     }
00302 #endif
00303 
00304     if (total_mass > 0) {
00305         real m_factor = total_mass/m_sum;
00306         for_all_daughters(node, b, bi)
00307             bi->set_mass(m_factor*bi->get_mass());
00308         b->set_mass(total_mass);
00309     }
00310 
00311     if(m_lower>=m_upper)
00312         sprintf(tmp_string,
00313                 "         %s mass function, total mass = %8.2f",
00314                 type_string(Equal_Mass), m_sum);
00315     else
00316         sprintf(tmp_string,
00317                 "       %s mass function, total mass = %8.2f Solar",
00318                 type_string(mf), m_sum);
00319     b->log_comment(tmp_string);
00320     cerr << "Mass function is " << type_string(mf) <<endl;
00321 
00322 }
00323 
00324 void main(int argc, char ** argv)
00325 {
00326     bool F_flag   = false;                        // Input mf via string
00327     mass_function mf = mf_Power_Law;             // Default = Power-law
00328     char *mfc = new char[64];
00329     real m_lower  = 1, m_upper = 1;             // Default = equal masses
00330     bool x_flag   = false;
00331     real exponent = -2.35;                      // Salpeter mass function
00332                                                 // If no exponent is given
00333                                                 // Miller & Scalo is used
00334     real m_total = -1;                          // Don't rescale
00335 
00336     bool renumber_stars = false;                // renumber stellar index
00337                                                 // from high to low mass
00338 
00339     int random_seed = 0;
00340     char seedlog[64];
00341 
00342     check_help();
00343 
00344     extern char *poptarg;
00345     int c;
00346     char* param_string = "E:e:iF:f:L:l:M:m:s:U:u:X:x:";
00347 
00348     while ((c = pgetopt(argc, argv, param_string)) != -1)
00349         switch(c) {
00350 
00351             case 'F': F_flag = true;
00352                       mfc = poptarg;
00353                       break;
00354             case 'f': mf = (mass_function)atoi(poptarg);
00355                       break;
00356             case 'i': renumber_stars = true;
00357                       break;
00358             case 'E':
00359             case 'e':
00360             case 'X':
00361             case 'x': x_flag = true;
00362                       exponent = atof(poptarg);
00363                       break;
00364             case 'L':
00365             case 'l': m_lower = atof(poptarg);
00366                       break;
00367             case 'M':
00368             case 'm': m_total = atof(poptarg);
00369                       break;
00370             case 's': random_seed = atoi(poptarg);
00371                       break;
00372             case 'U':
00373             case 'u': m_upper = atof(poptarg);
00374                       break;
00375             case '?': params_to_usage(cerr, argv[0], param_string);
00376                       get_help();
00377                       exit(1);
00378         }
00379 
00380     if (m_lower <= 0 ||
00381         m_upper <= 0 ||
00382         m_lower > m_upper)
00383         err_exit("mkmass: Illegal mass limits");
00384 
00385     if (F_flag)
00386       mf = extract_mass_function_type_string(mfc);
00387     delete mfc;
00388 
00389     node* b;
00390     b = get_node(cin);
00391 
00392     b->log_history(argc, argv);
00393 
00394     int actual_seed = srandinter(random_seed);
00395 
00396     sprintf(seedlog, "       random number generator seed = %d",actual_seed);
00397     b->log_comment(seedlog);
00398 
00399     mkmass(b, mf, m_lower, m_upper, exponent, m_total, renumber_stars); 
00400     real initial_mass = getrq(b->get_log_story(), "initial_mass");
00401 
00402     if (initial_mass > -VERY_LARGE_NUMBER)
00403         putrq(b->get_log_story(), "initial_mass", b->get_mass(),
00404               HIGH_PRECISION);
00405 
00406     // Add stellar conversion data, so add_star will work.
00407 
00408     real m_sum = b->get_mass();
00409     real old_mtot = b->get_starbase()->conv_m_dyn_to_star(1);
00410     if(old_mtot<=0) {
00411         real old_r_vir= b->get_starbase()->conv_r_star_to_dyn(1);
00412         real old_t_vir= b->get_starbase()->conv_t_star_to_dyn(1);
00413         b->get_starbase()->set_stellar_evolution_scaling(m_sum,
00414                                                          old_r_vir,
00415                                                          old_t_vir);
00416     }
00417     put_node(cout, *b);
00418     rmtree(b);
00419 
00420 }
00421 #endif
00422