Program Bubble Packing
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <math.h>
#include <vector>
#include <cstring>
#include <fstream>
#define pi 3.14159265358979323846264338327950288419716939937510
#define ntmax 50000
int every = 25;
using namespace std;
typedef struct{
double x;
double y;
} Vec_2D;
typedef struct{
Vec_2D pos;
Vec_2D force;
Vec_2D v;
double diameter;
int state;
}Buble;
typedef std::vector<Buble> Vbuble;
enum{fixbuble, freebuble};
double const c = 0.03; //edit
double const m = 1.7;
double const eps = 0.0034; //edit
double const diam_free_buble = 0.0126; //edit
double const diam_fix_buble = 0.0064; //edit
double const tol_dist = 1.16*diam_free_buble;
double const k0 = 1./(m)*pow(c/(2*0.7),2);
int Nbuble;
int Nfixbuble;
int Nfreebuble;
int nt = 1;
//-------------------------------------------------------------------------------------------------------
Buble* Initial_Probe_Position(int* n) {
Buble buble;
Vbuble bub_tmp;
int k = 0;
string line;
ifstream inputfile;
inputfile.open("initial_buble_position.txt");
char datanum[30];
while (!inputfile.eof()) {
getline(inputfile,line);
if (line.length() != 0 && line[0] != '#') {
strcpy(datanum,line.c_str());
sscanf(datanum, "%lf %lf %lf %lf %d", &buble.pos.x, &buble.pos.y, &buble.v.x, &buble.v.y, &buble.state);
}
if (buble.state == 0) {
buble.diameter = diam_fix_buble;
k++;
} else buble.diameter = diam_free_buble;
bub_tmp.push_back(buble);
}
inputfile.close();
Nbuble = bub_tmp.size();
Nfixbuble = k;
Nfreebuble = Nbuble - Nfixbuble;
(*n) = bub_tmp.size();
printf("#buble = %d\n",(*n));
Buble* buble_p = new Buble[(*n)];
for(int i = 0; i < (*n); i++){
buble_p[i] = bub_tmp[i];
}
bub_tmp.clear();
return buble_p;
}
double Buble_Distance(Buble b1, Buble b2){
return sqrt(pow(b1.pos.x-b2.pos.x,2) + pow(b1.pos.y-b2.pos.y,2));
}
Vec_2D Vec_Inbuble_Direction(Buble b1, Buble b2) {
double r = Buble_Distance(b1,b2);
Vec_2D ibf;
ibf.x = (b2.pos.x - b1.pos.x)/r;
ibf.y = (b2.pos.y - b1.pos.y)/r;
return ibf;
}
double Angle_of_Force(Buble b1, Buble b2) {
Vec_2D ibf = Vec_Inbuble_Direction(b1,b2);
double theta;
theta = atan2(ibf.y, ibf.x);
theta = (theta < 0 ? 2*pi + theta : theta);
return theta;
}
double Interbubble_Force(Buble b1, Buble b2) {
double r = Buble_Distance(b1,b2);
double r0 = 0.5*b1.diameter + 0.5*b2.diameter;
double w = r/r0;
double f = 0.;
if (w >= 0 && w <= 1.5)
f = 1.25*k0/r0*pow(w,3) - 2.375*k0/r0*pow(w,2) + 1.125*k0/r0;
else f = 0.;
return f;
}
void compute_f(const vector<double> &x_vector, vector<double> &f_vector, double f_ext){
//double x1 = x_vector[0];
double x2 = x_vector[1];
f_vector[0] = x2;
f_vector[1] = -c/m*x2 + 1./m*f_ext;
}
//4th order Runge-Kutta
Vec_2D Runge_Kutta_4(vector<double> &x_vector, double f_ext, double tf){
Vec_2D vecx;
double t0 = 0.;
double dt = 1.e-4;
//double t = t0;
int m = (tf-t0)/dt;
//------------------------------------------------------
int n = x_vector.size();
vector<double> x(n), f(n), k1(n), k2(n), k3(n), k4(n);
for (int i=1; i <= m+1; i++){
compute_f(x_vector,f,f_ext);
for(int i=0; i < n; i++)
k1[i] = dt*f[i];
for(int i=0; i < n; i++)
x[i] = x_vector[i] + 0.5*k1[i];
compute_f(x,f,f_ext);
for(int i=0; i < n; i++)
k2[i] = dt*f[i];
for(int i=0; i < n; i++)
x[i] = x_vector[i] + 0.5*k2[i];
compute_f(x,f,f_ext);
for(int i=0; i < n; i++)
k3[i] = dt*f[i];
for(int i=0; i < n; i++)
x[i] = x_vector[i] + k3[i];
compute_f(x,f,f_ext);
for(int i=0; i < n; i++)
k4[i] = dt*f[i];
for(int i=0; i < n; i++)
x_vector[i] += (k1[i] + 2*k2[i] + 2*k3[i] + k4[i])/6.;
}
vecx.x = x[0];
vecx.y = x[1];
return vecx;
}
void Sum_Buble_Force(Buble *buble, int Nb) {
double r0;
for (int i = 0; i < Nb; i++) {
if (buble[i].state == freebuble) {
//printf("i = %d free %d\n", i, buble[i].state);
double sumfx = 0.;
double sumfy = 0.;
double theta = 0.;
double r;
for (int j = 0; j < Nb; j++) {
if (i != j) {
r = Buble_Distance(buble[i], buble[j]);
r0 = 0.5*buble[i].diameter + 0.5*buble[j].diameter;
theta = Angle_of_Force(buble[i],buble[j]);
if (r/r0 > 1.) {
sumfx += fabs(Interbubble_Force(buble[i],buble[j]))*(cos(theta));
sumfy += fabs(Interbubble_Force(buble[i],buble[j]))*(sin(theta));
} else {
sumfx += -1*fabs(Interbubble_Force(buble[i],buble[j]))*(cos(theta));
sumfy += -1*fabs(Interbubble_Force(buble[i],buble[j]))*(sin(theta));
}
}
}
buble[i].force.x = sumfx;
buble[i].force.y = sumfy;
}
}
}
void New_Buble_Position(Buble *buble, double t, int Nb) {
vector<double> x(2);
vector<double> y(2);
Vec_2D fext;
Vec_2D vecx, vecy;
for (int i=0; i < Nb; i++) {
if (buble[i].state == freebuble) {
x[0] = buble[i].pos.x;
x[1] = buble[i].v.x;
fext.x = buble[i].force.x;
vecx = Runge_Kutta_4(x, fext.x, t);
buble[i].pos.x = vecx.x;
buble[i].v.x = vecx.y;
y[0] = buble[i].pos.y;
y[1] = buble[i].v.y;
fext.y = buble[i].force.y;
vecy = Runge_Kutta_4(y, fext.y, t);
}
}
}
int Stop_Searching(Buble const *buble, double eps, int Nb) {
int Nf = Nfreebuble;
int stop[Nf*(Nf-1)/2] = {0};
int keluar = 0;
int hit = 0;
int k = 0;
double f = 0.;
double d;
for (int i = Nfixbuble; i < Nb-1; i++) {
for (int j = i+1; j < Nb; j++) {
f = Interbubble_Force(buble[i], buble[j]);
d = Buble_Distance(buble[i],buble[j]);
if (((fabs(f) <= eps) && (fabs(f) >= 0.)) || ((d > 1.00*diam_free_buble) && (d < tol_dist))) {
stop[k] = 1;
printf("f%d-%d = %8.6lf stabil\n", i-Nfixbuble+1, j-Nfixbuble+1, fabs(f));
k++;
} else printf("f%d-%d = %8.6lf tidak stabil\n", i-Nfixbuble+1, j-Nfixbuble+1, fabs(f));
}
}
for (int i = 0; i < k; i++) {
if (stop[i] == 1) hit++;
//printf("%d ", stop[i]);
}
if (hit == Nf*(Nf-1)/2) keluar = 1;
printf("buble stabil = %d || target buble stabil = %d\n", hit, Nf*(Nf-1)/2);
return keluar;
}
void read_data_points(FILE* fp, double* x, double* y, int n){
for(int i = 0; i <= n-1; i++){
fscanf(fp, "%lf %lf", &(x[i]), &(y[i]));
}
}
void set_gnuplot_condition(FILE* gp)
{
fprintf(gp, "set xrange [%g:%g]\n", 0.33, 0.38);
fprintf(gp, "set yrange [%g:%g]\n", 0.17, 0.23);
fprintf(gp, "set xtics 0.01\n");
fprintf(gp, "unset mouse\n");
fprintf(gp, "set term x11\n");
fprintf(gp, "set size ratio -1\n");
}
void plot_buble(FILE* gp, Buble const* buble, double* x1, double* y1, int n1, int n, int Nb)
{
fprintf(gp, "set title \"iteration = %d\"\n", n);
fprintf(gp, "plot \"-\" with circles lc rgb \"blue\" fs transparent solid 0.1 border notitle, \"-\" w p lt 2 notitle\n");
for (int i=0; i < Nb; i++) {
//if (buble[i].state == freebuble)
fprintf(gp, "%f %f %f\n", buble[i].pos.x, buble[i].pos.y, 0.5*buble[i].diameter);
}
fprintf(gp, "e\n");
for (int i = 0; i <= n1-1; i++){
fprintf(gp, "%f %f\n", x1[i], y1[i]);
}
fflush(gp);
}
void Save_Buble(Buble const* buble, int Nb) {
FILE *file;
file=fopen("buble_position_optimal.txt","w");
fprintf(file, "#iteration = %d\n", nt-1);
fprintf(file, "#total buble = %d || fixbuble = %d || freebuble = %d\n", Nb, Nfixbuble, Nfreebuble);
fprintf(file, "#diameter fix buble = %10.6lf || diameter free buble = %10.6lf\n", diam_fix_buble, diam_free_buble);
fprintf(file,"#m = %4.2lf || c = %4.2lf || eps = %10.6lf || k0 = %10.6lf\n", m, c, eps, k0);
for (int i=0; i < Nb; i++) {
if (buble[i].state == freebuble) fprintf(file, "%8.6lf %8.6lf %8.6lf %8.6lf\n", buble[i].pos.x, buble[i].pos.y, 0.5*buble[i].diameter, 0.003);
}
fclose(file);
}
void save2file (Buble const* buble, int n, int Nb){
FILE *fp;
char filename[32];
sprintf(filename, "output%.4d.dat", n/every);
fp=fopen(filename,"w");
fprintf(fp, "#iteration = %d\n", n);
for (int i=0; i < Nb; i++) {
if (buble[i].state == freebuble)
fprintf(fp, "%f %f %f\n", buble[i].pos.x, buble[i].pos.y, 0.5*buble[i].diameter);
}
fclose(fp);
}
void save2gnuplot (int nt){
FILE *plot;
plot=fopen("buble_packing_animation.gp","w");
for (int n = 1; n <= nt/every; n++){
fprintf(plot,"reset\n");
fprintf(plot,"set terminal pngcairo background \"#ffffff\" enhanced fontscale 1.0 size 640, 480 dashlength 2\n");
fprintf(plot,"set output \"figure%.4d.png\"\n", n);
fprintf(plot,"set size ratio -1\n");
fprintf(plot,"set style line 1 linecolor rgb \"orange\" linewidth 2.000 dashtype 1 pointtype 2 pointsize default pointinterval 0\n");
fprintf(plot,"set style line 2 linecolor rgb \"blue\" linewidth 2.000 dashtype 1 pointtype 2 pointsize default pointinterval 0\n");
fprintf(plot, "set xrange [%g:%g]\n", 0.33, 0.38);
fprintf(plot, "set yrange [%g:%g]\n", 0.17, 0.23);
fprintf(plot, "set xtics 0.01\n");
fprintf(plot, "unset mouse\n");
fprintf(plot,"set xlabel \"x (m)\"\n");
fprintf(plot,"set ylabel \"y (m)\"\n");
fprintf(plot, "set title \"#iteration = %d\"\n", n*every);
fprintf(plot, "plot \"output%.4d.dat\" with circles lc rgb \"blue\" fs transparent solid 0.1 border notitle, \"\" with points notitle, \"lungcancer.txt\" w l ls 2 notitle\n",n);
fprintf(plot,"\n");
}
fclose(plot);
}
//----------------------------------------------------------------------------------------------------
int main(int argc, char* argv[])
{
FILE* gp = popen("gnuplot > /dev/null 2>&1", "w");
set_gnuplot_condition(gp);
//----------------------------------------------------------------
FILE* fp1 = fopen("lungcancer.txt", "r");
int const n1 = 1484;
double *x1 = new double[n1];
double *y1 = new double[n1];
read_data_points(fp1, x1, y1, n1);
//----------------------------------------------------------------
/*for (int i=0; i<n1; i++){
if (i % 20 == 0) printf("%7.6lf %7.6lf %d %d %d\n", x1[i], y1[i], 0, 0, 0);
}*/
int Nb = 0;
Buble *buble = Initial_Probe_Position(&Nb);
int stop = 0;
double dt = 1.e-3;
//----------------------------------------------------------------
plot_buble(gp, buble, x1, y1, n1, nt, Nb);
while ((stop == 0) && (nt <= ntmax)) {
printf("n = %d\n", nt);
Sum_Buble_Force(buble, Nb);
New_Buble_Position(buble, dt, Nb);
//plot_buble(gp, buble, x1, y1, n1, nt, Nb);
stop = Stop_Searching(buble, eps, Nb);
if (nt % every == 0) save2file(buble,nt,Nb);
nt++;
printf("--------------------------------------------------\n");
}
Save_Buble(buble, Nb);
save2gnuplot(nt-1);
delete[] buble;
return 0;
}