title: How code an integration formula using Python
body: <p>I have an integration equations to calculate key rate and need to convert it into Python. </p>
  
  <p><strong>The equation to calculate key rate is given by:</strong></p>
  
  <p><a href="https://i.stack.imgur.com/HLBeP.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/HLBeP.png" alt="enter image description here"></a></p>
  
  <p><strong>where R(n) is:</strong></p>
  
  <p><a href="https://i.stack.imgur.com/KaGyS.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/KaGyS.png" alt="enter image description here"></a></p>
  
  <p><strong>and p(n)dn is:</strong></p>
  
  <p><a href="https://i.stack.imgur.com/tAmUm.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/tAmUm.png" alt="enter image description here"></a></p>
  
  <p><strong>The key rate should be plotted like this:</strong></p>
  
  <p><a href="https://i.stack.imgur.com/QZ1BP.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/QZ1BP.png" alt="enter image description here"></a></p>
  
  <p>I have sucessfully plotted the static model of the graph using following equation:</p>
  
  <pre><code>import numpy as np
  import math
  from math import pi,e,log
  import matplotlib.pyplot as plt
  
  n1=np.arange(10, 55, 1)
  n=10**(-n1/10)
  
  Y0=1*(10**-5)
  nd=0.25
  ed=0.03
  nsys=nd*n
  QBER=((1/2*Y0)+(ed*nsys))/(Y0+nsys)
  H2=-QBER*np.log2(QBER)-(1-QBER)*np.log2(1-QBER)
  Rsp=np.log10((Y0+nsys)*(1-(2*H2)))
  print (Rsp)
  
  plt.plot(n1,Rsp)
  plt.xlabel('Loss (dB)')
  plt.ylabel('log10(Rate)')
  plt.show()
  </code></pre>
  
  <p>However, I failed to plot the R^ratewise model. This is my code:</p>
  
  <pre><code>import numpy as np
  import matplotlib.pyplot as plt
  
  def h2(x):
      return -x*np.log2(x)-(1-x)*np.log2(1-x)
  
  e0=0.5
  ed=0.03
  Y0=1e-5
  nd=0.25
  
  nt=np.linspace(0.1,0.00001,1000)
  y=np.zeros(np.size(nt))
  Rate=np.zeros(np.size(nt))
  eta_0=0.0015
  
  
  for (i,eta) in enumerate(nt):
      nsys=eta*nd
      sigma=0.9
      y[i]=1/(eta*sigma*np.sqrt(2*np.pi))*np.exp(-(np.log(eta/eta_0)+(1/2*sigma*sigma))**2/(2*sigma*sigma))
      Rate[i]=(max(0.0,(Y0+nsys)*(1-2*h2(min(0.5,(e0*Y0+ed*nsys)/(Y0+nsys))))))*y[i]
  
  
  plt.plot(nt,np.log10(Rate))
  plt.xlabel('eta')
  plt.ylabel('Rate')
  plt.show()
  </code></pre>
  
  <p>Hopefully that anyone can help me to code the key rate with integration p(n)dn as stated above. This is the paper for referrence:
  <a href="https://arxiv.org/pdf/1712.08949.pdf" rel="nofollow noreferrer">key rate</a></p>
  
  <p>Thank you.</p>