Anybody know MATLAB and Python? (Code conversion MATLAB>Python)
Question: I am trying to re-write this MATLAB program in Python. I haven't succeeded in
getting the same Python output, yet. But my attempt is given beneath the
MATLAB code. The code does not need any extra files/information to run. So
this should run OK on your MATLAB. And, if all goes to plan, Python also...

Summary of what the code does: Performs an integral taking in arguments `eV`
and `t`, for an array of variables `eV`. There are also more complicated
things considering a substitution for `E`. But, those familiar with both codes
should be able to follow.

Please feel free to ask if you have any questions, and many thanks for any
potential help/hints/solutions.

MATLAB CODE:

**Main.m**

    
    
    clear all %Remove items from MATLAB workspace and reset MuPAD engine
    clc %Clear command window
    clf %Clear figure window
    
    global d1 d2 T %Declare global variables
    
    T = 0.02; %Temperature value (K)
    d1 = 1; %Energy gap in electrode 1.
    d2 = 0.5; %Energy gap in electrode 2.
    
    small = 1e-9; 
    eV_values = linspace(0, 2.5, 2e3); %Row vector of 2e3 points linearly spaced between 0 and 0.25. These are the voltage values.
    
    current = zeros(size(eV_values)); %Zeros creates array all of zeros. size gives size of dataset array.
    
    tic %Start clock to measure performance
    for x = 1:numel(eV_values) %numel gives number of elements in array ev_values).
        eV = eV_values(x);
    
        clc %Clear command window
        disp(x) %Display array
    
        current(x) = quad(@(t)integrand(t, eV), -1 + small, 1 - small); 
    end
    toc %End clock
    
    clf %Clear figure window
    figure(1) %Create graphics object
    
    hold on %Retain current graph when adding new graphs/Delay evaluation.
    box on %Display the boundary of the current axes.
    plot(eV_values, real(current), 'b')
    plot(eV_values, imag(current), 'r')
    title('S-S')
    xlabel('eV/\Delta')
    ylabel('I(eV)')
    

**Integrand.m**

    
    
    function x = integrand(t, eV) %Declare function name and inputs
    
    global d1 d2 T %Declare global variable
    
    E = t./(1 - t.^2); %Variable substitution
    
    x = abs(E)./sqrt(E.^2 - d1^2).*abs(E + eV)./sqrt((E + eV).^2 - d2^2).*...
        (1./(1 + exp(E./T)) - 1./(1 + exp((E + eV)./T)));
    
    x = x.*heaviside(E.^2 - d1^2).*heaviside((E + eV).^2 - d2^2);
    x = x.*(1 + t.^2)./(1 - t.^2).^2;
    
    %heaviside step function
    

PYTHON CODE:

    
    
    from numpy import *
    import pylab as pl
    import array
    from scipy import integrate
    
    T = 0.02    # Global variable - Temperature (K)
    d1 = 1      # Global variable - Energy gap in electrode 1.
    d2 = 0.5    # Global variable - Energy gap in electrode 2.
    
    small = 1e-9 
    eV_values = linspace(0.0, 2.5, num=10)
    
    def heaviside( x ): 
        # Return 0 for x<0, 1 for x>0, 0.5 for x=0. #
        if x == 0:
            return 0.5
        return 0 if x < 0 else 1
    
    def integrand( t, eV ):
        #print(" t: %s, eV: %s" % (t, eV))
        E = t / ( 1 - t*t ) # E substitution.
        x1 = ( abs( E ) / sqrt( E*E - d1*d1 ) ) * ( abs( E + eV ) / ( sqrt( ( E + eV )**2 - d2*d2 ) ) ) * (1/(1 + exp(E/T)) - 1/(1 + exp((E + eV)/T)))
        x2 = x1*(heaviside( E*E - d1*d1 )*heaviside( (E + eV)**2 - d2*d2) )
        x = x2*( ( 1 + t*t ) / ( 1 - t*t )**2 )
        return x
    
    current = []
    for eV in eV_values:
        integral, err = integrate.quad( integrand, ( -1 + small ), ( 1 - small ), args=(eV, ) )
    
    #   print( eV, integral )
        print( eV, integral, err)
        current.append( integral )
    
    #print( 'current values')
    print( current )
    
    #pl.plot(eV_values,current,'b')
    #pl.plot(eV_values,imag(current),'r')
    #pl.title('S-S')
    #pl.xlabel(r'eV/$\Delta$')
    #pl.ylabel('I(eV)')
    #pl.show()
    

Notable problems:

  * The MATLAB code considers imaginary/real current values in quad. Python code currently doesn't, but ought to for obtaining the same output.
  * The current Python code outputs: Giving `nan` values for the `integral` and `err`. Again, this may be down to the program not considering imaginary and real values in the integral.

`In [3]: run IV.py IV.py:22: RuntimeWarning: invalid value encountered in sqrt
x1 = (abs(E)/sqrt(E*E - d1*d1)) * (abs(E + eV)/(sqrt((E + eV)**2 - d2*d2))) *
(1/(1 + exp(E/T)) - 1/(1 + exp((E + eV)/T))) IV.py:22: RuntimeWarning:
overflow encountered in exp x1 = (abs(E)/sqrt(E*E - d1*d1)) * (abs(E +
eV)/(sqrt((E + eV)**2 - d2*d2))) * (1/(1 + exp(E/T)) - 1/(1 + exp((E +
eV)/T))) (0.0, nan, nan) (0.27777777777777779, nan, nan) (0.55555555555555558,
nan, nan) (0.83333333333333337, nan, nan) (1.1111111111111112, nan, nan)
(1.3888888888888888, nan, nan) (1.6666666666666667, nan, nan)
(1.9444444444444446, nan, nan) (2.2222222222222223, nan, nan) (2.5, nan, nan)
[nan, nan, nan, nan, nan, nan, nan, nan, nan, nan]`


Answer: `numpy.sqrt` can't deal with minus values. Use `numpy.emath.sqrt` instead.

And because calculate with floating numbers is not precise, it's better to
return the real part of x:

    
    
    return x.real
    

In addition, calculations with numpy scalar is slow. It's better to use python
standard module `math` & `cmath` to do the calculations in `integrand()`.

