00001
00006
00007
00008
00009
00010
00011
00012
00013 #include "single_star.h"
00014 #include "sstar_to_dyn.h"
00015 #include "dyn.h"
00016
00017 #ifndef TOOLBOX
00018
00019 #define NNEBULAE 10
00020
00021 typedef struct nebulae {
00022 int id;
00023 real time;
00024 vector pos;
00025 };
00026 nebulae nebula[NNEBULAE];
00027
00028 #define NSN_REMNANT 4
00029 typedef struct sn_remnants {
00030 int id;
00031 bool bh;
00032 real time;
00033 vector pos;
00034 };
00035 sn_remnants sn_remnant[NSN_REMNANT];
00036
00037 #define NCOLLISIONS 10
00038 typedef struct collisions {
00039 int starid;
00040 char name[25];
00041 real time;
00042 };
00043 collisions collision[NCOLLISIONS];
00044
00045 enum filter_type {Visual=0, Radio, X_rays};
00046
00047 local filter_type get_filter_type(char fltr) {
00048
00049 switch(fltr) {
00050 case 'V': return Visual;
00051 break;
00052 case 'R': return Radio;
00053 break;
00054 case 'X': return X_rays;
00055 break;
00056 default: cerr << "No such filter"<<endl;
00057 exit(-1);
00058 };
00059
00060 }
00061
00062 void new_camera_position(vector &v_new,
00063 vector first_pos,
00064 vector last_pos,
00065 int nsteps,
00066 int nsnap) {
00067
00068 real fstep = (1.0*nsnap)/(1.0*nsteps);
00069 v_new[0] = first_pos[0] + fstep * (last_pos[0] - first_pos[0]);
00070 v_new[1] = first_pos[1] + fstep * (last_pos[1] - first_pos[1]);
00071 v_new[2] = first_pos[2] + fstep * (last_pos[2] - first_pos[2]);
00072
00073 }
00074
00075 #if 0
00076 void new_camera_position(vector &v_new,
00077 vector v_old,
00078 real theta_rotation,
00079 real phi_rotation,
00080 int nsteps,
00081 int nsnap) {
00082
00083
00084 real theta = nsnap*theta_rotation/nsteps;
00085 real phi = nsnap*phi_rotation/nsteps;
00086
00087 PRC(nsnap);PRC(theta);PRL(phi);
00088 real r = abs(v_old);
00089 v_new[0] = r * sin(theta);
00090 v_new[1] = r * sin(theta) * cos(phi);
00091 v_new[2] = r * cos(theta);
00092
00093 cerr << "vector = " << v_new << endl;
00094
00095 }
00096 #endif
00097
00098 local void print_camera_on_star(dyn* b, vector com_pos,
00099 int camera_on_star_id,
00100 real r_star,
00101 real aperture, int blur_samples) {
00102
00103 vector cam_view = b->get_pos() + b->get_vel();
00104 vector cam_pos = b->get_pos() + r_star*b->get_vel();
00105
00106
00107
00108 cam_view[0] = 0;
00109 cam_view[1] = 0;
00110 cam_view[2] = 0;
00111
00112 vector normal;
00113 normal[0] = sqrt(pow(abs(cam_pos), 2)
00114 - pow(cam_pos[1], 2)
00115 - pow(cam_pos[2], 2));
00116 normal[1] = 1;
00117 normal[2] = (-cam_pos[0] * normal[0] - cam_pos[1])/cam_pos[2];
00118
00119 cout << "// Normal to camera " << endl;
00120 cout << " #declare normal_to_camera = < " << normal[0] << ", "
00121 << normal[1] << ", "
00122 << normal[2] << " >" << endl;
00123
00124 cout << "// camera located on star #" << camera_on_star_id << endl;
00125 cout << "camera {" << endl
00126 << " location < " << cam_pos[0] << ", "
00127 << cam_pos[1] << ", "
00128 << cam_pos[2] << " >" << endl
00129 << " look_at < " << cam_view[0] << ", "
00130 << cam_view[1] << ", "
00131 << cam_view[2] << " >" << endl
00132 << " blur_samples " << blur_samples << endl;
00133 if (aperture>0)
00134 cout << " focal_point < " << cam_view[0] << ", "
00135 << cam_view[1] << ", "
00136 << cam_view[2] << " >" << endl
00137 << " aperture " << aperture << endl;
00138 cout << "}" << endl << endl;
00139
00140 }
00141
00142
00143 local bool print_camera_position_recursive(dyn* b, vector cam_pos,
00144 int camera_on_star_id,
00145 real r_star,
00146 real aperture, int blur_samples) {
00147
00148 if (b->get_oldest_daughter()) {
00149
00150 vector com_pos = b->get_pos() - cam_pos;
00151
00152
00153 for_all_daughters(dyn, b, bb)
00154 if (bb->n_leaves() >= 2)
00155 return print_camera_position_recursive(bb, com_pos,
00156 camera_on_star_id, r_star,
00157 aperture, blur_samples);
00158
00159 else if (bb->get_index()==camera_on_star_id) {
00160 print_camera_on_star(bb, com_pos,
00161 camera_on_star_id, r_star,
00162 aperture, blur_samples);
00163
00164 return true;
00165 }
00166
00167 }
00168 return false;
00169 }
00170
00171 local void print_filename_counter(int counter, ostream& s) {
00172
00173 if (counter>=10000) {
00174 cout << "\nToo namy filenames in print_povray_header()" << endl;
00175 exit(1);
00176 }
00177 else if (counter<10)
00178 s << "000" << counter;
00179 else if (counter<100)
00180 s << "00" << counter;
00181 else if (counter<1000)
00182 s << "0" << counter;
00183 else
00184 s << counter;
00185 }
00186
00187 local void print_pl_nebula(vector pos, real scale) {
00188
00189 cout << "object { Pl_Nebula scale " << scale
00190 << " translate < " << pos[0] << ", "
00191 << pos[1] << ", "
00192 << pos[2] << " > }"
00193 << endl;
00194 }
00195
00196 local void print_sn_nebula(vector pos, real scale) {
00197
00198 cout << "object { SN_Remnant scale " << scale
00199 << " translate < " << pos[0] << ", "
00200 << pos[1] << ", "
00201 << pos[2] << " > }"
00202 << endl;
00203 }
00204
00205 local void print_some_data(dyn *b) {
00206
00207 if (b->get_oldest_daughter()) {
00208
00209 real time = b->get_starbase()->conv_t_dyn_to_star(b->get_real_system_time());
00210
00211 int nts=0, nbs=0, nss=0;
00212 for_all_daughters(dyn, b, bb)
00213 if (bb->n_leaves() > 2)
00214 nts++;
00215 else if (bb->n_leaves() >= 2)
00216 nbs++;
00217 else
00218 nss++;
00219 }
00220 }
00221
00222 local void print_hertzsprung_Russell_diagram(dyn* b, vector cam_pos) {
00223
00224 if (b->get_oldest_daughter()) {
00225
00226 print_some_data(b);
00227
00228 int stype_s[no_of_star_type_summ];
00229 for (int i=0; i<no_of_star_type_summ; i++)
00230 stype_s[i]=0;
00231
00232 cout << "#declare Stellar_HRD = union {" << endl;
00233 cout << " object { Y_Axis } " << endl;
00234 cout << " object { X_Axis }\n " << endl;
00235
00236 for_all_leaves(dyn, b, bi) {
00237
00238 star_type_spec tpe_class = NAC;
00239 spectral_class star_class;
00240 stellar_type stype = NAS;
00241 stellar_type_summary sstype = ZAMS;
00242 real t_cur, t_rel, m_rel, m_env, m_core, co_core;
00243 real T_eff, L_eff, p_rot, b_fld;
00244
00245 if (bi->get_use_sstar()) {
00246 stype = bi->get_starbase()->get_element_type();
00247 sstype = summarize_stellar_type(stype);
00248 stype_s[dynamic_cast(int, sstype)]++;
00249
00250 T_eff = bi->get_starbase()->temperature();
00251 L_eff = bi->get_starbase()->get_luminosity();
00252 }
00253 else if (bi->get_star_story()) {
00254 extract_story_chapter(stype, t_cur, t_rel, m_rel, m_env,
00255 m_core, co_core,
00256 T_eff, L_eff, p_rot, b_fld,
00257 *bi->get_star_story());
00258 sstype = summarize_stellar_type(stype);
00259 stype_s[dynamic_cast(int, sstype)]++;
00260 star_class = get_spectral_class(T_eff);
00261 #if 0
00262 if (find_qmatch(bi->get_star_story(), "Type")) {
00263 cerr <<"Reading"<<endl;
00264 stype = extract_stellar_type_string(
00265 getsq(bi->get_star_story(), "Type"));
00266 sstype = summarize_stellar_type(stype);
00267 stype_s[dynamic_cast(int, sstype)]++;
00268 T_eff = getrq(bi->get_star_story(), "T_eff");
00269 star_class = get_spectral_class(T_eff);
00270 L_eff = getrq(bi->get_star_story(), "L_eff");
00271 #endif
00272 }
00273 else {
00274 cout << " No stellar information found for: ";
00275 bi->pretty_print_node(cout);
00276 return;
00277 }
00278
00279 real xt_pos = (4.78 - log10(T_eff))/4.78;
00280 real yl_pos = (log10(L_eff) + 1)/7.;
00281
00282 if (xt_pos>0&&xt_pos<1 && yl_pos>0&&yl_pos<1)
00283 cout << " object { Red_Sphere translate < "
00284 << xt_pos << ", "
00285 << yl_pos << ", "
00286 << "0 > }" << endl;
00287 }
00288
00289 cout << "\n} // End Stellar_HRD" << endl;
00290 cout << "object { Stellar_HRD translate < "
00291 << cam_pos[0]+0.5 << ", "
00292 << cam_pos[1]-1 << ", "
00293 << cam_pos[2] + 2
00294 << " > }\n" << endl;
00295
00296 }
00297 }
00298
00299 local void add_collision_effect(dyn *bi, vector pos, real time,
00300 real scale) {
00301
00302 for (int i=0; i<NCOLLISIONS; i++) {
00303 if(bi->get_name()) {
00304 if(!strcmp(bi->get_name(), collision[i].name) &&
00305 time >= collision[i].time && time < collision[i].time+0.125) {
00306 cerr << "Adding collision at time = " << time << endl;
00307 cerr << "object { OStar scale "
00308 << 4 * scale * pow(5 * (0.125 - (time-collision[i].time)), 3)
00309 << " translate < "
00310 << pos[0] << ", " << pos[1] << ", " << pos[2] << " > }"
00311 << endl;
00312 if(collision[i].starid == 0) {
00313 cout << "object { KStar scale ";
00314 }
00315 else {
00316 cout << "object { OStar scale ";
00317 }
00318 cout << 4 * scale * pow(5 * (0.125 - (time-collision[i].time)), 3)
00319 << " translate < "
00320 << pos[0] << ", " << pos[1] << ", " << pos[2] << " > }"
00321 << endl;
00322 }
00323 }
00324 }
00325 }
00326
00327 local void print_star(dyn *bi, vector pos,
00328 real scale_L, filter_type filter) {
00329
00330
00331
00332
00333
00334
00335
00336 real time = bi->get_starbase()->conv_t_dyn_to_star(bi->get_real_system_time());
00337
00338
00339 real Rsun_per_parsec = cnsts.parameters(solar_radius)
00340 / cnsts.parameters(parsec);
00341
00342
00343
00344
00345 star_type_spec tpe_class = NAC;
00346 spectral_class star_class;
00347 stellar_type stype = NAS;
00348 stellar_type_summary sstype = ZAMS;
00349 real t_cur, m_rel, m_env, m_core, co_core, T_eff, L_eff, p_rot, b_fld;
00350 real t_rel=0, R_eff=0;
00351 real M_tot, U, B, V, R, I;
00352
00353 if (bi->get_star_story()) {
00354
00355
00356
00357 stype = extract_stellar_type_string(
00358 getsq(bi->get_star_story(), "Type"));
00359 M_tot = getrq(bi->get_star_story(), "M_rel");
00360 T_eff = getrq(bi->get_star_story(), "T_eff");
00361 star_class = get_spectral_class(T_eff);
00362 L_eff = getrq(bi->get_star_story(), "L_eff");
00363
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373 sstype = summarize_stellar_type(stype);
00374 star_class = get_spectral_class(T_eff);
00375
00376
00377
00378 ltm_to_ubvri(log10(L_eff), log10(T_eff), M_tot,
00379 U, B, V, R, I);
00380
00381
00382 if (find_qmatch(bi->get_star_story(), "Class"))
00383 tpe_class = extract_stellar_spec_summary_string(
00384 getsq(bi->get_star_story(), "Class"));
00385 if (L_eff>0)
00386 R_eff = pow(T_eff/cnsts.parameters(solar_temperature), 2)
00387 / sqrt(L_eff);
00388 }
00389 else if (bi->get_use_sstar()) {
00390
00391 put_dyn(cerr, *bi);
00392 stype = bi->get_starbase()->get_element_type();
00393 M_tot = bi->get_starbase()->conv_m_dyn_to_star(bi->get_mass());
00394 t_cur = bi->get_starbase()->get_current_time();
00395 t_rel = bi->get_starbase()->get_relative_age();
00396 T_eff = bi->get_starbase()->temperature();
00397 L_eff = bi->get_starbase()->get_luminosity();
00398 star_class = get_spectral_class(T_eff);
00399 R_eff = bi->get_starbase()->get_effective_radius();
00400 ltm_to_ubvri(log10(L_eff), log10(T_eff), M_tot,
00401 U, B, V, R, I);
00402
00403
00404
00405
00406 }
00407 else {
00408 cout << " No stellar information found for: ";
00409 bi->pretty_print_node(cout);
00410 return;
00411 }
00412
00413 int known_at;
00414 bool is_known = false;
00415 bool should_be_known = false;
00416
00417 if (remnant(stype) && t_rel <= 1)
00418 should_be_known = true;
00419
00420 if (remnant(stype)) {
00421 for (int i=0; i<NNEBULAE; i++)
00422 if (bi->get_index() == nebula[i].id) {
00423 is_known = true;
00424 known_at = i;
00425 }
00426
00427 }
00428
00429 if (should_be_known) {
00430 if (!is_known)
00431 for (int i=0; i<NNEBULAE; i++)
00432 if (nebula[i].id < 0) {
00433 nebula[i].id = bi->get_index();
00434 nebula[i].pos = pos;
00435 nebula[i].time = t_rel;
00436 known_at = i;
00437 }
00438
00439 if (stype== Helium_Dwarf || Carbon_Dwarf || Oxygen_Dwarf)
00440 print_pl_nebula(nebula[known_at].pos, scale_L * (0.5 + t_rel));
00441 else {
00442 #if 0
00443 print_sn_nebula(nebula[known_at].pos, scale_L * (0.5 + t_rel));
00444 if (t_rel<0.01)
00445 cout << "object { single_star \n"
00446 << " finish { ambient <"
00447 << 0.4 << ", " << 0.6 << ", " << 0.7 << ">}\n"
00448 << " scale " << scale_L * pow(100 * (0.01-t_rel), 3)
00449 << " translate < "
00450 << pos[0] << ", " << pos[1] << ", " << pos[2]-0.5 << " > }"
00451 << endl;
00452 #endif
00453
00454 }
00455 }
00456 else if(is_known)
00457 nebula[known_at].id = -1;
00458
00459
00460 for (int i=0; i<NSN_REMNANT; i++)
00461 if(bi->get_index() == sn_remnant[i].id &&
00462 t_cur >= sn_remnant[i].time && t_cur < sn_remnant[i].time+0.025) {
00463 if (sn_remnant[i].bh)
00464 print_sn_nebula(pos, scale_L * (0.5 + t_rel));
00465 else
00466 print_pl_nebula(pos, scale_L * (0.5 + t_rel));
00467 if(t_cur >= sn_remnant[i].time && t_cur < sn_remnant[i].time+0.01) {
00468 cout << "object { single_star \n"
00469 << " finish { ambient <";
00470 if (sn_remnant[i].bh)
00471 cout << 0.4 << ", " << 0.6 << ", " << 0.7 << ">}\n";
00472 else
00473 cout << 0.56 << ", " << 0.14 << ", " << 0.14 << ">}\n";
00474 cout << " scale " << scale_L * pow(100 * (0.01-t_rel), 3)
00475 << " translate < "
00476 << pos[0] << ", " << pos[1] << ", " << pos[2]-0.5 << " > }"
00477 << endl;
00478 }
00479 }
00480
00481 add_collision_effect(bi, pos, time, scale_L);
00482 #if 0
00483
00484 cerr << "Add collisions"<<endl;
00485 for (int i=0; i<NCOLLISIONS; i++) {
00486 PRC(i);PRC(bi->get_index());PRL(collision[i].starid);
00487 if(bi->get_index() == collision[i].starid &&
00488 t_cur >= collision[i].time && t_cur < collision[i].time+0.125) {
00489 cout << "object { OStar scale "
00490 << scale_L * pow(5 * (0.125 - (t_cur-collision[i].time)), 3)
00491 << " translate < "
00492 << pos[0] << ", " << pos[1] << ", " << pos[2]-0.5 << " > }"
00493 << endl;
00494 }
00495 }
00496 #endif
00497
00498
00499 real logL_eff = scale_L * sqrt(log10(1 + L_eff));
00500 #if 0
00501 cout << "object { "
00502 << type_short_string(star_class) << "Star scale "
00503 << logL_eff << " translate < "
00504 << pos[0] << ", " << pos[1] << ", " << pos[2] << " > }"
00505 << endl;
00506
00507 if (tpe_class == Accreting)
00508 cout << "object { Accretion_Disc scale "
00509 << logL_eff << " translate < "
00510 << pos[0] << ", " << pos[1] << ", " << pos[2] << " > }"
00511 << endl;
00512 #endif
00513
00514 real G;
00515 real apparent_size = -1;
00516 if (filter == Visual) {
00517
00518
00519 real UM100 = -7.38974, UM1 = 5.17096, UM01 = 17.9608;
00520 real BM100 = -6.23550, BM1 = 5.06279, BM01 = 16.7472;
00521 real VM100 = -5.91088, VM1 = 4.46967, VM01 = 15.0682;
00522 real RM100 = -15.9099, RM1 = 4.13746, RM01 = 13.7278;
00523 real IM100 = -25.9089, IM1 = 3.81839, IM01 = 12.0133;
00524
00525
00526 RM100 = -6; RM01 = 11.4;
00527 R = max(0., min(1., (R-RM100)/(RM01-RM100)));
00528
00529
00530
00531 VM100 = -5.8; VM01 = 13;
00532 G = 1 - min(1., max(0., sqrt(abs((V-VM01)/(VM01-VM100)))));
00533
00534
00535
00536
00537
00538 BM100 = -7; BM01 = 15;
00539 B = pow(max(0., min(1., (BM01-B)/(BM01-BM100))), 2);
00540
00541
00542
00543
00544 apparent_size = 0.1 * scale_L * logL_eff;
00545 }
00546 else {
00547
00548 real VM100 = -5.91088, VM1 = 4.46967, VM01 = 15.0682;
00549 apparent_size = 0.1 *scale_L * M_tot;
00550
00551 switch(stype) {
00552 case Carbon_Star:
00553 case Helium_Star:
00554 case Helium_Giant: R = 0; B=0;
00555 G = 1 - min(1., max(0., sqrt(abs((V-VM1)/(VM01-VM100)))));
00556 apparent_size = 0.1 * scale_L * (VM01-V)/(VM01-VM100);
00557 break;
00558 case Carbon_Dwarf:
00559 case Helium_Dwarf:
00560 case Oxygen_Dwarf: G=0; B=0;
00561 R = 1 - min(1., max(0., sqrt(abs((V-VM1)/(VM01-VM100)))));
00562 apparent_size = 0.1 * scale_L * (VM01-V)/(VM01-VM100);
00563 break;
00564 case Thorn_Zytkow: R = 1; G=0; B=0;
00565 apparent_size = 0.1 * scale_L * (VM01-V)/(VM01-VM100);
00566 break;
00567 case Xray_Pulsar:
00568 case Radio_Pulsar:
00569 case Neutron_Star: R = 1; G=0; B=0;
00570 B = -0.3 * log10(min(1., p_rot));
00571 apparent_size = 0.1 * scale_L * B;
00572 break;
00573 case Black_Hole: R = 1; G=1; B=1;
00574 break;
00575 default:
00576 R = 0; G=0; B=0;
00577 apparent_size = -1;
00578 break;
00579
00580 };
00581 }
00582
00583 if(apparent_size>0) {
00584 cout << "object { single_star \n"
00585 << " finish { ambient <"
00586 << R << ", " << G << ", " << B << ">}\n"
00587 << " scale " << apparent_size << " translate < "
00588 << pos[0] << ", " << pos[1] << ", " << pos[2] << " > }"
00589 << endl;
00590 }
00591
00592 if (tpe_class == Accreting)
00593 cout << "object { Accretion_Disc scale "
00594 << 0.05*logL_eff << " translate < "
00595 << pos[0] << ", " << pos[1] << ", " << pos[2] << " > }"
00596 << endl;
00597
00598 }
00599
00600 local void print_node(dyn *bi, vector pos, real mass_scale,
00601 real mmax) {
00602
00603 real time = getrq(bi->get_root()->get_dyn_story(), "real_system_time");
00604
00605 real apparent_size = mass_scale * sqrt(bi->get_mass());
00606 real mmin = 0;
00607 real m = bi->get_mass();
00608 real R = sqrt((mmax-m)/(mmax-mmin));
00609 real G = 0.5;
00610 real B = sqrt(m/(mmax-mmin));
00611 if(m>0.0125) {
00612 R = 1;
00613 G = 0;
00614 B = 0;
00615 }
00616 else if(m<0.003125) {
00617 R = 0.8;
00618 G = 0.8;
00619 B = 0.8;
00620 }
00621 else {
00622 R = 0.6;
00623 G = 0.8;
00624 B = 0.196078;
00625
00626
00627
00628 }
00629 #if 0
00630 vector cam_pos;
00631 cam_pos[0] = 1;
00632 cam_pos[1] = 0;
00633 cam_pos[2] = 0;
00634 real V=100;
00635 real v = (V-pos[3])*tan(abs(cam_pos-pos)/pos[3]);
00636 R = 1;G=0;B=0;
00637 cout << "object { simple_star \n"
00638 << " finish { ambient <"
00639 << R << ", " << G << ", " << B << ">}\n";
00640 cout << " scale " << apparent_size;
00641 cout << " translate < "
00642 << pos[0]-v << ", " << pos[1] << ", " << pos[2] << " > }"
00643 << endl;
00644 R=0;G=0;B=1;
00645 cout << "object { simple_star \n"
00646 << " finish { ambient <"
00647 << R << ", " << G << ", " << B << ">}\n";
00648 cout << " scale " << apparent_size;
00649 cout << " translate < "
00650 << pos[0]+v << ", " << pos[1] << ", " << pos[2] << " > }"
00651 << endl;
00652 #endif
00653
00654
00655 R=G=B=1;
00656 cout << "object { simple_star \n"
00657 << " finish { ambient <"
00658 << R << ", " << G << ", " << B << ">}\n";
00659 if(time>=25)
00660 cout << " scale " << 1./(1+0.2*(25-20)) *apparent_size;
00661 else if(time>=20)
00662 cout << " scale " << 1./(1+0.2*(time-20)) *apparent_size;
00663 else
00664 cout << " scale " << apparent_size;
00665 cout << " translate < "
00666 << pos[0] << ", " << pos[1] << ", " << pos[2] << " > }"
00667 << endl;
00668
00669
00670 add_collision_effect(bi, pos, time, mass_scale);
00671
00672
00673 #if 0
00674
00675 for (int i=0; i<NCOLLISIONS; i++) {
00676
00677 if(bi->get_name()) {
00678
00679 if(!strcmp(bi->get_name(), collision[i].name) &&
00680 time >= collision[i].time && time < collision[i].time+0.125) {
00681 cerr << "Adding collision at time = " << time << endl;
00682 cerr << "object { OStar scale "
00683 << 2 *mass_scale * pow(5 * (0.125 - (time-collision[i].time)), 3)
00684 << " translate < "
00685 << pos[0] << ", " << pos[1] << ", " << pos[2]-apparent_size << " > }"
00686 << endl;
00687 cout << "object { OStar scale "
00688 << 2 *mass_scale * pow(5 * (0.125 - (time-collision[i].time)), 3)
00689 << " translate < "
00690 << pos[0] << ", " << pos[1] << ", " << pos[2]-apparent_size << " > }"
00691 << endl;
00692 }
00693 }
00694 }
00695 #endif
00696
00697
00698
00699
00700
00701
00702
00703 }
00704
00705 local int print_povray_binary_recursive(dyn *b,
00706 vector dc_pos,
00707 real mass_limit, real number_limit,
00708 bool povray, real scale_L,
00709 filter_type filter) {
00710
00711 int nb = 0;
00712 if (b->get_oldest_daughter()) {
00713
00714 vector r_com = b->get_pos() - dc_pos;
00715
00716 real m_tot = b->get_starbase()->conv_m_dyn_to_star(b->get_mass());
00717
00718 for_all_daughters(dyn, b, bb)
00719 if (bb->n_leaves() >= 2)
00720 nb += print_povray_binary_recursive(bb, dc_pos,
00721 mass_limit, number_limit,
00722 povray, scale_L, filter);
00723 else {
00724 print_star(bb, bb->get_pos()-r_com, scale_L, filter);
00725 }
00726
00727 }
00728 return nb;
00729 }
00730
00731 local void print_povtime(real time,
00732 vector pos,
00733 real scale=1,
00734 real depth=0.25) {
00735
00736 int p = cout.precision(LOW_PRECISION);
00737 cout << "text { ttf \"timrom.ttf\" ";
00738 if(time==0)
00739 cout << "\"Time = " << time << " \" ";
00740 else if(time>0)
00741 cout << "\"Time = " << time/73.387 << " Myr \" ";
00742 else
00743 cout << "\"Time = " << -time << " N-body \" ";
00744
00745 cout << depth << ", 0\n"
00746 << " pigment { Red }\n"
00747 << " translate < " << pos[0] << ", "
00748 << pos[1] << ", "
00749 << pos[2] << " >\n"
00750 << " scale " << scale
00751 << "}\n" << endl;
00752 cout.precision(p);
00753 }
00754
00755 local void print_start_text(vector cam_pos) {
00756 }
00757
00758 local void print_povray_stars(dyn *b, real mass_limit,
00759 real number_limit,
00760 bool povray,
00761 real scale_L,
00762 filter_type filter) {
00763
00764 cerr << "Print STARS"<<endl;
00765 for (int i=0; i<NNEBULAE; i++)
00766 nebula[i].id = -1;
00767
00768 bool cod = false;
00769
00770 vector dc_pos = 0;
00771 bool try_com = false;
00772 if(abs(dc_pos) == 0) {
00773 if (find_qmatch(b->get_dyn_story(), "density_center_pos")) {
00774
00775 if (getrq(b->get_dyn_story(), "density_center_time")
00776 != b->get_real_system_time()) {
00777 warning("mkpovfile: neglecting out-of-date density center");
00778 try_com = true;
00779 } else
00780 cod = true;
00781
00782 dc_pos = getvq(b->get_dyn_story(), "density_center_pos");
00783 }
00784
00785 if (try_com && find_qmatch(b->get_dyn_story(), "com_pos")) {
00786
00787 if (getrq(b->get_dyn_story(), "com_time")
00788 != b->get_real_system_time()) {
00789 warning("lagrad: neglecting out-of-date center of mass");
00790 } else
00791 dc_pos = getvq(b->get_dyn_story(), "com_pos");
00792 }
00793 }
00794
00795
00796 dc_pos = 0;
00797
00798 int ns=0, nb=0;
00799 for_all_daughters(dyn, b, bi)
00800 if (bi->is_leaf() &&
00801 ((bi->get_index()>0 && bi->get_index() <= number_limit) ||
00802 bi->get_starbase()->get_total_mass() >= mass_limit)) {
00803
00804 print_star(bi, bi->get_pos() - dc_pos, scale_L, filter);
00805 ns++;
00806 }
00807
00808 for_all_daughters(dyn, b, bi) {
00809 if (!bi->is_leaf())
00810 nb += print_povray_binary_recursive(bi, dc_pos,
00811 mass_limit, number_limit,
00812 povray, scale_L, filter);
00813 }
00814
00815 if (!povray && ns+nb==0)
00816 cout << " (none)\n";
00817 }
00818
00819 local void print_povray_bodies(dyn *b, real mass_limit,
00820 real number_limit, real mmax,
00821 real scale_M) {
00822
00823
00824
00825
00826
00827
00828
00829
00830
00831
00832
00833 int ns=0;
00834 for_all_leaves(dyn, b, bi)
00835 if(bi->get_mass() >= mass_limit) {
00836
00837 print_node(bi, bi->get_pos(), scale_M, mmax);
00838 ns++;
00839 }
00840
00841 if (ns==0)
00842 cout << "// (none)\n";
00843
00844 }
00845
00846 void print_povray_header(dyn* b, vector cam_pos,
00847 int camera_on_star_id,
00848 real aperture, real gamma,
00849 int blur_samples,
00850 int counter, real scale_L,
00851 int horizontal, int vertical,
00852 bool print_hrd) {
00853
00854 cout << "\n\n// End of POVRay scenary." << endl;
00855 cout << "\n\n// Start new POVRay scenary." << endl;
00856
00857 cout << "// filename = starlab_";
00858 print_filename_counter(counter, cout);
00859 cout << ".pov\n";
00860
00861 cout << "#include \"colors.inc\"" << endl
00862 << "#include \"shapes.inc\"" << endl
00863 << "#include \"textures.inc\"" << endl
00864 << "#include \"glass.inc\"" << endl;
00865 cout << "#include \"astro.inc\"" << endl;
00866
00867 cerr<< endl;
00868
00869 cout << "#version 3.0" << endl << endl;
00870
00871 cout << "global_settings { " << endl
00872 << " assumed_gamma " << gamma << endl
00873 << " max_trace_level 15" << endl
00874 << " ambient_light White" << endl
00875 << "}"
00876 << endl << endl;
00877
00878 if (camera_on_star_id<=0) {
00879
00880 vector normal;
00881 normal[1] = 1;
00882
00883
00884
00885
00886
00887 cout << "#declare camera_location = <0, 0, -1>" << endl;
00888 cout << "#declare normal_to_camera = -z" << endl;
00889 cout << "#declare camera_up = y" << endl;
00890
00891 cout << "camera {" << endl
00892 << " location < " << cam_pos[0] << ", "
00893 << cam_pos[1] << ", "
00894 << cam_pos[2] << " >" << endl
00895 << " look_at <0.0, 0.0, " << cam_pos[2]+5 << ">" << endl
00896 << " blur_samples " << blur_samples << endl;
00897 if (aperture>0)
00898 cout << " focal_point <0, 0, " << cam_pos[2]+5 << ">" << endl
00899 << " aperture " << aperture << endl;
00900 cout << "}" << endl << endl;
00901
00902
00903 cout << "light_source { < " << cam_pos[0] << ", "
00904 << cam_pos[1] << ", "
00905 << cam_pos[2] << " > White }" << endl;
00906
00907
00908
00909
00910
00911
00912 real time = -1*getrq(b->get_dyn_story(), "real_system_time");
00913
00914 vector time_pos;
00915 if(vertical<=400) {
00916 time_pos[0] = -10;
00917 time_pos[1] = 6;
00918 }
00919 else {
00920 time_pos[0] = -16;
00921 time_pos[1] = -11;
00922 }
00923
00924
00925
00926 time_pos[0] = -7;
00927 time_pos[1] = -5;
00928 time_pos[2] = 1;
00929
00930
00931
00932 real scale = 0.1;
00933 real depth = 0.125;
00934 print_povtime(time, time_pos, scale, depth);
00935
00936 if (print_hrd)
00937 print_hertzsprung_Russell_diagram(b, cam_pos);
00938
00939 }
00940 else {
00941 cam_pos = 0;
00942 if(!print_camera_position_recursive(b, cam_pos, camera_on_star_id, scale_L,
00943 aperture, blur_samples)) {
00944 cout << "No star id = " << camera_on_star_id << " found." << endl;
00945 cout << " in print_camera_position_recursive() " << endl;
00946 exit(-1);
00947 }
00948 }
00949
00950 cout << "// Here follow the definitions of the stars..." << endl
00951 << endl;
00952 }
00953
00954 void print_mpegplayer_param_file(ostream& s,
00955 int first_frame ,
00956 int last_frame ,
00957 int horizontal ,
00958 int vertical ,
00959 int GOP_size
00960 ) {
00961
00962
00963 s <<"# mpeg_encode parameter file\n"
00964 << "PATTERN IBBBPBBBBP\n"
00965 << "OUTPUT starlab.mpg\n"
00966 << endl;
00967
00968 s << "YUV_SIZE "<<vertical<<"x"<<horizontal<<"\n"
00969 << "BASE_FILE_FORMAT PPM\n"
00970 << "INPUT_CONVERT *\n"
00971 << "GOP_SIZE "<<GOP_size<<"\n"
00972 << "SLICES_PER_FRAME 1\n"
00973 << endl;
00974
00975 s << "INPUT_DIR .\n"
00976 << "INPUT\n"
00977 << "starlab_*.ppm [";
00978 print_filename_counter(first_frame, s);
00979 s << "-";
00980 print_filename_counter(last_frame-1, s);
00981 s << "]\n"
00982 << "END_INPUT\n" << endl;
00983
00984 s << "FORCE_ENCODE_LAST_FRAME\n"
00985
00986 << "PIXEL WHOLE\n"
00987 << "RANGE 10\n"
00988 << endl;
00989
00990 s << "PSEARCH_ALG LOGARITHMIC\n"
00991 << "BSEARCH_ALG CROSS2\n"
00992 << endl;
00993
00994 s << "IQSCALE 8\n"
00995 << "PQSCALE 10\n"
00996 << "BQSCALE 25\n"
00997 << endl;
00998
00999 s << "REFERENCE_FRAME ORIGINAL" << endl;
01000 }
01001
01002 void rdc_and_wrt_movie(dyn *b, bool povray, real scale_L, real mmax,
01003 char fltr) {
01004
01005 strcpy(collision[0].name, "255a+255b");
01006 collision[0].starid = 0;
01007 collision[0].time = 2.37956;
01008 strcpy(collision[1].name, "576a+576b");
01009 collision[1].starid = 0;
01010 collision[1].time = 3.79467;
01011 strcpy(collision[2].name, "463a+463b");
01012 collision[2].starid = 0;
01013 collision[2].time = 8.01511;
01014 strcpy(collision[3].name, "1747a+1747b");
01015 collision[3].starid = 0;
01016 collision[3].time = 8.97172;
01017 strcpy(collision[4].name, "1a+218b");
01018 collision[4].starid = 1;
01019 collision[4].time = 9.21207;
01020 strcpy(collision[5].name, "1b<+2>");
01021 collision[5].starid = 1;
01022 collision[5].time = 14.0345;
01023 strcpy(collision[6].name, "1a<+2>");
01024 collision[6].starid = 1;
01025 collision[6].time = 14.6999;
01026 strcpy(collision[7].name, "6a+6b");
01027 collision[7].starid = 0;
01028 collision[7].time = 15.9929;
01029 strcpy(collision[8].name, "1b<+3>");
01030 collision[8].starid = 1;
01031 collision[8].time = 22.9681;
01032 strcpy(collision[9].name, "14a+14b");
01033 collision[9].starid = 26;
01034 collision[9].time = 34.916;
01035
01036 #if 0
01037
01038 strcpy(collision[0].name, "1+167");
01039 collision[0].starid = 167;
01040 collision[0].time = 20.3479;
01041 strcpy(collision[1].name, "1<+2>");
01042 collision[1].starid = 69;
01043 collision[1].time = 21.1829;
01044 strcpy(collision[2].name, "30+64");
01045 collision[2].starid = 30;
01046 collision[2].time = 21.7115;
01047 strcpy(collision[3].name, "1<+3>");
01048 collision[3].starid = 4;
01049 collision[3].time = 23.6172;
01050 strcpy(collision[4].name, "1<+4>");
01051 collision[4].starid = 1435;
01052 collision[4].time = 25.2399;
01053 strcpy(collision[5].name, "155+938");
01054 collision[5].starid = 155;
01055 collision[5].time = 25.8208;
01056 strcpy(collision[6].name, "23+34");
01057 collision[6].starid = 23;
01058 collision[6].time = 26.3169;
01059 strcpy(collision[7].name, "6+403");
01060 collision[7].starid = 6;
01061 collision[7].time = 26.4171;
01062 strcpy(collision[8].name, "102+483");
01063 collision[8].starid = 102;
01064 collision[8].time = 26.901;
01065 strcpy(collision[9].name, "26+1007");
01066 collision[9].starid = 26;
01067 collision[9].time = 27.817;
01068 strcpy(collision[10].name, "29+173");
01069 collision[10].starid = 29;
01070 collision[10].time = 28.2324;
01071 strcpy(collision[11].name, "1<+6>");
01072 collision[11].starid = 2;
01073 collision[11].time = 28.234;
01074 strcpy(collision[12].name, "1<+7>");
01075 collision[12].starid = 1917;
01076 collision[12].time = 28.7764;
01077 strcpy(collision[13].name, "1<+8>");
01078 collision[13].starid = 46;
01079 collision[13].time = 29.3036;
01080 strcpy(collision[14].name, "1<+9>");
01081 collision[14].starid = 1562;
01082 collision[14].time = 29.6476;
01083 strcpy(collision[15].name, "1<+10>");
01084 collision[15].starid = 0;
01085 collision[15].time = 30.3179;
01086 strcpy(collision[16].name, "86+123");
01087 collision[16].starid = 0;
01088 collision[16].time = 31.8339;
01089 strcpy(collision[17].name, "1<+11>");
01090 collision[17].starid = 0;
01091 collision[17].time = 31.9796;
01092 strcpy(collision[18].name, "2+149");
01093 collision[18].starid = 0;
01094 collision[18].time = 32.3984;
01095 strcpy(collision[19].name, "1<+12>");
01096 collision[19].starid = 0;
01097 collision[19].time = 32.5874;
01098 strcpy(collision[20].name, "1<+13>");
01099 collision[20].starid = 0;
01100 collision[20].time = 32.7505;
01101 strcpy(collision[21].name, "15+19");
01102 collision[21].starid = 0;
01103 collision[21].time = 32.9756;
01104 strcpy(collision[22].name, "1<+14>");
01105 collision[22].starid = 0;
01106 collision[22].time = 33.4833;
01107 strcpy(collision[23].name, "1<+15>");
01108 collision[23].starid = 0;
01109 collision[23].time = 33.5243;
01110 strcpy(collision[24].name, "1<+16>");
01111 collision[24].starid = 0;
01112 collision[24].time = 33.7629;
01113 strcpy(collision[25].name, "1<+17>");
01114 collision[25].starid = 0;
01115 collision[25].time = 34.4345;
01116 strcpy(collision[26].name, "1<+18>");
01117 collision[26].starid = 0;
01118 collision[26].time = 34.4508;
01119 strcpy(collision[27].name, "1<+19>");
01120 collision[27].starid = 0;
01121 collision[27].time = 34.6628;
01122 strcpy(collision[28].name, "1<+20>");
01123 collision[28].starid = 0;
01124 collision[28].time = 34.7393;
01125 #endif
01126
01127
01128
01129
01130
01131
01132
01133
01134 #if 0
01135 sn_remnant[0].bh = true;
01136 sn_remnant[0].id = 4441;
01137 sn_remnant[0].time = 4.42;
01138 sn_remnant[1].bh = true;
01139 sn_remnant[1].id = 1;
01140 sn_remnant[1].time = 5.35;
01141 sn_remnant[1].bh = false;
01142 sn_remnant[2].id = 2056;
01143 sn_remnant[2].time = 6.14;
01144 sn_remnant[1].bh = false;
01145 sn_remnant[3].id = 3193;
01146 sn_remnant[3].time = 6.26;
01147 #endif
01148
01149 filter_type filter = get_filter_type(fltr);
01150
01151 real mass_limit = 0;
01152 real number_limit = VERY_LARGE_NUMBER;
01153 if(b->get_use_sstar())
01154 print_povray_stars(b, mass_limit, number_limit, povray, scale_L, filter);
01155 else
01156 print_povray_bodies(b, mass_limit, number_limit, mmax, scale_L);
01157
01158 }
01159
01160 #else
01161
01162
01163
01164
01165 enum SF_base_type {No_base=0, StTr_station, SW_ISD};
01166
01167 main(int argc, char ** argv)
01168 {
01169 bool c_flag = false;
01170 bool v_flag = false;
01171 bool B_flag = false;
01172 bool print_hrd = false;
01173 bool Stellar = true;
01174
01175 char filter = 'V';
01176 real aperture = -1;
01177 vector first_campos, last_campos, cam_pos;
01178 first_campos[2] = last_campos[2] = -10;
01179 real theta_rotate = 0;
01180 real phi_rotate = 0;
01181 int blur_samples = 10;
01182 int camera_on_star_id = -1;
01183 int nstart = 0;
01184
01185 char * mpeg_paramfile = "paramfile.mpeg";
01186 int horizontal = 90;
01187 int vertical = 120;
01188 int first_frame = 1;
01189 int last_frame = 0;
01190 real scale_L = 1;
01191
01192 real gamma = 2.2;
01193
01194 SF_base_type base_type = No_base;
01195 int nsteps = 32767;
01196 int nskip = 0;
01197 int nint_skip = 0;
01198
01199 char *comment;
01200 check_help();
01201
01202 extern char *poptarg;
01203 int c;
01204 char* param_string = "a:Bb:F:f:g:H:hI:L:N:n:W:X:x:Y:y:Z:z:P:T:S:s:c:";
01205
01206 while ((c = pgetopt(argc, argv, param_string)) != -1)
01207 switch(c)
01208 {
01209 case 'a': aperture = atof(poptarg);
01210 break;
01211
01212
01213 case 'B': Stellar = !Stellar;
01214 break;
01215 case 'b': blur_samples = atoi(poptarg);
01216 break;
01217 case 'F': filter = *poptarg;
01218 break;
01219 case 'f': mpeg_paramfile = poptarg;
01220 break;
01221 case 'g': gamma = atof(poptarg);
01222 break;
01223 case 'H': horizontal = atoi(poptarg);
01224 break;
01225 case 'h': print_hrd = true;
01226 break;
01227 case 'I': camera_on_star_id = atoi(poptarg);
01228 break;
01229 case 'W': vertical = atoi(poptarg);
01230 break;
01231 case 'L': scale_L *= atof(poptarg);
01232 break;
01233 case 'P': phi_rotate = atof(poptarg);
01234 phi_rotate *= cnsts.mathematics(pi)/180;
01235 break;
01236 case 'T': theta_rotate = atof(poptarg);
01237 theta_rotate *= cnsts.mathematics(pi)/180;
01238 break;
01239 case 'N': nsteps = atoi(poptarg);
01240 break;
01241 case 'n': nstart = atoi(poptarg);
01242 break;
01243 case 'X': first_campos[0] = atof(poptarg);
01244 break;
01245 case 'Y': first_campos[1] = atof(poptarg);
01246 break;
01247 case 'Z': first_campos[2] = atof(poptarg);
01248 break;
01249 case 'x': last_campos[0] = atof(poptarg);
01250 break;
01251 case 'y': last_campos[1] = atof(poptarg);
01252 break;
01253 case 'z': last_campos[2] = atof(poptarg);
01254 break;
01255 case 'S': nskip = atoi(poptarg);
01256 break;
01257 case 's': nint_skip = atoi(poptarg);
01258 break;
01259 case 'c': c_flag = true;
01260 comment = poptarg;
01261 break;
01262 case 'v': v_flag = true;
01263 break;
01264 case '?': params_to_usage(cout, argv[0], param_string);
01265 get_help();
01266 exit(1);
01267 }
01268
01269 dyn *b;
01270
01271 int GOP_size = min(1000, nsteps);
01272 remove(mpeg_paramfile);
01273 ofstream os(mpeg_paramfile, ios::app|ios::out);
01274 if (!os) cerr << "\nerror: couldn't create file "
01275 << mpeg_paramfile << endl;
01276
01277 print_mpegplayer_param_file(os, first_frame, GOP_size,
01278 horizontal, vertical, GOP_size);
01279 os.close();
01280
01281 int id;
01282 real time, mmax = -1;
01283 vector pos;
01284
01285 int nread = nstart;
01286 int j=0;
01287 for (int i = 0; i < nskip; i++){
01288 cerr << " skipping snapshot " << i << endl;
01289 if (! forget_node(cin)) exit(1);
01290 if(j==nint_skip) {
01291 nread++;
01292 j=0;
01293 }
01294 else
01295 j++;
01296 }
01297
01298 int nsnap = 0;
01299 while (b = get_dyn(cin)) {
01300
01301 if(nsnap==0)
01302 for_all_leaves(dyn, b, bi) {
01303 mmax = max(mmax, bi->get_mass());
01304 }
01305
01306 if(Stellar) {
01307 b->set_use_sstar(true);
01308 real T_start = 0;
01309 if(b->get_starbase()->get_stellar_evolution_scaling()) {
01310
01311
01312
01313
01314 }
01315 else {
01316 cerr << "No stellar evolution scaling present"<<endl;
01317 cerr << "Continue with node prescription"<<endl;
01318 }
01319
01320 }
01321
01322
01323
01324
01325 nsnap ++;
01326 nread ++;
01327 cerr << " reading snapshot " << nread << ", and "
01328 << nsnap << " done." << endl;
01329
01330 if (camera_on_star_id<=0)
01331 new_camera_position(cam_pos, first_campos, last_campos,
01332 nsteps, nsnap);
01333
01334 print_povray_header(b, cam_pos, camera_on_star_id,
01335 aperture, gamma, blur_samples, nread, scale_L,
01336 horizontal, vertical, print_hrd);
01337
01338 if (base_type == StTr_station) {
01339
01340 cout << "#include \"starbase.inc\"" << endl;
01341 cout << "object { Base_Station }\n" << endl;
01342 }
01343 else if (base_type == SW_ISD) {
01344 cout << "object { StarWars_Scene }\n" << endl;
01345 }
01346
01347 rdc_and_wrt_movie(b, true, scale_L, mmax, filter);
01348 last_frame++;
01349
01350 if (v_flag)
01351 put_dyn(cout, *b);
01352
01353 rmtree(b);
01354
01355 for (int i = 0; i < nint_skip; i++){
01356 nsnap ++;
01357 cerr << " skipping snapshot " << nsnap << endl;
01358 if (! forget_node(cin)) exit(1);
01359 }
01360
01361 if (nsnap==nsteps)
01362 break;
01363
01364 }
01365
01366 }
01367
01368 #endif
01369
01370
01371
