回复 2楼 麦香
程序代码：
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
#ifndef M_PI
#ifdef _PI
#define M_PI _PI
#else
#define M_PI 3.14159265358979323846
#endif
#endif
#define DIM 3
#define XX 0
#define YY 1
#define ZZ 2
double H , D , wall_radius , inner_radius ,x_min ,x_max , one_fill_H , two_fill_H;
int one_Num , two_Num;
int nx=0,ny=0, nz=0;
double lx=0.0, ly=0.0, lz=0.0;
double circle_radius = 0.0;
double x=0.0,y=0.0,z=0.0;
char temp[100];
int count = 0;
int i,j,k;
int main()
{
FILE *pout = fopen("para.ini", "r");
FILE *fout = fopen("x.xyz", "w");
FILE *fout1 = fopen("data.txt", "w");
if(!pout)
{
printf("Can not open file to read!\n");
exit(1);
}
if(!fout)
{
printf("Can not open file to write!\n");
exit(1);
}
if(!fout1)
{
printf("Can not open file1 to write!\n");
exit(1);
}
fscanf(pout, "%s%lf", temp, &H);
fscanf(pout, "%s%lf", temp, &D);
fscanf(pout, "%s%lf", temp, &x_min);
fscanf(pout, "%s%lf", temp, &one_fill_H);
fscanf(pout, "%s%lf", temp, &two_fill_H);
fscanf(pout, "%s%d", temp, &one_Num);
fscanf(pout, "%s%d", temp, &two_Num);
fscanf(pout, "%s%lf", temp, &wall_radius);
fscanf(pout, "%s%lf", temp, &inner_radius);
fclose(pout);
x_max = H + 2.0*wall_radius + x_min;
double box[DIM]={H+2.0*wall_radius+2.0*x_min, D+4.0*wall_radius, D+4.0*wall_radius}; //filling district : H,W,D
printf("box[x] = %lf\nbox[y] = %lf\nbox[z] = %lf\n",box[XX],box[YY],box[ZZ]);
//----------------------------------------------------------------------------------
////---------------------------------------------------------------------------
double fill_circle_radius = (D-4.0*inner_radius) / 2.0;
int one_fill_nx = (int)(one_fill_H / (2.0*inner_radius));
double fill_x_step = 2.0*inner_radius;
int fill_n_left_right = int(fill_circle_radius / inner_radius);
double fill_step_left = inner_radius;
double angle = 0.0;
int n = 0;
double tmp_x =0.0;
for(k=0;k<=one_fill_nx;k++)
{
x = x_min+4.0*inner_radius + inner_radius + k*fill_x_step;
if (x>tmp_x) tmp_x = x;
for(i=0;i<fill_n_left_right;i++)
{
double tmp_circle_radius = fill_circle_radius - (double)(i*fill_step_left);
double fill_sin_theta_2 = sin(inner_radius/tmp_circle_radius);
double fill_theta = 2.0 * asin(fill_sin_theta_2);
int fill_n_circle = int(2.0*M_PI / fill_theta) ;
for(j=0;j<fill_n_circle;j++)
{
angle = fill_theta*j;
y = tmp_circle_radius*cos(angle);
z = sin(angle)*tmp_circle_radius;
n ++;
if(n <= one_Num)
{
count ++;
fprintf(fout1, "%8d%12.4f%12.4f%12.4f%8.4f%8.4f%8.4f\n", count, x, y, z, 0.0, 0.0, 0.0);
fprintf(fout, "%s%12.4f%12.4f%12.4f\n", "one", x, y, z);
}
}
}
}
int one_fill_Num = count;
printf ("one_fill_Num = %d\n",one_fill_Num);
printf("n = %d\n",n);
if (n> one_Num) printf("one_fill_H is too large\n");
//-------------------------------------------------------------------
int m = 0;
int two_fill_nx = (int)(two_fill_H / (2.0*inner_radius));
for(k=0;k<=two_fill_nx;k++)
{
x = tmp_x + 2.0*inner_radius + k*fill_x_step;
for(i=0;i<fill_n_left_right;i++)
{
double tmp_circle_radius = fill_circle_radius - (double)(i*fill_step_left);
double fill_sin_theta_2 = sin(inner_radius/tmp_circle_radius);
double fill_theta = 2.0 * asin(fill_sin_theta_2);
int fill_n_circle = int(2.0*M_PI / fill_theta) ;
for(j=0;j<fill_n_circle;j++)
{
angle = fill_theta*j;
y = tmp_circle_radius*cos(angle);
z = sin(angle)*tmp_circle_radius;
m ++;
if(m <= two_Num)
{
count ++;
fprintf(fout1, "%8d%12.4f%12.4f%12.4f%8.4f%8.4f%8.4f\n", count, x, y, z, 0.0, 0.0, 0.0);
fprintf(fout, "%s%12.4f%12.4f%12.4f\n", "two", x, y, z);
}
}
}
}
int two_fill_Num = count - one_fill_Num;
int fill_Num = count;
printf ("two_fill_Num = %d\n",two_fill_Num);
printf("m = %d\n",m);
if (m> two_Num) printf("two_fill_H is too large\n");
printf ("fill_Num = %d\n",fill_Num);
///-------------------------------------------------------------------------------------
int nx = (int)(box[XX] / (2.0*wall_radius)) + 1;
double x_step = box[XX] / nx;
// compute n_circle
ly = box[YY] - 2.0*wall_radius;
circle_radius = ly / 2.0;
double sin_theta_2 = sin(wall_radius/circle_radius);
double theta = 2.0 * asin(sin_theta_2);
int n_circle = int(2.0*M_PI / theta) + 1;
theta = 2.0*M_PI / (double)n_circle;
// loop along the box[XX]
for(j=0;j<nx;j++)
{
x = x_step/2.0 + j*x_step;
// loop the circle
for(i=0;i<n_circle;i++)
{
count ++;
double angle = theta*i;
y = sin(angle)*circle_radius;
z = cos(angle)*circle_radius;
fprintf(fout1, "%8d%12.4f%12.4f%12.4f%8.4f%8.4f%8.4f\n", count, x, y, z, 0.0, 0.0, 0.0);
fprintf(fout, "%s%12.4f%12.4f%12.4f\n", "wall", x, y, z);
}
}
//----------------------------------------
int n_left_right = int(circle_radius / wall_radius) + 1;
double step_left = circle_radius / (double)n_left_right;
// ---------------------------------
x = x_min;
for(i=1;i<n_left_right;i++)
{
double tmp_circle_radius = circle_radius - (double)(i*step_left);
double sin_theta_2 = sin(wall_radius/tmp_circle_radius);
double theta = 2.0 * asin(sin_theta_2);
int n_circle = int(2.0*M_PI / theta) + 1;
theta = 2.0*M_PI / (double)n_circle;
for(j=0;j<n_circle;j++)
{
count ++;
double angle = theta*j;
y = sin(angle)*tmp_circle_radius;
z = cos(angle)*tmp_circle_radius;
fprintf(fout1, "%8d%12.4f%12.4f%12.4f%8.4f%8.4f%8.4f\n", count, x, y, z, 0.0, 0.0, 0.0);
fprintf(fout, "%s%12.4f%12.4f%12.4f\n", "wall", x, y, z);
}
}
// ---------------------------------------
x = x_max;
for(i=1;i<n_left_right;i++)
{
double tmp_circle_radius = circle_radius - (double)(i*step_left);
double sin_theta_2 = sin(wall_radius/tmp_circle_radius);
double theta = 2.0 * asin(sin_theta_2);
int n_circle = int(2.0*M_PI / theta) + 1;
theta = 2.0*M_PI / (double)n_circle;
for(j=0;j<n_circle;j++)
{
count ++;
double angle = theta*j;
y = sin(angle)*tmp_circle_radius;
z = cos(angle)*tmp_circle_radius;
fprintf(fout1, "%8d%12.4f%12.4f%12.4f%8.4f%8.4f%8.4f\n",count, x, y, z, 0.0, 0.0, 0.0);
fprintf(fout, "%s%12.4f%12.4f%12.4f\n", "wall", x, y, z);
}
}
int wall_Num = count - fill_Num;
printf("wall_Num = %d\n",wall_Num);
printf("total_Num = %d\n",count);
fclose(fout);
return 0;
}
