Code:
#include <fstream>
#include <iostream>
#include <stdlib.h>
#include <cstring>
#include <math.h>
#include <memory>
using namespace std;
//use 3layer no force 11 m/s trajectory to debug
#define Num_step 10000000
#define Dump_freq 5000
#define Num_part 3303
#define Num_bin 25
#define Num_coor 3
#define Lower_first_bin 10.000
#define Width_bin 1.0
#define Lx 36.96
#define Ly 34.676
#define Lambdax 0.8
#define Lambday 0.8
class Allocate_memory
{
protected:
int num_dump,
num_part,
num_coor;
public:
Allocate_memory (int i_dim_x, int i_dim_y, int i_dim_z):
num_dump (i_dim_x), num_part (i_dim_y), num_coor (i_dim_z) {};
double *** allocate_memory_3d_array (double *** array_3d)
{
array_3d=new double ** [num_dump];
for (int i=0;i<num_dump;i++)
{
array_3d[i]=new double * [num_part];
for (int j=0;j<num_part;j++)
{
array_3d[i][j]=new double [num_coor];
}
}
return array_3d;
}
double ** allocate_memory_2d_array (double **array_2d)
{
array_2d=new double * [num_dump];
for (int i=0;i<num_dump;i++)
{
array_2d[i]=new double [num_part];
}
return array_2d;
}
void free_memory_3d_array (double *** array_3d)
{
for (int i=0;i<num_dump;i++)
{
for (int j=0;j<num_part;j++)
{
delete [] array_3d[i][j];
array_3d[i][j]=NULL;
}
delete [] array_3d[i];
array_3d[i]=NULL;
}
delete [] array_3d;
array_3d=NULL;
}
void free_memory_2d_array (double ** array_2d)
{
for (int i=0;i<num_dump;i++)
{
delete [] array_2d[i];
array_2d[i]=NULL;
}
delete [] array_2d;
array_2d=NULL;
}
//virtual ~Allocate_memory ();
};
class Structure_factor:public Allocate_memory
{
private:
double *** coor,
// coor[Num_step/Dump_freq][Num_part][Num_coor],
coor_bin[Num_bin][Num_part][Num_coor],
time_averaged_factor_bin_3D[Num_bin][int (Lx/Lambdax)][int (Ly/Lambday)],
** mass,
* lower_bin,
* up_bin,
* kx,
* ky,
* num_atom_bin, //store number of atoms for each bin
lower_first_bin,
width_bin,
lx,
ly,
lambdax,
lambday,
no_time_averaged_factor_bin[Num_step/Dump_freq][int(Num_bin)][int(Lx/Lambdax-1)*int(Ly/Lambday-1)],
time_averaged_factor_bin_2D[Num_bin][int(Lx/Lambdax-1)*int(Ly/Lambday-1)];
const string name,
output_3d_array,
output_2d_array;
string tem,
mol,
type,
x_c,
y_c,
z_c,
i_x,
i_y,
i_z,
s_mass;
fstream lammpstrj;
int dump_freq,
num_bin,
num_kx,
num_ky;
public:
Structure_factor (const string i_name, const string i_output_3d_array, const string i_output_2d_array,
int i_dim_x, int i_dim_y, int i_dim_z, int i_dump_freq, int i_num_bin, double i_lower_first_bin,
double i_width_bin,double i_lx,double i_ly,double i_lambdax,double i_lambday):
Allocate_memory (i_dim_x,i_dim_y,i_dim_z),
dump_freq (i_dump_freq),
name (i_name),
output_3d_array (i_output_3d_array),
output_2d_array (i_output_2d_array),
num_bin (i_num_bin),
lower_first_bin (i_lower_first_bin),
width_bin (i_width_bin),
lx (i_lx),
ly (i_ly),
lambdax (i_lambdax),
lambday (i_lambday)
{
coor=allocate_memory_3d_array (coor);
mass=allocate_memory_2d_array (mass);
num_kx=int (lx/lambdax);
kx=new double [(num_kx-1)];
num_ky=floor (ly/lambday);
ky=new double [(num_ky-1)];
lower_bin=new double [(num_bin)];
up_bin=new double [int(num_bin)];
num_atom_bin=new double [(num_bin)];
}
void open_file ()
{
lammpstrj.open (name.c_str(),ios::in);
if (!lammpstrj)
{
cout << "Can not open the flow_solid_Couetteall.lammpstrj" << endl;
exit(1);
}
else
{
cout << "open the flow_solid_Couetteall.lammpstrj successfully" << endl;
read_data ();
}
lammpstrj.close ();
}
void read_data ()
{
for (int i=0;i<num_dump;i++)
{
for (int j=0;j<9;j++)
{
getline (lammpstrj,tem);
}
for (int k=0;k<num_part;k++)
{
getline (lammpstrj,mol,' ');
getline (lammpstrj,type,' ');
getline (lammpstrj,x_c,' ');
coor[i][k][0]=strtod(x_c.c_str(),NULL);
getline (lammpstrj,y_c,' ');
coor[i][k][1]=strtod(y_c.c_str(),NULL);
getline (lammpstrj,z_c,' ');
coor[i][k][2]=strtod(z_c.c_str(),NULL);
getline (lammpstrj,i_x,' ');
getline (lammpstrj,i_y,' ');
getline (lammpstrj,i_z,' ');
getline (lammpstrj,s_mass);
mass[i][k]=strtod(s_mass.c_str(),NULL);
}
}
}
~Structure_factor ()
{
free_memory_3d_array (coor);
free_memory_2d_array (mass);
delete [] lower_bin;
lower_bin=NULL;
delete [] up_bin;
up_bin=NULL;
delete [] num_atom_bin;
num_atom_bin=NULL;
delete [] kx;
kx=NULL;
delete [] ky;
ky=NULL;
}
};
int main ()
{
auto_ptr<Structure_factor> structure_factor (
new Structure_factor ("flow_solid_Couetteall.lammpstrj","coordinate","mass",
Num_step/Dump_freq,Num_part,Num_coor,Dump_freq,Num_bin,
Lower_first_bin,Width_bin,Lx,Ly,Lambdax,Lambday));
structure_factor->open_file ();
}
[