50 #include "spdeclarations.h" 57 #define PI 3.141592653589793238462643383279 66 const char gen_par_psfx[] =
".prj";
67 const char beg_coord_psfx[]=
".sco";
68 const char res_coord_psfx[]=
".rco";
69 const char inp_coord_psfx[]=
".pco";
70 const char res_data_psfx[] =
".rst";
74 char gen_par_name[80];
75 char beg_coord_name[80];
76 char res_coord_name[80];
77 char inp_coord_name[80];
78 char res_data_name[80];
81 FILE* gen_par_file= NULL;
82 FILE* beg_coord_file= NULL;
83 FILE* res_coord_file= NULL;
84 FILE* inp_coord_file= NULL;
85 FILE* res_data_file= NULL;
89 int No_cells_x, No_cells_y, No_cells_z;
90 int No_cells, Ncell_min;
91 long Max_steps= 100000, Ntau= 1638400;
92 double Epsilon= 0.6, Epsilon_scl;
93 double Din, Dout0, Dout, Dout2;
94 double Pactual, Pnom0= 0.6, Pnomin, Pactual_old;
96 double Box_x, Box_y, Box_z, Half_x, Half_y, Half_z;
98 double Diam_min, Diam_max;
102 int *Link_head, *Link_list;
103 int *Ncell_bound_x, *Ncell_bound_y, *Ncell_bound_z;
104 double *Pbc_x, *Pbc_y, *Pbc_z;
107 double *Force_x, *Force_y, *Force_z;
108 double *k_freq, *l_freq, *g_r;
113 double Lambda = 0.04;
133 int main(
int argc,
char * argv[])
135 strcpy(prj_name,
"Project01");
139 strcpy(prj_name, argv[1]);
166 void make_project(
void)
168 strcpy(gen_par_name, prj_name); strcat(gen_par_name, gen_par_psfx);
170 strcpy(beg_coord_name, prj_name); strcat(beg_coord_name, beg_coord_psfx);
171 strcpy(res_coord_name, prj_name); strcat(res_coord_name, res_coord_psfx);
172 strcpy(inp_coord_name, prj_name); strcat(inp_coord_name, inp_coord_psfx);
173 strcpy(res_data_name, prj_name); strcat(res_data_name, res_data_psfx);
175 if (( gen_par_file= fopen(gen_par_name,
"r")) == NULL)
177 printf(
"Keine Steuerparameterdatei %s.prj gefunden.\n", prj_name);
178 fehlerAusgabe(
"read_param", -3);
181 if (( beg_coord_file= fopen(beg_coord_name,
"w") ) == NULL)
183 printf(
"Konnte %s.sco nicht erzeugen.\n", prj_name);
184 fehlerAusgabe(
"read_param", -3);
187 if (( inp_coord_file= fopen(inp_coord_name,
"r") ) == NULL)
189 printf(
"Keine Datendatei %s.pco gefunden.\n", prj_name);
193 if (( res_data_file= fopen(res_data_name,
"w") ) == NULL)
195 printf(
"Konnte %s.rst nicht erzeugen.\n", prj_name);
196 fehlerAusgabe(
"read_param", -3);
209 int k= 0,
m= 0, l= 0;
214 fgets(line, 80, gen_par_file);
218 fgets(line, 80, gen_par_file);
221 if ((Flag_Random = atoi(line))== EOF)
223 printf(
"Fehler in Zeile %d beim Einlesen von Flag_Random.\n", l);
224 fehlerAusgabe(
"read_param", -1);
228 fgets(line, 80, gen_par_file);
232 fgets(line, 80, gen_par_file);
235 if (!(No_parts = atoi(line)))
237 printf(
"Fehler in Zeile %d.\n", l);
238 fehlerAusgabe(
"read_param", -1);
242 fgets(line, 80, gen_par_file);
246 fgets(line, 80, gen_par_file);
249 if (!(No_cells_x = atoi(line)))
251 printf(
"Fehler in Zeile %d.\n", l);
252 fehlerAusgabe(
"read_param", -1);
255 fgets(line, 80, gen_par_file);
258 if (!(No_cells_y = atoi(line)))
260 printf(
"Fehler in Zeile %d.\n", l);
261 fehlerAusgabe(
"read_param", -1);
264 fgets(line, 80, gen_par_file);
267 if (!(No_cells_z = atoi(line)))
269 printf(
"Fehler in Zeile %d.\n", l);
270 fehlerAusgabe(
"read_param", -1);
274 fgets(line, 80, gen_par_file);
278 fgets(line, 80, gen_par_file);
281 if (!(Epsilon = atof(line)))
283 printf(
"Fehler in Zeile %d.\n", l);
284 fehlerAusgabe(
"read_param", -1);
290 printf(
"Parameter Epsilon muss groesser Null sein.\n");
291 fehlerAusgabe(
"read_param", -1);
296 fgets(line, 80, gen_par_file);
300 fgets(line, 80, gen_par_file);
303 if (!(Ntau = atoi(line)))
305 printf(
"Fehler in Zeile %d.\n", l);
306 fehlerAusgabe(
"read_param", -1);
312 printf(
"Parameter Ntau muss groesser Null sein.\n");
313 fehlerAusgabe(
"read_param", -1);
318 fgets(line, 80, gen_par_file);
322 fgets(line, 80, gen_par_file);
325 if (!(Pnom0 = atof(line)))
327 printf(
"Fehler in Zeile %d.\n", l);
328 fehlerAusgabe(
"read_param", -1);
332 if(Pnom0 <= 0 || Pnom0 > 1.0)
334 printf(
"Parameter Pnom muss zwischen Null und Eins liegen.\n");
335 fehlerAusgabe(
"read_param", -1);
341 fgets(line, 80, gen_par_file);
345 fgets(line, 80, gen_par_file);
348 if (!(Diam_max = atof(line)))
350 printf(
"Fehler in Zeile %d.\n", l);
351 fehlerAusgabe(
"read_param", -1);
357 printf(
"Maximaler Durchmesser muss groesser Null sein.\n");
358 fehlerAusgabe(
"read_param", -1);
362 fgets(line, 80, gen_par_file);
365 if (!(Diam_min = atof(line)))
367 printf(
"Fehler in Zeile %d.\n", l);
368 fehlerAusgabe(
"read_param", -1);
372 if(Diam_min <= 0 || Diam_min > Diam_max)
374 printf(
"Minimaler Durchmesser muss groesser Null und kleiner als maximaler Durchmesser sein.\n");
375 fehlerAusgabe(
"read_param", -1);
380 fgets(line, 80, gen_par_file);
383 fgets(line, 80, gen_par_file);
386 if (!(Max_steps = atoi(line)))
388 printf(
"Fehler in Zeile %d.\n", l);
389 fehlerAusgabe(
"read_param", -1);
395 printf(
"Maximale Anzahl Verdichtungsschritte muss positiv sein.\n");
396 fehlerAusgabe(
"read_param", -1);
401 fgets(line, 80, gen_par_file);
404 fgets(line, 80, gen_par_file);
407 if ((Distrfunct = atoi(line))== EOF)
409 printf(
"Fehler in Zeile %d.\n", l);
410 fehlerAusgabe(
"read_param", -1);
416 printf(
"Nummer der Verteilungsfunktion muss positiv oder 0 sein.\nEs werden gleich gro�e Kugeln erzeugt.\n");
417 fehlerAusgabe(
"read_param", -1);
421 if (Flag_Random == 1)
425 fgets(line, 80, gen_par_file);
428 fgets(line, 80, gen_par_file);
431 if (!(Mu = atof(line)))
433 printf(
"Fehler in Zeile %d.\n", l);
434 fehlerAusgabe(
"read_param", -1);
440 printf(
"Mittelwert fuer Normal- oder Log-Normalvert. muss positiv sein.\n");
441 fehlerAusgabe(
"read_param", -1);
445 fgets(line, 80, gen_par_file);
448 if (!(GSigma = atof(line)))
450 printf(
"Fehler in Zeile %d.\n", l);
451 fehlerAusgabe(
"read_param", -1);
457 printf(
"Standardabweichung muss positiv sein.\n");
458 fehlerAusgabe(
"read_param", -1);
462 fgets(line, 80, gen_par_file);
465 if (!(Vv = atof(line)))
467 printf(
"Fehler in Zeile %d.\n", l);
468 fehlerAusgabe(
"read_param", -1);
474 printf(
"Volumenverh�ltnis muss positiv sein und zwischen 0 und 1 liegen.\n");
475 fehlerAusgabe(
"read_param", -1);
479 fgets(line, 80, gen_par_file);
482 if (!(Exp = atof(line)))
484 printf(
"Fehler in Zeile %d.\n", l);
485 fehlerAusgabe(
"read_param", -1);
489 fgets(line, 80, gen_par_file);
492 if (!(Perc_number = atoi(line)))
494 printf(
"Fehler in Zeile %d.\n", l);
495 fehlerAusgabe(
"read_param", -1);
501 printf(
"Anzahl Prozentbins muss groesser 1 sein.\n");
502 fehlerAusgabe(
"read_param", -1);
506 Percentage = (
double *) calloc(Perc_number+2,
sizeof(
double));
507 Perc_values = (
double *) calloc(Perc_number+1,
sizeof(
double));
513 for (k= 1; k<= Perc_number; k++)
515 fgets(line, 80, gen_par_file);
518 if (!(*(Percentage+k) = atof(line)))
520 printf(
"Fehler in Zeile %d.\n", l);
521 fehlerAusgabe(
"read_param", -1);
524 sum+= *(Percentage+k);
529 printf(
"\nProgrammausfuehrung fehlerhaft beendet.\nBitte geben Sie eine gueltige Prozentverteilung an!\n");
533 for (
m= 0;
m < Perc_number;
m++)
535 fgets(line, 80, gen_par_file);
538 if (!(*(Perc_values+
m) = atof(line)))
540 printf(
"Fehler in Zeile %d.\n", l);
541 fehlerAusgabe(
"read_param", -1);
546 fclose(gen_par_file);
560 int set_coordinates(
void)
562 int ipart= 0, cnt= 0, r;
563 double x,
y, z,
d, max_x= 0, max_y= 0, max_z= 0;
566 static long seed = 0;
570 if (Flag_Random != 0)
572 while(!feof(inp_coord_file))
574 if(fscanf(inp_coord_file,
"%lf %lf %lf %lf", &x, &y, &z, &d) > 1)
580 rewind(inp_coord_file);
583 Diam = (
double *) calloc(No_parts,
sizeof(
double));
584 X = (
double *) calloc(No_parts,
sizeof(
double));
585 Y = (
double *) calloc(No_parts,
sizeof(
double));
586 Z = (
double *) calloc(No_parts,
sizeof(
double));
588 Force_x = (
double *) calloc(No_parts,
sizeof(
double));
589 Force_y = (
double *) calloc(No_parts,
sizeof(
double));
590 Force_z = (
double *) calloc(No_parts,
sizeof(
double));
592 if (Flag_Random != 0)
594 r= fscanf(inp_coord_file,
"%lf %lf %lf %lf", &x, &y, &z, &d);
598 printf(
"Fehler in Koordinatendatei in Zeile 0.\n");
599 fehlerAusgabe(
"set_coordinates", -1);
604 while(ipart < No_parts)
606 if(max_x< x) max_x= x;
607 if(max_y< y) max_y=
y;
608 if(max_z< z) max_z= z;
615 if (*(X+ipart) < 0.0 || *(X+ipart) > No_cells_x ||
616 *(Y+ipart) < 0.0 || *(Y+ipart) > No_cells_y ||
617 *(Z+ipart) < 0.0 || *(Z+ipart) > No_cells_z )
619 fehlerAusgabe(
"set_coordinates", -1);
622 r= fscanf(inp_coord_file,
"%lf %lf %lf %lf", &x, &y, &z, &d);
630 if (max_x <= 1.0 && max_y <= 1.0 && max_z <= 1.0)
632 for (ipart = 0; ipart < No_parts; ipart++)
634 *(X+ipart) *= No_cells_x;
635 *(Y+ipart) *= No_cells_y;
636 *(Z+ipart) *= No_cells_z;
638 if (*(X+ipart) < 0.0 || *(X+ipart) > No_cells_x ||
639 *(Y+ipart) < 0.0 || *(Y+ipart) > No_cells_y ||
640 *(Z+ipart) < 0.0 || *(Z+ipart) > No_cells_z )
642 fehlerAusgabe(
"set_coordinates", -1);
648 fclose(inp_coord_file);
649 inp_coord_file = NULL;
653 for (ipart = 0; ipart < No_parts; ipart++)
655 *(X+ipart) = ((
double)rand()/(double)RAND_MAX) * (double)No_cells_x;
656 *(Y+ipart) = ((
double)rand()/(double)RAND_MAX) * (double)No_cells_y;
657 *(Z+ipart) = ((
double)rand()/(double)RAND_MAX) * (double)No_cells_z;
659 if (*(X+ipart) < 0.0 || *(X+ipart) > No_cells_x ||
660 *(Y+ipart) < 0.0 || *(Y+ipart) > No_cells_y ||
661 *(Z+ipart) < 0.0 || *(Z+ipart) > No_cells_z )
663 fehlerAusgabe(
"set_coordinates", -1);
686 int list_structure(
void)
689 int ndim_x, ndim_y, ndim_z;
690 double corr, Summe_Kubikdurchm;
691 double Power = 0.3333333333333333333333333333333;
692 double Sphere_vol = 1.9098593171027440292266051604702;
694 No_cells = No_cells_x * No_cells_y * No_cells_z;
695 Ncell_min = (No_cells_x < No_cells_y) ? No_cells_x : No_cells_y;
696 Ncell_min = (No_cells_z < Ncell_min) ? No_cells_z : Ncell_min;
697 Box_x = (double) No_cells_x;
698 Box_y = (double) No_cells_y;
699 Box_z = (double) No_cells_z;
700 Half_x = 0.5 * Box_x;
701 Half_y = 0.5 * Box_y;
702 Half_z = 0.5 * Box_z;
703 Diam_dens = No_cells * Sphere_vol;
705 Diam_max /= Diam_min;
711 if (Flag_Random == 0)
713 setDistribution(Distrfunct, Diam);
714 qsort ((
void *) Diam, No_parts,
sizeof(
double), dblcmp);
717 for (ipart = 0; ipart < No_parts; ipart++)
718 fprintf(beg_coord_file,
"%.6f\t%.6f\t%.6f\t%.6f\n", *(X+ipart), *(Y+ipart), *(Z+ipart), *(Diam+ipart));
720 fflush(beg_coord_file);
722 fclose(beg_coord_file);
723 beg_coord_file = NULL;
725 Summe_Kubikdurchm = 0.0;
727 for (ipart = 0; ipart < No_parts; ipart++)
728 Summe_Kubikdurchm += *(Diam+ipart) * *(Diam+ipart) * *(Diam+ipart);
730 Diam_dens /= Summe_Kubikdurchm;
733 Dout0 = pow (Diam_dens * Pnom0 , Power);
735 Relax = 0.5 * Dout / Ntau;
737 Epsilon_scl = Epsilon / Dout0;
739 Link_head = (
int *) calloc(No_cells,
sizeof(
int));
740 Link_list = (
int *) calloc(No_parts,
sizeof(
int));
742 ndim_x = 3 * No_cells_x - 2;
743 ndim_y = 3 * No_cells_y - 2;
744 ndim_z = 3 * No_cells_z - 2;
746 Ncell_bound_x = (
int *) calloc(ndim_x+1,
sizeof(
int));
747 if (Ncell_bound_x == NULL)
749 fehlerAusgabe(
"list_structure", -2);
752 Ncell_bound_y = (
int *) calloc(ndim_y+1,
sizeof(
int));
753 if (Ncell_bound_y == NULL)
755 fehlerAusgabe(
"list_structure", -2);
758 Ncell_bound_z = (
int *) calloc(ndim_z+1,
sizeof(
int));
759 if (Ncell_bound_z == NULL)
761 fehlerAusgabe(
"list_structure", -2);
764 Pbc_x = (
double *) calloc(ndim_x+1,
sizeof(
double));
767 fehlerAusgabe(
"list_structure", -2);
770 Pbc_y = (
double *) calloc(ndim_y+1,
sizeof(
double));
773 fehlerAusgabe(
"list_structure", -2);
776 Pbc_z = (
double *) calloc(ndim_z+1,
sizeof(
double));
779 fehlerAusgabe(
"list_structure", -2);
785 for (i= 1; i <= ndim_x; i++)
787 *(Ncell_bound_x+i) = imult * No_cells_y * No_cells_z;
792 if (imult >= No_cells_x)
802 for (i = 1; i <= ndim_y; i++)
804 *(Ncell_bound_y+i) = imult * No_cells_z;
809 if (imult >= No_cells_y)
819 for (i= 1; i <= ndim_z; i++)
821 *(Ncell_bound_z+i) = imult;
826 if (imult >= No_cells_z)
843 int packing_loop(
void)
845 date_time (res_data_file, (
long *) Time_out);
851 Pactual = Din*Din*Din / Diam_dens;
853 for (Step= 1; Step <= Max_steps && Flag_End == FALSE; Step++)
857 Epsilon_scl = Epsilon / Dout;
884 Pnomin = Dout * Dout * Dout / Diam_dens;
885 Pactual = Din * Din * Din / Diam_dens;
901 if (Din * *(Diam+No_parts-1) > Ncell_min)
902 Din = Ncell_min / *(Diam+No_parts - 1);
904 Pactual = Din * Din * Din / Diam_dens;
911 iprec= (int) -log10 (Pnomin - Pactual);
929 int interaction(
void)
936 for(icell= 0; icell < No_cells; icell++)
938 *(Link_head+icell)= -1;
941 for(ipart= 0; ipart < No_parts; ipart++)
943 *(Force_x+ipart)= 0.0;
944 *(Force_y+ipart)= 0.0;
945 *(Force_z+ipart)= 0.0;
948 for (ipart= 0; ipart < No_parts; ipart++)
950 Rcut= ((ipart == 0) ? 0.0 : 0.5) * Dout * (*(Diam + ipart - 1) + *(Diam + ipart));
952 if ((
int) (2.0 * Rcut) < Ncell_min - 1)
953 force_part(ipart, Rcut);
973 double x_buff, y_buff, z_buff;
975 for (i = 0; i < No_parts; i++)
977 x_buff = *(X+i)+(*(Force_x+i) * Epsilon_scl) / *(Diam+i);
978 y_buff = *(Y+i)+(*(Force_y+i) * Epsilon_scl) / *(Diam+i);
979 z_buff = *(Z+i)+(*(Force_z+i) * Epsilon_scl) / *(Diam+i);
982 x_buff+= (ceil(-(x_buff/Box_x))*Box_x);
985 x_buff-= (floor(x_buff/Box_x)*Box_x);
988 y_buff+= (ceil(-(y_buff/Box_y))*Box_y);
991 y_buff-= (floor(y_buff/Box_y)*Box_y);
994 z_buff+= (ceil(-(z_buff/Box_z))*Box_z);
997 z_buff-= (floor(z_buff/Box_z)*Box_z);
1012 void force_part(
int ipart_p,
double Rcut)
1014 int jpart= 0, ifirst= 0;
1015 int leap_x= 0, leap_y= 0, leap_z= 0;
1016 int icell_x= 0, icell_xy= 0, icell_y= 0, icell_z= 0, icell= 0;
1017 int low_cell_x, low_cell_y, low_cell_z;
1018 int lim_cell_x, lim_cell_y, lim_cell_z;
1019 int idif_x= 0, idif_y= 0, idif_z= 0, idif= 0;
1021 double pos_x, pos_y, pos_z, dist;
1022 double diam_i, dsum,
sigma;
1023 double dif_x, dif_y, dif_z, difpot;
1024 double f_x, f_y, f_z;
1026 pos_x = *(X+ipart_p);
1027 pos_y = *(Y+ipart_p);
1028 pos_z = *(Z+ipart_p);
1029 diam_i = *(Diam+ipart_p);
1031 low_cell_x = (int) (pos_x + No_cells_x - Rcut);
1032 low_cell_y = (int) (pos_y + No_cells_y - Rcut);
1033 low_cell_z = (int) (pos_z + No_cells_z - Rcut);
1034 lim_cell_x = (int) (pos_x + No_cells_x + Rcut);
1035 lim_cell_y = (int) (pos_y + No_cells_y + Rcut);
1036 lim_cell_z = (int) (pos_z + No_cells_z + Rcut);
1037 idif_x = (*(Ncell_bound_x+lim_cell_x) > *(Ncell_bound_x+low_cell_x)) ? 0 : 4;
1038 idif_y = (*(Ncell_bound_y+lim_cell_y) > *(Ncell_bound_y+low_cell_y)) ? 0 : 2;
1039 idif_z = (*(Ncell_bound_z+lim_cell_z) > *(Ncell_bound_z+low_cell_z)) ? 0 : 1;
1040 idif = idif_x + idif_y + idif_z;
1045 for (leap_x= low_cell_x; leap_x<= lim_cell_x; leap_x++)
1047 icell_x = *(Ncell_bound_x+leap_x);
1049 for (leap_y= low_cell_y; leap_y<= lim_cell_y; leap_y++)
1051 icell_xy = *(Ncell_bound_y+leap_y) + icell_x;
1053 for (leap_z= low_cell_z; leap_z<= lim_cell_z; leap_z++)
1055 icell = *(Ncell_bound_z+leap_z) + icell_xy;
1056 jpart = *(Link_head+icell);
1060 dif_x = *(X+jpart) - pos_x;
1061 dif_y = *(Y+jpart) - pos_y;
1062 dif_z = *(Z+jpart) - pos_z;
1064 dist = sqrt(dif_x * dif_x + dif_y * dif_y + dif_z * dif_z);
1065 dsum = 0.5 * (diam_i + *(Diam+jpart));
1066 sigma = Dout * dsum;
1070 if (dist < sqrt(Din)*dsum)
1071 Din = (dist*dist)/(dsum*dsum);
1073 difpot= (0.5 * Dout2 * diam_i * (*(Diam+jpart))) * (1.0 - (dist * dist) / (sigma *
sigma))/dist;
1075 f_x = difpot * dif_x;
1076 f_y = difpot * dif_y;
1077 f_z = difpot * dif_z;
1079 *(Force_x+ipart_p) -= f_x;
1080 *(Force_y+ipart_p) -= f_y;
1081 *(Force_z+ipart_p) -= f_z;
1082 *(Force_x+jpart) += f_x;
1083 *(Force_y+jpart) += f_y;
1084 *(Force_z+jpart) += f_z;
1086 jpart = *(Link_list+jpart);
1094 for (leap_x=low_cell_x; leap_x<=lim_cell_x; leap_x++)
1096 icell_x = *(Ncell_bound_x+leap_x);
1098 for (leap_y=low_cell_y; leap_y<=lim_cell_y; leap_y++)
1100 icell_xy = *(Ncell_bound_y+leap_y) + icell_x;
1102 for (leap_z=low_cell_z; leap_z<=lim_cell_z; leap_z++)
1104 icell = *(Ncell_bound_z+leap_z) + icell_xy;
1105 jpart = *(Link_head+icell);
1109 dif_x = *(X+jpart) - pos_x;
1110 dif_y = *(Y+jpart) - pos_y;
1111 dif_z = *(Z+jpart) - pos_z + *(Pbc_z+leap_z);
1113 dist = sqrt(dif_x * dif_x + dif_y * dif_y + dif_z * dif_z);
1114 dsum = 0.5 * (diam_i + *(Diam+jpart));
1115 sigma = Dout * dsum;
1119 if (dist < sqrt(Din)*dsum)
1120 Din = (dist*dist)/(dsum*dsum);
1122 difpot= (0.5 * Dout2 * diam_i * (*(Diam+jpart))) * (1.0 - (dist * dist) / (sigma *
sigma))/dist;
1124 f_x = difpot * dif_x;
1125 f_y = difpot * dif_y;
1126 f_z = difpot * dif_z;
1128 *(Force_x+ipart_p) -= f_x;
1129 *(Force_y+ipart_p) -= f_y;
1130 *(Force_z+ipart_p) -= f_z;
1131 *(Force_x+jpart) += f_x;
1132 *(Force_y+jpart) += f_y;
1133 *(Force_z+jpart) += f_z;
1135 jpart = *(Link_list+jpart);
1143 for (leap_x= low_cell_x; leap_x <= lim_cell_x; leap_x++)
1145 icell_x = *(Ncell_bound_x+leap_x);
1147 for (leap_y= low_cell_y; leap_y <= lim_cell_y; leap_y++)
1149 icell_xy = *(Ncell_bound_y+leap_y) + icell_x;
1151 for (leap_z= low_cell_z; leap_z <= lim_cell_z; leap_z++)
1153 icell = *(Ncell_bound_z+leap_z) + icell_xy;
1154 jpart = *(Link_head+icell);
1158 dif_x = *(X+jpart) - pos_x;
1159 dif_y = *(Y+jpart) - pos_y + *(Pbc_y+leap_y);
1160 dif_z = *(Z+jpart) - pos_z;
1162 dist = sqrt(dif_x * dif_x + dif_y * dif_y + dif_z * dif_z);
1163 dsum = 0.5 * (diam_i + *(Diam+jpart));
1164 sigma = Dout * dsum;
1168 if (dist < sqrt(Din)*dsum)
1169 Din = (dist*dist)/(dsum*dsum);
1171 difpot= (0.5 * Dout2 * diam_i * (*(Diam+jpart))) * (1.0 - (dist * dist) / (sigma *
sigma))/dist;
1173 f_x = difpot * dif_x;
1174 f_y = difpot * dif_y;
1175 f_z = difpot * dif_z;
1177 *(Force_x+ipart_p) -= f_x;
1178 *(Force_y+ipart_p) -= f_y;
1179 *(Force_z+ipart_p) -= f_z;
1180 *(Force_x+jpart) += f_x;
1181 *(Force_y+jpart) += f_y;
1182 *(Force_z+jpart) += f_z;
1184 jpart = *(Link_list+jpart);
1192 for (leap_x= low_cell_x; leap_x <= lim_cell_x; leap_x++)
1194 icell_x = *(Ncell_bound_x+leap_x);
1196 for (leap_y= low_cell_y; leap_y <= lim_cell_y; leap_y++)
1198 icell_xy = *(Ncell_bound_y+leap_y) + icell_x;
1200 for (leap_z= low_cell_z; leap_z <= lim_cell_z; leap_z++)
1202 icell = *(Ncell_bound_z+leap_z) + icell_xy;
1203 jpart = *(Link_head+icell);
1207 dif_x = *(X+jpart) - pos_x;
1208 dif_y = *(Y+jpart) - pos_y + *(Pbc_y+leap_y);
1209 dif_z = *(Z+jpart) - pos_z + *(Pbc_z+leap_z);
1211 dist = sqrt(dif_x * dif_x + dif_y * dif_y + dif_z * dif_z);
1212 dsum = 0.5 * (diam_i + *(Diam+jpart));
1213 sigma = Dout * dsum;
1217 if (dist < sqrt(Din)*dsum)
1218 Din = (dist*dist)/(dsum*dsum);
1220 difpot= (0.5 * Dout2 * diam_i * (*(Diam+jpart))) * (1.0 - (dist * dist) / (sigma *
sigma))/dist;
1222 f_x = difpot * dif_x;
1223 f_y = difpot * dif_y;
1224 f_z = difpot * dif_z;
1226 *(Force_x+ipart_p) -= f_x;
1227 *(Force_y+ipart_p) -= f_y;
1228 *(Force_z+ipart_p) -= f_z;
1229 *(Force_x+jpart) += f_x;
1230 *(Force_y+jpart) += f_y;
1231 *(Force_z+jpart) += f_z;
1233 jpart = *(Link_list+jpart);
1241 for (leap_x= low_cell_x; leap_x <= lim_cell_x; leap_x++)
1243 icell_x = *(Ncell_bound_x+leap_x);
1245 for (leap_y= low_cell_y; leap_y <= lim_cell_y; leap_y++)
1247 icell_xy = *(Ncell_bound_y+leap_y) + icell_x;
1249 for (leap_z= low_cell_z; leap_z <= lim_cell_z; leap_z++)
1251 icell = *(Ncell_bound_z+leap_z) + icell_xy;
1252 jpart = *(Link_head+icell);
1256 dif_x = *(X+jpart) - pos_x + *(Pbc_x+leap_x);
1257 dif_y = *(Y+jpart) - pos_y;
1258 dif_z = *(Z+jpart) - pos_z;
1260 dist = sqrt(dif_x * dif_x + dif_y * dif_y + dif_z * dif_z);
1261 dsum = 0.5 * (diam_i + *(Diam+jpart));
1262 sigma = Dout * dsum;
1266 if (dist < sqrt(Din)*dsum)
1267 Din = (dist*dist)/(dsum*dsum);
1269 difpot= (0.5 * Dout2 * diam_i * (*(Diam+jpart))) * (1.0 - (dist * dist) / (sigma *
sigma))/dist;
1271 f_x = difpot * dif_x;
1272 f_y = difpot * dif_y;
1273 f_z = difpot * dif_z;
1275 *(Force_x+ipart_p) -= f_x;
1276 *(Force_y+ipart_p) -= f_y;
1277 *(Force_z+ipart_p) -= f_z;
1278 *(Force_x+jpart) += f_x;
1279 *(Force_y+jpart) += f_y;
1280 *(Force_z+jpart) += f_z;
1282 jpart = *(Link_list+jpart);
1290 for (leap_x= low_cell_x; leap_x <= lim_cell_x; leap_x++)
1292 icell_x = *(Ncell_bound_x+leap_x);
1294 for (leap_y= low_cell_y; leap_y <= lim_cell_y; leap_y++)
1296 icell_xy = *(Ncell_bound_y+leap_y) + icell_x;
1298 for (leap_z= low_cell_z; leap_z <= lim_cell_z; leap_z++)
1300 icell = *(Ncell_bound_z+leap_z) + icell_xy;
1301 jpart = *(Link_head+icell);
1305 dif_x = *(X+jpart) - pos_x + *(Pbc_x+leap_x);
1306 dif_y = *(Y+jpart) - pos_y;
1307 dif_z = *(Z+jpart) - pos_z + *(Pbc_z+leap_z);
1309 dist = sqrt(dif_x * dif_x + dif_y * dif_y + dif_z * dif_z);
1310 dsum = 0.5 * (diam_i + *(Diam+jpart));
1311 sigma = Dout * dsum;
1315 if (dist < sqrt(Din)*dsum)
1316 Din = (dist*dist)/(dsum*dsum);
1318 difpot= (0.5 * Dout2 * diam_i * (*(Diam+jpart))) * (1.0 - (dist * dist) / (sigma *
sigma))/dist;
1320 f_x = difpot * dif_x;
1321 f_y = difpot * dif_y;
1322 f_z = difpot * dif_z;
1324 *(Force_x+ipart_p) -= f_x;
1325 *(Force_y+ipart_p) -= f_y;
1326 *(Force_z+ipart_p) -= f_z;
1327 *(Force_x+jpart) += f_x;
1328 *(Force_y+jpart) += f_y;
1329 *(Force_z+jpart) += f_z;
1331 jpart = *(Link_list+jpart);
1339 for (leap_x= low_cell_x; leap_x <= lim_cell_x; leap_x++)
1341 icell_x = *(Ncell_bound_x+leap_x);
1343 for (leap_y= low_cell_y; leap_y <= lim_cell_y; leap_y++)
1345 icell_xy = *(Ncell_bound_y+leap_y) + icell_x;
1347 for (leap_z= low_cell_z; leap_z <= lim_cell_z; leap_z++)
1349 icell = *(Ncell_bound_z+leap_z) + icell_xy;
1350 jpart = *(Link_head+icell);
1354 dif_x = *(X+jpart) - pos_x + *(Pbc_x+leap_x);
1355 dif_y = *(Y+jpart) - pos_y + *(Pbc_y+leap_y);
1356 dif_z = *(Z+jpart) - pos_z;
1358 dist = sqrt(dif_x * dif_x + dif_y * dif_y + dif_z * dif_z);
1359 dsum = 0.5 * (diam_i + *(Diam+jpart));
1360 sigma = Dout * dsum;
1364 if (dist < sqrt(Din)*dsum)
1365 Din = (dist*dist)/(dsum*dsum);
1367 difpot= (0.5 * Dout2 * diam_i * (*(Diam+jpart))) * (1.0 - (dist * dist) / (sigma *
sigma))/dist;
1369 f_x = difpot * dif_x;
1370 f_y = difpot * dif_y;
1371 f_z = difpot * dif_z;
1373 *(Force_x+ipart_p) -= f_x;
1374 *(Force_y+ipart_p) -= f_y;
1375 *(Force_z+ipart_p) -= f_z;
1376 *(Force_x+jpart) += f_x;
1377 *(Force_y+jpart) += f_y;
1378 *(Force_z+jpart) += f_z;
1380 jpart= *(Link_list+jpart);
1388 for (leap_x= low_cell_x; leap_x <= lim_cell_x; leap_x++)
1390 icell_x = *(Ncell_bound_x+leap_x);
1392 for (leap_y= low_cell_y; leap_y <= lim_cell_y; leap_y++)
1394 icell_xy = *(Ncell_bound_y+leap_y) + icell_x;
1396 for (leap_z= low_cell_z; leap_z <= lim_cell_z; leap_z++)
1398 icell = *(Ncell_bound_z+leap_z) + icell_xy;
1399 jpart = *(Link_head+icell);
1403 dif_x = *(X+jpart) - pos_x + *(Pbc_x+leap_x);
1404 dif_y = *(Y+jpart) - pos_y + *(Pbc_y+leap_y);
1405 dif_z = *(Z+jpart) - pos_z + *(Pbc_z+leap_z);
1407 dist = sqrt(dif_x * dif_x + dif_y * dif_y + dif_z * dif_z);
1408 dsum = 0.5 * (diam_i + *(Diam+jpart));
1409 sigma = Dout * dsum;
1413 if (dist < sqrt(Din)*dsum)
1414 Din = (dist*dist)/(dsum*dsum);
1416 difpot= (0.5 * Dout2 * diam_i * (*(Diam+jpart))) * (1.0 - (dist * dist) / (sigma *
sigma))/dist;
1418 f_x = difpot * dif_x;
1419 f_y = difpot * dif_y;
1420 f_z = difpot * dif_z;
1422 *(Force_x+ipart_p) -= f_x;
1423 *(Force_y+ipart_p) -= f_y;
1424 *(Force_z+ipart_p) -= f_z;
1425 *(Force_x+jpart) += f_x;
1426 *(Force_y+jpart) += f_y;
1427 *(Force_z+jpart) += f_z;
1429 jpart= *(Link_list+jpart);
1437 icell_x = (int) pos_x;
1438 icell_y = (int) pos_y;
1439 icell_z = (int) pos_z;
1440 icell = No_cells_z * (No_cells_y * icell_x + icell_y) + icell_z;
1441 ifirst = *(Link_head+icell);
1442 *(Link_list+ipart_p) = ifirst;
1443 *(Link_head+icell) = ipart_p;
1453 void force_all(
int ipart_p)
1455 int jpart= 0, ifirst= 0;
1456 int icell_x= 0, icell_y= 0, icell_z= 0, icell= 0;
1458 double pos_x= 0.0, pos_y= 0.0, pos_z= 0.0, dist= 0.0, sigma= 0.0;
1459 double diam_i= 0.0, dsum= 0.0;
1460 double dif_x= 0.0, dif_y= 0.0, dif_z= 0.0, difpot= 0.0;
1461 double f_x= 0, f_y= 0, f_z= 0;
1463 pos_x = *(X+ipart_p);
1464 pos_y = *(Y+ipart_p);
1465 pos_z = *(Z+ipart_p);
1466 diam_i = *(Diam+ipart_p);
1468 for (icell= 0; icell < No_cells; icell++)
1470 jpart = *(Link_head+icell);
1474 dif_x = *(X+jpart) - pos_x;
1475 dif_y = *(Y+jpart) - pos_y;
1476 dif_z = *(Z+jpart) - pos_z;
1478 if (dif_x < -Half_x)
1481 if (dif_y < -Half_y)
1484 if (dif_z < -Half_z)
1496 dist = sqrt(dif_x * dif_x + dif_y * dif_y + dif_z * dif_z);
1497 dsum = 0.5 * (diam_i + *(Diam+jpart));
1498 sigma = Dout * dsum;
1502 if (dist < sqrt(Din)*dsum)
1503 Din = (dist*dist)/(dsum*dsum);
1505 difpot= (0.5 * Dout2 * diam_i * (*(Diam+jpart))) * (1.0 - (dist * dist) / (sigma *
sigma))/dist;
1507 f_x = difpot * dif_x;
1508 f_y = difpot * dif_y;
1509 f_z = difpot * dif_z;
1511 *(Force_x+ipart_p) -= f_x;
1512 *(Force_y+ipart_p) -= f_y;
1513 *(Force_z+ipart_p) -= f_z;
1514 *(Force_x+jpart) += f_x;
1515 *(Force_y+jpart) += f_y;
1516 *(Force_z+jpart) += f_z;
1518 jpart = *(Link_list+jpart);
1522 icell_x = (int) pos_x;
1523 icell_y = (int) pos_y;
1524 icell_z = (int) pos_z;
1525 icell = No_cells_z * (No_cells_y * icell_x + icell_y) + icell_z;
1526 ifirst = *(Link_head+icell);
1527 *(Link_list+ipart_p) = ifirst;
1528 *(Link_head+icell) = ipart_p;
1539 void swap(
int i,
int j)
1542 temp = *(Diam+j), *(Diam+j) = *(Diam+i), *(Diam+i) =
temp;
1543 temp = *(X+j), *(X+j) = *(X+i), *(X+i) =
temp;
1544 temp = *(Y+j), *(Y+j) = *(Y+i), *(Y+i) =
temp;
1545 temp = *(Z+j), *(Z+j) = *(Z+i), *(Z+i) =
temp;
1558 for(gap = No_parts/2; gap > 0; gap /= 2)
1559 for(i = gap; i < No_parts; i++)
1560 for(j = i - gap; j >= 0 && *(Diam+j) > *(Diam+j+gap); j -= gap)
1571 void result(
int nstep)
1574 char* DistFunct[]= {
"constant",
"two-point",
"uniform",
"gaussian",
"log-normal",
"power",
"percentage",
"predefined (file)"};
1576 if(Distrfunct == 11)
1581 fprintf(res_data_file,
"----------------------------------------\n");
1582 if(Distrfunct == 11)
1583 fprintf(res_data_file,
"\nDistribution : %s",
"predefined (file)");
1585 fprintf(res_data_file,
"\nDistribution : %s", DistFunct[Distrfunct]);
1587 fprintf(res_data_file,
"\nCells in x-dir. : %d", No_cells_x);
1588 fprintf(res_data_file,
"\nCells in y-dir. : %d", No_cells_y);
1589 fprintf(res_data_file,
"\nCells in z-dir. : %d", No_cells_z);
1590 fprintf(res_data_file,
"\nNumber of Spheres: %d", No_parts);
1591 fprintf(res_data_file,
"\nDensity : %5.13f",Pactual);
1592 fprintf(res_data_file,
"\nCalculation steps: %d\n", nstep);
1593 fprintf(res_data_file,
"----------------------------------------\n");
1595 date_time (res_data_file, (
long *) Time_out);
1597 fflush(res_data_file);
1611 if (Link_head != NULL)
1614 if (Link_list != NULL)
1617 if (Force_x != NULL)
1619 if (Force_y != NULL)
1621 if (Force_z != NULL)
1624 if (Ncell_bound_x != NULL)
1625 free (Ncell_bound_x);
1626 if (Ncell_bound_y != NULL)
1627 free (Ncell_bound_y);
1628 if (Ncell_bound_z != NULL)
1629 free (Ncell_bound_z);
1647 if (Percentage != NULL)
1649 if (Perc_values != NULL)
1652 if (gen_par_file != NULL) fclose(gen_par_file);
1653 if (beg_coord_file != NULL) fclose(beg_coord_file);
1654 if (res_coord_file != NULL) fclose(res_coord_file);
1655 if (inp_coord_file != NULL) fclose(inp_coord_file);
1656 if (res_data_file != NULL) fclose(res_data_file);
1666 void write_sphere_data()
1669 char filename[80]=
"";
1671 strcpy(filename, res_coord_name);
1673 res_coord_file= fopen(filename,
"w");
1675 for (i= 0; i < No_parts; i++)
1677 fprintf(res_coord_file,
"%.14lf\t%.14lf\t%.14lf\t%.14lf\n", *(X+i), *(Y+i), *(Z+i), *(Diam+i)*Din);
1678 fflush(res_coord_file);
1691 int dblcmp(
const void *v1,
const void *v2)
1693 double t = *(
double *)v1 - *(
double *)v2;
1695 return ((t > 0) ? 1 : -1);
1705 void date_time(FILE * pfile_p,
long * time_p)
1709 *time_p =
time(&clock);
1710 fprintf(pfile_p,
"%s", ctime(&clock));
1720 int fehlerAusgabe(
char *funktion,
int fehlerNummer)
1722 int fehlerIndex= (-fehlerNummer) - 1;
1725 const char fehlerNachricht[FANZ][60+1]=
1727 " -1 Fehler in den Argumenten",
1728 " -2 Speicherplatz konnte nicht bereit gestellt werden",
1729 " -3 Datei konnte nicht zum Lesen geoeffnet werden",
1730 " -4 Datensatz konnte nicht gelesen werden",
1731 " -5 Datei konnte nicht zum Schreiben geoeffnet werden",
1732 " -6 Datensatz konnte nicht geschrieben werden",
1733 " -7 Datei konnte nicht geschlossen werden",
1736 if (fehlerIndex >= 0 && fehlerIndex < FANZ)
1738 printf(
"\nFehlermeldung: %s\n", fehlerNachricht[fehlerIndex]);
real(dp) temp
temperature (K)
real(dp), dimension(:), allocatable y
state variables
real(dp) sigma
interfacial tension
int main(int argc, char *argv[])
Reads parameters. Creates struts and walls. Saves foam morphology to a file.
m
Bubble Growth Application Recipe.
real(dp), dimension(:), allocatable d
diffusion coefficients (for each dissolved gas)