OortAB.C Source File

00001 
00002 
00003 
00004 
00005 
00006 
00007 
00008 
00009 
00010 
00011 
00012 
00013 #include "stdfunc.h"
00014 #ifndef TOOLBOX
00015 
00016 #else
00017 
00018 //local real Oort_AB(const real X, int nfit, const real Xarray[],
00019 //                 const real Yarray[]) {
00020 //}
00021 
00022 local real galaxy_mass_within_radius(real r) {
00023 
00024   real m_galaxy = 4.25E+6 * pow(r, 1.2); // Msun
00025 
00026   return m_galaxy;
00027 }
00028 
00029 
00030 local void print_Oort_constants(real r_gc, real m_cluster) {
00031 
00032   real m_galaxy = galaxy_mass_within_radius(r_gc);
00033 
00034   real r_tide = pow(m_cluster/(2*m_galaxy), ONE_THIRD) * r_gc;  // pc
00035 
00036   real v_c = 0.0657 * sqrt((m_galaxy+m_cluster)/r_gc);  // km/s
00037   real dvc_dr = 13.7 / pow(r_gc, 0.9);                  // km/s
00038 
00039   real Oort_A = (v_c/r_gc - dvc_dr)/2;                  // km/s/pc
00040   real Oort_B = -(v_c/r_gc + dvc_dr)/2;                 // km/s/pc
00041 
00042   real dr = 0.5;     // [pc]
00043   real rho_G = (galaxy_mass_within_radius(r_gc-dr)
00044                 - galaxy_mass_within_radius(r_gc+dr))
00045              / ((4*PI/3) * (pow(r_gc-dr, 3) - pow(r_gc+dr, 3)));
00046 
00047   cerr << "    r_gc= "<<r_gc<<" m_cluster= "<<m_cluster
00048        << " m_galaxy= "<<m_galaxy
00049        << " A= " <<Oort_A<< " B= " <<Oort_B<<" rho_G= " << rho_G << endl;
00050   cerr << " r_tide= "<<r_tide<<" v_c= "<<v_c<<" dvc_dr= "<< dvc_dr << endl;
00051 }
00052 
00053 
00054 main(int argc, char **argv) {
00055 
00056   bool m_flag = false;
00057   bool n_flag = false;
00058   bool r_flag = false;
00059 
00060   real r_min = 30;     // [pc]
00061   real r_max = 30;     // [pc]
00062   real m_cluster = 20000;  // [Msun]
00063   int n = 1;
00064 
00065   extern char *poptarg;
00066   int c;
00067   char* param_string = "M:m:N:n:R:r:";
00068 
00069     while ((c = pgetopt(argc, argv, param_string)) != -1)
00070         switch(c) {
00071             case 'R': r_max = atof(poptarg);
00072                       r_flag = true;
00073                       break;
00074             case 'r': r_min = atof(poptarg);
00075                       r_flag = true;
00076                       break;
00077             case 'N': 
00078             case 'n': n = atoi(poptarg);
00079                       n_flag = true;
00080                       break;
00081             case 'M':
00082             case 'm': m_cluster = atof(poptarg);
00083                       m_flag = true;
00084                       break;
00085             case '?': params_to_usage(cerr, argv[0], param_string);
00086                       exit(1);
00087         }
00088 
00089   if (r_min>r_max) {
00090     r_max = r_min;
00091     n=1;
00092   }
00093 
00094   cerr << "Calculate Oort constants. " << endl; 
00095    if (r_flag && m_flag) {
00096 
00097      real r_gc;
00098      real dr = (r_max-r_min)/n;
00099 
00100      for(int i=0; i<n; i++) { 
00101        r_gc = r_min + dr*i;
00102        print_Oort_constants(r_gc, m_cluster);
00103      }
00104    }
00105 }
00106 #endif