list_sdumb.C Source File

00001 
00005 
00006 
00007 #include "single_star.h"
00008 #include "sstar_to_dyn.h"
00009 #include "dyn.h"
00010 
00011 #ifdef TOOLBOX
00012 
00013 #if 0
00014 local void combine_ubvri(real Up, real Bp, real Vp, real Rp, real Ip,
00015                          real Us, real Bs, real Vs, real Rs, real Is,
00016                          real &U, real &B, real &V, real &R, real &I) {
00017 
00018     real ln10g25 = 2.3025851/2.5;
00019 
00020     U = -2.5 * log10(exp(-ln10g25*Up) + exp(-ln10g25*Us));
00021     B = -2.5 * log10(exp(-ln10g25*Bp) + exp(-ln10g25*Bs));
00022     V = -2.5 * log10(exp(-ln10g25*Vp) + exp(-ln10g25*Vs));
00023     R = -2.5 * log10(exp(-ln10g25*Rp) + exp(-ln10g25*Rs));
00024     I = -2.5 * log10(exp(-ln10g25*Ip) + exp(-ln10g25*Is));
00025 
00026 }
00027 
00028 #endif
00029 
00030 local void get_UBVRI_of_star(dyn *bi, vector pos,
00031         real &Us, real &Bs, real &Vs, real &Rs, real &Is) {
00032 
00033 
00034   // To solar radii
00035   //  vector pos = bi->get_pos() - dc_pos;
00036   pos[0] = bi->get_starbase()->conv_r_dyn_to_star(pos[0]);
00037   pos[1] = bi->get_starbase()->conv_r_dyn_to_star(pos[1]);
00038   pos[2] = bi->get_starbase()->conv_r_dyn_to_star(pos[2]);
00039 
00040   real time = bi->get_starbase()->conv_t_dyn_to_star(bi->get_system_time());
00041 
00042   // And now to parsec
00043   real Rsun_per_parsec = cnsts.parameters(solar_radius)
00044                        / cnsts.parameters(parsec);
00045   pos[0] *= Rsun_per_parsec;
00046   pos[1] *= Rsun_per_parsec;
00047   pos[2] *= Rsun_per_parsec;
00048   
00049 
00050      star_type_spec tpe_class = NAC;
00051      spectral_class star_class;
00052      stellar_type stype = NAS;
00053      stellar_type_summary sstype = ZAMS;
00054      real t_cur, m_rel, m_env, m_core, mco_core, T_eff, L_eff, p_rot, b_fld;
00055      real t_rel=0, R_eff=0;
00056      real M_tot;
00057      if (bi->get_use_sstar()) {
00058         stype = bi->get_starbase()->get_element_type();
00059         M_tot  = bi->get_starbase()->conv_m_dyn_to_star(bi->get_mass());
00060         t_cur = bi->get_starbase()->get_current_time();
00061         t_rel = bi->get_starbase()->get_relative_age();
00062         T_eff = bi->get_starbase()->temperature();
00063         L_eff = bi->get_starbase()->get_luminosity();
00064         star_class = get_spectral_class(T_eff);
00065         R_eff = bi->get_starbase()->get_effective_radius();
00066         ltm_to_ubvri(log10(L_eff), log10(T_eff), M_tot,
00067                      Us, Bs, Vs, Rs, Is);
00068 
00069      }
00070      else if (bi->get_star_story()) {
00071 
00072        extract_story_chapter(stype, t_cur, t_rel, 
00073                              m_rel, m_env, m_core, mco_core,
00074                              T_eff, L_eff, p_rot, b_fld,
00075                              *bi->get_star_story());
00076 
00077        M_tot = m_env + m_core;
00078        sstype = summarize_stellar_type(stype);
00079        star_class = get_spectral_class(T_eff);
00080        
00081        ltm_to_ubvri(log10(L_eff), log10(T_eff), M_tot,
00082                      Us, Bs, Vs, Rs, Is);
00083        
00084        if (find_qmatch(bi->get_star_story(), "Class"))
00085          tpe_class = extract_stellar_spec_summary_string(
00086              getsq(bi->get_star_story(), "Class"));
00087        if (L_eff>0)
00088           R_eff = pow(T_eff/cnsts.parameters(solar_temperature), 2)
00089                 * sqrt(L_eff);
00090      }
00091      else {
00092        cout << "    No stellar information found: " << endl;
00093        return;
00094      }
00095 }
00096 
00097 local void print_star(dyn *bi, vector pos, vector vel,
00098         real &Up, real &Bp, real &Vp, real &Rp, real &Ip) { 
00099 
00100 
00101   // To solar radii
00102   //  vector pos = bi->get_pos() - dc_pos;
00103   pos[0] = bi->get_starbase()->conv_r_dyn_to_star(pos[0]);
00104   pos[1] = bi->get_starbase()->conv_r_dyn_to_star(pos[1]);
00105   pos[2] = bi->get_starbase()->conv_r_dyn_to_star(pos[2]);
00106 
00107   real time = bi->get_starbase()->conv_t_dyn_to_star(bi->get_system_time());
00108 
00109   // And now to parsec
00110   real Rsun_per_parsec = cnsts.parameters(solar_radius)
00111                        / cnsts.parameters(parsec);
00112   pos[0] *= Rsun_per_parsec;
00113   pos[1] *= Rsun_per_parsec;
00114   pos[2] *= Rsun_per_parsec;
00115 
00116 //  real to_Rsun_Myr = cnsts.physics(km_per_s) * cnsts.physics(Myear)
00117 //               / cnsts.parameters(solar_radius);
00118                        
00119 //      real to_Rsun_Myr = cnsts.physics(km_per_s) * cnsts.physics(Myear)
00120 //                     / cnsts.parameters(solar_radius);
00121 //                       
00122 //      real to_dyn      = bi->get_starbase()->conv_r_star_to_dyn(1)
00123 //                       / bi->get_starbase()->conv_t_star_to_dyn(1);
00124 //      vel = vel/(to_Rsun_Myr * to_dyn);
00125                        
00126   real to_km_per_second = cnsts.parameters(solar_radius)
00127                         / (cnsts.physics(km_per_s) * cnsts.physics(Myear));
00128   real to_star      = bi->get_starbase()->conv_r_dyn_to_star(1)/
00129                       bi->get_starbase()->conv_t_dyn_to_star(1);
00130   to_km_per_second = to_star*to_km_per_second;
00131 
00132 //      PRC(to_star);PRL(to_km_per_second);
00133   vel[0] *= to_km_per_second;
00134   vel[1] *= to_km_per_second;
00135   vel[2] *= to_km_per_second;
00136 
00137      star_type_spec tpe_class = NAC;
00138      spectral_class star_class;
00139      stellar_type stype = NAS;
00140      stellar_type_summary sstype = ZAMS;
00141      real t_cur, m_rel, m_env, m_core, mco_core, T_eff, L_eff, p_rot, b_fld;
00142      real t_rel=0, R_eff=0;
00143      real M_tot, Us, Bs, Vs, Rs, Is;    
00144      if (bi->get_use_sstar()) {
00145         stype = bi->get_starbase()->get_element_type();
00146         M_tot  = bi->get_starbase()->conv_m_dyn_to_star(bi->get_mass());
00147         t_cur = bi->get_starbase()->get_current_time();
00148         t_rel = bi->get_starbase()->get_relative_age();
00149         T_eff = bi->get_starbase()->temperature();
00150         L_eff = bi->get_starbase()->get_luminosity();
00151         star_class = get_spectral_class(T_eff);
00152         R_eff = bi->get_starbase()->get_effective_radius();
00153         ltm_to_ubvri(log10(L_eff), log10(T_eff), M_tot,
00154                      Us, Bs, Vs, Rs, Is);
00155 
00156      }
00157      else if (bi->get_star_story()) {
00158 
00159        extract_story_chapter(stype, t_cur, t_rel, 
00160                              m_rel, m_env, m_core, mco_core,
00161                              T_eff, L_eff, p_rot, b_fld,
00162                              *bi->get_star_story());
00163 
00164        M_tot = m_env + m_core;
00165        sstype = summarize_stellar_type(stype);
00166        star_class = get_spectral_class(T_eff);
00167        
00168        ltm_to_ubvri(log10(L_eff), log10(T_eff), M_tot,
00169                      Us, Bs, Vs, Rs, Is);
00170        
00171        if (find_qmatch(bi->get_star_story(), "Class"))
00172          tpe_class = extract_stellar_spec_summary_string(
00173              getsq(bi->get_star_story(), "Class"));
00174        if (L_eff>0)
00175           R_eff = pow(T_eff/cnsts.parameters(solar_temperature), 2)
00176                 * sqrt(L_eff);
00177      }
00178      else {
00179        cout << "    No stellar information found for: ";
00180        bi->pretty_print_node(cout);
00181        return;
00182      }
00183 
00184      real U, B, V, R, I;
00185      combine_ubvri(Up, Bp, Vp, Rp, Ip,
00186                    Us, Bs, Vs, Rs, Is,
00187                    U, B, V, R, I);
00188 
00189      cout << " type= " << stype << " m= " << M_tot << " R= " << R_eff
00190           << " L= " << L_eff 
00191           << " T_eff= " << T_eff 
00192           << " r= "  << pos[0] << " " << pos[1] << " " << pos[2] 
00193           << " v= "  << vel[0] << " " << vel[1] << " " << vel[2] 
00194           << " ubvri= "<<  U << " " << B << " " << V << " " << R << " "
00195           << I  << " :: ";
00196 
00197      Up=U;
00198      Bp=B;
00199      Vp=V;
00200      Ip=I;
00201 }
00202 
00203 local int print_binary_recursive(dyn *b, vector dc_pos, vector dc_vel,
00204         real &U, real &B, real &V, real &R, real &I) { 
00205 
00206   int nb = 0;
00207   if (b->get_oldest_daughter()) {
00208     
00209     vector r_com  = b->get_pos() - dc_pos;
00210     vector v_com  = b->get_vel() - dc_vel;
00211 //    vector r_com  = dyn_something_relative_to_root(b, &dyn::get_pos) -dc_pos;
00212 //    vector v_com  = dyn_something_relative_to_root(b, &dyn::get_vel) -dc_vel;
00213 
00214     real m_tot = b->get_starbase()->conv_m_dyn_to_star(b->get_mass());
00215 
00216     for_all_daughters(dyn, b, bb)
00217       if (bb->n_leaves() >= 2) {
00218         cout << "\nBinary: ";
00219         U=B=V=R=I=VERY_LARGE_NUMBER;
00220         nb += print_binary_recursive(bb, dc_pos, dc_vel, U, B, V, R, I);
00221       }
00222       else {
00223         if (bb->get_parent()==bb->get_root()) {
00224           U=B=V=R=I=VERY_LARGE_NUMBER;
00225           cout << "\nStar:: ";
00226         }
00227         print_star(bb, bb->get_pos()-r_com, bb->get_vel()-v_com,
00228                    U, B, V, R, I);
00229 //      if (bb->get_parent()==bb->get_root())
00230 //         cerr << endl;
00231       }
00232 
00233   }
00234   return nb;
00235 }
00236 
00237 local void print_integrated_cluster(dyn *b, vector dc_pos,
00238                                     real r_min, real r_max,
00239                                     real v_min, real v_max,
00240         real &U, real &B, real &V, real &R, real &I,
00241         int project) { 
00242 
00243 
00244   int nb = 0;
00245   real rcom;
00246   vector r;
00247   real m_tot = 0;
00248   if (b->get_oldest_daughter()) {
00249 
00250     vector r_com  = b->get_pos() - dc_pos;
00251 
00252     real Uc, Bc, Vc, Rc, Ic;
00253     Uc=Bc=Vc=Rc=Ic=VERY_LARGE_NUMBER;
00254     real Us, Bs, Vs, Rs, Is;
00255     real U, B, V, R, I;
00256     U=B=V=R=I=VERY_LARGE_NUMBER;
00257     for_all_leaves(dyn, b, bb) {
00258         r = bb->get_pos()-dc_pos;
00259         if(project<=2)
00260           rcom = r[project];
00261         else
00262           rcom = abs(r);
00263       if (rcom>r_min && rcom<r_max) {
00264 
00265        m_tot += bb->get_mass();
00266        Us=Bs=Vs=Rs=Is=VERY_LARGE_NUMBER;
00267        get_UBVRI_of_star(bb, dc_pos, Us, Bs, Vs, Rs, Is);
00268        if (Vs<v_min && Vs>v_max) {
00269        combine_ubvri(Uc, Bc, Vc, Rc, Ic,
00270                      Us, Bs, Vs, Rs, Is,
00271                      U, B, V, R, I);
00272        Uc=U;
00273        Bc=B;
00274        Vc=V;
00275        Rc=R;
00276        Ic=I;
00277 }
00278 }
00279     }
00280 
00281     real L_tot=0, T_eff=0;
00282     cout << " Cluster: m= " << m_tot << " R= " << r_max << " L= " << L_tot
00283          << " T_eff= " << T_eff 
00284          << " r= "  << dc_pos[0] << " " << dc_pos[1] << " " << dc_pos[2] 
00285          << " ubvri= "<<  Uc << " " << Bc << " " << Vc << " " << Rc << " "
00286          << Ic  << " :: ";
00287         
00288   }
00289 }
00290 
00291 main(int argc, char ** argv)
00292 {
00293     check_help();
00294 
00295     bool C_flag = false;
00296 
00297     real r_max = VERY_LARGE_NUMBER;
00298     real r_min = 0;
00299     real v_max = -VERY_LARGE_NUMBER;
00300     real v_min = VERY_LARGE_NUMBER;
00301 
00302     int project = 0; // project along x, y, z or no projection [3]
00303 
00304     extern char *poptarg;
00305     int c;
00306     char* param_string = "CV:v:R:r:p:";
00307 
00308     while ((c = pgetopt(argc, argv, param_string)) != -1)
00309         switch(c)
00310             {
00311             case 'C': C_flag = true;
00312                       break;
00313             case 'V': v_max = atof(poptarg);
00314                       break;
00315             case 'v': v_min = atof(poptarg);
00316                       break;
00317             case 'R': r_max = atof(poptarg);
00318                       break;
00319             case 'r': r_min = atof(poptarg);
00320                       break;
00321             case 'p': project = atoi(poptarg);
00322                       break;
00323             case '?': params_to_usage(cerr, argv[0], param_string);
00324                       get_help();
00325                       exit(1);
00326             }
00327 
00328     dyn *b = NULL;
00329     int count = 0;
00330 
00331     bool cod, try_com = false;
00332     vector dc_pos = 0;
00333     vector dc_vel = 0;
00334     while (b = get_dyn(cin)) {
00335 
00336         real rd_min = b->get_starbase()->conv_r_star_to_dyn(r_min);
00337         real rd_max = b->get_starbase()->conv_r_star_to_dyn(r_max);
00338 
00339         cout << "\n\nTime = "
00340              << b->get_starbase()->conv_t_dyn_to_star(b->get_system_time())
00341              << endl;
00342 
00343         if (find_qmatch(b->get_dyn_story(), "density_center_pos")) {
00344       
00345            if (getrq(b->get_dyn_story(), "density_center_time")
00346                != b->get_system_time()) {
00347                warning("mkpovfile: neglecting out-of-date density center");
00348                try_com = true;
00349             } else
00350             cod = true;
00351       
00352             dc_pos = getvq(b->get_dyn_story(), "density_center_pos");
00353             dc_vel = getvq(b->get_dyn_story(), "density_center_vel");
00354          }
00355 
00356          if (try_com && find_qmatch(b->get_dyn_story(), "com_pos")) {
00357 
00358             if (getrq(b->get_dyn_story(), "com_time")
00359                != b->get_system_time()) {
00360 
00361                warning("lagrad: neglecting out-of-date center of mass");
00362             } else
00363                 dc_pos = getvq(b->get_dyn_story(), "com_pos");
00364                 dc_vel = getvq(b->get_dyn_story(), "com_vel");
00365          }
00366 
00367         real U,B,V,R,I;
00368         U=B=V=R=I=VERY_LARGE_NUMBER;
00369 
00370         if (!C_flag)
00371            print_binary_recursive(b, dc_pos, dc_vel, U, B, V, R, I);
00372         else
00373            print_integrated_cluster(b, dc_pos, rd_min, rd_max, v_min,
00374 v_max,
00375                                     U, B,
00376                                     V, R, I, project);
00377 
00378        rmtree(b);
00379     }
00380 }
00381 
00382 #endif