#!/usr/bin/python3

import numpy
import math
import matplotlib.pyplot as plt

##########################
Cx = 0.109e-9 # internal capacitance, F

C0 = 0.94e-9
C1 = 0.47e-9
C2 = 1.0e-9
C3 = 0.2e-9
C23 = 0.67e-9 #  added when both 2 and 3 is on

L = 35e-6
R = 10

fmin = 0.4e6 # plotting range, Hz
fmax = 1.0e6

##########################
C = numpy.zeros(8)
C[0] = C0
C[1] = C0 + C1
C[2] = C0 + C2
C[3] = C0 + C1 + C2
C[4] = C0 + C3
C[5] = C0 + C3 + C1
C[6] = C0 + C3 + C2 + C23
C[7] = C0 + C3 + C2 + C1 + C23

C = C + Cx

##########################
#F0 = 1/(2*math.pi*numpy.sqrt(C*L))

def plot_res(Cin, Cgr, Rin, col, t):

  f = numpy.linspace(fmin, fmax, 300)
  w = 2*math.pi*f;

  for i in range(8):
    Zin   = 1/(1j*w*Cin) + Rin;
    Zgr   = 1/(1j*w*Cgr)
    Zc    = 1/(1j*w*C[i])

    A = Zin / (Zin + Zc*(1j*w*L + R)/(1j*w*L + R + Zc) + Zgr);
    U = abs(1 - A);

    print(U)
    im = numpy.argmax(U)
    print(im)
    plt.text(f[im]/1e6, U[im]+0.001, "%d"%(i))

    if i!=0: t=''
    plt.plot(f/1e6, U, col, label=t)


#text(900, 1.1,...
#printf('L = %.1f mkH;  Cint = %.1f pF', L*1e6, Cin*1e12));

plot_res(1e-12, 4.7e-9, 0, 'r-', 'C_in = 1pF')
plot_res(2e-12, 4.7e-9, 0, 'b-', 'C_in = 2pF')

plt.xlabel('freq, MHz')
plt.ylabel('amplitude, a.u')
plt.legend()
plt.savefig("plot_nmr.png")