#!/usr/bin/env texlua

-- $Id: smithchart.lua 27 2024-07-26 18:01:43Z reinhard $

-- Copyright (C) 2020 Reinhard Kotucha <reinhard.kotucha@web.de>
-- 
-- You may freely use, modify, and/or distribute this file.

require('luagraphlib')
luagraphlib.import()

local fmt = string.format
local pi  = math.pi
local sin = math.sin
local cos = math.cos
local arc_tan  = math.atan
local arc_sin  = math.asin
local arc_cos  = math.acos
local pow = math.pow
local sqrt= math.sqrt
local deg = math.deg
local rad = math.rad

require('complex')
local Re  = complex.Re
local Im  = complex.Im
local conj = complex.conj
local add = complex.add
local sub = complex.sub
local mul = complex.mul
local div = complex.div
local abs = complex.abs
local cfmt   = complex.cfmt
local cprint = complex.cprint

local function abs_square (c)
  return Re(mul(c, conj(c)))
end

local opt = {}

local function Γ (z)
  -- reflection coefficient
  return div(add(z, -1), add(z, 1))
end

local function show_Gamma (x,y)
  gsave()
    setstrokecolor{0, 0.5, 0}
    moveto(x*size,y*size)
    circle(1)
    stroke()
  grestore()
end
  

local function markzero ()
  gsave()
    translate(0,0)
    setstrokecolor{.5, 0, 0}
    circle(pt(3))
    stroke()
  grestore()
end


local function arcspec (decoration)
  local arcs = {}
  if not decoration then decoration = 'simple' end
  
  arcs = {
    logo = {
      intersections = {
        { 0.37,  3 },
        { 1,     3 },
        { 3,     500 },
      },
      labels = {},
    },
    simple = {
      intersections = {
        { 0.2,  2 },
        { 0.5,  2,},
        { 1,    5 },
        { 2,    5 },
        { 5,   10 },
        { 10, 500 },
      },
      labels = { 0.2, 0.5, 1, 2, 5, 10 },
    },
    lores = {
      intersections = {
        { 0.05,  0.2 },
        { 0.1,  1 },
        { 0.15,  0.2 },
        { 0.2,  2 },
        { 0.3,  1 },
        { 0.4,  1 },
        { 0.6,  1,},
        { 0.7,  1,},
        { 0.8,  1,},
        { 0.9,  1,},
        { 1.2,  2,},
        { 1.4,  2,},
        { 1.6,  2,},
        { 1.8,  2,},
        
        { 0.5,  5,},
        
        { 1,   10 },
        { 2,   10 },
        { 3,    5 },
        { 4,    5 },
        { 5,   10 },
        { 10, 500 },
      },
      labels = { 0.2, 0.5, 1, 2, 5, 10 },
    },
  }
  return arcs
end


local function draw_X_arc (z1, z2, size)
  local x_center = 1
  local y_center = 1/Im(z1)
  local radius = abs(y_center)

  local X1 = Re(Γ(z1)) - 1
  local X2 = Re(Γ(z2)) - 1

  local Y1, Y2
  
  if Im(Γ(z1)) > 0 then 
    Y1 = Im(Γ(z1)) - radius
  else
    Y1 = Im(Γ(z1)) + radius
  end
  
  if Im(Γ(z2)) > 0 then 
    Y2 = Im(Γ(z2)) - radius
  else
    Y2 = Im(Γ(z2)) + radius
  end

  local _,𝜑1 = polar(X1, Y1)
  local _,𝜑2 = polar(X2, Y2)

  if 𝜑2 - 𝜑1 < 0 then 𝜑2 = 𝜑2 + 360 end
  if 𝜑1 < 0 then 𝜑1 = 𝜑1 + 360 end 

  local r = radius*size

  gsave()
    translate (x_center*size, y_center*size)
    gsave()
      moveto (cartesian(r, 𝜑1))
      arc (0, 0, r, 𝜑1, 𝜑2)
      stroke()
    grestore()
  grestore()
end


local function draw_R_arc (z1, z2, size)

  -- We have to set imaginary part of impedance to 0
  local radius = (1 - Re(Γ{Re(z1),0})) / 2
  local x_center = 1 - radius
  
  -- translate from Gamma space to circle space 
  local X1 = Re(Γ(z1)) - x_center
  local X2 = Re(Γ(z2)) - x_center
  
  local Y1 = Im(Γ(z1))
  local Y2 = Im(Γ(z2))
  
  local _, 𝜑1 = polar(X1, Y1)
  local _, 𝜑2 = polar(X2, Y2)

  if 𝜑2 - 𝜑1 < 0 then 𝜑2 = 𝜑2 + 360 end
  if 𝜑1 < 0 then 𝜑1 = 𝜑1 + 360 end 
  
  local r = radius*size

  gsave()
    translate ((1 - radius)*size, 0)
    moveto (cartesian(radius*size, 𝜑1))
    arc (0, 0, r, 𝜑1, 𝜑2)
    stroke()
  grestore()
end


function draw_smithchart (size)
  local default = default.smithchart

  local arclinewidth    = default.arclinewidth or pt(0.4)
  local circlecolor     = default.circlecolor or {0}
  local circlelinewidth = default.circlelinewidth or pt(0.6)
  
  setlinewidth(arclinewidth)
  setlinejoin('round')
  setlinecap('butt')

  gsave()
    setstrokecolor(circlecolor)
    setlinewidth(circlelinewidth)
    moveto(0,0)
    circle(size) -- outer circle -- Re(z) = 0
    stroke()
  grestore()

  local function draw_X_arcs (X, R, size)
    draw_X_arc ({0,  X}, {R,  X}, size)
    draw_X_arc ({R, -X}, {0, -X}, size)
  end

  local function draw_R_arcs (R, X, size)
    draw_R_arc ({R, X}, {R, -X}, size)
  end
  
  local reso = default.reso or 'simple'
  local arcs = arcspec()
  
  for _,v in ipairs(arcs[reso].intersections) do
    draw_X_arcs(v[1], v[2], size)
    draw_R_arcs(v[1], v[2], size)
  end
end


function smithchart_Z (attr)
  if not attr then attr={} end
  local default = default.smithchart
  local radius  = attr.radius or default.radius or 50
  local layer = attr.layer or default.layer or 'fg'
  local arccolor = default.Z[layer].arccolor
  local arclinewidth = attr.arclinewidth or default.arclinewidth or pt(0.4)
  
  gsave()
    setstrokecolor(arccolor)
    setlinewidth(arclinewidth)

    gsave() -- Im(z) = 0
      moveto (-radius, 0)
      lineto (radius, 0)
      stroke()
    grestore()

    draw_smithchart(radius)
  grestore()
end


function smithchart_Y (attr)
  if not attr then attr={} end
  local default = default.smithchart

  local radius = attr.radius or default.radius or 50
  local layer = attr.layer or default.layer or 'fg'
  local arccolor = default.Y[layer].arccolor
  local arclinewidth = attr.arclinewidth or default.arclinewidth or pt(0.4)
  
  gsave()
    setstrokecolor(arccolor)
    setlinewidth(arclinewidth)
    
    gsave() -- Im(z) = 0
      moveto (-radius, 0)
      lineto (radius, 0)
      stroke()
    grestore()

    rotate(180)
    draw_smithchart(radius)
  grestore()
end


function smithlabels_Z (attr)
  if not attr then attr={} end
  local default = default.smithchart
  
  local radius = attr.radius or default.radius or 40
  local fillcolor = attr.fillcolor or default.labels.fillcolor or {1}
  local strokecolor = attr.strokecolor or default.labels.strokecolor or {0}

  -- Booleans cannot be handled this way because nil is treated like false.
  local framed = true
  if attr.framed == false then
    framed = false
  elseif default.labels.framed == false then
    framed = false
  end

  local reso = 'simple'
  local arcs = arcspec()
  local labels = arcs[reso].labels 

  for _, R in ipairs (labels) do
    local x = Re(Γ{R, 0})

    gsave()
      translate(x*radius-1.3, 0)
      rotate(90)
      text(R,
           { sep = pt(2),
             scale = 1,
             fillcolor = {1},
             strokecolor = {1},
             linewidth = 0,
             framed = framed,
           }
      )
    grestore()
  end
  
  for _, X in ipairs (labels) do
    local x = Re(Γ{0, X})
    local y = Im(Γ{0, X})
    local _, phi = polar(x, y)

    gsave()
      translate(x*radius, y*radius)
      rotate(phi)
      translate(-1,1.2)
      text(X..'i',
           { sep = pt(2),
             scale = 1,
             align = 'right',
             fillcolor = {1},
             strokecolor = {1},
             linewidth = 0,
             framed = framed,
           }
      )
    grestore()
  end
  
  for _, X in ipairs (labels) do
    local x = Re(Γ{0, -X})
    local y = Im(Γ{0, -X})
    local _, phi = polar(x,y)

    gsave()
      translate(x*radius, y*radius)
      rotate(phi)
      translate(-1,-3.2)
      text('$-$'..X..'i',
           { sep = pt(2),
             scale = 1,
             align = 'right',
             fillcolor = {1},
             strokecolor = {1},
             linewidth = 0,
             framed = framed,
           }
      )
    grestore()
  end
end


function smithlabels_Y (size)
  if not attr then attr={} end
  local default = default.smithchart
  
  local radius = attr.radius or default.radius or 40
  local fillcolor = attr.fillcolor or default.labels.fillcolor or {1}
  local strokecolor = attr.strokecolor or default.labels.strokecolor or {0}

  -- Booleans cannot be handled this way because nil is treated like false.
  local framed = true
  if attr.framed == false then
    framed = false
  elseif default.labels.framed == false then
    framed = false
  end

  local reso = 'simple'
  local arcs = arcspec()
  local labels = arcs[reso].labels 

  for _, R in ipairs (labels) do
    local x = Re(Γ{R, 0})

    gsave()
      rotate(180)
      translate(x*radius-1.3, 0)
      rotate(90)
      text(R,
           { sep = pt(2),
             scale = 1,
             fillcolor = {1},
             strokecolor = {1},
             linewidth = 0,
             framed = true,
           }
      )
    grestore()
  end
  
  for _, X in ipairs (labels) do
    local x = Re(Γ{0, X})
    local y = Im(Γ{0, X})
    local _, phi = polar(x, y)

    gsave()
      rotate(180)
      translate(x*radius, y*radius)
      rotate(phi)
      translate(-1,1.2)
      text(X..'i',
           { sep = pt(2),
             scale = 1,
             align = 'right',
             fillcolor = {1},
             strokecolor = {1},
             linewidth = 0,
             framed = true
           }
      )
    grestore()
  end
  
  for _, X in ipairs (labels) do
    local x = Re(Γ{0, -X})
    local y = Im(Γ{0, -X})
    local _, phi = polar(x,y)

    gsave()
      rotate(180)
      translate(x*radius, y*radius)
      rotate(phi)
      translate(-1,-3.2)
      text('$-$'..X..'i',
           { sep = pt(2),
             scale = 1,
             align = 'right',
             fillcolor = {1},
             strokecolor = {1},
             linewidth = 0,
             framed = true
           }
      )
    grestore()
  end
end


function z_to_Gamma (re, im)
  local Gamma = Γ { re, im }
  return Re(Gamma), Im(Gamma)
end


function Z_to_Gamma (re, im, Z_0)
  local Gamma = Γ { re/Z_0, im/Z_0 }
  return Re(Gamma), Im(Gamma)
end


function smithdata (filename, ...)
  local attr = {}
  local fun
  
  for _,p in ipairs{...} do
    if type(p) == 'table' then
      attr = p
    end
    if type(p) == 'function' then
      fun = p
    end
  end

  local default = default.smithchart

  if type(fun) ~= 'function' then
    fun = default.lambda or function (a) return a[2], a[3] end
  end

  local radius = attr.radius or default.radius or 30
  local plotcolor = attr.plotcolor or default.plotcolor or {0, 1, 0}
  local plotlinewidth = attr.plotlinewidth or default.plotlinewidth or pt(1)
  
  gsave()
    to_path(readdatafile (filename, fun, attr))

    setstrokecolor(plotcolor)
    setlinewidth(plotlinewidth)
    setlinecap('round')
    setlinejoin('round')
    scalepath(radius,radius)
    stroke()
  grestore()
end


function smithbackground (attr)
  if not attr then attr={} end
  
  local radius  = default.smithchart.radius  or 40
  local bgcolor = attr.bgcolor or default.smithchart.bgcolor or {1}

  local inner = attr.inner or 0
  local outer = attr.outer or 1

  if attr.inner and attr.inner < 0 then
    inner = -math.pow(10,attr.inner/20)
    if attr.outer == 0 then
      attr.outer = -1e-12
    end
  end

  if attr.outer and attr.outer < 0 then
    outer = -math.pow(10,attr.outer/20)
  end

  gsave()
    setfillcolor(bgcolor)
    moveto(0,0)
    circle (radius*outer)
    moveto(0,0)
    circle(radius*inner)
    eofill()
  grestore()
end


function smithgamma (attr)
  -- print abs(Γ) circle[s]. attr.Gamma can be either a single value
  -- or a list.  In the latter case multiple circles are printed.  If
  -- values are negative they are interpreted as return loss in dB.
  
  if not attr then attr={} end
  
  local radius  = default.smithchart.radius  or 40

  local strokecolor = attr.strokecolor or
    default.smithchart.Gamma.strokecolor or {1}

  local linewidth = attr.linewidth or
    default.smithchart.Gamma.linewidth or {1}

  local Gamma = attr.Gamma or
    default.smithchart.Gamma.Gamma or .5

  local r = {} -- radius of circles

  if type(Gamma) == 'table' then
    r = Gamma
  elseif
    type(Gamma) == 'number' then
    r[1] = Gamma
  end

  for _,R in ipairs(r) do
    if R < 0 then
      R = -math.pow(10,R/20)
    end

    gsave()
      setlinewidth(linewidth)
      setstrokecolor(strokecolor)
      moveto(0,0)
      circle (radius*R)
      stroke()
    grestore()
  end
end


local function circle_from_3_points (c1, c2, c3)
  -- Return center and radius of a circle through 3 points.  Arguments
  -- are complex numbers.  Return values: center (complex), radius
  -- (real).
  -- 
  -- https://math.stackexchange.com/questions/213658/\
  --    get-the-equation-of-a-circle-when-given-3-points
  --
  -- by Prof. David Wright of Oklahoma State University

  --[[ Matlab
    function circle_from_3_points(z1, z2, z3)
    w = (z3 - z1)/(z2 - z1);
    if isreal(w)
      disp('Points are collinear');
    else
      c = (z2 - z1)*(w - abs(w)^2)/(2i*imag(w)) + z1;
      r = abs(z1 - c);
      disp(c);
      disp(r);
    end
  end

  circle_from_3_points(1 + 1i, 2 + 4i, 5 + 3i);
  3.0000 + 2.0000i
  2.2361
  --]]

  local w = div(sub(c3, c1), sub(c2, c1))

  if Im(w) == 0 then
    error ('Points are collinear')
  else
    local denom = { 0, 2*Im(w) }
    local center = add(div(mul(sub(c2 ,c1), sub(w, abs_square(w))), denom), c1)
    local radius = abs(sub(c1, center))
    return center, radius
  end
end


  
function smithQarcs (attr)
  if not attr then attr={} end
  local default = default.smithchart
  
  local smithradius = attr.radius or default.radius or 40
--  local fillcolor = attr.fillcolor or default.Qarc.fillcolor or {1}
  local fillcolor = attr.fillcolor  or {1}
--  local strokecolor = attr.strokecolor or default.Qarc.fillcolor or {0}
  local strokecolor = attr.strokecolor  or {0}
  local smithradius = default.radius or 40

  for _,Q in pairs (attr.Q) do
  
    local z_Q = {}
    -- In order to construct arcs of constant Q we have to determine
    -- the origin and radius of a circle from three points on its
    -- circumference.  One point is at Γ = -1 + 0i, another one is at
    -- Γ = 1 + 0i.  The third point should be close to the imaginary
    -- axis in order to avoid rounding errors.
    --
    -- Q = Im(z)/Re(z)
    --
    -- In order to stay close to the imaginary axis we keep Re(z)
    -- constant for Q < 1 and Im(z) otherwise.  Then Re(Γ) is always
    -- beween 0 and 0.2.
    
    if Q < 1 then
      z_Q = { 1, Q }
    else
      z_Q = { 1/Q, 1 }
    end
    
    local Gamma_Q = Γ(z_Q)
    local center, radius = circle_from_3_points({-1,0}, Gamma_Q, {1,0})
    local Qradius = smithradius * radius
    local Qx = smithradius * center[1]
    local Qy = smithradius * center[2]
    local alpha = deg(math.asin(1/radius))
  
    gsave()
      translate(0, Qy)
      moveto(smithradius, -Qy)
      setlinewidth(pt(0.4))
      arc(0, 0, Qradius, 90 - alpha, 90 + alpha)
      gsave()
      setstrokecolor(strokecolor)
      stroke()
      grestore()
      
      translate(0, Qradius+1)
      text('$Q='..Q..'$',
           { sep = pt(2),
             scale = 1,
             align = 'center',
             fillcolor = {1},
             strokecolor = {1},
             linewidth = 0,
             framed = false
           }
      )
    grestore()

  -- Same for negative imaginary values
    gsave()
      rotate(180)
      translate(0, Qy)
      moveto(smithradius, -Qy)
      setlinewidth(pt(0.4))
      arc(0, 0, Qradius, 90 - alpha, 90 + alpha)
      gsave()
        setstrokecolor(strokecolor)
        stroke()
      grestore()
    grestore()
  end
end

opt.debug = true
local function mark (c, smithradius, text)
  if opt.debug then
    gsave()
      translate(Re(c)*smithradius, Im(c)*smithradius)
      gsave()
        scale(.5,.5)
        tex.print([[\MTcircle]])
      grestore()
      scale(.5,.5)
      tex.print([[\hbox to0pt{\hss ]]..text..[[\hss}]])
    grestore()
  end
end


function y (z)
  return div(1, z)
end
  
function z (y)
  return div(1, y)
end


local function arc_from_3_points (c1, c2, c3)
  
  local center, radius = circle_from_3_points (c1, c2, c3)

  mark (center, 40, 'c')

  print('arc_from_3_points():  center = '..cfmt(center,3)..
        ', radius = '..fmt('%.3g', radius))
  
  -- c0 denotes the middle of the line segment c1--c3 

  local c0 = div(add(c3, c1), 2)
  print('c0 = '..cfmt(c0, 4)) 

  -- There are two rightangled triangles:
  --   c1--c0--center and c3--c0--center.
  --
  -- They are symmetric, hence we have to consider only one.
  --
  -- We are using the first in order to determine the angles.
  --
  -- The hypotenuse is equal to the radius.  The adjacent side
  -- corresponds to the line segment c0--center and the opposite side
  -- to c1--c0.
  --
  -- The angle is determined by the math.asin() function which expects
  -- real numbers.

  -- We also need the angle between center and c0 (midangle).

  local c_midang = sub(c0, center)
  
  local midangle = (deg(arc_tan(Im(c_midang), Re(c_midang))))
  print('midangle = '..fmt('%.2f°', midangle))

  local opposite = abs(sub(c1, c0)) / radius
  
  local halfangle = deg(arc_sin(opposite))
  print('halfangle = '..fmt('%.2g°', halfangle))
  
  print('opposite = '..fmt('%.2g', opposite))
  return Re(center), Im(center), radius, midangle, halfangle
end




function smithQarea (Q, attr)
  if not attr then attr={} end
  local default = default.smithchart
  
  local smithradius = attr.radius or default.radius or 40
--  local fillcolor = attr.fillcolor or default.Qarc.fillcolor or {1}
  local fillcolor = attr.fillcolor  or {1}
--  local strokecolor = attr.strokecolor or default.Qarc.fillcolor or {0}
  local strokecolor = attr.strokecolor  or {0}
  local smithradius = default.radius or 40

  local z_Q = {}
  -- In order to construct arcs of constant Q we have to determine the
  -- origin and radius of a circle from three points on it's
  -- circumference.  One point is at Γ = -1 + 0i, another one is at
  -- Γ = 1 + 0i.  The third point should be close to the imaginary
  -- axis in order to avoid rounding errors.
  --
  -- Q = Im(z)/Re(z)
  --
  -- In order to stay close to the imaginary axis we keep Re(z)
  -- constant for Q < 1 and Im(z) otherwise.  Then Re(Γ) is always
  -- beween 0 and 0.2.
  
  if Q < 1 then
    z_Q = { 1, Q }
  else
    z_Q = { 1/Q, 1 }
  end
  
  local Gamma_Q = Γ(z_Q)
  local center, radius = circle_from_3_points({-1,0}, Gamma_Q, {1,0})
  local Qradius = smithradius * radius
  local Qx = smithradius * center[1]
  local Qy = smithradius * center[2]
  local alpha = deg(math.asin(1/radius))

  local function draw ()
    gsave()
      translate(0, 0)
      moveto(smithradius, 0)
      arc(0, Qy, Qradius, 90 - alpha, 90 + alpha)
      -- print()
      -- print ('Qarc args: X = '..tostring(0)..
      --    ', Y = '.. 0 ..
      --    ', R = '.. Qy ..
      --    ', a1 = '..fmt('%.2f', 90 - alpha )..
      --    ', a2 = '..fmt('%.2f', 90 + alpha ))
   
      -- print()
      arcn(0, 0, 40, 180, 0)
      closepath()
      
      setlinewidth(0)
      setfillcolor{1, .9, .9}
      setstrokecolor(strokecolor)
      paint()
    grestore()
  end

  draw()
  gsave()
    scale(-1, -1)
    draw()
  grestore()
end


local function YinYang (...)
  local attr = {}
  local default = default.smithchart

  local radius = attr.radius or default.radius or 30
  local fillcolor = attr.forbiddencolor or default.forbiddencolor or {1, .9, .9}
  local linewidth = attr.forbiddenlinewidth or default.forbiddenlinewidth or pt(1)
  
  moveto(radius, 0)
  arc(0, 0, radius, 0, 180)
  arc(-radius/2, 0, radius/2, 180, 360)
  arcn(radius/2, 0, radius/2, 180, 0)
  closepath()
  setfillcolor(fillcolor)
  setlinewidth(linewidth)
  paint()
end


function smith_zC_yL (...)
  gsave()
    YinYang {...}
  grestore()
end


function smith_zL_yC (...)
  gsave()
    scale(1, -1)
    YinYang {...}
  grestore()
end


function smith_yC_zL (...)
  gsave()
    scale(-1, 1)
    YinYang {...}
  grestore()
end


function smith_yL_zC (...)
  gsave()
    scale(-1, -1)
    YinYang {...}
  grestore()
end


function smith_yL_zL_zL_yL (...)
  local attr = {}
  
  local default = default.smithchart

  local radius = attr.radius or default.radius or 30
  local plotcolor = attr.plotcolor or default.plotcolor or {0, 1, 0}
  local plotlinewidth = attr.plotlinewidth or default.plotlinewidth or pt(1)
  local fillcolor = {1, .9, .9}
  
  gsave()
    setlinewidth(mm(1))
  
    setfillcolor(fillcolor)
    moveto(radius,0)
    arc(0, 0, radius, 0, 180)
    arc(-radius/2, 0, radius/2, 180, 360)
    arc(radius/2, 0, radius/2, 180, 360)
    closepath()
    fill()
    grestore()
end


function smith_yC_zC_zC_yC (...)
  gsave()
    scale(-1,-1)
    smith_yL_zL_zL_yL (...)
  grestore()
end



-- Local Variables:
--  mode: Lua
--  TeX-engine: luatex
--  lua-indent-level: 2
--  indent-tabs-mode: nil
--  coding: utf-8-unix
-- End:
-- vim:set tabstop=2 expandtab:


-- function arcs ()
--   local r = 10
--   gsave()
--     setfillcolor(svg.PaleTurquoise)
--     setstrokecolor(dvips.MidnightBlue)
--     moveto(0, r)
--     arc(r, r, r, 180, 270)
--     lineto(100-r, 0)
--     arc(100-r, r, r, 270, 360)
--     lineto(100, 50-r)
--     arc(100-r, 50-r, r, 0, 90)
--     lineto(r, 50)
--     arc(r, 50-r, r, 90, 180)
--     closepath()
--     paint()
--   grestore()
-- end