nstarev.C Source File

00001 
00002 
00003 
00004 
00005 
00006 
00007 
00008 
00009 
00010 
00011 
00012 
00013 
00014 
00015 
00016 //++ Notes:
00017 //++  As in starev no kicks are applied in the supernova event.
00018 //++
00019 //++ Example of usage:      
00020 //++  mknode -n 10 | mkmass -u 100 -l 10 | addstar -T 1 | nstarev
00021 //++
00022 //++ See also: addstar
00023 //++           mkmass
00024 //++           mknode
00025 //++           starev
00026 
00027 
00029 
00030 //   version 1:  Februari 1993   Simon F. Portegies Zwart
00031 //                               spz@grape.c.u-tokyo.ac.jp
00032 
00033 //++ Due to overhull currently not operational.
00034 
00035 //#include "node.h"
00036 #include "dyn.h"
00037 #include "single_star.h"
00038 #include "main_sequence.h"
00039 //#include "sstar_to_dyn.h"
00040 
00041 #define EPSILON 1.e-10
00042 
00043 #ifdef TOOLBOX
00044 
00045 #define  SEED_STRING_LENGTH  60
00046 
00047 
00048 
00049 local void alter_star_from_random_time(dyn* b) {
00050 
00051     star* ss = ((star*)b->get_starbase());
00052     b->set_mass(b->get_starbase()
00053                  ->conv_m_star_to_dyn(b->get_starbase()
00054                                         ->get_total_mass()));
00055 
00056     ss->set_relative_age(b->get_starbase()
00057                            ->get_current_time());
00058     ss->set_current_time(b->get_starbase()
00059                            ->conv_t_dyn_to_star(b->get_system_time()));
00060     ss->set_current_time(0);
00061 }
00062 
00063 local void evolve_star_until_next_time(dyn* bi, const real out_time) {
00064 
00065           real current_time = ((star*)bi->get_starbase())->get_current_time();
00066           real time_step    =  bi->get_starbase()->get_evolve_timestep();
00067 
00068           while (out_time>current_time+time_step) {
00069              bi->get_starbase()->evolve_element(current_time+time_step);
00070              bi->get_starbase()->evolve_element(
00071                  min(current_time+time_step+EPSILON, out_time));
00072              current_time = ((star*)bi->get_starbase())->get_current_time();
00073              time_step    =  bi->get_starbase()->get_evolve_timestep();
00074 
00075              star_state ss(dynamic_cast(star*, bi->get_starbase()));
00076              put_state(ss, cerr);
00077 
00078 //             print_star(bi->get_starbase(), cerr);
00079              int p = cerr.precision(HIGH_PRECISION);
00080              bi->get_starbase()->dump(cerr, false);
00081              cerr.precision(p);
00082           }
00083          bi->get_starbase()->evolve_element(out_time);
00084 
00085          print_star(bi->get_starbase(), cerr);
00086 
00087 /*
00088 cerr<< "ov: " << bi->get_starbase()->get_element_type() 
00089         << " " <<bi->get_starbase()->get_total_mass() << " t= "
00090         << ((star*)bi->get_starbase())->get_current_time() << " "
00091         << ((star*)bi->get_starbase())->get_current_time()
00092            +bi->get_starbase()->get_evolve_timestep() << " -> "
00093         << out_time << endl;
00094 */
00095 }
00096 
00097 local void  evolve_the_stellar_system(dyn* b, real time) {
00098 
00099       dyn * bi;
00100 
00101       if (b->get_oldest_daughter() != NULL)
00102          for (bi=b->get_oldest_daughter(); bi != NULL;
00103                                            bi=bi->get_younger_sister())
00104          evolve_the_stellar_system(bi, time);
00105       else {
00106         evolve_star_until_next_time(bi, time);
00107 //         b->get_starbase()->evolve_element(time);
00108       }
00109    }
00110 
00111 void main(int argc, char ** argv) {
00112 
00113     bool s_flag = false;
00114     bool c_flag = false;
00115     bool m_flag = false;
00116     int  n = 1;
00117     int n_steps = 1;
00118 
00119     int  input_seed, actual_seed;
00120     int id=1;
00121     real m_tot=1,r_hm=1, t_hc=1;
00122     real m_rel=5, m_core=0.01;
00123     stellar_type type = Main_Sequence;
00124     real t_start = 0;
00125     real t_end   = 100;
00126     char  *comment;
00127     char  seedlog[SEED_STRING_LENGTH];
00128     extern char *poptarg;
00129     int c;
00130     char* param_string = "M:N:n:r:s:t:T:";
00131 
00132     check_help();
00133 
00134     while ((c = pgetopt(argc, argv, param_string)) != -1)
00135         switch(c) {
00136 
00137             case 'r': r_hm = atof(poptarg);
00138                       break;
00139             case 'M': m_flag = true;
00140                       m_tot = atof(poptarg);
00141                       break;
00142             case 'N': n = atoi(poptarg);
00143                       break;
00144             case 'n': n_steps = atoi(poptarg);
00145                       break;
00146             case 'T': t_hc = atof(poptarg);
00147                       break;
00148             case 't': t_end = atof(poptarg);
00149                       break;
00150             case 'S': type = (stellar_type)atoi(poptarg);
00151                       break;
00152             case 's': s_flag = TRUE;
00153                       input_seed = atoi(poptarg);
00154                       break;
00155             case 'c': c_flag = TRUE;
00156                       comment = poptarg;
00157                       break;
00158             case '?': params_to_usage(cerr, argv[0], param_string);
00159                       get_help();
00160                       exit(1);
00161         }
00162 
00163     cerr.precision(HIGH_PRECISION);
00164 
00165     if(!s_flag) input_seed = 0;
00166     actual_seed = srandinter(input_seed);
00167 
00168     // make flat tree 
00169     dyn* root = get_dyn(cin);
00170     //    node *root  = mknode(1);
00171     root->log_history(argc, argv);
00172 //    root->get_starbase()->set_stellar_evolution_scaling(m_tot, r_hm, t_hc);
00173 
00174 //    node *the_star = root->get_oldest_daughter();
00175 
00176     addstar(root, t_start, type);
00177 
00178     cerr.precision(STD_PRECISION);
00179 
00180     //    put_node(cout, *root);
00181 
00182     real dt, time = 0;
00183     real delta_t = t_end/((real)n_steps);
00184     real out_time, current_time; 
00185     real time1, time2;
00186     real previous_time;
00187     int  nstps=0;
00188     if(t_end>0)
00189     for_all_daughters(dyn, root, bi) {
00190        out_time = 0;
00191        do {
00192           out_time = min(out_time+delta_t, t_end);
00193           evolve_star_until_next_time(bi, out_time);
00194 //        bi->get_starbase()->dump(cerr, false);
00195        }
00196        while(out_time < t_end);
00197     }
00198     else if(t_end == 0) {
00199       real m, t_ZAMS;
00200       for_all_daughters(dyn, root, bi) {
00201         m = bi->get_starbase()->get_total_mass();
00202         t_ZAMS =  main_sequence_time(m);
00203         evolve_star_until_next_time(bi, randinter(0, t_ZAMS));
00204         bi->set_mass(bi->get_starbase()
00205                        ->conv_m_star_to_dyn(bi->get_starbase()
00206                                               ->get_total_mass()));
00207       }
00208     }
00209     else
00210       for_all_daughters(dyn, root, bi) {
00211         evolve_star_until_next_time(bi, randinter(0, -t_end));
00212 
00213         alter_star_from_random_time(bi);
00214 
00215       }
00216 #if 0
00217     for_all_daughters(dyn, root, bi) {
00218        cerr << "Time = " << bi->get_starbase()->get_current_time()
00219             << " [Myear],  mass = " << bi->get_starbase()->get_total_mass() 
00220             << " [Msun],  radius = " << bi->get_starbase()
00221                                           ->get_effective_radius() 
00222             << "   " << type_string(bi->get_starbase()->get_element_type())
00223             << endl;
00224     }
00225 #endif
00226 
00227     put_node(cout, *root);
00228     delete root;
00229 
00230   }
00231 
00232 #if 0
00233 
00235     if (s_flag == FALSE)
00236         input_seed = 0;                         /* default */
00237     actual_seed = srandinter(input_seed);
00238 
00239     node *b;
00240 
00241     b = get_node(cin);
00242     if (c_flag == TRUE)
00243         b->log_comment(comment);
00244     b->log_history(argc, argv);
00245     b->get_starbase()->set_stellar_evolution_scaling(m_tot, r_hm, t_hc);
00246 
00247     addstar(b, t_start, type);
00248 
00249     //put_node(cout, *b);
00250 
00251     real t_out, dt = t_end/n;
00252     for (int i=0; i<n; i++) {
00253       t_out = dt * (i+1);
00254 
00255       evolve_the_stellar_system(b, t_out);
00256       //      for_all_daughters(node, b, bi) 
00257       //        bi->get_starbase()->evolve_element(t_out);
00258 
00259       put_node(cout, *b);
00260     }
00261 }
00262 
00263 #endif
00264 #endif