#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <float.h>
#define CIRLEN 131072    //length of circle data arrays
   // global definitions
   const float pi=3.14159265;   //constant pi  
   float htpix,wdpix;           //picture ht and wd
   FILE *ab;     //file pointer for .ps
   FILE *cd;     //file pointer for text file
  	
	// function prototypes
   float rd();                       //returns random numbers from 0 to 1
   void placecirc_dem();             //fractal circle algorithm
   double zeta_hurw(float u,int N);  //Hurwitz zeta function
   double est_sum(float c,float n);  //auxiliary needed for Hurwitz zeta

    int main()  
   {
      float qqx,qqy;
   	
      char a1[8]="%!PS";
      char fn1[32]="circ_image.ps";         //Postscript data file name
      char fn2[32]="circ_data.txt";         //text data file name
   	
   	//====================== setup ============================
      srand(time(NULL));   //seeds the random number sequence based on clock time
   	
      htpix=11.;    //picture height and width
      wdpix=8.5;    //edit these for different aspect ratios
   	
      qqx=8.5/wdpix;
      qqy=11./htpix;
      ab=fopen(fn1,"w");    //open the file
      cd=fopen(fn2,"w");    //open the file
   		
      fprintf(ab,"%s\n",a1);     //header data for generic postscript
      fprintf(ab,"/in {72 mul} def\n");
      fprintf(ab,"%f %f scale\n",qqx,qqy);    	
      	
      // picture code starts ---------------------
         
      placecirc_dem();  
   	   
      // picture code ends -----------------------   
      	      			
      fprintf(ab,"showpage\n");    //required by Postscript
      fclose(ab);   //close it
      fclose(cd);   //close it
   	
      return 0;
   	 
   }	// end of main 
	
	// ---- start function definitions ----
	
    float rd()              //function creates random numbers
   {
      double rr;
      float x;
      rr=pow(2,31);         //for putting rand() into range 0-1
      x=(float)rand()/rr;   //random numbers distributed over interval 0-1
      return x;
   }
	
    void placecirc_dem()   
   {
   	//revised Oct 2010 for straight power law  
   	//revised June 2011 for free distribution   
      float *xc,*yc,*rc;
      long *nt;
      xc=malloc(sizeof(*xc)*CIRLEN);
      yc=malloc(sizeof(*yc)*CIRLEN);
      rc=malloc(sizeof(*rc)*CIRLEN);
      nt=malloc(sizeof(*nt)*CIRLEN);
   
      float rtst,rcir,rcir0,x1,y1,dx,dy;
      float fillmax,zval,acirc,gmean;
      float aratio,nNval,ffac;    
      float cval,exu,q1,q2,q3,q4,totarea,rsum,u,afill,ara;
      int Nval,tst1,k,m,j,nprint;
      long ncir,nitstot,itsmax,nmax,circmax,nits;
   	
      //----- setup data ----- edit as needed
      cval=1.28;            //exponent of power law for area
      Nval=2;               //first integer to be used in the sequence
      circmax=90000;        //number of circles placed (quits) -- can't exceed CIRLEN
      itsmax=400000;        //maximum iterations in a cycle (quits)
      nprint=5;             //prints data every nprint cycles
      fillmax=.99;          //maximum fill (quits)
   
      nmax=circmax+1;        
      nNval=Nval;          //integer conversion
      exu=.5*cval;         //exponent in power law for r
   	
   	//find area ratio for given cval
      zval=zeta_hurw(cval,nNval);
      acirc=htpix*wdpix/zval;
      aratio=1./zval;
      printf("zval %f  acirc %f  aratio %f\n",zval,acirc,aratio);
   	
      gmean=pow(htpix*wdpix,.5);    //setup
      q2=sqrt(aratio/pi);
      rcir=gmean*q2*pow(Nval,-exu);  //initialization
      ffac=gmean*q2;
      totarea=pi*pow(rcir,2);  //initialize area
      rcir0=rcir;
      rc[0]=rcir;
      xc[0]=rcir+rd()*(wdpix-2*rcir);   //place first circle
      yc[0]=rcir+rd()*(htpix-2*rcir);
      nt[0]=0;
   	   	
     // ------- start of computational loops ---------
      ncir=1;   //circle number init
      nitstot=0;
   	
      do   //loop on circle number
      {
         nits=0;   //iteration count init
         q1=(float)ncir+Nval;
         rtst=ffac*pow(q1,-exu);
         do    //loop on random search
         {
            x1=rtst+rd()*(wdpix-2*rtst);
            y1=rtst+rd()*(htpix-2*rtst);
            nits++;
            tst1=1;
            for(k=0;k<ncir;k++)   //loop over old placements
            {
               rsum=rtst+rc[k];  //sum of radii
               dx=xc[k]-x1;
               if(dx<0.) dx=-dx;
               if(dx<rsum)  //coarse test x
               {
                  dy=yc[k]-y1;
                  if(dy<0.) dy=-dy;
                  if(dy<rsum)  //coarse test y
                  {
                     q1=pow(dx,2)+pow(dy,2);
                     q2=sqrt(q1);
                     if(q2<rsum)  //fine test
                     {
                        tst1=0;
                        break;
                     }
                  }  //if dy<
               }  //if(dx<
            } //next k
         }while(tst1==0);	//repeat if too close to a circle
      
         m=fmod(ncir,nprint);
         if(m==0) printf("%i %7i %f %f %f \n",ncir,nits,u,rc[ncir-1],afill);
      
         nitstot=nitstot+nits;
         rc[ncir]=rtst;  //save data
         xc[ncir]=x1;
         yc[ncir]=y1;
         nt[ncir]=nitstot;
         u=(float)nitstot/ncir;
         ara=pi*pow(rtst,2);
         totarea=totarea+ara;
         afill=totarea/(wdpix*htpix);
         ncir++;
      }
      while(nits<itsmax & ncir<nmax & afill<fillmax);
   	
      printf("%i circles  c=%f  N=%i \n",ncir-1,cval,Nval);
      printf("extrapolated fill 1.0  actual fill = %f total its %i\n",afill,nitstot);
   	
      //create generic Postscript graphics file
      fprintf(ab,"1 1 1  setrgbcolor\n");
      fprintf(ab,"%7.4f in %7.4f in %7.4f in %7.4f in rectfill\n",0.,0.,wdpix,htpix);  //background
      for(j=0;j<=ncir-1;j++)     //circles with random color
      {
         fprintf(ab,"%6.3f %6.3f %6.3f  setrgbcolor\n",rd(),rd(),rd());
         fprintf(ab,"%f in %f in %f in 0 360 arc\n",xc[j],yc[j],rc[j]);
         fprintf(ab,"fill\n");
      }
   
   	//save circle data to file
      fprintf(cd,"w %f h %f %i\n",wdpix,htpix,ncir-1); 
      for(j=0;j<=ncir-1;j++)
         fprintf(cd,"%f  %f  %f  %i\n",xc[j],yc[j],rc[j],nt[j]);  
   		
      free(xc);
      free(yc);
      free(rc);
      free(nt);
   }
	
    double zeta_hurw(float u,int N)       //estimate Hurwitz zeta function
   {  //assumes pre-power factor is 1 for 1st term
      float nv;
      double v,sum;
      double j;
      sum=0.;
      for(j=N;j<=100000;j++)
      {
         nv=(float)j;
         sum=sum+pow(nv,-u);
      }
      v=sum+est_sum(u,100000.5);
      return v;
   }
	
    double est_sum(float c,float n)   //estimate sum of high-order terms
   {
      double v;
      v=(1./(c-1.))*pow(n-.5,1.-c);
      return v;
   }