#include "proj.h"
float Scape[NMax][NMax];
float scale[2*RANGE+1][2*RANGE+1];
Point pt[MAXPT];
int npt;
int cno;
int fdl;
float rmat[21][21] =
{{ 14.14213562, 13.45362405, 12.80624847, 12.20655562,
11.66190379, 11.18033989, 10.77032961, 10.44030651,
10.19803903, 10.04987562, 10 , 10.04987562,
10.19803903, 10.44030651, 10.77032961, 11.18033989,
11.66190379, 12.20655562, 12.80624847, 13.45362405,
14.14213562}, {13.45362405, 12.72792206, 12.04159458,
11.40175425, 10.81665383, 10.29563014, 9.848857802,
9.486832981, 9.219544457, 9.055385138, 9 ,
9.055385138, 9.219544457, 9.486832981, 9.848857802,
10.29563014, 10.81665383, 11.40175425, 12.04159458,
12.72792206, 13.45362405}, { 12.80624847, 12.04159458,
11.3137085 , 10.63014581, 10 , 9.433981132,
8.94427191 , 8.544003745, 8.246211251, 8.062257748,
8 , 8.062257748, 8.246211251, 8.544003745,
8.94427191 , 9.433981132, 10 , 10.63014581,
11.3137085 , 12.04159458, 12.80624847}, {12.20655562,
11.40175425, 10.63014581, 9.899494937, 9.219544457,
8.602325267, 8.062257748, 7.615773106, 7.280109889,
7.071067812, 7 , 7.071067812, 7.280109889,
7.615773106, 8.062257748, 8.602325267, 9.219544457,
9.899494937, 10.63014581, 11.40175425, 12.20655562},
{11.66190379, 10.81665383, 10 , 9.219544457,
8.485281374, 7.810249676, 7.211102551, 6.708203932,
6.32455532 , 6.08276253, 6 , 6.08276253,
6.32455532 , 6.708203932, 7.211102551, 7.810249676,
8.485281374, 9.219544457, 10 , 10.81665383,
11.66190379}, {11.18033989, 10.29563014, 9.433981132,
8.602325267, 7.810249676, 7.071067812, 6.403124237,
5.830951895, 5.385164807, 5.099019514, 5 ,
5.099019514, 5.385164807, 5.830951895, 6.403124237,
7.071067812, 7.810249676, 8.602325267, 9.433981132,
10.29563014, 11.18033989}, {10.77032961, 9.848857802,
8.94427191 , 8.062257748, 7.211102551, 6.403124237,
5.656854249, 5 , 4.472135955, 4.123105626,
4 , 4.123105626, 4.472135955, 5 ,
5.656854249, 6.403124237, 7.211102551, 8.062257748,
8.94427191 , 9.848857802, 10.77032961}, {10.44030651,
9.486832981, 8.544003745, 7.615773106, 6.708203932,
5.830951895, 5 , 4.242640687, 3.605551275,
3.16227766 , 3 , 3.16227766, 3.605551275,
4.242640687, 5 , 5.830951895, 6.708203932,
7.615773106, 8.544003745, 9.486832981, 10.44030651},
{10.19803903, 9.219544457, 8.246211251, 7.280109889,
6.32455532 , 5.385164807, 4.472135955, 3.605551275,
2.828427125, 2.236067977, 2 , 2.236067977,
2.828427125, 3.605551275, 4.472135955, 5.385164807,
6.32455532 , 7.280109889, 8.246211251, 9.219544457,
10.19803903}, {10.04987562, 9.055385138, 8.062257748,
7.071067812, 6.08276253, 5.099019514, 4.123105626,
3.16227766 , 2.236067977, 1.414213562, 1 ,
1.414213562, 2.236067977, 3.16227766, 4.123105626,
5.099019514, 6.08276253, 7.071067812, 8.062257748,
9.055385138, 10.04987562}, {10 , 9 ,
8 , 7 , 6 , 5 ,
4 , 3 , 2 , 1 ,
10 , 1 , 2 , 3 ,
4 , 5 , 6 , 7 ,
8 , 9 , 10 }, {10.04987562,
9.055385138, 8.062257748, 7.071067812, 6.08276253,
5.099019514, 4.123105626, 3.16227766, 2.236067977,
1.414213562, 1 , 1.414213562, 2.236067977,
3.16227766 , 4.123105626, 5.099019514, 6.08276253,
7.071067812, 8.062257748, 9.055385138, 10.04987562},
{10.19803903, 9.219544457, 8.246211251, 7.280109889,
6.32455532 , 5.385164807, 4.472135955, 3.605551275 ,
2.828427125, 2.236067977, 2 , 2.236067977,
2.828427125, 3.605551275, 4.472135955, 5.385164807,
6.32455532 , 7.280109889, 8.246211251, 9.219544457,
10.19803903}, {10.44030651, 9.486832981, 8.544003745,
7.615773106, 6.708203932, 5.830951895, 5 ,
4.242640687, 3.605551275, 3.16227766, 3 ,
3.16227766 , 3.605551275, 4.242640687, 5 ,
5.830951895, 6.708203932, 7.615773106, 8.544003745,
9.486832981, 10.44030651}, {10.77032961, 9.848857802,
8.94427191 , 8.062257748, 7.211102551, 6.403124237,
5.656854249, 5 , 4.472135955, 4.123105626,
4 , 4.123105626, 4.472135955, 5 ,
5.656854249, 6.403124237, 7.211102551, 8.062257748,
8.94427191 , 9.848857802, 10.77032961}, {11.18033989,
10.29563014, 9.433981132, 8.602325267, 7.810249676,
7.071067812, 6.403124237, 5.830951895, 5.385164807,
5.099019514, 5 , 5.099019514, 5.385164807,
5.830951895, 6.403124237, 7.071067812, 7.810249676,
8.602325267, 9.433981132, 10.29563014, 11.18033989},
{11.66190379, 10.81665383, 10 , 9.219544457,
8.485281374, 7.810249676, 7.211102551, 6.708203932,
6.32455532 , 6.08276253, 6 , 6.08276253,
6.32455532 , 6.708203932, 7.211102551, 7.810249676,
8.485281374, 9.219544457, 10 , 10.81665383,
11.66190379}, {12.20655562, 11.40175425, 10.63014581,
9.899494937, 9.219544457, 8.602325267, 8.062257748,
7.615773106, 7.280109889, 7.071067812, 7 ,
7.071067812, 7.280109889, 7.615773106, 8.062257748,
8.602325267, 9.219544457, 9.899494937, 10.63014581,
11.40175425, 12.20655562}, {12.80624847, 12.04159458,
11.3137085 , 10.63014581, 10 , 9.433981132,
8.94427191 , 8.544003745, 8.246211251, 8.062257748,
8 , 8.062257748, 8.246211251, 8.544003745,
8.94427191 , 9.433981132, 10 , 10.63014581,
11.3137085 , 12.04159458, 12.80624847}, {13.45362405,
12.72792206, 12.04159458, 11.40175425, 10.81665383,
10.29563014, 9.848857802, 9.486832981, 9.219544457,
9.055385138, 9 , 9.055385138, 9.219544457,
9.486832981, 9.848857802, 10.29563014, 10.81665383,
11.40175425, 12.04159458, 12.72792206, 13.45362405},
{14.14213562, 13.45362405, 12.80624847, 12.20655562,
11.66190379, 11.18033989, 10.77032961, 10.44030651,
10.19803903, 10.04987562, 10 , 10.04987562,
10.19803903, 10.44030651, 10.77032961, 11.18033989,
11.66190379, 12.20655562, 12.80624847, 13.45362405,
14.14213562}};
/* initialize the landscape */
void init()
{
int i,j,a,b;
for(i = 0;i < 2*RANGE+1;i++){
for(j = 0;j < 2*RANGE+1;j++){
scale[i][j] = 1 - (log(rmat[10-RANGE+i][10-RANGE+j])
- log(RANGE))/(log(0.5)-log(RANGE));
if(scale[i][j] < 0.)scale[i][j] = 0.;
}
}
for(i = 0;i < NMax;i++){
for(j = 0;j < NMax;j++){
Scape[i][j] = 1.;
}
}
}
void initwalker(Point *pt,int I,int J)
{
pt[npt].i = I;
pt[npt].j = J;
Scape[I][J] = OCCUPY;
npt++;
}
void main()
{
int k,ind,res,l,m,n,tmpi,tmpj,linei,tmpn,i,j,h;
float tmph,rg,x,y,xx,xy,a,b;
int i1,i2,j1,j2;
Point tmppt;
int fd,fdi,fdf,fdd,fdn,QUIT,maxnpt,nsites,inum;
char tmp[500];
float num[NMax*50][2];
fdd = open("sc.m",O_RDWR|O_CREAT|O_TRUNC,S_IRWXU);
sprintf(tmp,"function scape\n clear all\n cla\n hold on\n
x = 0:.1:10;\n y = 0:.1:10;\n global sc;\n sc =[");
write(fdd,tmp,strlen(tmp));
fdn = open("dim.m",O_RDWR|O_CREAT|O_TRUNC,S_IRWXU);
sprintf(tmp,"function scape\n clear all\n cla\n hold on\n
x = 1:1:50;\n global num\n num = [");
write(fdn,tmp,strlen(tmp));
fdl = open("list",O_RDWR|O_CREAT|O_TRUNC,S_IRWXU);
/* set starting walker */
init(); /* initialize the landscape*/
tmpn = 0;
h = (NMax-1)/2;
initwalker(pt,h,h);
cno = 4;
sprintf(tmp,"%d %d\n",h,h);
write(fdl,tmp,strlen(tmp));
ind = 0;
cno = 0;
QUIT = 0;
maxnpt = 1;
nsites = 1;
inum = 0;
while(1){
if(cno++%50000 == 0){
nsites = 0;
rg = 0.;
i1 = j1 = NMax;
i2 = j2 = 0;
for(i = 0;i < NMax;i++){
for(j = 0;j < NMax;j++){
if(Scape[i][j] != OCCUPY)continue;
nsites++;
if(rg < (1.*(i-h)*(i-h)+1.*(j-h)*(j-h)))
rg = (1.*(i-h)*(i-h)+1.*(j-h)*(j-h));
}
}
if(rg > (NMax)*(NMax)*.2){printf("exiting ...\n");break;}
}
ind = (int)((npt-1)*drand48());
res = step(ind);
if(npt == 0){printf("end npt %d\n",npt);break;}
}
for(i = 0;i < NMax;i++){
k = 10;
for(j = 0;j < NMax;j++){
tmph = Scape[i][j];
if(tmph == OCCUPY)tmph = 1.;
else tmph = 0.;
sprintf(tmp,"%.1f ",tmph);
write(fdd,tmp,strlen(tmp));
if(k == 0){
sprintf(tmp," ...\n");
write(fdd,tmp,strlen(tmp));
k = 10;
}
k--;
}
sprintf(tmp,";\n");
write(fdd,tmp,strlen(tmp));
}
sprintf(tmp,"];\n mesh(sc);\n axis equal;");
write(fdd,tmp,strlen(tmp));
close(fdd);
nsites = 0;
for(i = 0;i < NMax;i++){
for(j = 0;j < NMax;j++){
if(Scape[i][j] != OCCUPY)continue;
nsites++;
rg += (1.*(i-h)*(i-h)+1.*(j-h)*(j-h));
}
}
rg /= (nsites*1.);
rg = sqrt(rg);
for(k = 1;k <= (int)rg;k++){
nsites = 0;
for(i = h-k;i <= h+k ;i++){
for(j = h-k;j <= h+k;j++){
if(Scape[i][j] == OCCUPY)nsites++;
}
}
num[k][0] = log(2.*k+1.);
num[k][1] = log(nsites*1.);
}
inum = k-1;
printf("inum = %d...\n",inum);
k = 10;
x = y = xx = xy = 0.;
for(i = 1;i < inum;i++){
x += num[i][0];
y += num[i][1];
xx += num[i][0]*num[i][0];
xy += num[i][0]*num[i][1];
sprintf(tmp,"%.3f %.3f;",num[i][0],num[i][1]);
write(fdn,tmp,strlen(tmp));
if(k == 0){
sprintf(tmp," ...\n");
write(fdn,tmp,strlen(tmp));
k = 10;
}
k--;
}
sprintf(tmp, "];\nplot(num(:,1),num(:,2));");
write(fdn,tmp,strlen(tmp));
inum--;
a = (x*y-inum*xy)/(x*x-inum*xx);
b = (x*xy-xx*y)/(x*x-inum*xx);
sprintf(tmp,"\n legend('log(Area) = %.3f*log(r)+%.3f',4);\n",a,b);
write(fdn,tmp,strlen(tmp));
close(fdn);
printf("Fractal Dim %.3f\n",a);
close(fdl);
}