#include <stdio.h>
#include <stdlib.h>
#include <math.h>

/* Program to generate data for 2D plots by running many times my other program, bokeh.c

// v2: Writing the Fnumbers.dat file with F labels for L=0.3, 0.6, and 2m, and every lens.

The goal is to visualize your whole collection of lenses, so you can find the best bokeh lens
for any shooting conditions.

I used these sources:
http://en.wikipedia.org/wiki/Depth_of_field
http://en.wikipedia.org/wiki/Circle_of_confusion

The solution is exact for a thin lens.

File "lenses.txt" has to be present in the current directory. The format:
 - First line: crop factor of the camera;
 - Each prime lens (and each zoom of a zoom lens) has one line:

Lens_name  Focal_distance_in_mm  Fmin  Fmax

Here Fmin and Fmax are the minimum and maximum values of the aperture for the given lens and zoom.

*/

// Maximum number of lines in the lenses.txt file:
#define LENS_MAX 1000
// Plot dimensions:
#define NPLOT 300

struct lenses
{
  char name[20];
  int id;
  double focal;
  double Fmin;
  double Fmax;
  double F;
  double d;
  double c;
  float x0;
  float y0;
  int npix;
};

int main (int argc,char **argv)
{

  // The minimum circle of confusion for a full-frame camera (m):
  const double
    c_full = 3e-5,
    // Range of S (distance between the object and the background), in meters:
    S1 = 1,
    S2 = 100,
    // Range of L (frame height at the object's distance), in meters:
    L1 = 0.1,
    L2 = 3.16;

  struct lenses 
    lens[LENS_MAX], lens1, blens[2];

  double 
    L, DOF_min, S, crop, c0, l, B, Fmin, cmax, f;

  float
    bokeh[NPLOT][NPLOT],
    apert[NPLOT][NPLOT],
    dist[NPLOT][NPLOT],
    best[NPLOT][NPLOT];

  float
    al, jS1, F1, d1;

  int
    id[NPLOT][NPLOT];

  int
    i, Nmax, j, jmax, jS, jL, jLmax,jSmax, id1, id0, jLF[3], k, nF;

  char
    line[255];

  FILE
    *fp,
    *fp2;


  if (argc != 2)
    {
      (void) fprintf(stdout,"\nUsage: %s DOF_min\n",argv[0]);
      fprintf(stdout,"\nHere\n");
      fprintf(stdout,"  DOF_min: Minimum depth of field required (m) (can be zero if DOF is not important);\n");
      exit(0);
    }

  DOF_min = atof(argv[1]);

  // Reading the lenses database:
  i = -2;
  fp=fopen("lenses.txt","rt");
  while(fgets(line, 255, fp) != NULL)
    {
      i++;
      if (i == -1)
	{
	  sscanf (line, "%lf", &crop);
	} else
	{
	  if (i > LENS_MAX)
	    {
	      fprintf(stdout,"\nToo many lines in lenses.txt!\n");
	      exit(1);
	    }
	  if (line[0] == '#')
	    {
	      i--;
	    } else
	    {
	      sscanf (line, "%s %lf %lf %lf", lens[i].name, &lens[i].focal, &lens[i].Fmin, &lens[i].Fmax);
	      // Converting focal to meters:
	      lens[i].focal = lens[i].focal / 1000;
	      lens[i].id = i;
	    }
	}
    }
  fclose(fp);
  
  Nmax = i + 1;

  // Maximum circle of confusion for DOF computations (m):
  c0 = c_full / crop;
  // Sensor height (m):
  l = 0.024 / crop;

  // jL corresponding to 0.3m, 0.6m and 2m:
  L=0.3;  jLF[0] = (int) rint (log10(L/L1)/log10(L2/L1)*(NPLOT-1));
  L=0.6;  jLF[1] = (int) rint (log10(L/L1)/log10(L2/L1)*(NPLOT-1));
  L=2.0;  jLF[2] = (int) rint (log10(L/L1)/log10(L2/L1)*(NPLOT-1));
  //  printf("%d %d %d\n",jLF[0],jLF[1],jLF[2]);

  // 2D cycle for L and S:

  for (jL=0; jL<NPLOT; jL++)
    {
      L = L1*pow(10.0,(float)jL/(NPLOT-1.0)*log10(L2/L1));
      for (jS=0; jS<NPLOT; jS++)
	{
	  S = S1*pow(10.0,(float)jS/(NPLOT-1.0)*log10(S2/S1));

  // Cycle for all lenses and zooms:
  for (i=0; i<Nmax; i++)
    {
      // Distance from the lens to the foreground (m):
      lens[i].d = (L/l+1) * lens[i].focal;
      B = c0 * (lens[i].d-lens[i].focal) / (lens[i].focal*lens[i].focal);
      lens[i].x0 = 0.0;     lens[i].y0 = 0.0;  lens[i].npix = 0;  
      if (DOF_min == 0.0)
	{
	  lens[i].F = lens[i].Fmin;
	} else
	{
	  Fmin = (pow(lens[i].d*lens[i].d + DOF_min*DOF_min,0.5) - lens[i].d) / (B*DOF_min);
	  if (Fmin < lens[i].Fmin)
	    {
	      lens[i].F = lens[i].Fmin;
	    } else if (Fmin > lens[i].Fmax)
	    {
	      // Bad case:
	      lens[i].F = -1;
	      lens[i].c = -1;
	      lens[i].d = -1;
	      continue; 
	    } else
	    {
	      lens[i].F = Fmin;
	    }
	}
      lens[i].c = S/(lens[i].d+S) * lens[i].focal*lens[i].focal / lens[i].F / (lens[i].d-lens[i].focal) / c0;
    }

  // Sorting the structure array according to c (only for two best lenses):
  for (i=0; i<2; i++)
    {
      cmax = -2;
      for (j=i; j<Nmax; j++)
	{
	  if (lens[j].c > cmax)
	    {
	      cmax = lens[j].c;
	      jmax = j;
	    }
	}
      blens[i] = lens[jmax];
      lens[jmax].c = -1.0;
    }

  // Writing to the plot arrays:
  id[jL][jS] = blens[0].id;
  bokeh[jL][jS] = blens[0].c;
  apert[jL][jS] = blens[0].F;
  dist[jL][jS] = blens[0].d;
  if (blens[1].c > 0)
    {
      best[jL][jS] = blens[0].c/blens[1].c;
    } else
    {
      best[jL][jS] = -1.0;
    }

	} // S
    } // L

  // Writing the file with F labels:
  fp2 = fopen("Fnumbers.dat","wt");
  for (k=0; k<3; k++)
    {
      jL = jLF[k];
      for (jS=0; jS<NPLOT; jS++)
	{
	  if (jS == 0)
	    {
	      id1 = id[jL][jS];
	      F1 = apert[jL][jS];
	      d1 = dist[jL][jS];
	      jS1 = (float) jS;
	      nF = 1;
	    } else
	    {
	      id0 = id1;
	      id1 = id[jL][jS];
	      if (id1 == id0 && jS < NPLOT-1) 
		{
		  F1 = F1 + apert[jL][jS];
		  d1 = d1 + dist[jL][jS];
		  jS1 = jS1 + (float) jS;
		  nF = nF + 1;
		} else
		{
		  F1 = F1 / nF;
		  d1 = d1 / nF;
		  jS1 = jS1 / nF;
		  if (jS1<25) {jS1=25;}
		  if (jS1>NPLOT-25) {jS1=NPLOT-25;}
		  fprintf (fp2, "%7.1f %d \"f%.1f (%.1fm)\"\n", jS1,jL,F1,d1);
		  id1 = id[jL][jS];
		  F1 = apert[jL][jS];
		  d1 = dist[jL][jS];
		  jS1 = (float) jS;
		  nF = 1;
		}
	    }
	}
    }
  fclose (fp2);

  // Writing the data tables:

  fp=fopen("id.dat","wt");
  for (jL=0; jL<NPLOT; jL++)
    {
      for (jS=0; jS<NPLOT; jS++)
	{
	  if (bokeh[jL][jS] < 0)
	    {
	      id1 = -1;
	    } 
	  else
	    {
	      id1 = id[jL][jS];
	    }
	  fprintf (fp, "%4d", id1);
	}
      fprintf (fp, "\n");
    }
  fclose (fp);

  fp=fopen("bokeh.dat","wt");
  for (jL=0; jL<NPLOT; jL++)
    {
      for (jS=0; jS<NPLOT; jS++)
	{
	  fprintf (fp, "%9.3f", bokeh[jL][jS]);
	}
      fprintf (fp, "\n");
    }
  fclose (fp);

  fp=fopen("apert.dat","wt");
  for (jL=0; jL<NPLOT; jL++)
    {
      for (jS=0; jS<NPLOT; jS++)
	{
	  fprintf (fp, "%5.1f", apert[jL][jS]);
	}
      fprintf (fp, "\n");
    }
  fclose (fp);

  fp=fopen("dist.dat","wt");
  for (jL=0; jL<NPLOT; jL++)
    {
      for (jS=0; jS<NPLOT; jS++)
	{
	  fprintf (fp, "%7.2f", dist[jL][jS]);
	}
      fprintf (fp, "\n");
    }
  fclose (fp);

  fp=fopen("best.dat","wt");
  cmax = -1.0;
  for (jL=0; jL<NPLOT; jL++)
    {
      for (jS=0; jS<NPLOT; jS++)
	{
	  fprintf (fp, "%7.3f", best[jL][jS]);
	  if (best[jL][jS]>cmax)
	    {
	      cmax = best[jL][jS];
	      jLmax = jL;  jSmax = jS;
	    }
	}
      fprintf (fp, "\n");
    }
  fclose (fp);

  // Writing the file with lens labels:
  fp=fopen("labels.dat","wt");
  al = 1;
  for (jL=0; jL<NPLOT; jL++)
    {
      for (jS=0; jS<NPLOT; jS++)
	{
	  if (bokeh[jL][jS] > 0.0)
	    {
	      id1 = id[jL][jS];
	      lens[id1].x0 = lens[id1].x0 + pow((float)jS,al);
	      lens[id1].y0 = lens[id1].y0 + pow((float)jL,al);
	      lens[id1].npix = lens[id1].npix + 1;
	    }
	}
    }
  for (i=0; i<Nmax; i++)
    {
      if (lens[i].npix > 0)
	{
	  lens[i].x0 = pow(lens[i].x0/lens[i].npix,1.0/al);
	  lens[i].y0 = pow(lens[i].y0/lens[i].npix,1.0/al);
	  fprintf (fp, "%7.1f %7.1f %s\n",  lens[i].x0, lens[i].y0, lens[i].name);
	}
    }
  fclose (fp);


  L = L1*pow(10.0,(float)jLmax/(NPLOT-1.0)*log10(L2/L1));
  S = S1*pow(10.0,(float)jSmax/(NPLOT-1.0)*log10(S2/S1));
  f = 1000 * lens[id[jLmax][jSmax]].focal;
  printf ("\n  The best bokeh advantage =%7.3f for the lens %s (f=%3.0f) at L=%6.3f, S=%6.1f\n\n", 
	  best[jLmax][jSmax],lens[id[jLmax][jSmax]].name,f,L,S);

}