refine_mass.C Source File

00001 
00002        //=======================================================//    _\|/_
00003       //  __  _____           ___                    ___       //      /|\ ~
00004      //  /      |      ^     |   \  |         ^     |   \     //          _\|/_
00005     //   \__    |     / \    |___/  |        / \    |___/    //            /|\ ~
00006    //       \   |    /___\   |  \   |       /___\   |   \   // _\|/_
00007   //     ___/   |   /     \  |   \  |____  /     \  |___/  //   /|\ ~
00008  //                                                       //            _\|/_
00009 //=======================================================//              /|\ ~
00010 
00011 // refine_mass.C:  More careful determination of cluster mass.
00012 //                 Also flag escapers.
00013 //
00014 //                 Separate from refine_mass2 because the solution
00015 //                 for the Jacobi radius is analytical, and the
00016 //                 form of the potential is particularly simple.
00017 //
00018 // Externally visible functions:
00019 //
00020 //      void refine_cluster_mass
00021 
00022 #include "dyn.h"
00023 
00024 #define G_TOL 1.e-6
00025 
00026 local real solve3(real a, real b, real c)
00027 {
00028     // Iteratively solve
00029     //
00030     //          a z^3 + b z + c  =  0
00031     //
00032     // where we know a > 0, b > 0, c < 0.
00033 
00034     // Proceed by Newton-Raphson iteration.  Since a > 0, b > 0, so
00035     // there are are no turning points to complicate the logic.
00036 
00037     real z = 2*max(pow(abs(c/a), 1.0/3), sqrt(abs(b/a)));
00038     real g = c + z * (b + a*z*z);
00039     real gpr = b + 3*a*z*z;             // always > 0
00040 
00041     while (abs(g) > G_TOL) {
00042         z = z - g/gpr;
00043         g = c + z * (b + a*z*z);
00044         gpr = b + 3*a*z*z;
00045     }
00046 
00047     return z;
00048 }
00049 
00050 #define M_TOL 1.e-4
00051 #define M_ITER_MAX 20
00052 
00053 void refine_cluster_mass(dyn *b,
00054                          int verbose)           // default = 0
00055 {
00056     if (b->get_external_field() == 0) return;
00057 
00058     if (b->get_tidal_field() == 0) {
00059         refine_cluster_mass2(b, verbose);       // experimental
00060         return;
00061     }
00062 
00063     // Self-consistently determine the total mass within the outermost
00064     // closed zero-velocity surface.  Use a point-mass approximation for
00065     // the cluster potential and iterate until the actual mass within
00066     // the surface agrees with the mass used to generate the surface.
00067     // The total potential is
00068     //
00069     //          phi  =  -GM/r + (alpha1 x^2 + alpha3 z^2) / 2
00070     //
00071     // where we measure everything relative to the standard center (which
00072     // is the density center, if known, and the modified center of mass
00073     // otherwise).  The Jacobi radius for this field is
00074     //
00075     //          r_J  =  (-GM/alpha1)^{1/3}
00076     //
00077     // and the potential of the last closed surface is
00078     //
00079     //          phi_J  =  1.5 alpha1 r_J^2.
00080 
00081     vector center, vcenter;
00082     int which = get_std_center(b, center, vcenter);
00083 
00084     // which = 1 for density center, 2 for mcom.
00085 
00086     // Note from Steve (7/01):  The "center" should really be the
00087     // center of mass of particles within the Jacobi surface,
00088     // determined self-consistently.
00089     //
00090     // To do...  (See also check_and_remove_escapers().)
00091 
00092     real M_inside = total_mass(b), M = -1, M0 = M_inside;
00093     real r_J, r_x2, r_y2, r_z2, r_max_inside;
00094     int  N_inside, iter = 0;
00095 
00096     real phi_J;
00097 
00098     real M_J, M_x, M_y, M_z;
00099     int  N_J, N_x, N_y, N_z;
00100 
00101     if (verbose)
00102         cerr << endl << "  refine_cluster_mass: getting M by iteration"
00103              << endl << "  initial total system mass = " << M_inside
00104              << endl;
00105 
00106     while (iter++ < M_ITER_MAX
00107            && M_inside > 0
00108            && abs(M_inside/M - 1) > M_TOL) {
00109 
00110         M = M_inside;
00111         M_inside = 0;
00112         M_J = M_x = M_y = M_z = 0;
00113         N_J = N_x = N_y = N_z = 0;
00114 
00115         r_J = pow(-M/b->get_alpha1(), 0.3333333);
00116 
00117         phi_J = 1.5 * b->get_alpha1() * r_J * r_J;
00118         real r_max = 0;
00119 
00120         r_x2 = square(r_J);             // zero-velocity surface crosses x-axis
00121         r_y2 = square(-M/phi_J);        // zero-velocity surface crosses y-axis
00122 
00123         // zero-velocity surface crosses z-axis where
00124         //
00125         //      -GM/z + 0.5 alpha3 z^2  =  phi_J
00126         // so
00127         //      0.5 alpha3 z^3 -phi_J z -GM  =  0
00128 
00129         r_z2 = square(solve3(0.5*b->get_alpha3(), -phi_J, -M));
00130 
00131         N_inside = 0;
00132         r_max_inside = 0;
00133 
00134         for_all_daughters(dyn, b, bb) {
00135 
00136             vector dx = bb->get_pos() - center;
00137 
00138             real x = dx[0];
00139             real y = dx[1];
00140             real z = dx[2];
00141             real r = abs(dx);
00142 
00143             if (r < r_J) {
00144 
00145                 N_J++;
00146                 M_J += bb->get_mass();
00147 
00148                 if (r == 0
00149                     || -M/r + 0.5 * (b->get_alpha1()*x*x + b->get_alpha3()*z*z)
00150                             < phi_J) {
00151                     N_inside++;
00152                     M_inside += bb->get_mass();
00153                     r_max_inside = max(r, r_max_inside);
00154                 }
00155 
00156             }
00157             r_max = max(r, r_max);
00158 
00159             // Count projected masses and numbers.
00160 
00161             real xy = x*x + y*y;
00162             real xz = x*x + z*z;
00163             real yz = y*y + z*z;
00164 
00165             if (xy < r_x2) {            // z projection
00166                 M_x += bb->get_mass();
00167                 N_x++;
00168             }
00169             if (xz < r_x2) {            // y projection
00170                 M_x += bb->get_mass();
00171                 N_x++;
00172             }
00173 
00174             if (yz < r_y2) {            // x projection
00175                 M_y += bb->get_mass();
00176                 N_y++;
00177             }
00178             if (xy < r_y2) {            // z projection
00179                 M_y += bb->get_mass();
00180                 N_y++;
00181             }
00182 
00183             if (yz < r_z2) {            // x projection
00184                 M_z += bb->get_mass();
00185                 N_z++;
00186             }
00187             if (xz < r_z2) {            // y projection
00188                 M_z += bb->get_mass();
00189                 N_z++;
00190             }
00191         }
00192 
00193         M_x /= 2;
00194         N_x /= 2;
00195         M_y /= 2;
00196         N_y /= 2;
00197         M_z /= 2;
00198         N_z /= 2;
00199 
00200 #if 0
00201         fprintf(stderr, "    %2d  ", iter);
00202         PRC(M); PRC(r_J); PRL(r_max);
00203         PRI(8); PRC(N_inside); PRC(M_inside); PRL(r_max_inside);
00204 #endif
00205 
00206     }
00207 
00208     if (iter >= M_ITER_MAX) {
00209         if (verbose) PRI(2);
00210         cerr << "warning: refine_cluster_mass: too many iterations at time "
00211              << b->get_system_time() << endl;
00212     }
00213 
00214     if (verbose)
00215         cerr << "  within last zero-velocity surface:"
00216              << "  M = " << M_inside << "  N = " << N_inside
00217              << endl
00218              << "                                    "
00219              << "  r_J = " << r_J << "  r_max = " << r_max_inside
00220              << endl
00221              << "  within r_J:                       "
00222              << "  M = " << M_J << "  N = " << N_J
00223              << endl
00224              << "  within r_J (projected):           "
00225              << "  M = " << M_x << "  N = " << N_x
00226              << endl
00227              << "  within r_y (projected):           "
00228              << "  M = " << M_y << "  N = " << N_y << "  r_y = " << sqrt(r_y2)
00229              << endl
00230              << "  within r_z (projected):           "
00231              << "  M = " << M_z << "  N = " << N_z << "  r_z = " << sqrt(r_z2)
00232              << endl;
00233 
00234     // Repeat the inner loop above and flag stars as escapers or not.
00235     // Too many iterations probably means a limit cycle of some sort;
00236     // accept the results in that case.
00237 
00238     bool disrupted = (M_inside < 0.01*M0);
00239 
00240     for_all_daughters(dyn, b, bb) {
00241 
00242         vector dx = bb->get_pos() - center;
00243 
00244         real x = dx[0];
00245         real z = dx[2];
00246         real r = abs(dx);
00247 
00248         bool escaper = true;
00249         if (!disrupted) {
00250             if (r < r_J
00251                 && (r == 0 || -M/r + 0.5 * (b->get_alpha1()*x*x
00252                                              + b->get_alpha3()*z*z) < phi_J))
00253                 escaper = false;
00254         }
00255         putiq(bb->get_dyn_story(), "esc", escaper);
00256 
00257         // Print out t_esc if this is a new escaper.
00258 
00259         if (escaper && !find_qmatch(bb->get_dyn_story(), "t_esc"))
00260             putrq(bb->get_dyn_story(), "t_esc", (real)b->get_system_time());
00261 
00262         // Delete t_esc if this is not an escaper.
00263 
00264         if (!escaper && find_qmatch(bb->get_dyn_story(), "t_esc"))
00265             rmq(bb->get_dyn_story(), "t_esc");
00266     }
00267 }