title: Model I-V in Python
body: <p>Model I-V.</p>
  
  <p>Method:
  Perform an integral, as a function of E, which outputs Current for each Voltage value used. This is repeated for an array of v_values. The equation can be found below.</p>
  
  <p><img src="https://i.stack.imgur.com/5Dvz3.png" alt="enter image description here"></p>
  
  <p>Although the limits in this equation range from <code>-inf</code> to <code>inf</code>, the limits must be restricted so that (E+eV)^2-\Delta^2>0 and E^2-\Delta^2>0, to avoid poles. (\Delta_1 = \Delta_2). Therefore there are currently two integrals, with limits from <code>-inf</code> to <code>-gap-e*v</code>  and <code>gap</code> to <code>inf</code>. </p>
  
  <p>However, I keep returning a <code>math range error</code> although I believe I have excluded the troublesome E values by using the limits stated above. Pastie of errors: <a href="http://pastie.org/private/o3ugxtxai8zbktyxtxuvg" rel="nofollow noreferrer">http://pastie.org/private/o3ugxtxai8zbktyxtxuvg</a></p>
  
  <p>Apologies for the vagueness of this question. But, can anybody see obvious mistakes or code misuse? </p>
  
  <p>My attempt:</p>
  
  <pre><code>from scipy import integrate
  from numpy import *
  import scipy as sp
  import pylab as pl
  import numpy as np
  import math
  
  e = 1.60217646*10**(-19)
  r = 3000
  gap = 400*10**(-6)*e
  g = (gap)**2
  t = 0.02
  k = 1.3806503*10**(-23)
  kt = k*t
  
  v_values = np.arange(0,0.001,0.0001)
  
  I=[]
  for v in v_values:
      val, err = integrate.quad(lambda E:(1/(e*r))*(abs(E)/np.sqrt(abs(E**2-g)))*(abs(E+e*v)/(np.sqrt(abs((E+e*v)**2-g))))*((1/(1+math.exp((E+e*v)/kt)))-(1/(1+math.exp(E/k*t)))),-inf,(-gap-e*v)*0.9)
      I.append(val)
  I = array(I)
  
  I2=[]
  for v in v_values:
      val2, err = integrate.quad(lambda E:(1/(e*r))*(abs(E)/np.sqrt(abs(E**2-g)))*(abs(E+e*v)/(np.sqrt(abs((E+e*v)**2-g))))*((1/(1+math.exp((E+e*v)/kt)))-(1/(1+math.exp(E/k*t)))),gap*0.9,inf)
      I2.append(val2)
  I2 = array(I2)
  
  I[np.isnan(I)] = 0
  I[np.isnan(I2)] = 0
  
  pl.plot(v_values,I,'-b',v_values,I2,'-b')
  pl.show()
  </code></pre>