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

/* Program to find the lenses with the best (largest) bokeh (or "circle of confusion")
when the following parameters are given:

L: The height of the frame at the object's distance (m);
DOF_min: Minimum depth of field required (m); can be zero if DOF is not important;
S: The distance between the object and the background (set to a negative number if the background is in the infinity);
crop: The crop factor of the camera (1 for a full-frame camera; >1 for smaller frames).

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

struct lenses
{
  char name[20];
  double focal;
  double Fmin;
  double Fmax;
  double F;
  double d;
  double c;
};

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

  // The minimum circle of confusion for a full-frame camera (m):
  const double
    c_full = 3e-5;

  struct lenses 
    lens[LENS_MAX], lens1;

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

  int
    i, Nmax, j, jmax;

  char
    line[255];

  FILE
    *fp;

  if (argc != 4)
    {
      (void) fprintf(stdout,"\nUsage: %s L DOF_min S\n",argv[0]);
      fprintf(stdout,"\nHere\n");
      fprintf(stdout,"  L:       The height of the snapshot at the object's distance (m);\n");
      fprintf(stdout,"  DOF_min: Minimum depth of field required (m) (can be zero if DOF is not important);\n");
      fprintf(stdout,"  S:       The distance between the object and the background (set to a negative number if the background is in the infinity).\n\n");
      exit(0);
    }

  L = atof(argv[1]);
  DOF_min = atof(argv[2]);
  S = atof(argv[3]);
  if (S < 0.0)
    {
      S = 1e6;
    }

  // 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);
	    }
	  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;
	}
    }
  fclose(fp);
  
  Nmax = i + 1;

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

  // 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);
      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;
	      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:
  for (i=0; i<Nmax-1; i++)
    {
      cmax = -2;
      for (j=i; j<Nmax; j++)
	{
	  if (lens[j].c > cmax)
	    {
	      cmax = lens[j].c;
	      jmax = j;
	    }
	}
      lens1 = lens[i];
      lens[i] = lens[jmax];
      lens[jmax] = lens1;
    }

  // Printing the sorted output (maximum 50 lines):
  printf ("\n    Lenses sorted according to bokeh\n");
  printf ("  %15s%5s%8s%5s%8s\n\n","Name           ","f(mm)","d(m)","F ","c ");
  for (i=0; i<50; i++)
    {
      if (lens[i].F > 0)
	{
	  // Converting focal to mm:
	  f = lens[i].focal * 1000;
	  printf("%15s %5.0lf %8.2lf %5.1lf %8.2lf\n",lens[i].name,f,lens[i].d,lens[i].F,lens[i].c);	  
	}
    }


}