ClumpInput.h

Go to the documentation of this file.

// This file is part of the MercuryDPM project (https://www.mercurydpm.org).
// Copyright (c), The MercuryDPM Developers Team. All rights reserved.
// License: BSD 3-Clause License; see the LICENSE file in the root directory.
 
#ifndef CLUMP_INPUT_H
#define CLUMP_INPUT_H
 
#include <algorithm>
#include <dirent.h>
#include <sys/types.h>
#include <iostream>
#include <fstream>
#include <vector>
#include <string>
#include <cmath>
#include <sstream>
#include "Math/Matrix.h"
#include <CMakeDefinitions.h>
 
// Useful containers
typedef std::vector<double> DoubleVector;
typedef std::vector<DoubleVector> Double2DVector;
typedef std::vector<std::string> StringVector;
 
// Helper function for alphabetical sorting of clump names
bool CompareFunction (std::string a, std::string b) {return a < b;}
 
// Helper function to generate a random double in range (0, MAX)
double RandomDouble(double Max)
{return Max*((double) rand() / (RAND_MAX));}
 
struct ClumpData
{
public:
    
    std::string path;   // Path to MClump working directory
    StringVector clumpNames;    // Array of names of Clumps that will be used
    DoubleVector mass;      //  clump mass
    Double2DVector  pebblesX;   //  Pebbles geometry (outer index goes over Clumps, inner - over pebbles)
    Double2DVector  pebblesY;
    Double2DVector  pebblesZ;
    Double2DVector  pebblesR;
    
    Double2DVector  toi;        // Clump tensor of inertia (I11, I12, I13, I21..I33)
    Double2DVector  pd;     // Clump principal directions v1, v2, v3
    
};
 
 
void LoadConf(ClumpData &a)
{
    // Path to MCLump tool
    a.path = getMercuryDPMSourceDir() + "/Tools/MClump/Clumps/";
 
    struct dirent *entry;
    DIR *dir = opendir(a.path.c_str());
    while ((entry = readdir(dir)) != NULL) {
        if (entry->d_type == DT_DIR){
            a.clumpNames.push_back(entry->d_name);
        }
    }
    closedir(dir);
 
    // remove "." and ".." from the list of dirs (position of those in the list is OS-sensitive, hence this code)
 
    std::sort(a.clumpNames.begin(), a.clumpNames.end(), CompareFunction);//sort the vector
    a.clumpNames.erase(a.clumpNames.begin());
    a.clumpNames.erase(a.clumpNames.begin());
 
    // Show the names of available Clumps
    for (int i = 0; i<a.clumpNames.size(); i++) std::cout << a.clumpNames[i] << std::endl;
 
}
 
void LoadPebbles(ClumpData &a)
{
    logger(INFO, "Loading clump pebbles...");
 
    a.pebblesX.resize(a.clumpNames.size());
    a.pebblesY.resize(a.clumpNames.size());
    a.pebblesZ.resize(a.clumpNames.size());
    a.pebblesR.resize(a.clumpNames.size());
 
    for (int i = 0; i < a.clumpNames.size(); i++ ){
        std::ifstream infile((a.path + a.clumpNames[i] + "/Clump/Clump.txt").c_str(), std::ios::in | std::ios::binary);
        std::string line{};
        while (std::getline(infile, line)) {
            std::istringstream iss(line);
            std::string substring{};
            StringVector val;
            while (std::getline(iss, substring, ',')) val.push_back(substring);
            a.pebblesX[i].push_back(std::stof(val[0]));
            a.pebblesY[i].push_back(std::stof(val[1]));
            a.pebblesZ[i].push_back(std::stof(val[2]));
            a.pebblesR[i].push_back(std::stof(val[3]));
        }
        infile.close();
    }
    logger(INFO, "\t OK\n");
}
 
void LoadMass(ClumpData &a)
{
    logger(INFO, "Loading clump masses...");
 
    a.mass.resize(a.clumpNames.size());
    
    for (int i = 0; i < a.clumpNames.size(); i++ ){
        std::ifstream infile((a.path + a.clumpNames[i] + "/Inertia/Mass.txt").c_str(), std::ios::in | std::ios::binary);
        std::string mass;
        infile >> mass;
        a.mass[i] = std::stof(mass);
        infile.close();
    }
    logger(INFO, "\t OK\n");
}
 
void LoadTOI(ClumpData &a)
{
    logger(INFO, "Loading clump TOI..");
    a.toi.resize(a.clumpNames.size());
 
    for (int i = 0; i < a.clumpNames.size(); i++ ){
        std::ifstream infile((a.path + a.clumpNames[i] + "/Inertia/TOI.txt").c_str(), std::ios::in | std::ios::binary);
        std::string line{};
        while (std::getline(infile, line)) {
            std::istringstream iss(line);
            std::string substring{};
            StringVector val;
            while (std::getline(iss, substring, ',')) val.push_back(substring);
            for (int k = 0; k < 3; k++) a.toi[i].push_back(std::stof(val[k]));
        }
        infile.close();
    }
    logger(INFO, "\t OK\n");
}
 
void LoadPD(ClumpData &a)
{
    logger(INFO, "Loading clump PD..");
    a.pd.resize(a.clumpNames.size());
 
    for (int i = 0; i < a.clumpNames.size(); i++ ){
        std::ifstream infile((a.path + a.clumpNames[i] + "/Inertia/PD.txt").c_str(), std::ios::in | std::ios::binary);
        std::string line{};
        while (std::getline(infile, line)) {
            std::istringstream iss(line);
            std::string substring{};
            StringVector val;
            while (std::getline(iss, substring, ',')) val.push_back(substring);
            for (int k = 0; k < 3; k++) a.pd[i].push_back(std::stof(val[k]));
        }
        infile.close();
    }
    logger(INFO, "\t OK\n");
}
 
 
void LoadClumps(ClumpData &data, bool VERBOSE = false)
{
 
    logger(INFO, "LOAD CLUMP DATA");
    LoadConf(data);
    LoadPebbles(data);
    LoadMass(data);
    LoadTOI(data);
    LoadPD(data);
    if (VERBOSE) {
        std::cout<<"LOADED CLUMPS"<<std::endl;
        for (int i = 0; i < data.pebblesX.size(); i++) {
            std::cout << data.clumpNames[i] << " mass:" << data.mass[i] << std::endl;
            std::cout << data.clumpNames[i] << " list of pebbles:" << std::endl;
            for (int j = 0; j < data.pebblesX[i].size(); j++) {
                std::cout << "Pebble " << j << ": (" << data.pebblesX[i][j] << "," << data.pebblesY[i][j] << ","
                          << data.pebblesZ[i][j] << ")," << data.pebblesR[i][j] << std::endl;
            }
 
            std::cout << data.clumpNames[i] << " TOI:" << std::endl;
            std::cout << data.toi[i][0] << "," << data.toi[i][1] << "," << data.toi[i][2] << std::endl;
            std::cout << data.toi[i][3] << "," << data.toi[i][4] << "," << data.toi[i][5] << std::endl;
            std::cout << data.toi[i][6] << "," << data.toi[i][7] << "," << data.toi[i][8] << std::endl;
 
            std::cout << data.clumpNames[i] << " Principal directions:" << std::endl;
            std::cout << data.pd[i][0] << "," << data.pd[i][1] << "," << data.pd[i][2] << std::endl;
            std::cout << data.pd[i][3] << "," << data.pd[i][4] << "," << data.pd[i][5] << std::endl;
            std::cout << data.pd[i][6] << "," << data.pd[i][7] << "," << data.pd[i][8] << std::endl;
        }
    }
 
}
 
ClumpData RotateClump(ClumpData data, int clump_index, DoubleVector new_pd)
{
    ClumpData new_data = data;
    new_data.pd[clump_index] = new_pd;
    return new_data;
}
 
 
 
 
 
 
DoubleVector UniformRandomPDs(){
 
    Vec3D n1, n2, n3, ref;
 
    // basis vector n1
    double r1 = RandomDouble(2) - 1.0;
    double r2 = RandomDouble(1);
 
    double theta = acos(r1); // Note that for isotropy of n1 theta is NOT uniform!
    double phi = 2 * M_PI * r2;
 
    n1 = Vec3D(sin(theta) * cos(phi), sin(theta) * sin(phi), cos(theta));
 
    // basis vector n2
    r1 = RandomDouble(1);
    r2 = RandomDouble(1);
 
    theta = acos(2 * r1 - 1); // Note that for isotropy of n1 theta is NOT uniform!
    phi = 2 * M_PI * r2;
 
    ref = Vec3D(sin(theta) * cos(phi), sin(theta) * sin(phi), cos(theta));
    n2 = Vec3D::cross(ref, n1); n2.normalise();
    n3 = Vec3D::cross(n1, n2); n3.normalise();
    return DoubleVector{n1.X, n1.Y, n1.Z, n2.X, n2.Y, n2.Z, n3.X, n3.Y, n3.Z};
 
}
 
 
#endif // CLUMP_INPUT_H