34 return makeUnexpected<Storage>(result.
error());
41 std::cout <<
"Processing pkdgrav file ... " << std::endl;
44 std::cout <<
"Invalid file: " << storage.
error() << std::endl;
52 logRadiiSfd.
write(p.value,
" ", p.count);
58 std::cout <<
"Processing pkdgrav file ... " << std::endl;
61 std::cout <<
"Invalid file: " << storage.
error() << std::endl;
77 omega.
begin(), omega.
end(), [](
Vector& v1,
Vector& v2) { return getLength(v1) > getLength(v2); });
87 std::cout <<
"Processing pkdgrav file ... " << std::endl;
90 std::cout <<
"Invalid file: " << storage.
error() << std::endl;
96 std::cout <<
"Moon count = " << moonCnt << std::endl;
101 std::cout <<
"Processing SPH file ... " << std::endl;
105 Outcome outcome = input->load(filePath, storage, stats);
107 std::cout <<
"Cannot load particle data, " << outcome.
error() << std::endl;
116 logSfd.
write(p.value,
" ", p.count);
122 const Path& omegaPath,
123 const Path& omegaDPath,
124 const Path& omegaDirPath) {
125 std::cout <<
"Processing SPH file ... " << std::endl;
129 Outcome outcome = input.
load(filePath, storage, stats);
131 std::cout <<
"Cannot load particle data, " << outcome.
error() << std::endl;
142 const Float m_total = std::accumulate(m.begin(), m.end(), 0._f);
154 logOmega.
write(p.value,
" ", p.count);
163 logOmegaDir.
write(p.value,
" ", p.count);
170 for (
Size i = 0; i < m.size(); ++i) {
176 for (
Size i = 0; i < m.size(); ++i) {
177 if (m[i] > 3._f * params.
massCutoff * m_total) {
179 logOmegaD.
write(2._f * h[i],
" ",
getLength(w[i]) * 3600._f * 24._f / (2._f *
PI));
188 std::cout <<
"Processing SPH file ... " << std::endl;
192 Outcome outcome = input->load(filePath, storage, stats);
194 std::cout <<
"Cannot load particle data, " << outcome.
error() << std::endl;
206 for (
Size i : idxs) {
214 for (
Size i = 0; i < r.size(); ++i) {
226 logSfd.
write(p.value,
" ", p.count);
233 std::cout <<
"Processing SPH file ... " << std::endl;
237 Outcome outcome = input.
load(filePath, storage, stats);
239 std::cout <<
"Cannot load particle data, " << outcome.
error() << std::endl;
251 for (
Size i : idxs) {
259 for (
Size i = 0; i < r.size(); ++i) {
266 for (
Size i = 0; i < v.size(); ++i) {
280 logSfd.
write(p.value,
" ", p.count);
296 ifs >> number >> dummy;
300 for (
int i = 0; i < 6; ++i) {
306 for (
int i = 0; i < 5; ++i) {
312 for (
int i = 0; i < 10; ++i) {
317 fromString<Size>(number),
319 fromString<Float>(
radius),
320 fromString<Float>(period),
335 bool firstLine =
true;
337 while (std::getline(ifs, line)) {
338 if (line.empty() || line[0] ==
'#') {
348 std::stringstream ss(line);
369 if (ast1.number && ast2.number && ast1.number.value() == ast2.number.value()) {
379 if (iter == catalog.
end()) {
382 if (iter->period && iter->radius) {
391 const Path& outputPath,
395 std::ifstream ifs(familyData.
native());
402 if (
auto harAst = findInHarris(famAst, catalog)) {
403 found.
push(harAst.value());
404 rangePeriod.
extend(harAst->period.value());
405 rangeR.
extend(harAst->radius.value());
408 if (found.
size() < 10) {
414 std::ofstream ofs(outputPath.
native());
422 auto periodToOmega = [](
const Float P) {
return 1._f / (P / 24._f); };
426 const Float omega = periodToOmega(ast.period.value());
427 if (&ast != &*largestRemnant) {
428 std::string printedName = ast.name;
430 printedName =
"(" + std::to_string(ast.number.value()) +
") " + printedName;
432 ofs << ast.radius.value() <<
" " << omega <<
" " << printedName << std::endl;
437 rangeOmega.
extend(periodToOmega(rangePeriod.
lower()));
438 rangeOmega.
extend(periodToOmega(rangePeriod.
upper()));
444 if (points.
size() > 36) {
447 std::ofstream histofs(histPath.
native());
456 histofs << p.value <<
" " << p.count << std::endl;
460 outPoints = std::move(points);
464 static std::string elliptize(
const std::string& s,
const Size maxSize) {
466 if (s.size() < maxSize) {
469 const Size cutSize = (maxSize - 3) / 2;
470 return s.substr(0, cutSize) +
"..." + s.substr(s.size() - cutSize);
474 Path harrisPath(
"/home/pavel/projects/astro/asteroids/grant3/harris.out");
475 std::ifstream ifs(harrisPath.
native());
479 Size below3 = 0, below7 = 0, below12 = 0, total = 0;
480 Float avgPeriod = 0._f;
482 for (
auto& h : harris) {
487 if (h.period.value() < 3) {
490 if (h.period.value() < 7) {
493 if (h.period.value() < 12) {
497 if (h.radius.valueOr(0._f) > 100._f) {
498 avgPeriod += h.period.value();
503 std::cout <<
"Below 3h: " << (100.f * below3) / total <<
"%" << std::endl;
504 std::cout <<
"Below 7h: " << (100.f * below7) / total <<
"%" << std::endl;
505 std::cout <<
"Below 12h: " << (100.f * below12) / total <<
"%" << std::endl;
508 std::cout <<
"Average period for R>100km: " << avgPeriod / count << std::endl;
515 Path(
"/home/pavel/projects/astro/asteroids/families.txt"),
Path(
"LR_D_omega.txt"), harris, points);
517 const Path parentPath(
"/home/pavel/projects/astro/asteroids/families");
519 std::mutex printMutex;
520 std::mutex plotMutex;
522 std::ofstream kss(
"KS.txt");
523 kss <<
"# name D_ks probability" << std::endl;
525 std::ofstream alls(
"D_omega_all.txt");
530 std::unique_lock<std::mutex> lock(printMutex);
531 std::cout << paths[index].native() << std::endl;
533 const Path name = paths[index].fileName();
534 const Path targetPath =
Path(
"D_omega") /
Path(name.native() +
".txt");
535 const Path ksPath =
Path(
"KS") /
Path(name.native() +
".txt");
539 std::unique_lock<std::mutex> lock(plotMutex);
544 kss << std::left << std::setw(35) << elliptize(
Path(name).removeExtension().native(), 30)
545 << std::setw(15) << ks->D << std::setw(15) << ks->prob << std::endl;
548 alls << p.x <<
" " << p.y <<
" " << index << std::endl;
553 Process gnuplot(
Path(
"/bin/gnuplot"), {
"doplot.plt" });
558 if (points.
size() > 25) {
561 Process gnuplot2(
Path(
"/bin/gnuplot"), {
"dohistogram.plt" });
565 Path(
"hist.png"), uniquePaths.
getPath(
Path(histPath).replaceExtension(
"png")));
570 Process gnuplot(
Path(
"/bin/gnuplot"), {
"doplot_all.plt" });
580 const Float x = rng() * a * 10._f;
581 const Float y = rng() / a;
583 if (maxwellBoltzmann(x, a) > y) {
627 std::ifstream ifs(filePath.
native());
631 while (std::getline(ifs, line)) {
632 const Float d = std::stof(line);
633 rs.
push(d / 2 * 1000);
636 std::ofstream yarko(
"yarko.in");
637 yarko << rs.
size() << std::endl;
639 yarko << r <<
" 2860.0 1500.0 0.0010 680.0 0.10 0.90" << std::endl;
642 std::ofstream spin(
"spin.in");
643 spin << rs.
size() <<
"\n-1\n1\n";
645 for (
Size i = 0; i < rs.
size(); ++i) {
646 const Float phi = rng() * 2._f *
PI;
647 const Float cosTheta = rng() * 2._f - 1._f;
653 std::ofstream yorp(
"yorp.in");
654 for (
Size i = 1; i <= rs.
size(); ++i) {
655 yorp << i <<
" " <<
clamp(
int(rng() * 200), 0, 199) << std::endl;
663 Outcome res1 = input.
load(lastDumpPath, lastDump, stats);
664 Outcome res2 = input.
load(firstDumpPath, firstDump, stats);
665 if (!res1 || !res2) {
666 throw IoError((res1 || res2).error());
680 output.
dump(firstDump, stats);
687 Outcome res = input.
load(inputPath, storage, stats);
698 const Size componentCnt =
700 std::cout <<
"Component cnt = " << componentCnt << std::endl;
703 for (
Size i = 0; i < components.
size(); ++i) {
704 if (components[i] != 0) {
717 output.
dump(storage, stats);
723 for (
Size i = 0; i < m.size(); ++i) {
724 volume += m[i] / rho[i];
727 std::cout <<
"eq. diameter = " <<
Sph::cbrt(3._f * volume / (4._f *
PI)) * 2._f / 1000._f <<
"km"
732 std::cout <<
"period = " << 2._f *
PI / omega / 3600._f <<
"h" << std::endl;
746 SPH_ASSERT(a > 0._f && b > 0._f && c > 0._f, a, b, c);
747 std::cout <<
"a/b = " << a / b << std::endl;
748 std::cout <<
"b/c = " << b / c << std::endl;
752 std::cout <<
"Expected usage: post mode [parameters]" << std::endl
753 <<
" where 'mode' is one of:" << std::endl
754 <<
" - pkdgravToSfd - computes the cumulative SFD from pkdgrav output file" << std::endl
755 <<
" - pkdgravToOmega - computes the spin rate distribution from pkdgrav output file"
757 <<
" - pkdgravToMoons - finds satellites of the largest remnant (fragment) from pkdgrav "
760 <<
"- ssfToSfd - computes the cumulative SFD from SPH output file" << std::endl
761 <<
"- ssfToVelocity - computes the velocity distribution from SPH output file" << std::endl
762 <<
"- harris - TODO" << std::endl
763 <<
"- stats - prints ejected mass and the period of the largest remnant" << std::endl
764 <<
"- swift - makes yarko.in, yorp.in and spin.in input file for swift" << std::endl;
767 int main(
int argc,
char** argv) {
774 std::string mode(argv[1]);
775 if (mode ==
"pkdgravToSfd") {
777 std::cout <<
"Expected parameters: post pkdgravToSfd ss.50000.bt sfd.txt";
781 }
else if (mode ==
"pkdgravToOmega") {
783 std::cout <<
"Expected parameters: post pkdgravToOmega ss.50000.bt omega.txt";
787 }
else if (mode ==
"pkdgravToMoons") {
789 std::cout <<
"Expected parameters: post pkdgravToMoons ss.50000.bt 0.1";
792 const Float limit = std::atof(argv[3]);
794 }
else if (mode ==
"ssfToSfd") {
796 std::cout <<
"Expected parameters: post ssfToSfd [--components] output.ssf sfd.txt";
799 if (std::string(argv[2]) ==
"--components") {
804 }
else if (mode ==
"ssfToVelocity") {
806 }
else if (mode ==
"ssfToOmega") {
808 std::cout <<
"Expected parameters: post ssfToOmega output.ssf omega.txt omega_D.txt "
813 }
else if (mode ==
"ssfToVelDir") {
815 std::cout <<
"Expected parameters: post ssfToVelDir output.ssf veldir.txt" << std::endl;
819 }
else if (mode ==
"harris") {
821 }
else if (mode ==
"swift") {
823 std::cout <<
"Expected parameters: post maketp D.dat";
827 }
else if (mode ==
"origComponents") {
829 std::cout <<
"Expected parameters: post origComponents lastDump.ssf firstDump.ssf "
834 }
else if (mode ==
"extractLr") {
836 std::cout <<
"Expected parameters: post extractLr input.ssf lr.ssf" << std::endl;
844 }
catch (
const std::exception& e) {
845 std::cout <<
"ERROR: " << e.what() << std::endl;
Various function for interpretation of the results of a simulation.
INLINE bool isReal(const AntisymmetricTensor &t)
#define SPH_ASSERT(x,...)
const EmptyFlags EMPTY_FLAGS
uint32_t Size
Integral type used to index arrays (by default).
double Float
Precision used withing the code. Use Float instead of float or double where precision is important.
Vector moveToCenterOfMassSystem(ArrayView< const Float > m, ArrayView< Vector > r)
Modifies particle positions so that their center of mass lies at the origin.
Generating initial conditions of SPH particles.
Logging routines of the run.
constexpr INLINE T clamp(const T &f, const T &f1, const T &f2)
constexpr INLINE T sqr(const T &f) noexcept
Return a squared value.
INLINE Float cbrt(const Float f)
Returns a cubed root of a value.
INLINE T sqrt(const T f)
Return a squared root of a value.
constexpr INLINE Float pow< 3 >(const Float v)
INLINE T atan2(const T y, const T x)
constexpr Float RAD_TO_DEG
constexpr Float PI
Mathematical constants.
const NothingType NOTHING
Saving and loading particle data.
Simple thread pool with fixed number of threads.
Creating and managing processes.
@ FLAG
ID of original body, used to implement discontinuities between bodies in SPH.
@ POSITION
Positions (velocities, accelerations) of particles, always a vector quantity,.
@ DENSITY
Density, always a scalar quantity.
@ MASS
Paricles masses, always a scalar quantity.
Holder of quantity values and their temporal derivatives.
@ ZERO
Quantity without derivatives, or "zero order" of quantity.
INLINE void parallelFor(IScheduler &scheduler, const Size from, const Size to, TFunctor &&functor)
Executes a functor concurrently from all available threads.
StaticArray< T0 &, sizeof...(TArgs)+1 > tie(T0 &t0, TArgs &... rest)
Creates a static array from a list of l-value references.
Statistics gathered and periodically displayed during the run.
Container for storing particle quantities and materials.
Eigen eigenDecomposition(const SymmetricTensor &t)
Computes eigenvectors and corresponding eigenvalues of symmetric matrix.
INLINE Float getLength(const Vector &v)
Returns the length of the vector. Enabled only for vectors of floating-point precision.
INLINE Float getSqrLength(const Vector &v)
INLINE Vector sphericalToCartesian(const Float r, const Float theta, const Float phi)
Constructs a vector from spherical coordinates.
int main(int argc, char *argv[])
Object providing safe access to continuous memory of data.
INLINE Iterator< StorageType > begin()
INLINE Iterator< StorageType > end()
INLINE Iterator< StorageType > end() noexcept
INLINE void push(U &&u)
Adds new element to the end of the array, resizing the array if necessary.
INLINE T & back() noexcept
INLINE TCounter size() const noexcept
INLINE Iterator< StorageType > begin() noexcept
Array clone() const
Performs a deep copy of all elements of the array.
Wrapper of pointer that deletes the resource from destructor.
INLINE const TError & error() const
Returns the error message.
Output saving data to binary data without loss of precision.
virtual Expected< Path > dump(const Storage &storage, const Statistics &stats) override
Saves data from particle storage into the file.
Wrapper of type that either contains a value of given type, or an error message.
const Error & error() const
Returns the error message.
Type & value()
Returns the reference to expected value.
void write(TArgs &&... args)
Creates and logs a message by concatenating arguments.
Object representing a 1D interval of real numbers.
INLINE Float lower() const
Returns lower bound of the interval.
INLINE void extend(const Float &value)
Extends the interval to contain given value.
INLINE Float upper() const
Returns upper bound of the interval.
Exception thrown when file cannot be read, it has invalid format, etc.
Simple (forward) iterator over continuous array of objects of type T.
INLINE Type & value()
Returns the reference to the stored value.
Object representing a path on a filesystem.
std::string native() const
Returns the native version of the path.
Path fileName() const
Returns the filename of the path.
Path parentPath() const
Returns the parent directory. If the path is empty or root, return empty path.
Holds a handle to a created process.
void wait()
Blocks the calling thread until the managed process exits. The function may block indefinitely.
Object holding various statistics about current run.
Container storing all quantities used within the simulations.
Quantity & insert(const QuantityId key, const OrderEnum order, const TValue &defaultValue)
Creates a quantity in the storage, given its key, value type and order.
StaticArray< Array< TValue > &, 3 > getAll(const QuantityId key)
Retrieves quantity buffers from the storage, given its key and value type.
Array< TValue > & getDt(const QuantityId key)
Retrieves a quantity derivative from the storage, given its key and value type.
Size getParticleCnt() const
Returns the number of particles.
@ INDICES_SORTED
Use if the given array is already sorted (optimization)
void remove(ArrayView< const Size > idxs, const Flags< IndicesFlag > flags=EMPTY_FLAGS)
Removes specified particles from the storage.
bool has(const QuantityId key) const
Checks if the storage contains quantity with given key.
Array< TValue > & getValue(const QuantityId key)
Retrieves a quantity values from the storage, given its key and value type.
Symmetric tensor of 2nd order.
Thread pool capable of executing tasks concurrently.
static SharedPtr< ThreadPool > getGlobalInstance()
Returns the global instance of the thread pool.
Object generating unique paths.
Path getPath(const Path &expected)
Generates a unique path, based on given input path.
Creating code components based on values from settings.
AutoPtr< IInput > getInput(const Path &path)
bool pathExists(const Path &path)
Checks if a file or directory exists (or more precisely, if a file or directory is accessible).
Array< Path > getFilesInDirectory(const Path &directory)
Alternatitve to iterateDirectory, returning all files in directory in an array.
Size fileSize(const Path &path)
Returns the size of a file.
Outcome createDirectory(const Path &path, const Flags< CreateDirectoryFlag > flags=CreateDirectoryFlag::ALLOW_EXISTING)
Creates a directory with given path. Creates all parent directories as well.
Outcome copyFile(const Path &from, const Path &to)
Copies a file on given path to a different path.
Array< MoonEnum > findMoons(const Storage &storage, const Float radius=1._f, const Float limit=0._f)
Find a potential satellites of the largest body.
Array< HistPoint > getDifferentialHistogram(ArrayView< const Float > values, const HistogramParams ¶ms)
Computes the differential histogram from given values.
HistogramSource
Source data used to construct the histogram.
@ COMPONENTS
Equivalent radii of connected chunks of particles (SPH framework)
@ PARTICLES
Radii of individual particles, considering particles as spheres (N-body framework)
KsFunction getUniformKsFunction(Interval rangeX, Interval rangeY)
Array< Size > findLargestComponent(const Storage &storage, const Float particleRadius, const Flags< ComponentFlag > flags)
Returns the indices of particles belonging to the largest remnant.
Vector getAngularFrequency(ArrayView< const Float > m, ArrayView< const Vector > r, ArrayView< const Vector > v, const Vector &r0, const Vector &v0, ArrayView< const Size > idxs=nullptr)
Computes the immediate vector of angular frequency of a rigid body.
KsResult kolmogorovSmirnovTest(ArrayView< const Float > data, const Function< Float(Float)> &expectedCdf)
One-dimensional Kolmogorov-Smirnov test with given CDF of expected probability distribution.
Size findComponents(const Storage &storage, const Float particleRadius, const Flags< ComponentFlag > flags, Array< Size > &indices)
Finds and marks connected components (a.k.a. separated bodies) in the array of vertices.
SymmetricTensor getInertiaTensor(ArrayView< const Float > m, ArrayView< const Vector > r, const Vector &r0, ArrayView< const Size > idxs=nullptr)
Computes the total inertia tensor of particles with respect to given center.
@ MOON
Body is on elliptical trajectory, it is a potential sattelite.
Array< HistPoint > getCumulativeHistogram(const Storage &storage, const ExtHistogramId id, const HistogramSource source, const HistogramParams ¶ms)
Computes cumulative (integral) histogram of particles in the storage.
@ VELOCITIES
Particle velocities.
@ ROTATIONAL_AXIS
Distribution of axis directions, from -pi to pi.
@ ROTATIONAL_FREQUENCY
Rotational frequency in revs/day.
void origComponents(const Path &lastDumpPath, const Path &firstDumpPath, const Path &colorizedDumpPath)
int pkdgravToOmega(const Path &filePath, const Path &omegaPath)
int ssfToSfd(const Post::HistogramSource source, const Path &filePath, const Path &sfdPath)
void ssfToVelDir(const Path &filePath, const Path &outPath)
int ssfToVelocity(const Path &filePath, const Path &outPath)
int pkdgravToSfd(const Path &filePath, const Path &sfdPath)
void extractLr(const Path &inputPath, const Path &outputPath)
int ssfToOmega(const Path &filePath, const Path &omegaPath, const Path &omegaDPath, const Path &omegaDirPath)
int pkdgravToMoons(const Path &filePath, const Float limit)
void makeSwift(const Path &filePath)
Vector values
Eigenvalues.
Optional< std::string > name
Parameters of the histogram.
Function< bool(Size index)> validator
Function used for inclusiong/exclusion of values in the histogram.
Float massCutoff
Cutoff value (lower bound) of particle mass for inclusion in the histogram.
Interval range
Range of values from which the histogram is constructed.
Size binCnt
Number of histogram bins.