TransimpedanceAmplifierCalibration.swift

import simd

protocol LinearSystem {
  associatedtype Matrix
  var matrix: Matrix { get }
  
  associatedtype Solution
  func solve(matrix: Matrix) -> Solution
}

// MARK: - Three Equations

struct ThreeEquationsCase: LinearSystem {
  var Rf: Float = .zero
  var Cf: Float = .zero
  var C2: Float = .zero
  var C3: Float = .zero
  
  var matrix: simd_float3x3 {
    let column1 = SIMD3<Float>(
      0.0,
      4.45 / (Rf * Cf),
      3.42 / (Rf * Cf))
    let column2 = SIMD3<Float>(
      0.0,
      0.0,
      3.42 / (Rf * Cf))
    let column3 = SIMD3<Float>(
      5.59 / (Rf * Cf),
      0.0,
      0.0)
    
    return simd_float3x3(
      column1, column2, column3)
  }
  
  struct Solution {
    var R1: Float = .zero
    var R2: Float = .zero
    var R3: Float = .zero
  }
  
  func solve(matrix: simd_float3x3) -> Solution {
    // Construct the right-hand side.
    let rhs = SIMD3<Float>(
      1.0 / C2,
      1.0 / C3,
      1.0 / C3)
    
    // Solve for the unknowns.
    let x = matrix.inverse * rhs
    
    // Convert to the output format.
    var output = Solution()
    output.R1 = x[0]
    output.R2 = x[1]
    output.R3 = x[2]
    return output
  }
}

// MARK: - Four Equations

struct FourEquationsCase: LinearSystem {
  var Rf: Float = .zero
  var Cf: Float = .zero
  var C2: Float = .zero
  var zeroPoleFrequencyRatio: Float = 1.00
  
  var matrix: simd_float4x4 {
    let column1 = SIMD4<Float>(
      0.0,
      4.45 / (Rf * Cf),
      3.42 / (Rf * Cf),
      0.0)
    let column2 = SIMD4<Float>(
      0.0,
      0.0,
      3.42 / (Rf * Cf),
      zeroPoleFrequencyRatio / (Rf * Cf))
    let column3 = SIMD4<Float>(
      5.59 / (Rf * Cf),
      0.0,
      0.0,
      zeroPoleFrequencyRatio / (Rf * Cf))
    let column4 = SIMD4<Float>(
      0.0,
      -1.0,
      -1.0,
      0.0)
    
    return simd_float4x4(
      column1, column2, column3, column4)
  }
  
  struct Solution {
    var R1: Float = .zero
    var R2: Float = .zero
    var R3: Float = .zero
    var C3: Float = .zero
  }
  
  func solve(matrix: simd_float4x4) -> Solution {
    // Construct the right-hand side.
    let rhs = SIMD4<Float>(
      1.0 / C2,
      0.0,
      0.0,
      1.0 / C2)
    
    // Solve for the unknowns.
    let x = matrix.inverse * rhs
    
    // Convert to the output format.
    var output = Solution()
    output.R1 = x[0]
    output.R2 = x[1]
    output.R3 = x[2]
    output.C3 = 1 / x[3]
    return output
  }
}

// MARK: - Scripting

// The chief independent variable controlling the others.
let C2: Float = 3e-9

struct Combinations {
  let Rf: [Float] = [330e6 * 0.95, 330e6 * 1.00, 330e6 * 1.05]
  let Cf: [Float] = [80e-15, 150e-15, 220e-15]
  let zeroPoleFrequencyRatio: [Float] = [1.000, 1.075]
}
let combinations = Combinations()

struct Distribution {
  var minimum: Float = .greatestFiniteMagnitude
  var maximum: Float = -.greatestFiniteMagnitude
  
  var sum: Float = .zero
  var populationSize: Int = 0
  var average: Float {
    sum / Float(populationSize)
  }
  
  mutating func accumulate(_ input: Float) {
    if input < minimum {
      minimum = input
    }
    if input > maximum {
      maximum = input
    }
    
    sum += input
    populationSize += 1
  }
  
  func repr() -> String {
    """
    minimum = \(minimum)
    average = \(average)
    maximum = \(maximum)
    """
  }
}
var distributionR1 = Distribution()
var distributionR2 = Distribution()
var distributionR3 = Distribution()
var distributionC3 = Distribution()

for Rf in combinations.Rf {
  for Cf in combinations.Cf {
    for zeroPoleFrequencyRatio in combinations.zeroPoleFrequencyRatio {
      var system = FourEquationsCase()
      system.Rf = Rf
      system.Cf = Cf
      system.C2 = C2
      system.zeroPoleFrequencyRatio = zeroPoleFrequencyRatio
      
      //        print()
      //        print("Rf =", system.Rf / 1e6, "MΩ")
      //        print("Cf =", system.Cf / 1e-12, "pF")
      //        print("C2 =", system.C2 / 1e-12, "pF")
      //        print("zero-pole ratio =", system.zeroPoleFrequencyRatio)
      
      let matrix = system.matrix
      let solution = system.solve(matrix: matrix)
      
      //        print("----------")
      //        print("R1 =", solution.R1)
      //        print("R2 =", solution.R2)
      //        print("R3 =", solution.R3)
      //        print("C3 =", solution.C3 / 1e-12, "pF")
      
      distributionR1.accumulate(solution.R1)
      distributionR2.accumulate(solution.R2)
      distributionR3.accumulate(solution.R3)
      distributionC3.accumulate(solution.C3)
    }
  }
}

print()
print("R1")
print(distributionR1.repr())

print()
print("R2")
print(distributionR2.repr())

print()
print("R3")
print(distributionR3.repr())

print()
print("C3")
print(distributionC3.repr())

// Constraints on bandwidth to ensure loop stability
do {
  let lowFrequencyCriterion = 1 / (2 * Float.pi * C2 * 5.1e3)
  let highFrequencyCriterion = 1 / (2 * Float.pi * C2 * 61e3)
  print()
  print(highFrequencyCriterion, "< R3 <<", lowFrequencyCriterion)
}

// Possible maximum currents demanded by compensation network
do {
  let amplifierMaximumSlewRate: Float = 17e6
  let carrierWaveMaximumSlewRate: Float = 0.6e6 // 10 kHz, 10 V
  
  let amplifierCurrent = amplifierMaximumSlewRate * C2
  let carrierWaveCurrent = carrierWaveMaximumSlewRate * C2
  let resistorCurrent = 10.0 / distributionR3.minimum
  print()
  print("i = \(amplifierCurrent / 1e-3) mA")
  print("i = \(carrierWaveCurrent / 1e-3) mA")
  print("i = \(resistorCurrent / 1e-3) mA")
}