# Description --------------------------------------------------------------------------------------------------------------
# Simulation of simple diodes in the time and frequency domains.
#
# Name:        diode_i(v)_time_domain.py
# Author:      W. Hetherington, Physics, Oregon State University
# Created:     2012/01/28
# Modified:    2012
# Copyright:   (c) 2012
# License:     No restrictions
#---------------------------------------------------------------------------------------------------------------------------

from scipy import *
from matplotlib import pyplot as plt
import string as s
import os
from numpy.random import normal
from numpy import *

# Graph I(t) given a sinusoidal V(t) or graph the power spectra of V and I.

def PowerSpectrum(x, y) :
    samples = x.size    # 2500 for scope
    halfway = samples/2
    # Multiplying a waveform by a window function is important when the waveform is not exactly periodic
    # over the range, that is when point[0] != point[2499].  Numpy has hanning, blackman, bartlett, hamming and kaiser.
    window = hanning(samples)   # 2500 for scope
    fty = fft.fft(y*window)
    #fty = fft.fft(window)
    n = fty.size
    ftyabs = abs(fty[0:halfway])
    ftydb = 20.0*log10(ftyabs)
    xstep = x[1] - x[0] #float(raw_input('Enter the time per division in seconds, such as 5e-06 : '))   #scope.chan_refs[0].xinc
    print n, xstep
    ftx = fft.fftfreq(n, xstep)[0:halfway]  #/1.0e03    # in kHz
    print ('Note that the frequency axis will be in Hz.')   #kHz.')
    ftmaxindex = ftyabs.argmax()
    ftmax = ftyabs[ftmaxindex]
    print ftmax, ftmaxindex
    #t = [ftx, ftydb]    # not ftyabs
    return ftx, ftydb

def Test() :
    # Construct time array of arbitrary unit
    points = 10000
    times = arange(points, dtype=float)
    
    # Define operational parameters
    isat = 1.e-09		# ampere
    k = 8.6173324e-05		# in eV/kelvin
    #e = 1.6e-19		# coulomb
    q = 1.0			# in e
    t = 300
    alpha =q/k/t
    time_inc = 1.0e-6	# Assume that the time increment is 1 us.
    
    # There are two choices for graphing the behavior: time domain or frequency domain.
    # Choose the type of graph by setting the domain variable below.
    
    domain = 'time'
    #domain = 'freq'
    
    print('Working ...')
    if domain == 'time' :
	nup = 10.0/points
	nu = nup/time_inc
	vin = 0.41 * sin(2.0 * pi * nup * times)
	current = isat*(exp(alpha * vin) -1) * 1000.0	# in milliampere
	plt.plot(times, vin)
	plt.plot(times, current)
	plt.xlabel('Time (' + str(time_inc) +')')
	plt.ylabel('Potential (V) or Current (mA)')
	plt.title('Diode: Isat = ' + str(isat) + ', q = ' + str(q) + ' e, frequency = ' + str(nu) + ' Hz')
    elif domain == 'freq' :
	nup = 1000.0/points
	nu = nup/time_inc
	vin = 0.41 * sin(2.0 * pi * nup * times)
	current = isat*(exp(alpha * vin) -1) * 1000.0	# in milliampere
	freqs, power_spec = PowerSpectrum(times * time_inc, vin)
	plt.plot(freqs, power_spec)
	freqs, power_spec = PowerSpectrum(times * time_inc, current)
	plt.plot(freqs, power_spec)
	plt.xlabel('Frequency (Hz)')
	plt.ylabel('Power (dB)')
	plt.title('Diode: Isat = ' + str(isat) + ', q = ' + str(q) + ' e, frequency = ' + str(nu) + ' Hz')
    else :
	print('Select time or frequency domain graphs in the program.')
    
    legend = ['Potential', 'Current']
    plt.legend(legend, 'upper right')
    plt.grid()
    plt.show()
	
#-----------------------------------------------------------------------------------

if __name__ == '__main__' :
	Test()