adddouble.C Source File

00001 
00002        //=======================================================//    _\|/_
00003       //  __  _____           ___                    ___       //      /|\ ~
00004      //  /      |      ^     |   \  |         ^     |   \     //          _\|/_
00005     //   \__    |     / \    |___/  |        / \    |___/    //            /|\ ~
00006    //       \   |    /___\   |  \   |       /___\   |   \   // _\|/_
00007   //     ___/   |   /     \  |   \  |____  /     \  |___/  //   /|\ ~
00008  //                                                       //            _\|/_
00009 //=======================================================//              /|\ ~
00010 
00028 //++ Notes:
00029 //++  Online arguments over rule parameters from input node.
00030 //++
00031 //++  If no time unit is explicitly specified it is derived
00032 //++  from the other units.
00033 //++
00034 //++  The mass of the system can be set independedly from the total
00035 //++  mass of the input N-body system (e.g. when modeling a large
00036 //++  cluster by a smaller N-body system).
00037 //++
00038 //++ Example of usage:
00039 //++  adddouble
00040 //++
00041 //++ See also: addstar
00042 
00044 //                           spz@grape.c.u-tokyo.ac.jp
00045 
00046 #include "double_star.h"
00047 #include "util_io.h"
00048 #include "dstar_to_dyn.h"
00049 #include "kepler.h"
00050 
00051 #ifndef TOOLBOX
00052 
00053 //#define MASS_UPDATE_LIMIT 0.0001
00054 //#define SMA_UPDATE_LIMIT 0.0001
00055 
00056 #define REPORT_ADD_DOUBLE        false
00057 #define REPORT_DEL_DOUBLE        false
00058 #define REPORT_EVOLVE_DOUBLE     false
00059 
00060 //This is a challenge, does it improve bugginess.
00061 
00062 local void randomize_created_binary(dyn *the_binary_node)
00063 {
00064      dyn *d1 = the_binary_node->get_oldest_daughter();
00065      dyn *d2 = d1->get_younger_sister();
00066      real m1 = d1->get_mass();
00067      real m2 = d2->get_mass();
00068      real m_total = m1 + m2;
00069      real sma = the_binary_node->get_starbase()->get_semi();
00070      real ecc = the_binary_node->get_starbase()->get_eccentricity();
00071 
00072      if (REPORT_ADD_DOUBLE) {
00073        cerr << "randomize_created_binary: ";
00074        PRC(m1);PRC(m2);PRC(sma);PRL(ecc);
00075      }
00076 
00077      if (sma<0 || ecc<0 || ecc>1) {
00078         cerr << "randomize_created_binary(): This is not a binary!"<<endl;
00079         return;
00080      }
00081 
00082      kepler johannes;
00083 
00084      real peri = 1; // Default value (unimportant unless ecc = 1).
00085 
00086      // For now, binary phase is random.
00087 
00088      real mean_anomaly = randinter(-PI, PI);
00089 
00090      make_standard_kepler(johannes, 0, m_total, -0.5 * m_total / sma, ecc,
00091                           peri, mean_anomaly);
00092      set_random_orientation(johannes);
00093      johannes.print_all();
00094 }
00095 
00096 
00097 void update_dyn_from_binary_evolution(dyn* the_binary_node,
00098                                       dyn* d1, real dm1_fast,
00099                                                real dm2_fast)
00100 {
00101 
00102     dyn *d2 = NULL;
00103     if (d1==NULL) {
00104         d1 = the_binary_node->get_oldest_daughter();
00105         d2 = d1->get_younger_sister();
00106     }
00107     else
00108         d2 = d1->get_binary_sister();
00109 
00110     real dyn_m1 = get_total_mass(d1) - dm1_fast;
00111     real dyn_m2 = get_total_mass(d2) - dm2_fast;
00112     real dyn_mtot = dyn_m1 + dyn_m2;
00113 
00114     // kepler and binary exists.
00115     // As there is a kepler it should have been initialized.
00116 
00117     //          Needed to think about more carefully.
00118     kepler *johannes = the_binary_node->get_oldest_daughter()->get_kepler();
00119     d1->set_pos(-dyn_m2 * johannes->get_rel_pos() / dyn_mtot);
00120     d1->set_vel(-dyn_m2 * johannes->get_rel_vel() / dyn_mtot);
00121 
00122     d2->set_pos(dyn_m1 * johannes->get_rel_pos() / dyn_mtot);
00123     d2->set_vel(dyn_m1 * johannes->get_rel_vel() / dyn_mtot);
00124     d1->set_mass(dyn_m1);
00125     d2->set_mass(dyn_m2);
00126 
00127 
00128 }
00129 
00130 
00131 /*-----------------------------------------------------------------------------
00132  *  adddouble  -- for all particles, add a star part using "new double()".
00133  *
00134  *  Adds a double_star to a parent node of two leaves, i.e.; the
00135  *  double_star is attached to the center-of-mass node of the
00136  *  two stars.
00137  *  Note that the kepler pointer is attached to the leave nodes.
00138  *  The function checks (too extensively) wheter or not a double_star
00139  *  is really needed. The function is fully self-supporting:
00140  *  call adddouble from any root-node and it finds itself what it needs,
00141  *  (unlike delete_double see dstar_to_kira.C).
00142  *-----------------------------------------------------------------------------
00143  */
00144 
00145 local void do_the_adding(hdyn *bi, real stellar_time,   // bi is binary CM
00146                          binary_type type,
00147                          real sma, real ecc)
00148 {
00149     hdyn* od = bi->get_oldest_daughter();
00150 
00151     if (stellar_time <= 0) {
00152         //cerr << "adding at Time=zero"<<endl;
00153         // Should be the individual binaries time!!
00154         stellar_time = bi->get_starbase()
00155                          ->conv_t_dyn_to_star(od->get_time());
00156                                       // 3/99 was get_system_time()
00157     }
00158 
00159     if (REPORT_ADD_DOUBLE) {
00160         int p = cerr.precision(HIGH_PRECISION);
00161         cerr << "Before adddouble: at " << stellar_time << " to: ";
00162         bi->pretty_print_tree(cerr);
00163         cerr << " a=" << sma << " e=" << ecc << endl;
00164         od->get_kepler()->print_all(cerr);
00165         cerr.precision(p);
00166     }
00167 
00168     // addstar(bi, stellar_time); // Should not be needed.
00169 
00170     int id;
00171     if((od->is_low_level_node() &&
00172         od->get_younger_sister()->is_low_level_node()) &&
00173        (od->get_elder_sister() == NULL) &&
00174        bi->n_leaves()==2 ) {
00175 
00176         if (!has_dstar(od)) {
00177             //bi->get_parent()->get_starbase()->get_element_type()!=Double) {
00178 
00179             //    if (bi->is_parent() && !bi->is_root())
00180             story * old_story = bi->get_starbase()->get_star_story();
00181             bi->get_starbase()->set_star_story(NULL);
00182 
00183             id = bi->get_index();
00184 
00185             // cerr << "Adding binary to "<< id << " at time = "
00186             //      << stellar_time << endl;
00187 
00188             double_star* new_double
00189                 = new_double_star(bi, sma, ecc, stellar_time, id, type);
00190             // synchronize_binary_components(dynamic_cast(double_star*,
00191             //                                            bi->get_starbase()));
00192 
00193             if (REPORT_ADD_DOUBLE) {
00194                 cerr<<"Adddouble: id, t, a, e: "
00195                     <<new_double->get_identity()<<" "
00196                     <<new_double->get_current_time() << " "
00197                     <<new_double->get_semi()<<" "
00198                     <<new_double->get_eccentricity()<<endl;
00199                 put_state(make_state(new_double));
00200             }
00201 
00202             // Give the new binary the old star_story.
00203 
00204             new_double->set_star_story(old_story);
00205 
00206             if (REPORT_ADD_DOUBLE) {
00207                 int p = cerr.precision(HIGH_PRECISION);
00208                 cerr<<"After adddouble: at "<<stellar_time<< " to: ";
00209                 new_double->dump(cerr);
00210                 od->get_kepler()->print_all(cerr);
00211                 put_state(make_state(new_double));
00212                 cerr << endl;
00213                 cerr.precision(p);
00214             }
00215         }
00216         else {
00217             if (REPORT_ADD_DOUBLE)
00218                 cerr << "No double_star to node needed in adddouble"<<endl;
00219         }
00220     }
00221     else {
00222         if (REPORT_ADD_DOUBLE)
00223             cerr << "Sorry, no binary node in adddouble"<<endl;
00224     }
00225 }
00226 
00227 void adddouble(hdyn *b,                                 // b is binary CM
00228                real dyn_time, binary_type type,
00229                bool random_initialization,
00230                real a_min, real a_max,
00231                real e_min, real e_max)
00232 {
00233     if (REPORT_ADD_DOUBLE) {
00234         int p = cerr.precision(HIGH_PRECISION);
00235         cerr<<"adddouble: "<<b<<" at T="<<dyn_time<<endl;
00236         cerr.precision(p);
00237     }
00238 
00239     real stellar_time = b->get_starbase()->conv_t_dyn_to_star(dyn_time);
00240     addstar(b, stellar_time);
00241 
00242     real ecc, sma;
00243     real binary_age=0;
00244 //      binary_type local_type = Unknown_Binary_Type;
00245     real a_const = log(a_max) - log(a_min);
00246     int id;
00247 //      kepler johannes;
00248 
00249     for_all_nodes(hdyn, b, bi) {
00250         if (bi->is_parent() && !bi->is_root()) {
00251 
00252             // Note that b's star_story is DELETED here.
00253             // Is this intentional?
00254 
00255             story * old_story = b->get_starbase()->get_star_story();
00256 //            extract_story_chapter(local_type, sma, ecc, *old_story);
00257             b->get_starbase()->set_star_story(NULL);
00258 //            delete old_story;
00259 
00260             binary_age = max(bi->get_oldest_daughter()->get_starbase()
00261                                ->get_current_time(),
00262                              bi->get_oldest_daughter()->get_binary_sister()
00263                                ->get_starbase()->get_current_time());
00264             // cerr << "times: " << binary_age<<" "<<stellar_time<<endl;
00265             id = bi->get_index();
00266 
00267             if (random_initialization) {
00268 
00269                sma = a_min*exp(randinter(0., a_const));
00270                do {
00271                   ecc = sqrt(randinter(0., 1.));
00272                }
00273                while (ecc<e_min && ecc>=e_max);
00274 
00275             } else { //if (local_type==Unknown_Binary_Type)
00276 
00277               // cerr << "Non-random creation"<<endl;
00278 //               type = local_type;
00279 //            else {
00280 
00281                 if (!bi->get_oldest_daughter()->get_kepler())
00282                     new_child_kepler(bi, dyn_time, MAX_CIRC_ECC);
00283                 else {
00284                     bi->get_oldest_daughter()->get_kepler()->     // Probably
00285                         set_circular_binary_limit(MAX_CIRC_ECC);  // unnecessary
00286                     dyn_to_child_kepler(bi, dyn_time);
00287                 }
00288                 
00289                 sma = b->get_starbase()->
00290                                 conv_r_dyn_to_star(bi->get_oldest_daughter()->
00291                                 get_kepler()->get_semi_major_axis());
00292                 ecc = bi->get_oldest_daughter()->get_kepler()->
00293                                 get_eccentricity();
00294             }
00295 
00296             do_the_adding(bi, binary_age, type, sma, ecc);
00297 
00298             if (random_initialization) {
00299                 randomize_created_binary(bi);
00300                 //            update_dynamical_part_of_binary(bi);
00301                 // The above function no longer exist
00302                 // I hope the following does the expected
00303                 // thing.... (JM)
00304                 cerr << "try update dyn here\n";
00305                 update_dyn_from_binary_evolution(bi);
00306                 cerr << "return update dyn here\n";
00307             }
00308         }
00309     }
00310 }
00311 
00312 #else
00313 
00314 /*-----------------------------------------------------------------------------
00315  *  main  --
00316  *-----------------------------------------------------------------------------
00317  */
00318 main(int argc, char ** argv)
00319     {
00320     int  c;
00321     bool  A_flag = FALSE;
00322     bool  a_flag = FALSE;
00323     bool  E_flag = FALSE;
00324     bool  e_flag = FALSE;
00325     bool  t_flag = FALSE;
00326     bool  M_flag = FALSE;
00327     bool  R_flag = FALSE;
00328     bool  Q_flag = FALSE;
00329     bool  T_flag = FALSE;
00330     bool  p_flag = FALSE;
00331     bool  S_flag = FALSE;
00332     bool  c_flag = FALSE;
00333     bool  random_initialization = false;
00334     real  m_tot = 1;
00335     real  r_vir = 1;
00336     real  t_hc = 1;
00337     real  Q_vir = 0.5;
00338     real  a_min = 1;
00339     real  a_max = 1.e+6;
00340     real  e_min = 0;
00341     real  e_max = 1;
00342     real  t_start = 0;
00343     binary_type type = Detached;
00344 
00345     int random_seed = 0;
00346     char seedlog[64];
00347     char  *comment;
00348     extern char *poptarg;
00349     int  pgetopt(int, char **, char *);
00350 
00351     char * param_string = "A:a:E:e:M:R:Q:T:t:Ss:c:";
00352     check_help();
00353 
00354     while ((c = pgetopt(argc, argv, param_string)) != -1)
00355         switch(c)
00356             {
00357             case 'A': A_flag = true;
00358                       a_max = atof(poptarg);
00359                       break;
00360             case 'a': a_flag = true;
00361                       a_min = atof(poptarg);
00362                       break;
00363             case 'E': E_flag = true;
00364                       e_max = atof(poptarg);
00365                       break;
00366             case 'e': e_flag = true;
00367                       e_min = atof(poptarg);
00368                       break;
00369             case 'M': M_flag = true;
00370                       m_tot = atof(poptarg);
00371                       break;
00372             case 'R': R_flag = true;
00373                       r_vir = atof(poptarg);
00374                       break;
00375             case 'Q': Q_flag = true;
00376                       Q_vir = atof(poptarg);
00377                       break;
00378             case 'T': T_flag = true;
00379                       t_hc = atof(poptarg);
00380                       break;
00381             case 't': t_start = atof(poptarg);
00382                       break;
00383             case 's': random_seed = atoi(poptarg);
00384                       break;
00385             case 'S': S_flag = true;
00386                       break;
00387             case 'c': c_flag = true;
00388                       comment = poptarg;
00389                       break;
00390             case '?': params_to_usage(cerr, argv[0], param_string);
00391                       get_help();
00392                       exit(1);
00393             }
00394 
00395     int actual_seed = srandinter(random_seed);
00396     sprintf(seedlog, "       random number generator seed = %d",actual_seed);
00397 
00398        dyn* b = get_dyn(cin);
00399 
00400        b->log_history(argc, argv);
00401        b->log_comment(seedlog);
00402 
00403        if (!T_flag) t_hc = 21.0 * sqrt(Q_vir/m_tot) * pow(r_vir, 1.5);
00404 
00405        if (S_flag)
00406           if (M_flag && R_flag && T_flag)
00407              b->get_starbase() ->set_stellar_evolution_scaling(m_tot, r_vir,
00408                                                                t_hc);
00409            else {
00410                cerr << "Proper scaling not set." << endl; cerr
00411                     << "Use the -M -R an -T options to set the scaling for "
00412                     << "the mass, \nthe size and the time" << endl; exit(1);
00413            }
00414 
00415           if (A_flag || a_flag || E_flag || e_flag)
00416            random_initialization=true;
00417        adddouble(b, t_start, type, random_initialization,
00418                  a_min, a_max, e_min, e_max);
00419        put_dyn(cout, *b);       
00420        delete b;
00421     }
00422 
00423 #endif
00424 
00425 /* endof: adddouble.c */