/*
 * Calculate depth of field and depth of focus
 *
 * $Id: dof.c,v 1.17 2020/12/11 00:00:02 grog Exp grog $
 */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <sys/file.h>
#include <sys/limits.h>
#include <sys/stat.h>
#include <sys/types.h>
#include <unistd.h>
#include "photos.h"

float focal_length;                             /* focal length of lens, in mm */
float f;                                        /* focal length of lens, in m */
float u = 0.0;					/* distance of subject, metres */
float end_distance = 0.0;			/* end distance for focus stacking */
float v;                                        /* distance from lens to focal plane, metres */
float aperture;                                 /* aperture of lens, f/ */
float absolute_aperture;			/* aperture of lens, mm */
float hyperfocal;                               /* hyperfocal distance, metres */
float magnification;                            /* magnification of image, calculated */
float mag;					/* magnification parameter */
float confusion = 0.000008;                     /* circle of confusion */
float dof;                                      /* for macros only XXX */
float near_limit;                               /* DoF limits */
float nearest_limit;				/* near limit in first iteration */
float far_limit;
float exposure_compensation;			/* for macros */
int steps = 1;					/* number of focus steps */
int verbose = 0;                                /* show uninteresting stuff */
int virtual = 0;				/* show distances beyond ∞ if set */
float EV;					/* log(10) of 1.4 (1 EV) */
char *me;
int i;						/* loop counter */

void usage ()
{
  fprintf (stderr,
	   "Usage:\n"
	   "%s  [--] [-c confusion in µm] [-e end distance] [-m mag] [-s sensor] [-S steps] [-v] focal-length aperture [distance]\n"
	   "\tdistance defaults to hyperfocal distance\n"
	   "\tmag is magnification, overrides distance\n"
	   "\tdefault confusion is %1.2f µm\n"
	   "\t-S\tcalculate <steps> focus steps\n\n"
	   "\tIf distance < 0, it sepecifies the focal plane distance\n"
	   "sensor can be:\n"
           "     F\tFull frame\n"
           "     H\tAPS H\n"
           "     C\tAPS C (Nikon)\n"
           "     Canon\tAPS C (Canon)\n"
           "     S\tFoveon (Sigma)\n",
	   me,
	   confusion * 1000000);
  exit (1);
}

void missparm (char *flag)
{
  fprintf (stderr, "Missing arg for %s\n", flag);
  usage ();
  }

int main (int argc, char *argv [])
{
  char *fl_format;				/* format in which to display focal length */
  char *ap_format;				/* format in which to display aperture */
  int argcm1 = argc - 1;
  int parmcnt = 0;				/* positional paramter count */
  int focussteps = 0;				/* number of steps in focus stack */
  EV = log (M_SQRT2);				/* log(10) of 1.4 (1 EV) */
  struct sensor *s = &fourthirds;		/* sensor details */

  /* Parse arguments */
  me = argv [0];				/* global, for error messages */
  for (i = 1; i < argc; i++)
    {
    if (argv [i] [0] == '-')			/* flag */
      {
      if (strlen (argv [i]) != 2)		/* invalid parm */
	usage ();
      switch (argv [i] [1])
	{
	case '-':
	  virtual = 1;				/* show the virtual distance beyond ∞ */
	  break;

	case 'c':				/* circle of confusion */
	  i++;					/* point at arg */
	  if (i == argc)
	    missparm (argv [i]);
	  confusion = atof (argv [i]) / 1000000.0; /* in micro-m */
	  break;

	case 'm':
	  i++;					/* point at arg */
	  if (i == argc)
	    missparm (argv [i]);
	  mag = atof (argv [i]);
	  if (mag <= 0.0)			/* invalid, */
	    {
	    fprintf (stderr, "Magnification must be > 0, not %f\n", mag);
	    exit (1);
	    }
	  break;

	case 'e':
	  i++;					/* point at arg */
	  if (i == argc)
	    missparm (argv [i]);
	  end_distance = atof (argv [i]);
	  far_limit = end_distance;
	  if (mag <= 0.0)			/* invalid, */
	    {
	    fprintf (stderr, "Magnification must be > 0, not %f\n", mag);
	    exit (1);
	    }
	  break;
	case 's':
	  i++;					/* point to arg */
	  if (i == argc)
	    usage ();
	  else if (strcmp (argv [i], "F") == 0)
	    s = &fullframe;
	  else if (strcmp (argv [i], "H") == 0)
	    s = &aps_h;
	  else if (strcmp (argv [i], "C") == 0)
	    s = &aps_c;
	  else if (strcmp (argv [i], "Canon") == 0)
	    s = &aps_canon;
	  else if (strcmp (argv [i], "S") == 0)
	    s = &foveon;
	  else
	    usage ();
	  break;

	case 'S':
	  i++;					/* point at arg */
	  if (i == argc)
	    missparm (argv [i]);
	  if (! (steps = atoi (argv [i])))
	    {
	    fprintf (stderr, "Invalid step value: %s\n", argv [i]);
	    exit (1);
	    }
	  break;

	case 'v':
	  verbose = 1;
	  break;
	}
      }
    else					/* positional parameter */
      {
      switch (parmcnt)
	{
	case 0:					/* focal length in mm */
	  focal_length = atof (argv [i]);
	  f = focal_length / 1000.0;		/* convert to metres */
	  break;

	case 1:					/* aperture */
	  aperture = atof (argv [i]);
	  break;

	case 2:					/* distance (u) */
	  u = atof (argv [i]);
	  break;

	default:
	  fprintf (stderr, "Too many positional parameters\n");
	  usage ();
	}
      parmcnt++;
      }
    }

  if (parmcnt < 2)
    {
    if (parmcnt)				/* just usage with no parameters */
      fprintf (stderr, "Not enough parameters (%d)\n", parmcnt);
    usage ();
    }

  absolute_aperture = focal_length / aperture;
  hyperfocal = f * f / (aperture * confusion);
  if (! u)
    u = hyperfocal;				/* Distance defaults to hyperfocal */

  if (fabs (u) < f)
    {
    fprintf (stderr, "Distance %f m too close.  Minimum distance is focal length %f m\n", fabs (u), f);
    exit (1);
    }

  /* set output format prescisions based on focal length and aperture */
  if (focal_length < 10)
    fl_format = "Lens focal length:\t%5.1f mm\n";
  else
    fl_format = "Lens focal length:\t%4.0f mm\n";
  if (aperture < 1)
    ap_format = "Aperture:\t\t f/%4.2f\n";
  else if (aperture < 10)
    ap_format = "Aperture:\t\t f/%3.1f\n";
  else
    ap_format = "Aperture:\t\tf/%4.1f\n";
  printf (fl_format, focal_length);
  printf (ap_format, aperture);
  printf ("Absolute aperture:\t  %4.2f mm\n", absolute_aperture);
  printf ("Hyperfocal distance:\t%8.3f m\n"
	  "Circle of confusion:\t%7.2f µm\n"
	  "\n",
          hyperfocal,
          confusion * 1000000);
  printf ("Subject        Focal plane     Magnification  Exposure       Near        Far         Depth of\n"
          "distance (m)   distance (mm)                    comp  (EV)   limit (m)   limit (m)  field (m)\n");

  while ((steps > 0) || (far_limit < end_distance))
    {
    steps--;
    focussteps++;				/* one more focus step */
    if (mag)					/* specifying magnification, not distance */
       {
       magnification = mag;			/* use what's gven */
       v = f * (mag + 1);
       u = 1 / (1 / f - 1 / v);
       mag = 0;					/* continue with u value if we have more steps */
       }
    else if (u > 0)				/* normal distance spec */
      v = 1 / (1 / f - 1 / u);
    else					/* negative: u is really v */
      {
      v = -u;
      u = 1 / (1 / f - 1 / v);
      }
    exposure_compensation = log (v / f) / EV;

    magnification = v / u;
    near_limit = hyperfocal * u / (hyperfocal + u);
    if (! nearest_limit)
      nearest_limit = near_limit;
    /* This part of the format goes up to near limit */
    if (u < 1.5)				/* we call this a macro, different precision */
      fl_format = "%9.4f       %9.4f    %8.3f            %2.1f      %9.6f";
    else
      fl_format = "%7.2f         %8.3f     %7.2f             %2.1f      %8.3f";
    printf (fl_format,
	    u,
	    v * 1000,
	    magnification,
	    exposure_compensation,
	    near_limit );
    /*
     * If we're less than the hyperfocal distance, the far limit is less than ∞.
     * We calculate it only then.  Otherwise it's left at 0, implying ∞,
     * *unless* we're really interested in how far beyond ∞ it goes.
     */
    if ((u < hyperfocal) || virtual)
      {
      far_limit = hyperfocal * u / (hyperfocal - u);
      /*
       * Set precision based on subject distance (in metres).
       * This includes far limit and depth of field.
       */
      if (u < 0.2)
	fl_format ="    %9.6f %9.6f\n";
      else if (u < 1.5)
	fl_format ="  %9.4f %9.4f\n";
      else
	fl_format ="   %8.3f  %8.3f\n";
      printf (fl_format,
	      far_limit,
	      far_limit - near_limit );
      u = far_limit;
      if (far_limit < 0)			/* wrapped around, */
	end_distance = far_limit;		/* don't try any further */
      }
    else					/* far limit ∞ or further */
      {
      printf ("      infinity   infinity\n");
      steps = 0;				/* go no further */
      }
    }

  if (focussteps > 1)				/* print this only if we really have steps */
    printf ("\n%d focus steps", focussteps);

  /*
   * There's no point trying to print a "depth of field" if the far limit is
   * infinity (which here is indicated by the opposite, 0).
   */
  if (far_limit > 0)				/* not infinity or wrapped around */
    {
    if (nearest_limit < 0.1)
      {
      if (far_limit < 1)				/* do it in mm */
	printf ("\nTotal depth of field: %0.3f mm (%0.3f mm - %0.3f mm)\n",
		(far_limit - nearest_limit) * 1000,
		nearest_limit * 1000,
		far_limit * 1000);
      else
	printf ("\nTotal depth of field: %0.6f (%0.6f m - %0.6f m)\n",
		far_limit - nearest_limit,
		nearest_limit,
		far_limit);
      }
    else if (nearest_limit < 1)
      printf ("\nTotal depth of field: %0.4f (%0.4f m - %0.4f m)\n",
	      far_limit - nearest_limit,
	      nearest_limit,
	      far_limit);
    else
      printf ("\nTotal depth of field: %0.2f (%0.2f m - %0.2f m)\n",
	      far_limit - nearest_limit,
	      nearest_limit,
	      far_limit);
    }
  if (magnification > 0.4)			/* mm appropriate */
    printf ("Field of view: %0.1f x %0.1f mm\n",
	    s->width / magnification,
	    s->height / magnification );
  else if (magnification > 0.04)		/* whole mm appropriate */
    printf ("Field of view: %0.0f x %0.0f mm\n",
	    s->width / magnification,
	    s->height / magnification );
  else						/* metres */
    printf ("Field of view: %0.2f x %0.2f m\n",
	    s->width / magnification / 1000,
	    s->height / magnification / 1000 );


  return 0;
}