"""
Created on Tue Mar 26 22:53:41 2019

@author: 150193878
"""
import numpy as np
import matplotlib.pyplot as plt
import scipy.constants as constants

T=300 #Kelvin
q=1.6e-19 #Coloumbs
kb=constants.k
beta=q/(kb*T)  #Coloumbs/Joule
V_oc=0.5 #Volts

def Power(I,V):
    if np.isclose(V,V_oc-(1/beta)*np.log(1+beta*V),atol=1e-3):
        return V*I


for I in np.linspace(0,2,500):
    for V in np.linspace(0,(1/beta)*np.log(1+beta*V_oc),1000):
        IV=(Power(I,V))

print(IV)

None