sdyn3_ev.C Source File

00001 
00002 // sdyn3_ev.C
00003 
00004 // Member functions for the sdyn3 class that are related to orbit integration.
00005 // Also include some useful general functions (move elsewhere?...)
00006 
00007 #include "sdyn3.h"
00008 
00009 real potential_energy(sdyn3 * b)
00010 {
00011     real u = 0;
00012 
00013     for_all_daughters(sdyn3, b, bi) {
00014         for (sdyn3 * bj = bi->get_younger_sister(); bj != NULL;
00015              bj = bj->get_younger_sister())
00016           u -= bi->get_mass() * bj->get_mass()
00017                / abs(bi->get_pos() - bj->get_pos());
00018     }
00019 
00020     return u;
00021 }
00022 
00023 real energy(sdyn3 * b)
00024 {
00025     real k = 0, u = 0, dissipation = 0;
00026 
00027     for_all_daughters(sdyn3, b, bi) {
00028 
00029         // PRC(bi), PRL(bi->get_pos());
00030 
00031         k += bi->get_mass() * bi->get_vel() * bi->get_vel();
00032         dissipation += bi->get_energy_dissipation();
00033 
00034         for (sdyn3 * bj = bi->get_younger_sister(); bj != NULL;
00035              bj = bj->get_younger_sister()) {
00036 
00037             // PRC(bj), PRL(bj->get_pos());
00038 
00039             u -= bi->get_mass() * bj->get_mass()
00040                   / abs(bi->get_pos() - bj->get_pos());
00041         }
00042     }
00043 
00044     return 0.5*k + u + dissipation;
00045 }
00046 
00047 real total_angular_momentum(sdyn3* b)
00048 {
00049     vector am = vector(0);
00050     for_all_daughters(sdyn3, b, bb)
00051         am += bb->get_mass() * (bb->get_pos() ^ bb->get_vel());
00052     return abs(am);
00053 }
00054 
00055 void sdyn3::accumulate_new_acc_and_jerk_from_new(
00056                         sdyn3* bj,                  // n-body system pointer
00057                         real eps2,                  // softening length squared
00058                         int no_ssd_flag,
00059                         int& collision_flag)
00060 {
00061     if (bj->get_oldest_daughter() != NULL)
00062         for_all_daughters(sdyn3, bj, bb)
00063             accumulate_new_acc_and_jerk_from_new(bb, eps2,
00064                                                  no_ssd_flag, collision_flag);
00065     else
00066         if (this != bj) {
00067 
00068             vector d_pos = new_pos - bj->get_new_pos();
00069             vector d_vel = new_vel - bj->get_new_vel();
00070             real r2 = d_pos*d_pos;
00071 
00072             if (r2 < (radius + bj->get_radius()) * (radius + bj->get_radius()))
00073                 collision_flag = 1;
00074 
00075             real r2inv = 1.0/(r2 + eps2);
00076             real a3 = -3.0*(d_pos*d_vel)*r2inv;
00077             real rinv = sqrt(r2inv);
00078             real mrinv = bj->get_mass() * rinv;
00079             new_pot -= mrinv;
00080             real mr3inv = mrinv * r2inv;
00081             new_acc -= mr3inv*d_pos;
00082             new_jerk -= mr3inv*(d_vel + a3*d_pos);
00083 
00084             if (!no_ssd_flag) {
00085 
00086                 // Update nn and oscillation data.
00087 
00088                 if (r2 < nn_dr2) {
00089                     nn_dr2 = r2;
00090                     nn_label = bj->get_index();
00091                     nn = bj;
00092                 }
00093                 if (r2 < min_nn_dr2) {
00094                     min_nn_dr2 = r2;
00095                     min_nn_label = bj->get_index();
00096                 }
00097 
00098                 ((sdyn3 *) parent)->inc_ssd(0.5*r2);
00099             }                                // correct for double counting
00100 
00101             // Determine characteristic time scales:
00102 
00103             real encounter_time_squared = 1 / (r2inv * d_vel * d_vel);
00104             if (min_encounter_time_sq > encounter_time_squared)
00105                 min_encounter_time_sq = encounter_time_squared;
00106 
00107             real inverse_free_fall_time_squared = 
00108                 r2inv * rinv * (mass + bj->get_mass());
00109             if (min_free_fall_time_sq * inverse_free_fall_time_squared > 1)
00110                 min_free_fall_time_sq = 1/inverse_free_fall_time_squared;
00111         }
00112 }
00113 
00114 void sdyn3::calculate_new_acc_and_jerk_from_new(sdyn3* b,
00115                                                 real eps_squared,
00116                                                 int  no_ssd_flag,
00117                                                 int& collision_flag)
00118 {
00119     if (oldest_daughter != NULL) {
00120 
00121         collision_flag = 0;
00122 
00123         if (!no_ssd_flag){
00124             ssd = 0;
00125             nn_label = 0;
00126             for_all_daughters(sdyn3, this, c)
00127                 c->set_nn_dr2(VERY_LARGE_NUMBER);
00128         }
00129 
00130         for_all_daughters(sdyn3, this, bb)
00131             bb->calculate_new_acc_and_jerk_from_new(b, eps_squared,
00132                                                     no_ssd_flag,
00133                                                     collision_flag);
00134         if (!no_ssd_flag)
00135             for_all_daughters(sdyn3, this, d) {
00136                 if (d->get_init_nn_label() <= 0)
00137                     d->set_init_nn_label(d->get_nn_label());
00138                 if (d->get_init_nn_label() != d->get_nn_label())
00139                         d->set_nn_change_flag(1);
00140             }
00141 
00142         if (!no_ssd_flag) {
00143 
00144             if (min_min_ssd > ssd) min_min_ssd = ssd;
00145             
00146             if (ssd_ingoing_flag) {
00147                 if (min_ssd > ssd) min_ssd = ssd;
00148                 if (ssd > min_ssd * SSD_HYSTERESIS_FACTOR)
00149                     ssd_ingoing_flag = 0, max_ssd = ssd, n_ssd_osc++;
00150             } else {
00151                 if (max_ssd < ssd) max_ssd = ssd;
00152                 if (ssd < max_ssd / SSD_HYSTERESIS_FACTOR)
00153                     ssd_ingoing_flag = 1, min_ssd = ssd;
00154             }
00155             
00156             real e_tot = 0;
00157             for_all_daughters(sdyn3, this, bbb)
00158                 e_tot += 0.5 * bbb->get_mass() * (bbb->get_new_pot()
00159                                     + bbb->get_new_vel() * bbb->get_new_vel());
00160             if (de_tot_abs_max < abs(e_tot - e_tot_init))
00161                 de_tot_abs_max = abs(e_tot - e_tot_init);
00162         }
00163 
00164     } else {
00165 
00166         new_pot = 0;
00167         new_acc = new_jerk = 0;
00168         min_encounter_time_sq = min_free_fall_time_sq = VERY_LARGE_NUMBER;
00169         accumulate_new_acc_and_jerk_from_new(b, eps_squared, no_ssd_flag,
00170                                              collision_flag);
00171     }
00172 }
00173 
00174 void  sdyn3::taylor_pred_new_pos_and_vel(const real dt)
00175 {
00176     taylor_pred_new_pos(dt);
00177     taylor_pred_new_vel(dt);
00178 }
00179 
00180 void  sdyn3::correct_new_acc_and_jerk(const real new_dt, const real prev_new_dt)
00181 {
00182     if (new_dt == prev_new_dt) return;
00183 
00184     real dt_off = new_dt - 0.5 * prev_new_dt;
00185                                     // offset from midpoint of prev_new_dt step
00186     real theta = 0.5 * prev_new_dt;
00187     real tau = dt_off / theta;          // equals 1 if new_dt = prev_new_dt
00188 
00189     real inv_theta = 1 / theta;
00190     real tau2 = tau * tau;
00191     real tau3 = tau2 * tau;
00192     vector prev_new_acc = new_acc;
00193     vector prev_new_jerk = new_jerk;
00194 
00195     new_acc = 0.25 * (
00196                       acc * (2 - 3 * tau + tau3)
00197                       + prev_new_acc * (2 + 3 * tau - tau3)
00198                       + jerk * theta * (1 - tau - tau2 + tau3)
00199                       + prev_new_jerk * theta * (-1 - tau + tau2 + tau3)
00200                       );
00201 
00202     new_jerk = 0.25 * (
00203                        acc * inv_theta * 3 * (-1 + tau2)
00204                        + prev_new_acc * inv_theta * 3 * (1 - tau2)
00205                        + jerk * (-1 - 2*tau + 3*tau2)
00206                        + prev_new_jerk * (-1 + 2*tau + 3*tau2)
00207                        );
00208 }
00209 
00210 void  sdyn3::correct_new_pos_and_vel(const real new_dt)  // timestep
00211 {
00212     real new_dt2 = new_dt * new_dt;
00213     real new_dt3 = new_dt2 * new_dt;
00214 
00215     new_vel = vel + new_dt * (new_acc + acc)/2
00216                   - new_dt2 * (new_jerk - jerk)/12;
00217 
00218     new_pos = pos + new_dt * (new_vel + vel)/2 
00219                   - new_dt2 * (new_acc - acc)/12;          // alpha = -5/3
00220 }
00221 
00222 //----------------------------------------------------------------------
00223 
00224 // Useful functions (used by low_n_evolve and scatter3):
00225 
00226 void set_kepler_from_sdyn3(kepler& k, sdyn3* b1, sdyn3* b2)
00227 {
00228     if (b1->get_time() != b2->get_time()) 
00229         err_exit("set_kepler_from_sdyn3: inconsistent times");
00230 
00231     k.set_time(b1->get_time());
00232     k.set_total_mass( b1->get_mass() + b2->get_mass() );
00233     k.set_rel_pos( b2->get_pos() - b1->get_pos() );
00234     k.set_rel_vel( b2->get_vel() - b1->get_vel() );
00235 
00236     k.initialize_from_pos_and_vel();
00237 }
00238 
00239 void kepler_pair_to_triple(kepler & k1, // Inner binary (b1 + b2)
00240                            kepler & k2, // Outer binary
00241                            sdyn3 * b1,
00242                            sdyn3 * b2,
00243                            sdyn3 * b3)
00244 {
00245     // Construct positions and velocities from the specified kepler structures.
00246 
00247     // Assume that labels, masses, times, and the center of mass are handled
00248     // elsewhere.
00249 
00250     // Set up the outer orbit first.  Note that the center of mass of
00251     // the three-body system is taken to be at rest at the origin.
00252 
00253     // Make no assumptions about masses!
00254 
00255     real m1 = b1->get_mass();
00256     real m2 = b2->get_mass();
00257     real m3 = b3->get_mass();
00258 
00259     real m12 = m1 + m2;
00260     real m123 = m12 + m3;
00261 
00262     // Relative position and velocity from (1,2) to 3
00263 
00264     b3->set_pos(m12 * k2.get_rel_pos() / m123);
00265     b3->set_vel(m12 * k2.get_rel_vel() / m123);
00266     
00267     b1->set_pos(-m3 * k2.get_rel_pos() / m123);
00268     b1->set_vel(-m3 * k2.get_rel_vel() / m123);
00269 
00270     b2->set_pos(b1->get_pos());
00271     b2->set_vel(b1->get_vel());
00272 
00273     // Then set up the inner binary.
00274 
00275     b1->inc_pos(-m2 * k1.get_rel_pos() / m12);      // rel_pos is from 1 to 2
00276     b1->inc_vel(-m2 * k1.get_rel_vel() / m12);
00277 
00278     b2->inc_pos( m1 * k1.get_rel_pos() / m12);
00279     b2->inc_vel( m1 * k1.get_rel_vel() / m12);
00280 
00281 }