source: proiecte/NBody/NBody 2.0/OutputData.cpp @ 34

#include "NBody.hpp"

extern FILE *snapshotsFile, *diagnosticsFile;
extern void* me;                                                // Pointing to NBody class; 

/////////////////////////////////// NBody::put_snapshot() Func \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
//
// Writes a single snapshot on the output snapshot file.
//
/////////////////////////////////// NBody::put_snapshot() Func \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

void NBody::put_snapshot()
{
        if (snapshotsFile == NULL){
                std::cout << "Error: snapshot file not opened." << std::endl;
                ((NBody*)me)->cleanup();
                exit(1);
        }

        // cout.precision(16);                                          // Full double precision;
        fprintf(snapshotsFile, "%d\n", numParticles);           // Total particle number;
        fprintf(snapshotsFile, "%f\n", curr_time_step);         // Current time;

        for (int i = 0; i < numParticles ; i++){
                
                fprintf(snapshotsFile, "%f\t", pos[i*4+3]);             // Mass of particle i;
                        
        for (int k = 0; k < 3; k++)
                        fprintf(snapshotsFile, "%f\t", pos[i*4+k]);     // Position of particle i;
                        
        for (int k = 0; k < 3; k++)
                fprintf(snapshotsFile, "%f\t", vel[i*4+k]);     // Velocity of particle i;
        
                fprintf(snapshotsFile, "\n");
        }
        fprintf(snapshotsFile, "\n");
}


/////////////////////////////////// NBody::epot_particle() Func \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
//
// Computes the potential energy of a single particle;
//
/////////////////////////////////// NBody::epot_particle() Func \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

float NBody::epot_particle(int currPart)
{   
        float ePot = 0.0;
        for (int i=0; i<numParticles; i++){
                
                if (currPart != i){

                        int idx_curr = currPart*4;      // The current particle's position index in the "pos" vector;
                        int idx_other = i*4;            // The other particle's position index in the "pos" vector;

                        float r[3], r2= 0.0;
                        for (int k = 0; k < 3 ; k++){
                                r[k]= pos[idx_other+k] - pos[idx_curr+k];
                                r2+= r[k]*r[k];
                        }
                                
                        // To avoid Division by Zero:
                        if (r2 == 0) r2= 0.1e-50;

                        // ePot = -M1*M2/r;
                        ePot += -pos[idx_curr+3]*pos[idx_other+3]/sqrt(r2);
                }
        }

        return ePot;
}

/////////////////////////////////// NBody::write_diagnostics() Func \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
//
// Writes diagnostics on the diagnostics file:
// current time; number of integration steps so far;
// kinetic, potential, and total energy; absolute and
// relative energy errors since the start of the run.
//
/////////////////////////////////// NBody::write_diagnostics() Func \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

void NBody::write_diagnostics()
{

    float Ekin = 0;                                             // The Kinetic energy of the n-body system;
        for (int i = 0; i < numParticles ; i++){
                
                int idx = i*4;

                for (int k = 0; k < 3 ; k++){
                        // pos[idx+3] = mass;   
                        // Ekin = 1/2*m*(v^2)

            Ekin += 0.5 * pos[idx+3] * vel[idx+k] * vel[idx+k];
                }
        }

    float Epot= 0;                                                      // The Potential energy of the n-body system;
        for (int i=0; i<numParticles; i++)
        Epot += epot_particle(i);
    
        Epot *= 0.5;                                                    // Against double counting;

    float Etot =  Ekin+Epot;                            // The total energy of the actual n-body system;

        if (initFlag){                          // At first pass, the initial energy is equal to the
                Etot_init = Etot;                                       // total energy of the system. We will keep this initial
                                                                                        // energy for error detections;
        }                                                                               

        if (!initFlag){                                                 // At first pass we don't write any diagnostics;
                if (diagnosticsFile == NULL){
                        std::cout << "Error: diagnostics file not opened." << std::endl;
                        ((NBody*)me)->cleanup();
                        exit(1);
                }

                // cout.precision(16);                          // Full double precision;
                fprintf(diagnosticsFile, "At time t = %f, after %ld steps:\n", curr_time_step, taken_steps);
                fprintf(diagnosticsFile, " E_init = %f,\n E_tot = %f,\n E_kin = %f,\n E_pot = %f\n", Etot_init, Etot, Ekin, Epot);
                fprintf(diagnosticsFile, "                Absolute energy error: E_tot - E_init = %f\n", Etot - Etot_init);
                fprintf(diagnosticsFile, "                Relative energy error: (E_tot - E_init) / E_init = %f\n", (Etot - Etot_init) / Etot_init);
        }else{
                initFlag = false;
        }
}