// Uses Frink's ability to render accurately-sized graphics to draw Fresnel
// zone plates.  You can print these out accurately on a laser printer with
// transparencies.  A Fresnel Zone plate is almost like a printable lens.
// it takes advantage of quantum diffraction of light to focus light onto
// a target. They are freaking awesome.
//
// The main limitation of creating your own zone plates is that you have to
// be able to print at very high resolution to be able to focus visible light
// with a short focal length.
//
// http://zoneplate.lbl.gov/theory
// http://en.wikipedia.org/wiki/Zone_plate

// Draw a (printable, exactly-scaled) zone plate.
//   * f is the focal length,
//   * lambda is the target wavelength of light you want to focus
//   * phase is 0 or 1.  If phase=0, the center will be clear.  If phase=1,
//     the center will be dark.  Surprisingly, both will focus the light at
//     the geometric center of the target, even if the center is geometrically
//     in shadow!
//   * r is the radius of the zone plate to be drawn.
//   * res is the limiting resolution of your printer.  (This is at best half
//     the actual resolution of your printer.)  This affects the "max effective
//     zones" that your printer can meaningfully render.
drawZonePlate[f=27 feet, lambda=700 nm, phase=0, r=3 in, res=600/in] :=
{
   g = new graphics

   // maxn is the number of zones that will be drawn
   maxn = floor[-2 f lambda^-1 + 2 lambda^-2 (f^2 lambda^2 + lambda^2 r^2)^(1/2)]
   if (maxn > 10000)
      maxn=10000

   // maxr is the radius that will actually be drawn
   maxr = sqrt[maxn lambda f + maxn^2 lambda^2 / 4]
   n = maxn
   while n > 0
   {
      r = sqrt[n lambda f + n^2 lambda^2 / 4]
      c = (n+phase) mod 2 
      g.color[c,c,c]
      g.fillEllipseCenter[0 in, 0 in, 2 r, 2r]
      n = n - 1
   }

   // maxEffectiveN is the maximum number of zones that can be meaningfully
   // rendered (due to your printer's resolution limits)
   maxEffectiveN = floor[(lambda (4 f + lambda) res^2 (4 - lambda^2 res^2) - 2 sqrt[(lambda^2 res^2 (-2 + lambda res) (2 + lambda res)(-4 + (-16 f^2 + lambda^2) res^2))]) / (2 lambda^2 res^2 (-4 + lambda^2 res^2))]

   // maxEffectiveR is the maximum radius of the zone plate that will be
   // meaningfully rendered (due to your printer's resolution limits)
   maxEffectiveR = sqrt[maxEffectiveN lambda f + maxEffectiveN^2 lambda^2 / 4]

   // maxEffectiveArea is the maximum area of the zone plate that will be
   // meaningfully rendered (due to your printer's resolution limits)
   maxEffectiveArea = pi maxEffectiveR^2

   maxPower = maxEffectiveArea * 500 W/m^2
   
   str = ""
   g.color[0,0,0]
   g.font["SansSerif", .5 cm]
   str = str + "wavelength =" + format[lambda, "nm", 2] + "\n"
   str = str + "focal length =" + format[f, "m", 3] + "\n"
   str = str + "$maxn zones drawn" + "\n"
   str = str + "diameter =" + format[2 maxr,"in",3] + "\n"
   str = str + "Max effective zones=$maxEffectiveN\n" 
   str = str + "Max effective diameter=" + format[2 maxEffectiveR, "cm", 4]  + "\n"
   str = str + "Max effective area=" + format[maxEffectiveArea, "m^2", 6] + "\n"
   str = str + "Max power=" + format[maxPower, "W", 3]

   g.text[str, 0 cm, maxr, "center", "top"]
   
   g.write["zoneplate.svg", 1200, 1200]
   g.show[]
}


// Alter this line to match the focal length and wavelength you desire.
drawZonePlate[50 feet, 700. nm, 0, 3.5 in, 600/in]