diff --git a/file.dat b/file.dat new file mode 100644 index 0000000..6ba8741 --- /dev/null +++ b/file.dat @@ -0,0 +1,77 @@ +# subhalo 411321 +# snap, t [Gyr], M [MSun], rh [kpc] +25 2.145 1.842577e+11 2.742022e+01 +26 2.238 1.892209e+11 2.671429e+01 +27 2.384 1.931019e+11 2.774201e+01 +28 2.539 2.084169e+11 2.707294e+01 +29 2.685 2.154326e+11 2.746115e+01 +30 2.839 2.560636e+11 2.403104e+01 +31 2.981 3.377362e+11 2.154270e+01 +32 3.129 3.875323e+11 2.139171e+01 +33 3.285 4.428067e+11 2.109969e+01 +34 3.447 4.814599e+11 2.066491e+01 +35 3.593 5.435782e+11 1.958987e+01 +36 3.744 5.832241e+11 1.588560e+01 +37 3.902 5.910308e+11 1.318808e+01 +38 4.038 6.372445e+11 1.241970e+01 +39 4.206 6.685910e+11 1.169071e+01 +40 4.293 6.742185e+11 1.183042e+01 +41 4.502 6.817118e+11 1.224346e+01 +42 4.657 7.010794e+11 1.133662e+01 +43 4.816 7.106476e+11 1.083693e+01 +44 4.980 7.341724e+11 1.081872e+01 +45 5.115 7.524877e+11 1.087563e+01 +46 5.289 8.035750e+11 1.090464e+01 +47 5.431 8.380413e+11 1.074618e+01 +48 5.577 8.604913e+11 1.094197e+01 +49 5.726 9.188854e+11 1.095124e+01 +50 5.878 9.939528e+11 1.101752e+01 +51 6.073 1.050354e+12 1.165214e+01 +52 6.193 1.106218e+12 1.186527e+01 +53 6.356 1.173352e+12 1.247234e+01 +54 6.522 1.255338e+12 1.303894e+01 +55 6.692 1.383761e+12 1.599427e+01 +56 6.822 1.373275e+12 1.645093e+01 +57 6.998 1.422295e+12 1.514309e+01 +58 7.132 1.462187e+12 1.522601e+01 +59 7.314 1.523746e+12 1.505336e+01 +60 7.453 1.625764e+12 1.554669e+01 +61 7.642 1.715243e+12 1.789185e+01 +62 7.786 1.705854e+12 1.600661e+01 +63 7.932 1.748553e+12 1.663414e+01 +64 8.079 1.708907e+12 1.684767e+01 +65 8.280 1.682277e+12 1.765823e+01 +66 8.432 1.674851e+12 1.764632e+01 +67 8.587 1.679583e+12 1.768747e+01 +68 8.743 1.670373e+12 1.757611e+01 +69 8.902 1.685143e+12 1.758471e+01 +70 9.062 1.693599e+12 1.784006e+01 +71 9.225 1.707019e+12 1.780380e+01 +72 9.389 1.700025e+12 1.720548e+01 +73 9.556 1.716124e+12 1.685487e+01 +74 9.724 1.698637e+12 1.614731e+01 +75 9.837 1.690106e+12 1.639938e+01 +76 10.009 1.713624e+12 1.581026e+01 +77 10.182 1.725170e+12 1.540258e+01 +78 10.299 1.725714e+12 1.529604e+01 +79 10.535 1.741246e+12 1.508748e+01 +80 10.654 1.786631e+12 1.489346e+01 +81 10.834 1.775003e+12 1.485126e+01 +82 11.016 1.864654e+12 1.507335e+01 +83 11.138 1.871185e+12 1.536455e+01 +84 11.323 1.885768e+12 1.537435e+01 +85 11.509 1.896120e+12 1.494602e+01 +86 11.635 1.917398e+12 1.486861e+01 +87 11.824 1.948357e+12 1.470634e+01 +88 11.951 1.953813e+12 1.459721e+01 +89 12.143 2.008370e+12 1.469724e+01 +90 12.337 2.004639e+12 1.439044e+01 +91 12.467 2.029410e+12 1.444626e+01 +92 12.663 2.029940e+12 1.426860e+01 +93 12.795 2.056368e+12 1.414617e+01 +94 12.993 2.067279e+12 1.379220e+01 +95 13.127 2.078900e+12 1.382175e+01 +96 13.328 2.182388e+12 1.354795e+01 +97 13.463 2.223654e+12 1.352882e+01 +98 13.667 2.275025e+12 1.345206e+01 +99 13.803 2.306498e+12 1.335429e+01 diff --git a/loadtxt.cpp b/loadtxt.cpp new file mode 100644 index 0000000..939e0f9 --- /dev/null +++ b/loadtxt.cpp @@ -0,0 +1,88 @@ +#include +#include +#include +#include + +//TODO if cols is an empty vector, get all columns from the file +//TODO error checking + +class Loadtxt { +public: + Loadtxt(std::string file_name, std::vector cols) + { + std::sort(cols.begin(), cols.end()); + n_cols = cols.size(); + const int tmp_number_of_rows = 16384; + int n_rows_alloc = tmp_number_of_rows; + buffer = (double*)malloc(n_cols * sizeof(double) * n_rows_alloc); + std::ifstream file(file_name); + std::string line; + int row = -1; + while (getline(file, line)) { + if (line[line.find_first_not_of(whitespaces)]=='#') continue; + if (++row >= n_rows_alloc) { + n_rows_alloc += tmp_number_of_rows; + buffer = (double*)realloc(buffer, n_cols * sizeof(double) * n_rows_alloc); + } + line_to_buf(cols, line, buffer + row*n_cols); + } + file.close(); + buffer = (double*)realloc(buffer, n_cols * sizeof(double) * (++row)); + n_rows = row; + } + ~Loadtxt() + { + free(buffer); + } + std::vector> get_cols() + { + std::vector> data(n_cols); + for (int col=0; col(n_rows); + for (int row=0; row cols, std::string line, double *buffer) + { + int n_cols = cols.size(); + line = line.substr(line.find_first_not_of(whitespaces)); + auto pos = line.find_first_of(whitespaces, 1); + int col=0, idx=0; + std::vector data; + while (pos != std::string::npos) { + std::string num_as_string; + num_as_string = line.substr(0, pos); + pos = line.find_first_not_of(whitespaces, pos); + line = line.substr(pos, std::string::npos); + pos = line.find_first_of(whitespaces, 1); + if (col++ == cols[idx]) buffer[idx++] = std::stod(num_as_string); + } + if (col++ == cols[idx]) buffer[idx++] = std::stod(line); + if (idx < n_cols) throw; + } + double *buffer; + int n_rows, n_cols; +}; + +// Below is a deomonstration. The file has multiple columns but we are only +// interested in the second and fourth. We pass the file name and the column +// vector {1, 3} (since numbering starts at zero) and immediately call the +// get_cols() member, which returns a vector of vectors. We then manually assign +// the data members into named vectors. +// +// #include +// int main() +// { + +// auto data = Loadtxt("file.dat", {1, 3}).get_cols(); +// auto& time = data[0]; +// auto& value = data[1]; +// for (size_t i=0; i t, M_bulge_data, b_bulge_data; + std::vector t_data, M_halo_data, b_halo_data; - t.push_back(0); - t.push_back(1); - t.push_back(2); - - M_bulge_data.push_back(9e9); - M_bulge_data.push_back(1e10); - M_bulge_data.push_back(2e10); - - interp_M_bulge = new Interp(t, M_bulge_data); - std::cout << "exiting..." << std::endl; - exit(0); + std::ifstream file(file_name); + std::string line; + while (std::getline(file, line)) { + auto pos = line.find('#'); + if (pos != std::string::npos) line = line.substr(0, pos); + pos = line.find_first_not_of(" \t"); + if (pos == std::string::npos) continue; + double data[3]; + sscanf(line.c_str(), "%*s %lf %lf %lf", &data[0], &data[1], &data[2]); + t_data.push_back(data[0]); + M_halo_data.push_back(data[1]); + b_halo_data.push_back(data[2]); // note, this is not half-mass radius + } + interp_M_halo = new Interp(t_data, M_halo_data); + interp_b_halo = new Interp(t_data, b_halo_data); } int func(double t, const double y[], double f[], void *params) { - double M_bulge = (*interp_M_bulge)(t); - Plummer plummer(M_bulge, 0); + double M_halo = (*interp_M_halo)(t); + double b_halo = (*interp_b_halo)(t); + /*printf("xxxxxxxxx %e, %e msun\n", t, M_halo); + printf("xxxxxxxxx %e, %e kpc\n", t, b_halo); + exit(0);*/ + Plummer plummer(M_halo, b_halo); f[0] = y[3]; // vx -> x' f[1] = y[4]; // vy -> y' f[2] = y[5]; // vz -> z' @@ -80,8 +88,8 @@ public: } private: - Interp *interp_M_bulge; - Interp *interp_b_bulge; + Interp *interp_M_halo; + Interp *interp_b_halo; } galaxy("file.dat"); // Not very nice to have it as a global variable but GSL will have problem otherwise. @@ -99,7 +107,7 @@ int func(double t, const double y[], double f[], void *params) extern "C" int integrate(const double y0[], const double t_max, double y[]) { - double t = 0; + double t = 2.145; constexpr double h = 1./4096.; constexpr double epsabs = 1e-7; constexpr double epsrel = 0; diff --git a/plot.py b/plot.py index d2df0db..e498eb7 100644 --- a/plot.py +++ b/plot.py @@ -20,5 +20,5 @@ gsl_success = libmain.gsl_success() #t_array = linspace(plot_tmin, plot_tmax, plot_points) #r_array = empty_like(t_array) -res = integrate([10,0,0,0,200,0], 0.1) +res = integrate([10,0,0,0,200,0], 2.245) print(res, gsl_success) \ No newline at end of file