#!/usr/bin/env python
  
  import numpy
  import matplotlib.pyplot as plt
  
  L = 14e-2
  c = 299792458
  R = 0.99991
  
  fsr = c/(2*L)
  print("fsr = %f GHz"%(fsr/1e9))
  
  finesse = numpy.pi*R**0.25/(1-R**0.5)
  print("finesse = %i"%finesse)
  
  print("finesse/fsr = %f s"%(finesse/fsr))
  
  fc = fsr*130000
  df = 20e6
  #print(fc)
  
  f = 2*numpy.pi*numpy.linspace(fc-df, fc+df, 100000) # frequency array
  i = complex(0,1)
  
  fpt = ((1.-R)*numpy.exp(-2.*L*i*f/c))/(1.-R*numpy.exp(-2.*L*i*f/c))
  fpr = R**0.5*(1.-numpy.exp(-2.*L*i*f/c))/(1.-R*numpy.exp(-2.*L*i*f/c))
  
  plt.subplot(2, 1, 1)
  plt.plot(f, abs(fpr))
  plt.title('FP cavity')
  plt.ylabel('Reflective intensity')
  plt.grid()
  
  
  fm = 23e6
  ram = 0
  
  fm1 = R**0.5*(1.-numpy.exp(-2.*L*i*(f+fm)/c))/(1.-R*numpy.exp(-2.*L*i*(f+fm)/c)) + ram*i
  fm2 = R**0.5*(1.-numpy.exp(-2.*L*i*(f-fm)/c))/(1.-R*numpy.exp(-2.*L*i*(f-fm)/c)) + ram*i
  
  ff = fpr*fm1.conjugate()-fpr.conjugate()*fm2
  iff = ff.imag
  
  
  plt.subplot(2, 1, 2)
  plt.plot(f, iff)
  plt.xlabel('Frequency')
  plt.ylabel('PDH')
  plt.grid()
  
  
  plt.show()