#!/usr/bin/python
  # -*- coding: utf-8 -*-
  
  '''TF of FP cavity in reflection
  
                                     ----------
   X(s)          -----               | -tau⋅s |
  -------------->|1-R|----->(X)----->|e       |------
         |       -----       ^       ----------     |
         |                   |                      |
         |                   |                      V
       ----                ----                   ----
       |-R|                |-R|                   |-R|
       ----                ----                   ----
         |                   ^                      |
         |                   |       ----------     |
   Y(s)  V       -----       |       | -tau⋅s |     |
  <-----(X)<-----|1-R|<--------------|e       |<-----
                 -----               ----------
  
  '''
  
  '''import matplotlib as mpl
  mpl.use('pgf')
  
  def figsize(scale):
      fig_width_pt = 469.755                          # Get this from LaTeX using \the\textwidth
      inches_per_pt = 1.0/72.27                       # Convert pt to inch
      golden_mean = (5.0**0.5-1.0)/2.0            # Aesthetic ratio (you could change this)
      fig_width = fig_width_pt*inches_per_pt*scale    # width in inches
      fig_height = fig_width*golden_mean              # height in inches
      fig_size = [fig_width,fig_height]
      return fig_size
  
  pgf_with_latex = {                      # setup matplotlib to use latex for output
      "pgf.texsystem": "pdflatex",        # change this if using xetex or lautex
      "text.usetex": True,                # use LaTeX to write all text
      "font.family": "serif",
      "font.serif": [],                   # blank entries should cause plots to inherit fonts from the document
      "font.sans-serif": [],
      "font.monospace": [],
      "axes.labelsize": 10,               # LaTeX default is 10pt font.
      "text.fontsize": 10,
      "legend.fontsize": 8,               # Make the legend/label fonts a little smaller
      "xtick.labelsize": 8,
      "ytick.labelsize": 8,
      "figure.figsize": figsize(0.9),     # default fig size of 0.9 textwidth
      "pgf.preamble": [
          r"\usepackage[utf8x]{inputenc}",    # use utf8 fonts becasue your computer can handle it :)
          r"\usepackage[T1]{fontenc}",        # plots will be generated using this preamble
          ]
      }
  mpl.rcParams.update(pgf_with_latex)
  '''
  from numpy import *
  import matplotlib.pyplot as plt
  
  r = 0.99998
  L = 140e-3
  c = 299792458
  
  tau = L/c
  fsr = 1/(2*tau)
  finesse = pi*r/(1-r**2)
  print(' Finesse = %f
   FSR = %f
   FSR/Finesse = %f'%(finesse, fsr, fsr/finesse))
  
  f0 = (193e12//fsr)*fsr
  f = linspace(f0-0.1*fsr, f0+0.1*fsr, int(1e5))
  
  fm = 30e6
  
  ff = linspace(1, 1e6, int(1e5))
  
  #TF of cavity in reflection
  def Fr(w):
      s = 1j*w
      return r*(exp(-2*tau*s)-1)/(1-r**2*exp(-2*tau*s))
  
  #TF of cavity frequency-to-signal
  def Ff2s(w):
      s = 1j*w
      return (1-exp(-2*tau*s))/(2*tau*s)*(1-r**2)/(1-r**2*exp(-2*tau*s))
  
  
  G = Fr(2*pi*f)
  Glsb = Fr(2*pi*(f-fm))
  Gusb = Fr(2*pi*(f+fm))
  
  Gf2s = Ff2s(2*pi*ff)
  
  #plot setup
  fig, axarr = plt.subplots(3, 3, sharex='col')
  
  #plot bode of F(w)
  axarr[0, 1].set_ylabel('Magnitude')
  #axarr[0, 1].set_xlabel('Frequency (Hz)')
  axarr[0, 1].plot(f, abs(G))
  axarr[0, 1].grid()
  axarr[1, 1].set_ylabel('Phase (deg)')
  #axarr[1, 1].set_xlabel('Frequency (Hz)')
  axarr[1, 1].plot(f, angle(G, deg = True))
  axarr[1, 1].grid()
  
  #plot bode of Ff2s(w)
  axarr[1, 2].set_ylabel('Magnitude')
  #axarr[1, 2].set_xlabel('Frequency (Hz)')
  axarr[1, 2].set_xscale('log')
  axarr[1, 2].plot(ff, 20*log10(abs(Gf2s)))
  axarr[1, 2].grid()
  axarr[2, 2].set_ylabel('Phase (deg)')
  axarr[2, 2].set_xlabel('Frequency (Hz)')
  axarr[2, 2].set_xscale('log')
  axarr[2, 2].xaxis.set_label_coords(0.5, -0.3)
  axarr[2, 2].plot(ff, angle(Gf2s, deg = True))
  axarr[2, 2].grid()
  
  
  #PDH signal
  iG = (G*Gusb.conjugate()-G.conjugate()*Glsb).imag
  axarr[2, 1].set_ylabel('epsillon normalized')
  axarr[2, 1].set_xlabel('Frequency (Hz)')
  axarr[2, 1].xaxis.set_label_coords(0.5, -0.3)
  axarr[2, 1].plot(f, iG)
  axarr[2, 1].grid()
  
  
  #nyquist
  axarr[2, 0].set_ylabel('$j\omega$')
  axarr[2, 0].set_xlabel('$\sigma$')
  axarr[2, 0].plot(real(G), imag(G))
  axarr[2, 0].grid()
  axarr[2, 0].axis('equal')
  axarr[2, 0].set_xlim([-1.1, 1.1])
  
  #layout
  axarr[0, 0].axis('off')
  axarr[1, 0].axis('off')
  axarr[0, 2].axis('off')
  fig.subplots_adjust(wspace=0.5)
  plt.tight_layout()
  #plt.savefig('fig.pgf')
  plt.show()
  
  '''results
  
           ⎛      -2.0⋅ⅈ⋅L⋅ω ⎞
           ⎜      ───────────⎟
           ⎜           c     ⎟
         r⋅⎝-1 + ℯ           ⎠
  F(w) = ─────────────────────
               -2.0⋅ⅈ⋅L⋅ω
               ───────────
            2       c
         - r ⋅ℯ            + 1
  
  '''