--[[focus_brkt_frames_for_X-T3_R16 
To compute the number of frames to use in the focus bracketing 
menu of the X-T3 when the magnification is smaller than 1. 
  Set the parameters below and press Run on the command bar 
above.  The result prints in the window down the screen. 
  In the case of the XF80mm and starting at the minimum distance 
add a small number like 0.00001m to the start distance. 
  There are four usage scenarios depending on whether one 
measures 
  1. the lens-to-frame distance for the first and last frame 
     (Df and Dl). 
  2. the sensor-to-frame distance (Dsf and Dsl), 
  3. the height of the first frame and the distance between the 
     first and the last frame (firstFrameHeight if T.f then and zoneDepth), 
  4. the frame height at the mid-point between the first and 
     the last frame (middleFrameHeight and zoneDepth). 
  "veryFar" will be used for Dl, Dsl and zoneDepth when the 
parameter is set to "nil" (like Dl = nil). A scenario is 
evaluated when the first parameter is a number and not nil. 
 --]] 

veryFar        = 1.0e20 
sensorHeight   = 0.0156 -- meter. X-Trans-4. 
sensorHPixels  = 4160  -- X-Trans-4 

f              = 0.014 -- Lens focal length. 
Af             = 10   -- Aperture number. 
cameraStep     = 7    -- step setting from focus brkt menu. 

Df = 0.175-f; Dl = 0.54-f; -- Lens-to-frame distance for first and last frame. 
Dsf = 0.175; Dsl = 0.54; -- Sensor-to-frame distance. 
firstFrameHeight = 0.16; zoneDepth = 0.35; 
middleFrameHeight = nil; 

-- Circle of confusion for the stack. Thanks to SVJIM. 
coc = 0.2 * cameraStep * sensorHeight / sensorHPixels 

--[[Focus stack for Fujifilm X. Estimate the number of frames (ENF). 
   * Copyright 2020 Richard Lemieux 
     Free under the Free Software Foundation GPL-V3 license. 
   * Exact computation for the thin lens geometric optic model. 
     Algebra done with MapleV. 
   * All distances as well as focal lengths are in meters. 
   * A frame here is the canvas centered at the focus point on 
     the subject side of the lens. The image of the frame fits 
     exactly the sensor. 
   * The thickness of each frame is the DOF. The same value of 
     the  circle of confusion is used in the DOF formula for all 
     the frames in the stack. 
   * The distance between frames is taken from center to center. 
   * The value of the coc used to compute the sequence 
     here is NOT the actual coc of the pictures as shot by the 
     camera. It is the coc the camera uses to advance the focus 
     element at each step of the stack. Check the actual formula 
     used here by looking at the parameter section above. 
       The coc used in the stack is the same as the coc of the 
     full resolution picture when cameraStep is 5. And it is 
     the same as the coc of half resolution pictures when 
     cameraStep is 10. 
   * The sequence of frames is computed so the frames just 
     touch each other and do not overlap using the coc computed 
     here. The actual overlap comes from the actual coc of the 
     pictures. My assumption is that people at Fujifilm use the 
     same sequence :-) . Test results look good 0 far. However 
     discrepancies appeared after I took into account the 
     fixed sensor setup ... D2Ds() and Ds2D(). 
   * The lens-to-frame distance (D) is used internally for the 
     the computations in the framework of the thin lens model 
     used here. 
       "D" is measured between the center of the thin lens and 
     the center of the frame in focus. 
       Alternatively Fuji cameras have a mark on the top plate 
     (a line with a circle in the middle) at the position of the 
     sensor. Function Ds2D() below will compute D if Ds is 
     provided. 
       "Ds" denotes a distance measured from the sensor to the 
     frame that is in focus. 
       At small magnifications (sensor size / frame size) as 
     occur in people photography Ds ~= D + f. 
   * Distances xf, xl and xc refer to the position of the focus 
     point on the sensor side of the lens. xf, xl and xc are 
     measured relative to the Infinity focus point position. 
     For instance xf + f is the distance between the sensor and 
     the lens when the focus is on the first frame. 
   * All formulas here result from symbolic manipulations done 
     in MapleV starting from the very basic geometry of thin 
     lens geometric optics as teached in most High-School 
     programs. Any errors are mine. 
   * The programs included here come with no warranty of any kind. 
   * Remember that Ds can't be smaller than four times the focal 
     length :-) . Look at the Ds2D function. 
--]] 

if cls then cls() end 

local math = require("math") 
local abs = abs or math.abs 
local floor = floor or math.floor 
local log10 = log10 or function (n1) return math.log(n1)/math.log(10) end 
local max = max or math.max 

ENF_L = function (xf,xl,nf) 
  -- With the lens fixed, step the sensor from its position when 
  -- focused on the first frame at the close end (xf) to its 
  -- position at the far end (xl) while counting the number of 
  -- steps/frames (nf). 
  if (not xf) then return "Please set frame height or distance" 
  elseif (xf < xl) then return nf 
  else local  xStepSize = 2*Af*coc*(xf+f)/(Af*coc+f) --No gap. 
       return ENF_L(xf-xStepSize,xl,nf+1) 
  end 
end 

gap = function (SS,xc) 
  -- The gap computed here is the distance between the back 
  -- of the DOF range of the current frame and the front of the 
  -- DOF range of the next frame. That's on the subject side 
  -- of the lens. 
  -- "SS" is the distance considered for the next move of the 
  -- lens. "xc-SS" would become the next position of the lens. 
  -- "xc" is the distance between the sensor and the focal 
  -- point of the lens with the focus point set on the current 
  -- frame. "xc+f" is the distance between the sensor and 
  -- the thin lens. 
  local gapN = SS*(Af*coc+f) - 2*Af*coc*(xc+f) 
  local gapD = (xc-SS+Af*coc)*(xc-Af*coc) 
  return f*gapN/gapD 
end 

-- First derivative of gap (d gap / d SS) 
gap_ = function (SS,xc) 
  return f*(f-Af*coc)/(xc-SS+Af*coc)^2 
end 

-- Second derivative of gap (d2 gap / d2 SS) 
gap__ = function(SS,xc) 
  return -2*gap_(SS,xc)/(xc-SS+Af*coc) 
end 

ENF_S1 = function (xf,xl,nf) 
  -- With the sensor fixed, step the lens from the first frame 
  -- to the last frame. Do it like ENF_L does it but with a 
  -- small gap between the DOF zones to account for the lens 
  -- movement. 
  --   When the sensor is fixed and the lens moves by the 
  -- distance SS, we want gap(SS,x) = SS. 
  --   Think of each step as being done in two movements. (1) 
  -- With the lens fixed move the sensor forward by distance SS 
  -- until the gap with the next frame DOF is also SS. (2) Move 
  -- the whole camera backward by the distance SS. The gap is 
  -- gone and the sensor has returned to where it was initially, 
  --   A numerical 2-step iteration is used here to compute SS. 
  -- The case when magnification is 1 is handled with a hack. 
  if (not xf) then return "Please set frame height or distance" 
  elseif (xf < xl) then return nf 
  else 
    local xStepSize 
    local SS0 = 2*Af*coc*(xf+f)/(Af*coc+f) --candidate xStepSize 
    local DSS0,SS1 
    local G_0 = gap_(SS0,xf) - 1 
    if (abs(G_0) < 0.1) then 
      DSSh0 = (7*max(1,f/0.08))*SS0 --  This occurs at magnification 1 
      local SSh1 = SS0 + DSSh0 
      local DSSh1 = (SSh1 - gap(SSh1,xf))/(gap_(SSh1,xf) - 1) 
      SS1 = SSh1 + DSSh1 
    else 
      DSS0 = SS0/G_0 -- gap is 0 here. 
      SS1 = SS0 + DSS0 
    end 
    local DSS1 = (SS1 - gap(SS1,xf))/(gap_(SS1,xf) - 1) 
    local SS2 = SS1 + DSS1 
    xStepSize = SS2 
    --print(string.format( 
    --   "x %.3g, M %.3g, SS %.3g, G() %.3g, G_() %.3g", 
    --   xf, xf/f, xStepSize, gap(xStepSize,xf)-xStepSize, 
    --        gap(xStepSize,xf)-1)) 
    return ENF_S1(xf-xStepSize,xl,nf+1) 
  end 
end 

-- Select the setup: either the lens or the sensor is fixed. 
ENF = ENF_S1 

-- Now compute the lens position when the focus is on the frames 
-- at the close end and at the far end of the stack. 
FFH2X = function (firstFrameHeight,zoneDepth) 
  -- Sets xf and xl given the height of the first frame and 
  -- the distance between the first and the last frame. 
  local M = sensorHeight / firstFrameHeight -- Magnification of first frame 
  local xf = f*M 
  local Df = f + f/M;  -- M = f / (Df - f), Df: Distance of first frame 
  local Dsf = D2Ds(Df) 
  local Dsl = Dsf + zoneDepth 
  local Dl = Ds2D(Dsl) 
  local xl = f*f/(Dl - f) 
  return xf, xl, Dsf, Dsl 
end 
DD2X = function (firstFrameDistance,lastFrameDistance) 
  -- Same as above except we are given the lens-to-frame distance 
  -- in the framework of the thin lens model. Use the function 
  -- Ds2D() given below to convert when the distance is 
  -- measured from the sensor and not from the lens. 
  local xf = f*f/(firstFrameDistance - f) 
  local xl = f*f/(lastFrameDistance - f) 
  return xf, xl 
end 
MFH2X = function (middleFrameHeight,zoneDepth) 
  -- Sets xf and xl given the height of the middle frame and 
  -- the distance between the first and the last frame. 
  local M = sensorHeight / middleFrameHeight -- Magnification of middle frame 
  local Dm = f + f/M;  -- M = f / (Dm - f), Dm: Distance of middle frame 
  local Dsm = D2Ds(Dm) 
  local Dsf = Dsm - zoneDepth/2 
  local Dsl = Dsm + zoneDepth/2 
  local Df = Ds2D(Dsf) 
  local Dl = Ds2D(Dsl) 
  local xf = f*f/(Df - f) 
  local xl = f*f/(Dl - f) 
  return xf, xl, Dsf, Dsl, Dsm 
end 

-- Convert between distance-to-sensor Ds and distance-to-lens D 
-- for a thin single glass lens. 
Ds2D = function (Ds) 
  -- This root applies at magnifications smaller than 1. 
  local D = 0.5 * (Ds + (Ds*Ds - 4*Ds*f)^0.5) 
  return D 
end 
D2Ds = function(D) 
  local Ds = D*D/(D - f) 
  return Ds 
end 

function n2t(v1) 
  -- Transform a number to text to "sd" significant digits. 
  -- This code does not handle negative numbers. 
  if (v1 < 0) then return V1 end 
  if (v1 >= veryFar) then return "veryFar" end 
  local sd = 3  -- Number of significant digits to print. 
  local fl10 = floor(log10(v1)) 
  local roundfactor = 10^(sd-1-fl10) 
  local v1rounded = floor(v1*roundfactor+0.5)/roundfactor; 
  local decimals = max(0,sd-fl10-1) ; 
  return string.format("%0."..decimals.."f",v1rounded) 
end 

-- Scenarios 
sD = function (Df,Dl)   -- lens-to-frame distance. 
  if (Df < 2*f) then print("Df is too small. Setting Df to 2*f:",n2t(2*f)) 
                     Df =  2*f end 
  local xf,xl = DD2X(Df,Dl or veryFar) 
  local frames = string.format("%0d", ENF(xf,xl,1)) 
  local Dsf = D2Ds(Df) 
  local Dsl = D2Ds(Dl or veryFar) 
  print( "Using lens-to-frame distance "..frames.." frames".. 
    " (Dsf "..n2t(Dsf).."m, Dsl "..n2t(Dsl).."m)") 
end 
sDs = function (Dsf,Dsl)  -- sensor-to-frame fistance. 
  if (Dsf < 4*f) then print("Dsf is too small. Setting Dsf to 4*f:",n2t(4*f)) 
                     Dsf =  4*f end 
  local Df = Ds2D(Dsf) 
  local Dl = Ds2D(Dsl or veryFar) 
  local xf,xl = DD2X(Df,Dl) 
  local frames = string.format("%0d", ENF(xf,xl,1)) 
  print("Using sensor-to-frame distance "..frames.." frames") 
end 
sFFH = function(firstFrameHeight,zoneDepth) 
  local xf,xl,Dsf,Dsl = FFH2X(firstFrameHeight, zoneDepth or veryFar) 
  local frames = string.format("%0d", ENF(xf,xl,1)) 
  print("Using first frame height "..frames.." frames".. 
    " (Dsf "..n2t(Dsf).."m, Dsl "..n2t(Dsl).."m)") 
end 
sMFH = function (middleFrameHeight, zoneDepth) 
  if zoneDepth then 
    local xf,xl,Dsf,Dsl,Dsm = MFH2X(middleFrameHeight, zoneDepth) 
    if (xf > f or xf < 0) then xf = f; Dsf = 4*f; 
      print("zoneDepth too large for this frame size. Clipping.") 
    end 
    local frames = string.format("%0d", ENF(xf,xl,1)) 
    print("Using middle frame height "..frames.." frames".. 
      " (Dsf "..n2t(Dsf).."m, Dsl "..n2t(Dsl).. 
      "m, Dsm "..n2t(Dsm).."m)") 
  else 
    print( "Using middle frame height:  Please set zoneDepth.") 
  end 
end 

GO = function () 
  local c = 0 -- Counts the scenarios that activate. 
  if Df then c=c+1; sD(Df,Dl); end 
  if Dsf then c=c+1; sDs(Dsf,Dsl); end 
  if firstFrameHeight then c=c+1; sFFH(firstFrameHeight,zoneDepth); end 
  if middleFrameHeight then c=c+1; sMFH(middleFrameHeight,zoneDepth); end 
  if (c == 0) then 
    print("Please set either firstFrameHeight or Df or Dsf") 
  end 
end 

GO() 

testP=false 
if testP then 
  print("2*f "..n2t(Ds2D(D2Ds(2*f)))) -- OK down to lens-to-frame distance = 2*f. 
  print("Df "..Df..", Dl "..Dl.. 
    ", Dsf "..Dsf..", Dsl "..Dsl) -- Didn't change. 
  print( 4*f ) -- Reminder. Ds may not be set smaller than that. 
  print("xf",xf,"xl",xl,"Dsm",Dsm) -- Shoudn't be global. 
end 

-- Plot 
if false then 
  cls() 
  local G = function(SS,xc) return gap(SS,xc) - SS end 
  local gSS, gG = {},{} 
  local i1,i2,i3 = 1,21,1 
  local SS1,SS2 = -2e-3, 2e-3 
  -- (SS - SS1)/(SS2-SS1) = (i-i1)/(i2-i1) 
  local i2ss=function(i) return SS1+(SS2-SS1)*(i-i1)/(i2-i1) end 
  for i4=i1,i2,i3 do gSS[i4]=i2ss(i4); gG[i4]=G(gSS[i4],f) end 
  p = plot.new() 
  plot.add(p, gSS, gG) 
  plot.set(p, "title", "(gap - SS) vs SS for x = f") 
  plot.set(p, 1, "size", 2) -- Curve #1 
  plot.set(p, 1, "style", "-") 
end