/*7:*/
#line 758 "./boxcount.w"

/*8:*/
#line 789 "./boxcount.w"

#include <math.h> 
#include <time.h> 
#include <stdio.h> 
#include <stdlib.h> 
#include <string.h> 
#include <ctype.h> 

/*:8*/
#line 759 "./boxcount.w"

/*9:*/
#line 810 "./boxcount.w"

#define VERSION "1.6"
#define COPYRIGHT "Copyright (C) 2006-2011, Fredrik Jonsson"
#define SUCCESS (0)
#define FAILURE (1)
#define NCHMAX (256)
#define APPEND 1
#define OVERWRITE 2

/*:9*/
#line 760 "./boxcount.w"

/*10:*/
#line 825 "./boxcount.w"

extern char*optarg;
char*progname;

/*:10*/
#line 761 "./boxcount.w"

/*23:*/
#line 1306 "./boxcount.w"

/*24:*/
#line 1320 "./boxcount.w"

void moment(double data[],int nnmin,int nnmax,double*ave,double*adev,
double*sdev,double*var,double*skew,double*curt){
int j;
double ep= 0.0,s,p;
if(nnmax-nnmin<=1){
fprintf(stderr,"%s: Error in routine moment()! ",progname);
fprintf(stderr,"(nnmax-nnmin>1 is a requirement)\n");
exit(FAILURE);
}
s= 0.0;
for(j= nnmin;j<=nnmax;j++)s+= data[j];
*ave= s/(nnmax-nnmin+1);
*adev= (*var)= (*skew)= (*curt)= 0.0;
for(j= nnmin;j<=nnmax;j++){
*adev+= fabs(s= data[j]-(*ave));
ep+= s;
ep+= s;
*var+= (p= s*s);
*skew+= (p*= s);
*curt+= (p*= s);
}
*adev/= (nnmax-nnmin+1);
*var= (*var-ep*ep/(nnmax-nnmin+1))/(nnmax-nnmin);
*sdev= sqrt(*var);
if(*var){
*skew/= ((nnmax-nnmin+1)*(*var)*(*sdev));
*curt= (*curt)/((nnmax-nnmin+1)*(*var)*(*var))-3.0;
}else{
fprintf(stderr,"%s: Error in routine moment()! ",progname);
fprintf(stderr,"No skew/kurtosis for zero variance\n");
exit(FAILURE);
}
}

/*:24*/
#line 1307 "./boxcount.w"

/*25:*/
#line 1361 "./boxcount.w"

long int num_coordinate_pairs(FILE*trajectory_file){
double tmp;
int tmpch;
long int mm= 0;
fseek(trajectory_file,0L,SEEK_SET);
while((tmpch= getc(trajectory_file))!=EOF){
ungetc(tmpch,trajectory_file);
fscanf(trajectory_file,"%lf",&tmp);
fscanf(trajectory_file,"%lf",&tmp);
mm++;
tmpch= getc(trajectory_file);
while((tmpch!=EOF)&&(!isdigit(tmpch)))tmpch= getc(trajectory_file);
if(tmpch!=EOF)ungetc(tmpch,trajectory_file);
}
fseek(trajectory_file,0L,SEEK_SET);
return(mm);
}


/*:25*/
#line 1308 "./boxcount.w"

/*26:*/
#line 1390 "./boxcount.w"

/*27:*/
#line 1399 "./boxcount.w"

short pathcharacter(int ch){
return(isalnum(ch)||(ch=='.')||(ch=='/')||(ch=='\\')||(ch=='_')||
(ch=='-')||(ch=='+'));
}

/*:27*/
#line 1391 "./boxcount.w"

/*28:*/
#line 1411 "./boxcount.w"

char*strip_away_path(char filename[]){
int j,k= 0;
while(pathcharacter(filename[k]))k++;
j= (--k);
while(isalnum((int)(filename[j])))j--;
j++;
return(&filename[j]);
}

/*:28*/
#line 1392 "./boxcount.w"


/*:26*/
#line 1309 "./boxcount.w"

/*29:*/
#line 1426 "./boxcount.w"

/*30:*/
#line 1437 "./boxcount.w"

double*dvector(long nl,long nh){
double*v;
v= (double*)malloc((size_t)((nh-nl+2)*sizeof(double)));
if(!v){
fprintf(stderr,"Error: Allocation failure in dvector()\n");
exit(FAILURE);
}
return v-nl+1;
}

/*:30*/
#line 1427 "./boxcount.w"

/*33:*/
#line 1473 "./boxcount.w"

void free_ivector(int*v,long nl,long nh){
free((char*)(v+nl-1));
}

/*:33*//*35:*/
#line 1495 "./boxcount.w"

void free_livector(long int*v,long nl,long nh){
free((char*)(v+nl-1));
}

/*:35*/
#line 1428 "./boxcount.w"

/*32:*/
#line 1459 "./boxcount.w"

int*ivector(long nl,long nh){
int*v;
v= (int*)malloc((size_t)((nh-nl+2)*sizeof(int)));
if(!v){
fprintf(stderr,"Error: Allocation failure in ivector()\n");
exit(FAILURE);
}
return v-nl+1;
}

/*:32*//*34:*/
#line 1481 "./boxcount.w"

long int*livector(long nl,long nh){
long int*v;
v= (long int*)malloc((size_t)((nh-nl+2)*sizeof(long int)));
if(!v){
fprintf(stderr,"Error: Allocation failure in livector()\n");
exit(FAILURE);
}
return v-nl+1;
}

/*:34*/
#line 1429 "./boxcount.w"

/*31:*/
#line 1451 "./boxcount.w"

void free_dvector(double*v,long nl,long nh){
free((char*)(v+nl-1));
}

/*:31*/
#line 1430 "./boxcount.w"

/*36:*/
#line 1504 "./boxcount.w"

short int**simatrix(long nrl,long nrh,long ncl,long nch)
{
long i,nrow= nrh-nrl+1,ncol= nch-ncl+1;
short int**m;
m= (short int**)malloc((size_t)((nrow+1)*sizeof(short int*)));
if(!m){
fprintf(stderr,"%s: Allocation failure 1 in simatrix()\n",progname);
exit(FAILURE);
}
m+= 1;
m-= nrl;
m[nrl]= (short int*)malloc((size_t)((nrow*ncol+1)*sizeof(short int)));
if(!m[nrl]){
fprintf(stderr,"%s: Allocation failure 2 in simatrix()\n",progname);
exit(FAILURE);
}
m[nrl]+= 1;
m[nrl]-= ncl;
for(i= nrl+1;i<=nrh;i++)m[i]= m[i-1]+ncol;
return m;
}

/*:36*/
#line 1431 "./boxcount.w"

/*37:*/
#line 1530 "./boxcount.w"

void free_simatrix(short int**m,long nrl,long nrh,long ncl,long nch){
free((char*)(m[nrl]+ncl-1));
free((char*)(m+nrl-1));
}

/*:37*/
#line 1432 "./boxcount.w"


/*:29*/
#line 1310 "./boxcount.w"

/*38:*/
#line 1538 "./boxcount.w"

void showsomehelp(void){
fprintf(stderr,"Usage: %s [options]\n",progname);
fprintf(stderr,"Options:\n");
fprintf(stderr," -h, --help\n");
fprintf(stderr,"    Display this help message and exit clean.\n");
fprintf(stderr," -v, --verbose\n");
fprintf(stderr,"    Toggle verbose mode on/off.\n");
fprintf(stderr," -o, --outputfile <str>\n");
fprintf(stderr,"    Specifies the base name of the output files where the program\n");
fprintf(stderr,"    is to save the calculated data. If the --outputfile or -o\n");
fprintf(stderr,"    option is followed by 'append' the estimate for the fractal\n");
fprintf(stderr,"    dimension will be appended to the file named <str>.dat, which\n");
fprintf(stderr,"    will be created if it does not exist. If the following word\n");
fprintf(stderr,"    instead is `overwrite' the file will instead be overwritten.\n");
fprintf(stderr," -i, --trajectoryfile\n");
fprintf(stderr,"    Specifies the input trajectory of the graph whose fractal\n");
fprintf(stderr,"    dimension is to be estimated.\n");
fprintf(stderr," -w, --computationwindow llx <num> lly <num> urx <num> ury <num>\n");
fprintf(stderr,"    This option explicitly specifies the domain over which the\n");
fprintf(stderr,"    box-counting algorithm will be applied, in terms of the\n");
fprintf(stderr,"    lower-left and upper-right corners (llx,lly) and (urx,ury),\n");
fprintf(stderr,"    respectively. By specifying this option, any automatic\n");
fprintf(stderr,"    calculation of computational window will be neglected.\n");
fprintf(stderr," -m, --boxmaps\n");
fprintf(stderr,"    If this option is present, the program will generate\n");
fprintf(stderr,"    MetaPost code for maps of the distribution of boxes.\n");
fprintf(stderr,"    In doing so, also the input trajectory is included in\n");
fprintf(stderr,"    the graphs. The convention for the naming of the output\n");
fprintf(stderr,"    map files is that they are saved as <str>.<N>.dat,\n");
fprintf(stderr,"    where <str> is the base filename as specified using the\n");
fprintf(stderr,"    -o or --outputfile option, <N> is the automatically appended\n");
fprintf(stderr,"    current level of resolution refinement in the box-counting,\n");
fprintf(stderr,"    and where '.dat' is the file suffix as automatically appended\n");
fprintf(stderr,"    by the program.\n");
fprintf(stderr," --graphsize <width in mm> <height in mm>\n");
fprintf(stderr,"    If the -m or --boxmaps option is present at the command line,\n");
fprintf(stderr,"    then the --graphsize option will override the default graph\n");
fprintf(stderr,"    size of the generated boxmaps. (Default graph size is 80 mm\n");
fprintf(stderr,"    width and 56 mm height.)\n");
fprintf(stderr," -Nmin, --minlevel <num>\n");
fprintf(stderr,"    Specifies the minimum level Nmin of grid refinement that \n");
fprintf(stderr,"    will be used in the evaluation of the estimate of the fractal\n");
fprintf(stderr,"    dimension. With this option specified, the coarsest level\n");
fprintf(stderr,"    used in the box-counting will be a [(2^Nmin)x(2^Nmin)]-grid\n");
fprintf(stderr,"    of boxes.\n");
fprintf(stderr," -Nmax, --maxlevel <num>\n");
fprintf(stderr,"    This option is similar to the --minlevel or -Nmin options,\n");
fprintf(stderr,"    with the difference that it instead specifies the maximum\n");
fprintf(stderr,"    level Nmax of grid refinement that will be used in the\n");
fprintf(stderr,"    evaluation of the estimate of the fractal dimension.\n");
}

/*:38*/
#line 1311 "./boxcount.w"

/*39:*/
#line 1647 "./boxcount.w"

short int lines_intersect(double p1x,double p1y,double q1x,double q1y,
double p2x,double p2y,double q2x,double q2y){
double a,b,c,d,e,f,det,tmp1,tmp2;
short int intersect;
a= q1x-p1x;
b= p2x-q2x;
c= q1y-p1y;
d= p2y-q2y;
e= p2x-p1x;
f= p2y-p1y;
det= a*d-b*c;
tmp1= e*d-b*f;
tmp2= a*f-e*c;
intersect= 0;
if(det> 0){
if(((0.0<=tmp1)&&(tmp1<=det))&&((0.0<=tmp2)&&(tmp2<=det)))intersect= 1;
}else if(det<0){
if(((det<=tmp1)&&(tmp1<=0.0))&&((det<=tmp2)&&(tmp2<=0.0)))intersect= 1;
}
return(intersect);
}

/*:39*/
#line 1312 "./boxcount.w"

/*40:*/
#line 1692 "./boxcount.w"

short int box_intersection(double px,double py,double qx,double qy,
double llx,double lly,double urx,double ury){
if(lines_intersect(px,py,qx,qy,llx,lly,urx,lly))
return(1);
if(lines_intersect(px,py,qx,qy,urx,lly,urx,ury))
return(1);
if(lines_intersect(px,py,qx,qy,urx,ury,llx,ury))
return(1);
if(lines_intersect(px,py,qx,qy,llx,ury,llx,lly))
return(1);
return(0);
}

/*:40*/
#line 1313 "./boxcount.w"

/*41:*/
#line 1714 "./boxcount.w"

/*42:*/
#line 1762 "./boxcount.w"

long unsigned int get_num_covering_boxes_with_boxmaps(
double*x_traj,double*y_traj,
long int mm,int resolution,double global_llx,double global_lly,
double global_urx,double global_ury,char boxgraph_filename[],
double boxgraph_width,double boxgraph_height,
char trajectory_filename[]){
short int**box;
long unsigned int i,j,m,nn,istart,jstart,istop,jstop,iincr,jincr,num_boxes;
double*x_box,*y_box;

double px,py,qx,qy;
FILE*boxgraph_file;

/*43:*/
#line 1805 "./boxcount.w"

{
nn= 1;
for(i= 1;i<=resolution;i++)nn= 2*nn;
box= simatrix(1,nn,1,nn);
for(i= 1;i<=nn;i++)for(j= 1;j<=nn;j++)box[i][j]= 0;
x_box= dvector(1,nn+1);
y_box= dvector(1,nn+1);
for(m= 1;m<=nn+1;m++){
x_box[m]= global_llx+((double)(m-1))*(global_urx-global_llx)/((double)(nn));
y_box[m]= global_lly+((double)(m-1))*(global_ury-global_lly)/((double)(nn));
}
}

/*:43*/
#line 1776 "./boxcount.w"

/*44:*/
#line 1822 "./boxcount.w"

{
istart= 0;
jstart= 0;
for(m= 1;m<=nn;m++){
if((x_box[m]<=x_traj[1])&&(x_traj[1]<=x_box[m+1]))istart= m;
if((y_box[m]<=y_traj[1])&&(y_traj[1]<=y_box[m+1]))jstart= m;
}
if((istart==0)||(jstart==0)){
fprintf(stderr,"%s: Error! Cannot find box indices of 1st coordinate!\n",
progname);
fprintf(stderr,"%s: Please check data for input trajetory.\n",progname);
exit(FAILURE);
}
}


/*:44*/
#line 1777 "./boxcount.w"

for(m= 1;m<=mm-1;m++){
/*45:*/
#line 1842 "./boxcount.w"
{
px= x_traj[m];
py= y_traj[m];
qx= x_traj[m+1];
qy= y_traj[m+1];
istop= istart;
jstop= jstart;
if(px<qx){
iincr= 1;
while(x_box[istop+1]<qx)istop++;
}else{
iincr= -1;
while(x_box[istop]> qx)istop--;
}
if(py<qy){
jincr= 1;
while(y_box[jstop+1]<qy)jstop++;
}else{
jincr= -1;
while(y_box[jstop]> qy)jstop--;
}
if(0==1){
fprintf(stdout,"%s: Endpoint box indices: (i=%ld,j=%ld)\n",
progname,istop,jstop);
}
}

/*:45*/
#line 1779 "./boxcount.w"

/*46:*/
#line 1878 "./boxcount.w"

{
for(i= istart;i!=(istop+iincr);i+= iincr){
for(j= jstart;j!=(jstop+jincr);j+= jincr){
if(box_intersection(px,py,qx,qy,
x_box[i],y_box[j],x_box[i+1],y_box[j+1])){
box[i][j]= 1;
}
}
}
istart= istop;
jstart= jstop;
}

/*:46*/
#line 1780 "./boxcount.w"

}
/*47:*/
#line 1898 "./boxcount.w"
{
if(!strcmp(boxgraph_filename,"")){
boxgraph_file= NULL;
}else{
if((boxgraph_file= fopen(boxgraph_filename,"w"))==NULL){
fprintf(stderr,"%s: Could not open %s for box graphs!\n",
progname,boxgraph_filename);
exit(FAILURE);
}
fseek(boxgraph_file,0L,SEEK_SET);
fprintf(boxgraph_file,"input graph;\n");
fprintf(boxgraph_file,"def box(expr i,j)=\n");
fprintf(boxgraph_file,"begingroup\n");
fprintf(boxgraph_file,"llx:=%2.8f+(i-1)*%2.8f;\n",
global_llx,(global_urx-global_llx)/(nn));
fprintf(boxgraph_file,"lly:=%2.8f+(j-1)*%2.8f;\n",
global_lly,(global_ury-global_lly)/(nn));
fprintf(boxgraph_file,"urx:=%2.8f+(i)*%2.8f;\n",
global_llx,(global_urx-global_llx)/(nn));
fprintf(boxgraph_file,"ury:=%2.8f+(j)*%2.8f;\n",
global_lly,(global_ury-global_lly)/(nn));
fprintf(boxgraph_file,"gdraw (llx,lly)--(urx,lly);\n");
fprintf(boxgraph_file,"gdraw (urx,lly)--(urx,ury);\n");
fprintf(boxgraph_file,"gdraw (urx,ury)--(llx,ury);\n");
fprintf(boxgraph_file,"gdraw (llx,ury)--(llx,lly);\n");
fprintf(boxgraph_file,"endgroup\n");
fprintf(boxgraph_file,"enddef;\n");
fprintf(boxgraph_file,"beginfig(1);\n");
fprintf(boxgraph_file,"w:=%2.4fmm; h:=%2.4fmm;\n",
boxgraph_width,boxgraph_height);
fprintf(boxgraph_file,"draw begingraph(w,h);\n");
fprintf(boxgraph_file,"pickup pencircle scaled .3pt;\n");
fprintf(boxgraph_file,"setrange(%2.6f,%2.6f,%2.6f,%2.6f);\n",
global_llx,global_lly,global_urx,global_ury);
}
}

/*:47*/
#line 1782 "./boxcount.w"

/*48:*/
#line 1964 "./boxcount.w"
{
if(boxgraph_file!=NULL){
fprintf(boxgraph_file,"pickup pencircle scaled .5pt;\n");
i= 0;
for(m= 1;m<=mm;m++){
if(i==0){
if(m==1){
fprintf(boxgraph_file,"gdraw (%2.4f,%2.4f)",
x_traj[m],y_traj[m]);
}else if(2<mm){
fprintf(boxgraph_file,"gdraw (%2.4f,%2.4f)--(%2.4f,%2.4f)",
x_traj[m-1],y_traj[m-1],x_traj[m],y_traj[m]);
}
}else{
fprintf(boxgraph_file,"--(%2.4f,%2.4f)",x_traj[m],y_traj[m]);
}
i++;
if((i==5)||(m==mm)){
fprintf(boxgraph_file,";\n");
i= 0;
}
}
}
}

/*:48*/
#line 1783 "./boxcount.w"

/*49:*/
#line 1999 "./boxcount.w"

{
num_boxes= 0;
for(i= 1;i<=nn;i++){
for(j= 1;j<=nn;j++){
if(box[i][j]==1){
num_boxes++;
if(boxgraph_file!=NULL){
fprintf(boxgraph_file,"box(%ld,%ld);\n",i,j);
}
}
}
}
free_simatrix(box,1,nn,1,nn);
}

/*:49*/
#line 1784 "./boxcount.w"

/*50:*/
#line 2020 "./boxcount.w"
{
if(boxgraph_file!=NULL){
fprintf(boxgraph_file,"autogrid(itick bot,itick lft);\n");
fprintf(boxgraph_file,"glabel.bot(btex $x$ etex,OUT);\n");
fprintf(boxgraph_file,"glabel.lft(btex $y$ etex,OUT);\n");
fprintf(boxgraph_file,"endgraph;\n");
fprintf(boxgraph_file,"endfig;\n");
fprintf(boxgraph_file,"end\n");
}
}

/*:50*/
#line 1785 "./boxcount.w"

/*51:*/
#line 2033 "./boxcount.w"
{
if(boxgraph_file!=NULL){
fprintf(stdout,
"%s: MetaPost output box distribution graph written to %s.\n",
progname,boxgraph_filename);
fclose(boxgraph_file);
}
}

/*:51*/
#line 1786 "./boxcount.w"

return(num_boxes);
}

/*:42*/
#line 1715 "./boxcount.w"

/*52:*/
#line 2078 "./boxcount.w"

long unsigned int get_num_covering_boxes(double*x_traj,double*y_traj,long int mm,
int i,double global_llx,double global_lly,
double global_urx,double global_ury){
return(get_num_covering_boxes_with_boxmaps(x_traj,y_traj,mm,i,
global_llx,global_lly,global_urx,global_ury,"",0.0,0.0,""));
}

/*:52*/
#line 1716 "./boxcount.w"


/*:41*/
#line 1314 "./boxcount.w"



/*:23*/
#line 762 "./boxcount.w"


int main(int argc,char*argv[])
{
/*11:*/
#line 904 "./boxcount.w"

short verbose,user_set_compwin,output_mode,make_boxmaps;
long int*num_boxes;
time_t initime;
time_t now;
long mm,nn,nnmin,nnmax;
double global_llx,global_lly,global_urx,global_ury;
double*x_traj,*y_traj,*x,*y;
FILE*trajectory_file,*frac_estimate_file,*boxmap_file;
char trajectory_filename[NCHMAX],output_filename[NCHMAX],frac_estimate_filename[NCHMAX];
char boxmap_filename[NCHMAX];
double boxgraph_width,boxgraph_height;
int no_arg;
long i;
double*fracdimen_estimates,ave,adev,sdev,var,skew,curt;

/*:11*/
#line 766 "./boxcount.w"

/*12:*/
#line 922 "./boxcount.w"

verbose= 0;
user_set_compwin= 0;
output_mode= OVERWRITE;
make_boxmaps= 0;
nnmin= 0;
nnmax= 10;
strcpy(output_filename,"out");
strcpy(trajectory_filename,"");
boxgraph_width= 80.0;
boxgraph_height= 56.0;
trajectory_file= NULL;
frac_estimate_file= NULL;
boxmap_file= NULL;
initime= time(NULL);

/*:12*/
#line 767 "./boxcount.w"

/*13:*/
#line 944 "./boxcount.w"

{
progname= strip_away_path(argv[0]);
fprintf(stdout,"This is %s v.%s. %s\n",progname,VERSION,COPYRIGHT);
no_arg= argc;
while(--argc){
if(!strcmp(argv[no_arg-argc],"-o")||
!strcmp(argv[no_arg-argc],"--outputfile")){
--argc;
if(!strcmp(argv[no_arg-argc],"append")||
!strcmp(argv[no_arg-argc],"a")){
output_mode= APPEND;
}else if(!strcmp(argv[no_arg-argc],"overwrite")||
!strcmp(argv[no_arg-argc],"o")){
output_mode= OVERWRITE;
}else{
fprintf(stderr,"%s: Error in '-o' or '--outputfile' option!",
progname);
exit(FAILURE);
}
--argc;
strcpy(output_filename,argv[no_arg-argc]);
}else if(!strcmp(argv[no_arg-argc],"-w")||
!strcmp(argv[no_arg-argc],"--computationwindow")){
user_set_compwin= 1;
--argc;
if(!strcmp(argv[no_arg-argc],"llx")){
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&global_llx)){
fprintf(stderr,"%s: Error in parsing lower-left x-value.\n",
progname);
exit(FAILURE);
}
}else{
fprintf(stderr,"%s: Error in computation window option\n",
progname);
fprintf(stderr,"%s: Expecting 'llx'\n",progname);
exit(FAILURE);
}
--argc;
if(!strcmp(argv[no_arg-argc],"lly")){
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&global_lly)){
fprintf(stderr,"%s: Error in parsing lower-left y-value.\n",
progname);
exit(FAILURE);
}
}else{
fprintf(stderr,"%s: Error in computation window option\n",
progname);
fprintf(stderr,"%s: Expecting 'lly'\n",progname);
exit(FAILURE);
}
--argc;
if(!strcmp(argv[no_arg-argc],"urx")){
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&global_urx)){
fprintf(stderr,"%s: Error in parsing lower-left x-value.\n",
progname);
exit(FAILURE);
}
}else{
fprintf(stderr,"%s: Error in computation window option\n",
progname);
fprintf(stderr,"%s: Expecting 'urx'\n",progname);
exit(FAILURE);
}
--argc;
if(!strcmp(argv[no_arg-argc],"ury")){
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&global_ury)){
fprintf(stderr,"%s: Error in parsing lower-left y-value.\n",
progname);
exit(FAILURE);
}
}else{
fprintf(stderr,"%s: Error in computation window option\n",
progname);
fprintf(stderr,"%s: Expecting 'ury'\n",progname);
exit(FAILURE);
}
}else if(!strcmp(argv[no_arg-argc],"-i")||
!strcmp(argv[no_arg-argc],"--trajectoryfile")){
--argc;
strcpy(trajectory_filename,argv[no_arg-argc]);
}else if(!strcmp(argv[no_arg-argc],"-v")||
!strcmp(argv[no_arg-argc],"--verbose")){
verbose= (verbose?0:1);
}else if(!strcmp(argv[no_arg-argc],"-h")||
!strcmp(argv[no_arg-argc],"--help")){
showsomehelp();
exit(SUCCESS);
}else if(!strcmp(argv[no_arg-argc],"-m")||
!strcmp(argv[no_arg-argc],"--boxmaps")){
make_boxmaps= 1;
}else if(!strcmp(argv[no_arg-argc],"--graphsize")){
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&boxgraph_width)){
fprintf(stderr,"%s: Error in width of '--graphsize' option.\n",
progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&boxgraph_height)){
fprintf(stderr,"%s: Error in height of '--graphsize' option.\n",
progname);
exit(FAILURE);
}
}else if(!strcmp(argv[no_arg-argc],"--minlevel")||
!strcmp(argv[no_arg-argc],"-Nmin")){
--argc;
if(!sscanf(argv[no_arg-argc],"%ld",&nnmin)){
fprintf(stderr,"%s: Error in '--minlevel' or '-Nmin' option.\n",
progname);
exit(FAILURE);
}
}else if(!strcmp(argv[no_arg-argc],"--maxlevel")||
!strcmp(argv[no_arg-argc],"-Nmax")){
--argc;
if(!sscanf(argv[no_arg-argc],"%ld",&nnmax)){
fprintf(stderr,"%s: Error in '--maxlevel' or '-Nmax' option.\n",
progname);
exit(FAILURE);
}
}else{
fprintf(stderr,"%s: Specified option '%s' invalid!\n",
progname,argv[no_arg-argc]);
showsomehelp();
exit(FAILURE);
}
}
}

/*:13*/
#line 768 "./boxcount.w"

/*14:*/
#line 1079 "./boxcount.w"

{
fprintf(stdout,"%s: Program execution started %s",progname,ctime(&initime));
}

/*:14*/
#line 769 "./boxcount.w"

/*15:*/
#line 1089 "./boxcount.w"
{
if(!strcmp(trajectory_filename,"")){
fprintf(stderr,"%s: No input trajectory specified!\n",progname);
fprintf(stderr,"%s: Please use the '--trajectoryfile' option.\n",
progname);
fprintf(stderr,"%s: Use '--help' option to display a help message.\n",
progname);
exit(FAILURE);
}
if((trajectory_file= fopen(trajectory_filename,"r"))==NULL){
fprintf(stderr,"%s: Could not open %s for loading trajectory!\n",
progname,trajectory_filename);
exit(FAILURE);
}
mm= num_coordinate_pairs(trajectory_file);
x_traj= dvector(1,mm);
y_traj= dvector(1,mm);
for(i= 1;i<=mm;i++){
fscanf(trajectory_file,"%lf",&x_traj[i]);
fscanf(trajectory_file,"%lf",&y_traj[i]);
}
fclose(trajectory_file);
if(verbose){
fprintf(stdout,"%s: Loaded %ld coordinate pairs from file '%s'.\n",
progname,mm,trajectory_filename);
}
}

/*:15*/
#line 770 "./boxcount.w"

/*16:*/
#line 1119 "./boxcount.w"
{
if(!strcmp(output_filename,"")){
fprintf(stderr,"%s: No output base name specified!\n",progname);
fprintf(stderr,"%s: Please use the '--outputfile' option.\n",
progname);
exit(FAILURE);
}
sprintf(frac_estimate_filename,"%s.dat",output_filename);
if(output_mode==APPEND){
if((frac_estimate_file= fopen(frac_estimate_filename,"a"))==NULL){
fprintf(stderr,"%s: Could not open %s for output!\n",
progname,frac_estimate_filename);
exit(FAILURE);
}
fseek(frac_estimate_file,0L,SEEK_END);

}else if(output_mode==OVERWRITE){
if((frac_estimate_file= fopen(frac_estimate_filename,"w"))==NULL){
fprintf(stderr,"%s: Could not open %s for loading trajectory!\n",
progname,frac_estimate_filename);
exit(FAILURE);
}
fseek(frac_estimate_file,0L,SEEK_SET);

}else{
fprintf(stderr,"%s: Error: Output mode (output_mode) undefined!",progname);
exit(FAILURE);
}
}

/*:16*/
#line 771 "./boxcount.w"

/*17:*/
#line 1158 "./boxcount.w"

{
if(!user_set_compwin){
global_llx= x_traj[1];
global_lly= y_traj[1];
global_urx= global_llx;
global_ury= global_lly;
for(i= 1;i<=mm;i++){
if(x_traj[i]> global_urx)global_urx= x_traj[i];
if(y_traj[i]> global_ury)global_ury= y_traj[i];
if(x_traj[i]<global_llx)global_llx= x_traj[i];
if(y_traj[i]<global_lly)global_lly= y_traj[i];
}
if(verbose){
fprintf(stdout,"%s: Global box-counting window:\n",progname);
fprintf(stdout,"%s:   (llx,lly)=(%2.8f,%2.8f)\n",progname,
global_llx,global_lly);
fprintf(stdout,"%s:   (urx,ury)=(%2.8f,%2.8f)\n",progname,
global_urx,global_ury);
}
}
}

/*:17*/
#line 772 "./boxcount.w"

/*18:*/
#line 1195 "./boxcount.w"
{
num_boxes= livector(1,nnmax);
nn= 1;
for(i= 1;i<=nnmin-1;i++)nn= 2*nn;
for(i= nnmin;i<=nnmax;i++){
nn= 2*nn;
if(make_boxmaps){

sprintf(boxmap_filename,"%s-%02ld.mp",output_filename,i);
num_boxes[i]= get_num_covering_boxes_with_boxmaps(x_traj,y_traj,mm,i,
global_llx,global_lly,global_urx,global_ury,boxmap_filename,
boxgraph_width,boxgraph_height,trajectory_filename);
}else{

num_boxes[i]= get_num_covering_boxes(x_traj,y_traj,mm,i,
global_llx,global_lly,global_urx,global_ury);
}
if(verbose){
fprintf(stdout,"%s: N=%ld (%ldx%ld-grid of boxes): ",progname,i,nn,nn);
fprintf(stdout,"Trajectory covered by %ld boxes\n",num_boxes[i]);
}
}
}

/*:18*/
#line 773 "./boxcount.w"

/*19:*/
#line 1245 "./boxcount.w"
{
x= dvector(nnmin,nnmax);
y= dvector(nnmin,nnmax);
fracdimen_estimates= dvector(nnmin,nnmax);
nn= 1;
for(i= 1;i<=nnmax;i++){
nn= 2*nn;
if(nnmin<=i)x[i]= log((double)nn);
}
for(i= nnmin;i<=nnmax;i++)y[i]= log(num_boxes[i]);
for(i= nnmin;i<=nnmax;i++)fracdimen_estimates[i]= y[i]/x[i];
if(verbose){
for(i= nnmin;i<=nnmax;i++){
fprintf(stdout,"%s: N=%ld, fractal dimension estimate=%f\n",progname,
i,fracdimen_estimates[i]);
}
}
moment(fracdimen_estimates,nnmin,nnmax,&ave,&adev,&sdev,&var,&skew,&curt);
if(verbose){
fprintf(stdout,"%s: Estimate of fractal dimension: %f\n",progname,ave);
fprintf(stdout,"%s: Average deviation: %f\n",progname,adev);
fprintf(stdout,"%s: Standard deviation: %f\n",progname,sdev);
fprintf(stdout,"%s: Skewness: %f\n",progname,skew);
}
free_livector(num_boxes,1,nnmax);
free_dvector(fracdimen_estimates,nnmin,nnmax);
free_dvector(x,nnmin,nnmax);
free_dvector(y,nnmin,nnmax);
}

/*:19*/
#line 774 "./boxcount.w"

/*20:*/
#line 1277 "./boxcount.w"
{
if(output_mode==APPEND){
fseek(frac_estimate_file,0L,SEEK_END);
}else if(output_mode==OVERWRITE){
fseek(frac_estimate_file,0L,SEEK_SET);
}
fprintf(frac_estimate_file,"%f %f %f\n",ave,sdev,skew);
}


/*:20*/
#line 775 "./boxcount.w"

/*21:*/
#line 1289 "./boxcount.w"
{
fclose(frac_estimate_file);
}

/*:21*/
#line 776 "./boxcount.w"

/*22:*/
#line 1295 "./boxcount.w"
{
now= time(NULL);
if(verbose)
fprintf(stdout,"%s: Total execution time: %d s\n",progname,
((int)difftime(now,initime)));
fprintf(stdout,"%s: Program execution closed %s",progname,ctime(&now));
}

/*:22*/
#line 777 "./boxcount.w"

return(SUCCESS);
}

/*:7*/