makedisk.C Source File

00001 
00023 
00024 #define G       6.67e-8         // cgs units
00025 #define PC      3.086e18        // cm
00026 #define MSUN    1.989e33        // g
00027 #define YR      3.156e7         // s
00028 
00029 #include "dyn.h"
00030 
00031 #ifdef TOOLBOX
00032 
00033 local void mkdisk(dyn* b, int n,
00034                   real m_central, real m_disk,
00035                   real r_inner, real r_outer, real radius, real v_disp)
00036 {
00037     real pmass = m_disk / n;
00038     v_disp *= sqrt(3.0);        // so a uniform distribution has proper v_disp;
00039 
00040     vector cmpos = 0, cmvel = 0;
00041 
00042     for_all_daughters(dyn, b, bi) {
00043 
00044         if (bi->get_elder_sister() == NULL) {
00045 
00046             // Central object.
00047 
00048             bi->set_mass(m_central);
00049             bi->set_pos(0);
00050             bi->set_vel(0);
00051 
00052             // Particle 0 is a massive black hole.  We can include
00053             // tidal destruction within the kira collision framework
00054             // by giving the hole a radius such that it has the same
00055             // "density" as the other particles in the system.  This
00056             // works so long as the radius of the hole thus defined is
00057             // much greater than the radius of any other particle in
00058             // the system.
00059 
00060             putrq(bi->get_dyn_story(), "R_eff",
00061                   radius * pow(m_central/pmass, 1.0/3));
00062 
00063             // Note that R_eff is meaningful only to hdyn tools.
00064 
00065         } else {
00066 
00067             // Disk particle; disk lies in the (x-y) plane, mean motion
00068             // in the positive sense.
00069 
00070             bi->set_mass(pmass);        
00071 
00072             real r = sqrt(r_inner*r_inner
00073                              + randinter(0,1)
00074                                   * (r_outer*r_outer - r_inner*r_inner));
00075             real theta = randinter(0, 2*M_PI);
00076             bi->set_pos(vector(r*cos(theta), r*sin(theta), 0));
00077 
00078             real v_orb = sqrt(m_central/r);
00079             bi->set_vel(vector(-v_orb*sin(theta) + v_disp*randinter(-1,1),
00080                                 v_orb*cos(theta) + v_disp*randinter(-1,1),
00081                                                    v_disp*randinter(-1,1)));
00082 
00083             // All particles have the same mass and radius, and hence
00084             // the same density.
00085 
00086             putrq(bi->get_dyn_story(), "R_eff", radius);
00087 
00088             cmpos += pmass * bi->get_pos();
00089             cmvel += pmass * bi->get_vel();
00090         }
00091     }
00092 
00093     b->set_mass(m_central + m_disk);
00094 
00095     // Force the system CM to be at rest at the origin.
00096 
00097     cmpos /= (m_central + m_disk);
00098     cmvel /= (m_central + m_disk);
00099 
00100     for_all_daughters(dyn, b, bi) {
00101         bi->inc_pos(-cmpos);
00102         bi->inc_vel(-cmvel);
00103     }
00104 }
00105 
00106 #define  SEED_STRING_LENGTH  60
00107 
00108 main(int argc, char ** argv)
00109 {
00110     bool c_flag = false;
00111     bool n_flag = false;
00112     bool s_flag = false;
00113     bool i_flag = false;
00114     bool l_flag = false;
00115     bool o_flag = false;
00116     bool verbose = false;
00117 
00118     int  i;
00119     int  n = 0;
00120     int  input_seed, actual_seed;
00121 
00122     real m_central = 0, m_disk = 0;
00123     real r_inner = 0, r_outer = 0, radius = 0;
00124     real v_disp = 0;
00125 
00126     char  *comment;
00127     char  seedlog[SEED_STRING_LENGTH];
00128 
00129     check_help();
00130 
00131     extern char *poptarg;
00132     int c;
00133     char* param_string = "c:il:m:M:n:or:R:s:v:V";
00134 
00135     while ((c = pgetopt(argc, argv, param_string)) != -1)
00136         switch(c)
00137             {
00138             case 'c': c_flag = true;
00139                       comment = poptarg;
00140                       break;
00141             case 'i': i_flag = true;
00142                       break;
00143             case 'l': radius = atof(poptarg);
00144                       break;
00145             case 'm': m_disk = atof(poptarg);
00146                       break;
00147             case 'M': m_central = atof(poptarg);
00148                       break;
00149             case 'n': n = atoi(poptarg);
00150                       break;
00151             case 'o': o_flag = true;
00152                       break;
00153             case 'r': r_inner = atof(poptarg);
00154                       break;
00155             case 'R': r_outer = atof(poptarg);
00156                       break;
00157             case 's': s_flag = true;
00158                       input_seed = atoi(poptarg);
00159                       break;
00160             case 'v': v_disp = atof(poptarg);
00161                       break;
00162             case 'V': verbose = true;
00163                       break;
00164             case '?': params_to_usage(cerr, argv[0], param_string);
00165                       get_help();
00166                       exit(1);
00167             }            
00168     
00169     if (n <= 0) {
00170         cerr << "mkdisk: must specify number of particles with -n#\n";
00171         exit(1);
00172     }
00173     
00174     if (m_disk <= 0) {
00175         cerr << "mkdisk: must specify disk mass with -m\n";
00176         exit(1);
00177     }
00178     
00179     if (m_central <= 0) {
00180         cerr << "mkdisk: must specify central mass with -M\n";
00181         exit(1);
00182     }
00183     
00184     if (r_inner <= 0) {
00185         cerr << "mkdisk: must specify inner disk radius with -r\n";
00186         exit(1);
00187     }
00188     
00189     if (r_outer <= 0) {
00190         cerr << "mkdisk: must specify outer disk radius with -R\n";
00191         exit(1);
00192     }
00193     
00194     if (r_outer <= r_inner) {
00195         cerr << "mkdisk: must specify r_outer > r_inner\n";
00196         exit(1);
00197     }
00198 
00199     if (verbose) {
00200 
00201         real t_unit = sqrt(pow(PC,3)/(G*1.e6*MSUN));            // unit: s
00202         cerr << "[M] = 1.e6 solar masses" << endl;
00203         cerr << "[L] = 1 pc" << endl;
00204         cerr << "[T] = " << t_unit/YR << " yr" << endl;
00205         real v_unit = 1.e-5*PC/t_unit;
00206         cerr << "[v] = " << v_unit << " km/s" << endl;
00207 
00208         cerr << endl;
00209 
00210         PRL(n);
00211         PRC(m_central), PRL(m_disk);
00212         PRC(r_inner), PRC(r_outer), PRL(radius);
00213         cerr << "v_disp = " << v_disp
00214              << " = " << v_disp*v_unit << " km/s" << endl;
00215 
00216         cerr << endl;
00217 
00218         real p_inner = 2*M_PI*sqrt(pow(r_inner,3)/m_central);
00219         cerr << "orbital period at disk inner edge = " << p_inner
00220              << " = " << p_inner*t_unit/YR << " yr" << endl;
00221 
00222         real p_outer = 2*M_PI*sqrt(pow(r_outer,3)/m_central);
00223         cerr << "orbital period at disk outer edge = " << p_outer
00224              << " = " << p_outer*t_unit/YR << " yr" << endl;
00225 
00226         cerr << endl;
00227 
00228         real v_inner = sqrt(m_central/r_inner);
00229         cerr << "orbital speed at disk inner edge = " << v_inner
00230              << " = " << v_inner*v_unit << " km/s" << endl;
00231 
00232         real v_outer = sqrt(m_central/r_outer);
00233         cerr << "orbital speed at disk outer edge = " << v_outer
00234              << " = " << v_outer*v_unit << " km/s" << endl;
00235 
00236         if (radius > 0) {
00237             real v_esc = sqrt(2*(m_disk/n)/radius);
00238             cerr << "cloud escape speed = " << v_esc
00239                  << " = " << v_esc*v_unit << " km/s" << endl;
00240         }
00241     }
00242     
00243     dyn *b, *by, *bo;
00244     b = new dyn();                              // root node
00245 
00246     bo = new dyn();                             // central mass
00247     if (i_flag) 
00248         bo->set_label(0);
00249     b->set_oldest_daughter(bo);
00250     bo->set_parent(b);
00251 
00252     for (i = 1; i <= n; i++) {                  // disk particles
00253         by = new dyn();
00254         if (i_flag) by->set_label(i);
00255         by->set_parent(b);
00256         bo->set_younger_sister(by);
00257         by->set_elder_sister(bo);
00258         bo = by;
00259     }
00260 
00261     if (c_flag == true) b->log_comment(comment);
00262 
00263     b->log_history(argc, argv);
00264 
00265     if (s_flag == false) input_seed = 0;
00266     actual_seed = srandinter(input_seed);
00267 
00268     if (o_flag) cerr << "mkdisk: random seed = " << actual_seed << endl;
00269     sprintf(seedlog, "       random number generator seed = %d",actual_seed);
00270     b->log_comment(seedlog);
00271 
00272     mkdisk(b, n,
00273            m_central, m_disk,
00274            r_inner, r_outer, radius, v_disp);
00275 
00276     put_node(cout, *b);
00277 }
00278 
00279 #endif