#!/usr/bin/env python
# model shading effects on PV panels

# Copyright Andrew Tridgell 2010
# Released under GNU GPL v3 or later

import math
import numpy, pylab, scipy.optimize, sys
from optparse import OptionParser

# some panel parameters
panels = {
    # Isc = short circuit current at 25C
    # Voc = open circuit voltage at 25C
    # Imp = current at maximum power point at 25C
    # Vmp = voltage at maximum power point at 25C
    # Vt  = Voltage change per degree C (V/C)
    # It  = Current change per degree C (A/C)
    #
    #            Isc   Voc   Imp   Vmp       Vt     It
    "SPR225" : (5.87, 48.5, 5.49, 41.0, -0.1325,  3.50e-3),
    "STP180" : (5.40, 44.4, 5.05, 35.6, -0.1510,  3.50e-3),
    "BP175"  : (5.45, 43.6, 4.94, 35.4, -0.1570,  9.27e-4), 
    }

(panel_Isc, panel_Voc, panel_Imp, panel_Vmp, panel_VT, panel_IT) = panels["SPR225"]

# some inverters
inverters = {
    # only the MPP voltage range so far
    #                    MPPmin  MPPmax
    "SMC5000A"         : (246,   480),
    "SB5000TL"         : (175,   440),
    "XantrexGT5"       : (240,   550),
    "FroniusIGPlus120" : (230,   500),
    "None"             : (0,     1000),
    }

inverter_range = inverters["SMC5000A"]


def find_best(err, lower, upper):
    '''wrapper around scipy optimisation function'''
    if upper == lower:
        return lower
    if upper <= lower:
        raise Exception('lower<upper  %g %g' % (lower, upper))
    def err1(x):
        return err(x[0])
    (ret, n, err) = scipy.optimize.fmin_tnc(err1, [float(0.5*(lower+upper))],
                                            bounds=[(float(lower), float(upper))],
                                            approx_grad=True, messages=0, xtol=-1,
                                            maxfun=1000)
    return ret[0]


def find_best2(err, lower, upper):
    '''fmin_tnc is not much good for global optimums'''
    nsteps = 20
    ret = []
    stepsize = (upper-lower)/nsteps
    best_ret = -1
    for i in range(0, nsteps):
        r = find_best(err, lower + stepsize*i, lower+(stepsize*(i+1)))
        if best_ret == -1 or err(r) < best_err:
            best_ret = r
            best_err = err(r)
    return best_ret


def round(v, m):
    '''round v to a multiple of m'''
    return m*math.ceil(0.5+(v/m))

def round_down(v, m):
    '''round v down to a multiple of m'''
    return m*math.floor(v/m)


def temperature_adjust(t):
    '''adjustment for temperature sensitive coefficients. This is very rough'''
    V_addition = cell_VT  * (t-25)
    I_addition = panel_IT * (t-25)
    return (V_addition, I_addition)


def irradiance_adjust(w):
    '''adjustment for irradiance sensitive coefficients. This is based on squinting
       at the graphs in the Sunpower-225-WHT spec sheet. It should be in
       the panel type table.'''
    V_multiplier = math.log(1+(1.0e6*w))/math.log(1e6*1000)
    I_multiplier = w/1000.0
    return (V_multiplier, I_multiplier)


# a cache hack ... we evaluate the same parms a lot of times
# assumes that return only depends on function args
cacheargs_ret = {}
def cacheargs(function):
    def cacheargs_wrapper(*args):
        argsf = (args, str(function))
        if argsf in cacheargs_ret:
            return cacheargs_ret[argsf]
        else:
            ret = function(*args)
            cacheargs_ret[argsf] = ret
            return ret
    return cacheargs_wrapper

# cacheing for PV object functions
# assumes that return only depends on Isc, Voc and function args
cacheobj_ret = {}
def cacheobj(function):
    def cacheobj_wrapper(*args):
        obj=args[0]
        argsf = (args[1:], getattr(obj, 'Isc', 0), getattr(obj, 'Voc', 0), str(function))
        if argsf in cacheobj_ret:
            return cacheobj_ret[argsf]
        else:
            ret = function(*args)
            cacheobj_ret[argsf] = ret
            return ret
    return cacheobj_wrapper

@cacheargs
def Imodel(V, Voc, Isc, a0, a1, w, t):
    '''the core of the cell model for current. Very rough'''
    # adjust for temperature and irradiance
    (V_addition, I_addition) = temperature_adjust(t)
    (V_multiplier, I_multiplier) = irradiance_adjust(w)
    Voc += V_addition
    Isc += I_addition
    Voc *= V_multiplier
    Isc *= I_multiplier

    V   *= cell_scale
    Voc *= cell_scale


    if V > Voc:
        # reverse current
        return -(V-Voc)/opts.panel_Roc
    if V < 0:
        # over current
        return Isc - (V/opts.panel_Rsc)
    if V/a1 > 100:
        V = a1*100
    return Isc - (a0/1e7) * math.exp(V/a1)

@cacheargs
def Vmodel(I, Voc, Isc, a0, a1, w, t):
    '''the core of the cell model for voltage. Very rough'''
    # adjust for temperature and irradiance
    (V_addition, I_addition) = temperature_adjust(t)
    (V_multiplier, I_multiplier) = irradiance_adjust(w)
    Voc += V_addition
    Isc += I_addition
    Voc *= V_multiplier
    Isc *= I_multiplier
    
    Voc *= cell_scale

    if I >= Isc:
        # over current
        return (- (I - Isc) * opts.panel_Rsc)/cell_scale
    if I < 0:
        # reverse current
        return (Voc - (I * opts.panel_Roc))/cell_scale
    return (math.log((Isc - I)/(a0/1e7))*a1)/cell_scale


@cacheargs
def find_cell_params(Vmp, Voc, Isc, Imp, w, t):
    '''find the cell variables given the key parameters'''
    # we have 6 data points to work with, and we try
    # to fit a simple exponential based on a diode model
    #
    # build the error function for the cell model
    # by combining the current error and the power
    # error. This produces much better parameters
    # than just fitting by current

    # mp params are not adjusted yet
    (V_addition, I_addition) = temperature_adjust(t)
    (V_multiplier, I_multiplier) = irradiance_adjust(w)
    Vmp += V_addition
    Voc += V_addition
    Imp += I_addition
    Isc += I_addition
    Vmp *= V_multiplier
    Voc *= V_multiplier
    Imp *= I_multiplier
    Isc *= I_multiplier

    dataI = [ (0, Isc), (Vmp, Imp), (Voc, 0) ]
    dataP = [ (0, 0), (Vmp, Vmp*Imp), (Voc, 0) ]
    def err(a):
        ret = 0
        for (x,y) in dataI:
            I = Imodel(x, Voc, Isc, a[0], a[1], 1000, 25)
            ret += (I-y)*(I-y)
        for (x,y) in dataP:
            I = Imodel(x, Voc, Isc, a[0], a[1], 1000, 25)
            e = (I*x) - y
            ret += e*e
        return ret
    (ret, x1, x2) = scipy.optimize.fmin_tnc(err, [1, 1],
                                            bounds=[(float(0.01), float(100)),
                                                    (float(0.01), float(100))],
                                            messages=0,
                                            maxfun=10000,
                                            approx_grad=True,
                                            xtol=1e-30)
    return ret


def sumV(m, I):
    '''add up a set of PV like elements in series'''
    s = m[0]
    lastV = sumV = s.V(I)
    lastIn = (s.Voc, s.Isc)
    for i in range(1, len(m)):
        s = m[i]
        input = (s.Voc, s.Isc)
        if input != lastIn:
            lastV = s.V(I)
            lastIn = input
        sumV += lastV
    return sumV

class cell:
    '''model a single PV cell'''
    def __init__(self):
        self.a = [ 1, 1 ]
        self.temperature = opts.temperature
        self.irradiance = opts.irradiance
        self._recalc()
        
    def _recalc(self):
        self.a = find_cell_params(cell_Vmp, cell_Voc, panel_Isc, panel_Imp, self.irradiance, self.temperature)
        self.Voc = self.V(0)
        self.Isc = self.I(0)

    def I(self, V):
        return Imodel(V, cell_Voc, panel_Isc, self.a[0], self.a[1], self.irradiance, self.temperature)

    def V(self, I):
        return Vmodel(I, cell_Voc, panel_Isc, self.a[0], self.a[1], self.irradiance, self.temperature)

    def set_temperature(self, t):
        self.temperature = t
        self._recalc()

    def set_irradiance(self, w):
        self.irradiance = w
        self._recalc()

    def P(self, V):
        I = self.I(V)
        return I * V


class subpanel:
    '''model a subpanel, which is normally 24 cells'''
    def __init__(self):
        self.cells = []
        for i in range(0, opts.cells_per_subpanel):
            self.cells.append(cell())
        self._recalc()
        
    def _recalc(self, recursive=False):
        if recursive:
            for c in self.cells:
                c._recalc()
        # subpanel Voc is sum of cells
        self.Voc = 0
        for c in self.cells:
            self.Voc += c.V(0)                
            
        # subpanel Isc is min of cells
        self.Isc = self.cells[0].I(0)
        for c in self.cells:
            if c.I(0) < self.Isc:
                self.Isc = c.I(0)

    def V(self, I):
        if I >= self.Isc:
            # this implements the bypass diode
            return - (I - self.Isc) * subpanel_Rbp
        ret = sumV(self.cells, I)
        return ret

    @cacheobj
    def I(self, V):
        if V < 0:
            # bypass diode activated
            return self.Isc - (V/subpanel_Rbp)
        elif V > self.Voc:
            # reverse current
            lower = -1000*self.Isc
            upper = 0
        else:
            # normal range
            lower = 0
            upper = self.Isc
        def Verr(x):
            r = abs(V - self.V(x))
            return r
        return find_best(Verr, lower, upper)

    def P(self, V):
        I = self.I(V)
        return I * V

    def set_temperature(self, t):
        for c in self.cells:
            c.set_temperature(t)
        self._recalc()

    def set_irradiance(self, w):
        for c in self.cells:
            c.set_irradiance(w)
        self._recalc()


class panel:
    '''model a panel, which is normally 3 subpanels'''
    def __init__(self):
        self.subpanels = []
        for i in range(0, opts.num_subpanels):
            self.subpanels.append(subpanel())
        self._recalc()

    def _recalc(self, recursive=False):
        if recursive:
            for s in self.subpanels:
                s._recalc(recursive)
        # panel Voc is sum of subpanels
        self.Voc = 0
        for c in self.subpanels:
            self.Voc += c.V(0)                
            
        # panel Isc is max of cells. It is the maximum
        # not the minimum due to the effect of the bypass diodes
        # on the subpanels
        self.Isc = self.subpanels[0].I(0)
        for c in self.subpanels:
            if c.I(0) > self.Isc:
                self.Isc = c.I(0)

    def V(self, I):
        return sumV(self.subpanels, I)

    @cacheobj
    def I(self, V):
        if V < 0:
            lower = self.Isc
            upper = self.Isc*1000
        elif V < self.Voc:
            lower = 0
            upper = self.Isc
        else:
            lower = -self.Isc*1000
            upper = 0
            
        def Verr(x):
            return abs(V - self.V(x))
        return find_best(Verr, lower, upper)

    def P(self, V):
        I = self.I(V)
        return I * V
    
    def set_temperature(self, t):
        for s in self.subpanels:
            s.set_temperature(t)
        self._recalc()

    def set_irradiance(self, w):
        for s in self.subpanels:
            s.set_irradiance(w)
        self._recalc()


class string:
    '''model a string'''
    def __init__(self, numpanels):
        self.panels = []
        for i in range(0, numpanels):
            self.panels.append(panel())
        self._recalc()

    def _recalc(self, recursive=False):
        if recursive:
            for p in self.panels:
                p._recalc(recursive)
        # string Isc is max of the panel currents
        # it is max, not min, as any underperforming
        # panels will go into bypass
        self.Isc = self.panels[0].I(0)
        for c in self.panels:
            if c.I(0) > self.Isc:
                self.Isc = c.I(0)

        # string Voc is sum of panels
        self.Voc = 0
        for c in self.panels:
            self.Voc += c.V(0)                

    def V(self, I):
        if I <= 0 and opts.blocking_diodes:
            return self.Voc
        return sumV(self.panels, I)

    @cacheobj
    def I(self, V):
        if V >= self.Voc and opts.blocking_diodes:
            return 0
        if V < 0:
            lower = self.Isc
            upper = self.Isc*1000
        elif V < self.Voc:
            lower = 0
            upper = self.Isc
        else:
            lower = -self.Isc*1000
            upper = 0
        def Verr(x):
            return abs(V - self.V(x))
        return find_best(Verr, lower, upper)

    def P(self, V):
        I = self.I(V)
        return I * V

    def set_temperature(self, t):
        for p in self.panels:
            p.set_temperature(t)
        self._recalc()

    def set_irradiance(self, w):
        for p in self.panels:
            p.set_irradiance(w)
        self._recalc()


class pvarray:
    '''model an array with 1 MPPT'''
    def __init__(self, stringsizes):
        self.strings = []
        for s in stringsizes:
            self.strings.append(string(s))
        self._recalc()

    def _recalc(self, recursive=False):
        if recursive:
            for s in self.strings:
                s._recalc(recursive)
        self.Isc = self.I(0)
        self.Voc = self.V(0)

    def I(self, V):
        sumI = 0
        for s in self.strings:
            sumI += s.I(V)
        return sumI

    @cacheobj
    def V(self, I):
        def Ierr(x):
            ret = abs(I - self.I(x))
            return ret
        if I < 0:
            lower = 0
            upper = panel_Voc*len(self.strings[0].panels)
        elif I > self.Isc:
            lower = -panel_Voc*opts.string_length*10
            upper = 0
        else:
            lower = 0
            upper = panel_Voc*opts.string_length
        print "aV I=%g lower=%g upper=%g" % (I, lower, upper)
        return find_best2(Ierr, lower, upper)

    def P(self, V):
        if V < inverter_range[0] or V > inverter_range[1]:
            return 0
        I = self.I(V)
        return I * V

    def Ppos(self, V):
        p = self.P(V)
        if (p < 0):
            return 0
        return p

    def set_temperature(self, t):
        for s in self.strings:
            s.set_temperature(t)
        self._recalc()

    def set_irradiance(self, w):
        for s in self.strings:
            s.set_irradiance(w)
        self._recalc()

    def MPPV(self):
        def Perr(x):
            p = self.P(x)
            if p <= 0:
                return 1e100
            return 1.0/p
        return find_best2(Perr, 0, inverter_range[1])

    def MPPP(self):
        return self.P(self.MPPV())


class pvsystem:
    '''model a set of PV arrays, each with 1 MPPT'''
    def __init__(self, arraysizes):
        self.arrays = []
        for a in arraysizes:
            self.arrays.append(pvarray(a))
        self._recalc()

    def _recalc(self, recursive=False):
        if recursive:
            for a in self.arrays:
                a._recalc(recursive)

    def set_temperature(self, t):
        for a in self.arrays:
            a.set_temperature(t)
        self._recalc()

    def set_irradiance(self, w):
        for a in self.arrays:
            a.set_irradiance(w)
        self._recalc()

    def MPPP(self):
        ret = 0
        for a in self.arrays:
            ret += a.MPPP()
        return ret


def plotit(fn, low, high, label='', loc=None, plot_steps=None):
    '''convenience plotting function'''
    if plot_steps is None:
        plot_steps = opts.plot_steps
    def array_func(fn, a):
        ret = []
        for i in range(0, len(a)):
            ret.append(fn(a[i]))
        return ret
    xrange = numpy.arange(low, high, (high-low)/plot_steps)
    y = array_func(fn, xrange)
    pylab.ion()
    pylab.figure(num=1, figsize=(12,6))
    pylab.plot(xrange, y, label=label)
    pylab.draw()
    if loc is not None:
        pylab.legend(loc=loc)
        pylab.draw()


def shaded_1mpp(a):
    '''plot shading effects with 1 MPPT'''
    maxv = round(inverter_range[1], 200)
    minv = round_down(inverter_range[0]/2, 100)
    plotit(a.P, minv, maxv, label='Unshaded', loc='upper left')
    MPPV=a.MPPV()
    max_MPPP = a.P(MPPV)
    print("MPP=%.2fW at %.2fV" % (max_MPPP, MPPV))

    pylab.xlabel('Voltage (V)')
    pylab.ylabel('Power (W)')
    pylab.title("%d strings, 1 MPPT" % len(a.strings))

    for i in range(0, opts.string_length):
        a.strings[0].panels[i].set_irradiance(opts.irradiance * opts.shade_factor)
        a._recalc(True)
        MPPV=a.MPPV()
        MPPP=a.P(MPPV)
        print("MPP=%.2fW at %.2fV (%d shaded lost %4.1f%%)" % (MPPP, MPPV, i+1, 100*(max_MPPP-MPPP)/max_MPPP))
        plotit(a.P, minv, maxv, label='%d shaded' % (i+1), loc='upper left')

    pylab.show()

def shaded_1mpp_1string():
    '1MPPT with 1 string'
    a=pvarray([opts.string_length])
    shaded_1mpp(a)

def shaded_1mpp_2string():
    '1MPPT with 2 strings'
    a=pvarray([opts.string_length, opts.string_length])
    shaded_1mpp(a)

def shaded_1mpp_3string():
    '1MPPT with 3 strings'
    a=pvarray([opts.string_length, opts.string_length, opts.string_length])
    shaded_1mpp(a)

def shaded_1mpp_4string():
    '1MPPT with 4 strings'
    a=pvarray([opts.string_length, opts.string_length, opts.string_length, opts.string_length])
    shaded_1mpp(a)

def shaded_1mpp_half():
    '''1MPPT, 2 strings, 1 string half shaded'''
    a=pvarray([opts.string_length, opts.string_length])
    maxv = round(inverter_range[1], 200)
    minv = round_down(inverter_range[0], 100)
    plotit(a.P, minv, maxv, label='Unshaded', loc='upper left')
    MPPV=a.MPPV()
    max_MPPP = a.P(MPPV)
    print("MPP=%.2fW at %.2fV" % (max_MPPP, MPPV))

    pylab.xlabel('Voltage (V)')
    pylab.ylabel('Power (W)')
    pylab.title("%d strings, 1 MPPT" % len(a.strings))

    for i in range(0, opts.string_length/2):
        a.strings[0].panels[i].set_irradiance(opts.irradiance * opts.shade_factor)
    a._recalc(True)
        
    MPPV=a.MPPV()
    MPPP=a.P(MPPV)
    print("MPP=%.2fW at %.2fV (%d shaded lost %4.1f%%)" % (MPPP, MPPV, i+1, 100*(max_MPPP-MPPP)/max_MPPP))
    plotit(a.P, minv, maxv, label='%d shaded' % (i+1), loc='upper left')

    pylab.show()

def shaded_2mpp():
    '''plot shading effects with 2 MPPTs'''
    a1=pvarray([opts.string_length])
    a2=pvarray([opts.string_length])
    Pa1 = a1.MPPP()
    Pa2 = a2.MPPP()
    print("A1=%.2f A2=%.2f MPPP=%.2fW" % (Pa1, Pa2, Pa1+Pa2))

    for i in range(0, opts.string_length):
        a1.strings[0].panels[i].set_irradiance(opts.irradiance * opts.shade_factor)
        a1._recalc(True)
        Pa1 = a1.MPPP()
        Pa2 = a2.MPPP()
        print("A1=%.2f A2=%.2f MPPP=%.2fW (%d shaded)" % (Pa1, Pa2, Pa1+Pa2, i+1))


def panel_shading():
    '''plot shading effects on 1 panel'''
    p = panel()
    maxv = p.Voc
    minv = 0
    plotit(p.P, minv, maxv, label='Unshaded', loc='upper left')

    pylab.xlabel('Voltage (V)')
    pylab.ylabel('Power (W)')
    pylab.title("1 Panel")

    for i in range(0, opts.num_subpanels):
        p.subpanels[i].set_irradiance(opts.irradiance * opts.shade_factor)
        p._recalc(True)
        plotit(p.P, minv, maxv, label='%d shaded' % (i+1), loc='upper left')
    pylab.show()

def panel_shading_IV():
    '''plot shading effects on 1 panel'''
    p = panel()
    maxv = p.Voc
    minv = 0
    plotit(p.I, minv, maxv, label='Unshaded', loc='upper left')

    pylab.xlabel('Voltage (V)')
    pylab.ylabel('Current (A)')
    pylab.title("1 Panel")

    for i in range(0, opts.num_subpanels):
        p.subpanels[i].set_irradiance(opts.irradiance * opts.shade_factor)
        p._recalc(True)
        plotit(p.I, minv, maxv, label='%d shaded' % (i+1), loc='upper left')
    pylab.show()

def panel_shading_VI():
    '''plot shading effects on 1 panel'''
    p = panel()
    maxi = p.Isc
    mini = 0
    plotit(p.V, mini, maxi, label='Unshaded', loc='upper left')

    pylab.xlabel('Current (A)')
    pylab.ylabel('Voltage (V)')
    pylab.title("1 Panel")

    for i in range(0, opts.num_subpanels):
        p.subpanels[i].set_irradiance(opts.irradiance * opts.shade_factor)
        p._recalc(True)
        plotit(p.V, mini, maxi, label='%d shaded' % (i+1), loc='upper left')
    pylab.show()

def mpp_loss():
    '''show loss from not having separate MPPTs per string'''
    n=opts.string_length
    s1=pvsystem([[n, n]])   # 2 strings on one MPPT
    s2=pvsystem([[n],[n]])  # 1 string on two MPPTs
    mppt1 = s1.MPPP()
    mppt2 = s2.MPPP()
    loss = 100.0*(mppt2-mppt1)/mppt2
    print("2mppt=%.2fW 1mppt=%.2f loss=%4.0f%%" % (mppt2, mppt1, loss))

    for i in range(0, n):
        s1.arrays[0].strings[0].panels[i].set_irradiance(opts.irradiance * opts.shade_factor)
        s2.arrays[0].strings[0].panels[i].set_irradiance(opts.irradiance * opts.shade_factor)
        s1._recalc(True)
        s2._recalc(True)
        mppt1 = s1.MPPP()
        mppt2 = s2.MPPP()
        loss = 100.0*(mppt2-mppt1)/mppt2
        print("2mppt=%.2fW 1mppt=%.2f loss=%4.1f%% (%d shaded)" % (mppt2, mppt1, loss, i+1))


def mppt_sfr_old():
    '''model stephens choice of 1 or 2 MPPTs'''
    n = opts.string_length
    s1=pvsystem([[n, n, n]])   # he currently has 3 strings on 1 MPPT
    s2=pvsystem([[n], [n, n]]) # he is considering moving one string to another MPPT
    mppt1 = s1.MPPP()
    mppt2 = s2.MPPP()
    loss = 100.0*(mppt2-mppt1)/mppt2
    print("2mppt=%.2fW 1mppt=%.2f loss=%4.0f%%" % (mppt2, mppt1, loss))

    for i in range(0, n):
        s1.arrays[0].strings[0].panels[i].set_irradiance(opts.irradiance * opts.shade_factor)
        s2.arrays[0].strings[0].panels[i].set_irradiance(opts.irradiance * opts.shade_factor)
        s1._recalc(True)
        s2._recalc(True)
        mppt1 = s1.MPPP()
        mppt2 = s2.MPPP()
        loss = 100.0*(mppt2-mppt1)/mppt2
        print("2mppt=%.2fW 1mppt=%.2f loss=%4.1f%% (%d shaded)" % (mppt2, mppt1, loss, i+1))

def mppt_sfr():
    '''model stephens choice of 1 or 2 MPPTs'''
    n = opts.string_length

    # on one system, shade subpanels along one string at a time
    def vertical_stringing(nshaded):
        s.set_irradiance(opts.irradiance)
        anum = snum = pnum = bnum = 0
        for i in range(0, nshaded):
            anum = (i / n_array)
            snum = (i / n_string) % len(s.arrays[0].strings)
            pnum = (i / n_panel) % len(s.arrays[0].strings[0].panels)
            bnum = i % len(s.arrays[0].strings[0].panels[0].subpanels)
            s.arrays[anum].strings[snum].panels[pnum].subpanels[bnum].set_irradiance(opts.irradiance * opts.shade_factor)
        print nshaded, anum, snum, pnum, bnum
        s._recalc(True)
        return s.MPPP()

    # 1 MPPTs, 3 strings
    s=pvsystem([[n,n,n]])
    n_panel  = len(s.arrays[0].strings[0].panels[0].subpanels)
    n_string = n_panel  * len(s.arrays[0].strings[0].panels)
    n_array  = n_string * len(s.arrays[0].strings)
    n_total  = n_array * len(s.arrays)

    pylab.xlabel('Number of subpanels shaded')
    pylab.ylabel('Total power (W)')
    pylab.title("MPPT choice and shading")
    pylab.draw()
    plotit(vertical_stringing, 0, n*opts.num_subpanels, label='1 MPPT',
           plot_steps=n*opts.num_subpanels, loc='upper right')    

    # 2 MPPTs, 1 strings on each
    s=pvsystem([[n], [n,n]])
    n_panel  = len(s.arrays[0].strings[0].panels[0].subpanels)
    n_string = n_panel  * len(s.arrays[0].strings[0].panels)
    n_array  = n_string * len(s.arrays[0].strings)
    n_total  = n_array * len(s.arrays)

    plotit(vertical_stringing, 0, n*opts.num_subpanels, label='2 MPPTs',
           plot_steps=n*opts.num_subpanels, loc='upper right')    

    pylab.show()


def pv_82():
    '''model question from itshot'''
    n = opts.string_length

    s2 = pvarray([2])
    s8 = pvarray([8])
    s82 = pvarray([8,2])

    pylab.xlabel('Voltage')
    pylab.ylabel('Power (W)')
    pylab.title("PV for 8-2 combination")
    pylab.draw()
    plotit(s2.Ppos, 0, 400, label='2 Panels', loc='upper left')    
    plotit(s8.Ppos, 0, 400, label='8 Panels', loc='upper left')    
    plotit(s82.Ppos, 0, 400, label='8+2 parallel', loc='upper left')    
    pylab.show()


def stringing_tridge():
    '''model my choice of stringing approaches'''
    n = opts.string_length

    # on one system, shade subpanels along one string at a time
    def vertical_stringing(nshaded):
        s.set_irradiance(opts.irradiance)
        anum = snum = pnum = bnum = 0
        for i in range(0, nshaded):
            anum = (i / n_array)
            snum = (i / n_string) % len(s.arrays[0].strings)
            pnum = (i / n_panel) % len(s.arrays[0].strings[0].panels)
            bnum = i % len(s.arrays[0].strings[0].panels[0].subpanels)
            s.arrays[anum].strings[snum].panels[pnum].subpanels[bnum].set_irradiance(opts.irradiance * opts.shade_factor)
        print nshaded, anum, snum, pnum, bnum
        s._recalc(True)
        return s.MPPP()

    # on the other system, shade subpanels spread over the strings
    def horiz_stringing(nshaded):
        s.set_irradiance(opts.irradiance)
        anum = snum = pnum = bnum = 0
        for i in range(0, nshaded):
            anum = i % len(s.arrays)
            snum = (i / len(s.arrays)) % len(s.arrays[0].strings)
            pnum = (i / (len(s.arrays)*len(s.arrays[0].strings))) % len(s.arrays[0].strings[0].panels)
            bnum = (i / (len(s.arrays)*len(s.arrays[0].strings)*len(s.arrays[0].strings[0].panels))) % len(s.arrays[0].strings[0].panels[0].subpanels)
            s.arrays[anum].strings[snum].panels[pnum].subpanels[bnum].set_irradiance(opts.irradiance * opts.shade_factor)
        print nshaded, anum, snum, pnum, bnum
        s._recalc(True)
        return s.MPPP()

    print "This may take some time ...."

    # 12 MPPTs, 1 strings on each
    s=pvsystem([[n],[n],[n],[n],[n],[n],[n],[n],[n],[n],[n],[n]]) 
    n_panel  = len(s.arrays[0].strings[0].panels[0].subpanels)
    n_string = n_panel  * len(s.arrays[0].strings[0].panels)
    n_array  = n_string * len(s.arrays[0].strings)
    n_total  = n_array * len(s.arrays)

    plotit(horiz_stringing, 0, n_total, label='Horizontal, 12 inverters',
           plot_steps=n_total, loc='upper right')
    pylab.xlabel('Number of subpanels shaded')
    pylab.ylabel('Total power (W)')
    pylab.title("Stringing effect on shading")
    pylab.draw()
    plotit(vertical_stringing, 0, n_total, label='Vertical, 12 inverters',
           plot_steps=n_total, loc='upper right')    

    # 6 MPPTs, 2 strings on each
    s=pvsystem([[n,n], [n,n], [n,n], [n,n], [n,n], [n,n]]) 
    n_panel  = len(s.arrays[0].strings[0].panels[0].subpanels)
    n_string = n_panel  * len(s.arrays[0].strings[0].panels)
    n_array  = n_string * len(s.arrays[0].strings)
    n_total  = n_array * len(s.arrays)

    plotit(horiz_stringing, 0, n_total, label='Horizontal, 6 inverters',
           plot_steps=n_total, loc='upper right')
    pylab.xlabel('Number of subpanels shaded')
    pylab.ylabel('Total power (W)')
    pylab.title("Stringing effect on shading")
    pylab.draw()
    plotit(vertical_stringing, 0, n_total, label='Vertical, 6 inverters',
           plot_steps=n_total, loc='upper right')    

    pylab.show()


def string_IV():
    '''string IV graph'''
    s=string(opts.string_length)
    plotit(s.I, 0, round(s.V(0),2), label='IV for 1 string')
    pylab.xlabel('Voltage (V)')
    pylab.ylabel('Current (A)')
    pylab.title("String IV")
    pylab.show()

def string_VI():
    '''string VI graph'''
    s=string(opts.string_length)
    plotit(s.V, 0, round(s.I(0),2), label='VI for 1 string')
    pylab.ylabel('Voltage (V)')
    pylab.xlabel('Current (A)')
    pylab.title("String VI")
    pylab.show()

def string_P():
    '''power vs voltage for one string'''
    a=pvarray([opts.string_length])
    maxv = round(inverter_range[1], 100)
    minv = round_down(inverter_range[0], 100)
    plotit(a.P, minv, maxv, label='Unshaded', loc='upper left')
    MPPV=a.MPPV()
    print("MPP=%.2fW at %.2fV" % (a.P(MPPV), MPPV))

    pylab.xlabel('Voltage (V)')
    pylab.ylabel('Power (W)')
    pylab.title("Single string shading")

    for i in range(0, opts.string_length):
        a.strings[0].panels[i].set_irradiance(opts.irradiance * opts.shade_factor)
        a._recalc(True)
        MPPV=a.MPPV()
        print("MPP=%.2fW at %.2fV (%d shaded)" % (a.P(MPPV), MPPV, i+1))
        plotit(a.P, minv, maxv, label='%d shaded' % (i+1), loc='upper left')

    pylab.show()

def array_IV():
    '''array IV graph'''
    a=pvarray([opts.string_length, opts.string_length])
    plotit(a.I, 0, round(a.V(0),2), label='IV for 2 string array')
    pylab.xlabel('Voltage (V)')
    pylab.ylabel('Current (A)')
    pylab.title("Array IV (2 strings)")
    pylab.show()

def array_VI():
    '''array VI graph'''
    a=pvarray([opts.string_length, opts.string_length])
    plotit(a.V, 0, round(a.I(0),2), label='IV for 2 string array')
    pylab.ylabel('Voltage (V)')
    pylab.xlabel('Current (A)')
    pylab.title("Array VI (2 strings)")
    pylab.show()
    

def cell_IV():
    '''one cell IV graph'''
    c=cell()
    plotit(c.I, -1, round(c.V(0), 2), label='IV for 1 cell')
    pylab.xlabel('Voltage (V)')
    pylab.ylabel('Current (A)')
    pylab.title("Single cell IV")
    pylab.show()

def cell_VI():
    '''one cell VI graph'''
    c=cell()
    plotit(c.V, -0.1*c.Isc, 1.1*c.I(0), label='VI for 1 cell')
    pylab.ylabel('Voltage (V)')
    pylab.xlabel('Current (A)')
    pylab.title("Single cell VI")
    pylab.show()

def cell_P():
    '''one cell P graph'''
    c=cell()
    plotit(c.P, 0, round(c.V(0), 2))
    pylab.xlabel('Voltage (V)')
    pylab.ylabel('Power (W)')
    pylab.title("Cell power")
    pylab.show()

def panel_IV():
    '''one panel IV graph'''
    c=panel()
    plotit(c.I, 0, round(c.V(0),2), label='IV for 1 panel')
    pylab.xlabel('Voltage (V)')
    pylab.ylabel('Current (A)')
    pylab.title("Single panel IV")
    pylab.show()

def panel_VI():
    '''one panel VI graph'''
    c=panel()
    plotit(c.V, 0, round(c.I(0), 1), label='VI for 1 panel')
    pylab.ylabel('Voltage (V)')
    pylab.xlabel('Current (A)')
    pylab.title("Single panel VI")
    pylab.show()

def panel_P():
    '''one panel P graph'''
    p=panel()
    plotit(p.P, 0, round(p.V(0), 20))
    pylab.xlabel('Voltage (V)')
    pylab.ylabel('Power (W)')
    pylab.title("Panel power")
    pylab.show()


def subpanel_IV():
    '''one subpanel IV graph'''
    c=subpanel()
    plotit(c.I, 0, round(c.V(0), 20), label='IV for 1 subpanel')
    pylab.xlabel('Voltage (V)')
    pylab.ylabel('Current (A)')
    pylab.title("Single subpanel IV")
    pylab.show()

def subpanel_VI():
    '''one subpanel VI graph'''
    c=subpanel()
    plotit(c.V, 0, round(c.I(0), 1), label='VI for 1 subpanel')
    pylab.ylabel('Voltage (V)')
    pylab.xlabel('Current (A)')
    pylab.title("Single subpanel VI")
    pylab.show()

def subpanel_P():
    '''one subpanel P graph'''
    s=subpanel()
    plotit(s.P, 0, round(s.V(0), 20))
    pylab.xlabel('Voltage (V)')
    pylab.ylabel('Power (W)')
    pylab.title("Subpanel power")
    pylab.show()

def cell_T():
    '''cell temperature sensitivity'''
    c=cell()
    maxv = round(c.V(0), 2)
    for t in range(-5, 75, 10):
        c.set_temperature(t)
        plotit(c.P, 0, maxv, label='%d C' % t, loc='upper left')
        pylab.xlabel('Voltage (V)')
        pylab.ylabel('Power (W)')
        pylab.title("Cell temperature effect")
    pylab.show()

def cell_L():
    '''cell power light sensitivity'''
    c=cell()
    for w in range(1500, 0, -250):
        c.set_irradiance(w)
        plotit(c.P, 0, round(c.V(0), 5), label='%d W/m2' % w, loc='upper left')
    pylab.xlabel('Voltage (V)')
    pylab.ylabel('Power (W)')
    pylab.title("Cell light effect")
    pylab.show()

def panel_T():
    '''panel power temperature sensitivity'''
    p=panel()
    maxv = round(p.V(0), 20)
    for t in range(-5, 75, 10):
        p.set_temperature(t)
        plotit(p.P, 0, maxv, label='%d C' % t, loc='upper left')
        pylab.xlabel('Voltage (V)')
        pylab.ylabel('Power (W)')
        pylab.title("Panel power temperature effect")
    pylab.show()

def panel_IV_T():
    '''panel IV temperature sensitivity'''
    p=panel()
    maxv = round(p.V(0), 20)
    for t in range(-5, 65, 10):
        p.set_temperature(t)
        plotit(p.I, 0, maxv, label='%d C' % t, loc='upper right')
        pylab.xlabel('Voltage (V)')
        pylab.ylabel('Current (A)')
        pylab.title("Panel IV temperature effect")
    pylab.show()

def panel_L():
    '''panel power light sensitivity'''
    p=panel()
    for w in range(1500, 0, -250):
        p.set_irradiance(w)
        plotit(p.P, 0, round(p.V(0), 20), label='%d W/m2' % w, loc='upper left')
    pylab.xlabel('Voltage (V)')
    pylab.ylabel('Power (W)')
    pylab.title("Panel light effect")
    pylab.show()

def panel_IV_L():
    '''panel VI light sensitivity'''
    p=panel()
    for w in range(1500, 0, -250):
        p.set_irradiance(w)
        plotit(p.I, 0, round(p.V(0), 25), label='%d W/m2' % w, loc='upper right')
    pylab.xlabel('Voltage (V)')
    pylab.ylabel('Current (A)')
    pylab.title("Panel VI light effect")
    pylab.show()


def selftest():
    '''check model internal behaviour'''
    def check(tag, v1, v2, tolerance=1.0):
        err = abs((v1-v2)/float(v1))
        tolerance /= 100.0
        if err > tolerance:
            raise Exception("%s: v1=%f v2=%g" % (tag, v1, v2))

    def check_IV(tag, m, tolerance=10.0):
        irange = numpy.arange(-m.Isc, m.Isc*2, (3*m.Isc/100))
        for i in irange:
            check('%s: I=%g Isc=%g' % (tag, i, m.Isc), m.I(m.V(i)), i,
                  tolerance=tolerance)
        vrange = numpy.arange(-m.Voc, m.Voc*2, (3*m.Voc/100))
        for v in vrange:
            check('%s: V=%g I=%g Voc=%g' % (tag, v, m.I(v), m.Voc), m.V(m.I(v)), v,
                  tolerance=tolerance)

    print "Checking cell behaviour"
    c=cell()
    check("Voc", cell_Voc, c.V(0))
    check("Isc", panel_Isc, c.I(0))
    check("Vmp", cell_Vmp, c.V(panel_Imp))
    check("Imp", panel_Imp, c.I(cell_Vmp))
    check("P",   panel_Imp*cell_Vmp, c.P(cell_Vmp))
    check("Imp2", panel_Imp, c.I(c.V(panel_Imp)))
    check("Vmp2", cell_Vmp, c.V(c.I(cell_Vmp)))

    check_IV("cell IV", c)
    c.set_irradiance(500)
    check_IV("cell IV 500W/m2", c)

    c.set_irradiance(500)
    p1 = c.P(cell_Vmp)
    c.set_irradiance(1000)
    p2 = c.P(cell_Vmp)
    check("irradiance", p1, p2/2, tolerance=10)

    check("leakage", c.I(c.Voc+1), -1.0/(opts.panel_Roc/cell_scale), tolerance=10)

    
    

    print "Checking subpanel behaviour"
    s=subpanel()
    check("Voc", cell_Voc*opts.cells_per_subpanel, s.V(0))
    check("Isc", panel_Isc, s.I(0))
    check("Vmp", cell_Vmp*opts.cells_per_subpanel, s.V(panel_Imp))
    check("Imp", panel_Imp, s.I(panel_Vmp/opts.num_subpanels))

    check_IV("subpanel IV", s)

    s2=subpanel()
    s2.set_irradiance(500)
    check_IV("subpanel IV 500W/m2", s2, tolerance=2.0)
    Vmp = panel_Vmp/opts.num_subpanels
    check("irradiance", s2.P(Vmp), s.P(Vmp)/2, tolerance=10)

    check("bypass1", s.I(-1.0) - s.Isc, 1.0/subpanel_Rbp)
    check("bypass2", s.I(-2.0) - s.Isc, 2.0/subpanel_Rbp)
    check("bypass3", s.V(s.Isc + 1.0/subpanel_Rbp), -1)

    print "Checking panel behaviour"
    p=panel()
    Vmp = panel_Vmp
    check("Voc", panel_Voc, p.V(0))
    check("Isc", panel_Isc, p.I(0))
    check("Vmp", panel_Vmp, p.V(panel_Imp))
    check("Imp", panel_Imp, p.I(Vmp))
    check("bypass", p.I(-1.0) - p.Isc, 1.0/(subpanel_Rbp*opts.num_subpanels))

    check_IV("panel IV", p)

    p2=subpanel()
    p2.set_irradiance(500)
    check_IV("panel IV 500W/m2", p2, tolerance=2.0)
    Vmp = panel_Vmp/opts.num_subpanels
    check("irradiance", p2.P(Vmp), p.P(Vmp)/2, tolerance=15)

    print "Checking string behaviour"
    s=string(opts.string_length)
    check("Voc", panel_Voc*opts.string_length, s.V(0))
    check("Isc", panel_Isc, s.I(0))
    check("Vmp", panel_Vmp*opts.string_length, s.V(panel_Imp))
    check("Imp", panel_Imp, s.I(panel_Vmp*opts.string_length))
    check("bypass", s.I(-1.0*opts.string_length) - s.Isc, 1.0/(subpanel_Rbp*opts.num_subpanels))

    check_IV("string IV", s, tolerance=5.0)

    s2=string(opts.string_length)
    s2.set_irradiance(500)

    check_IV("string IV 500W/m2", s2, tolerance=5.0)

    check("irradiance", s2.P(panel_Vmp*opts.string_length),
                             s.P(panel_Vmp*opts.string_length)/2, tolerance=10)

    print "Checking array behaviour"
    a=pvarray([opts.string_length, opts.string_length])
    check("Voc", panel_Voc*opts.string_length, a.V(0))
    check("Isc", panel_Isc*2, a.I(0))
#    check("Vmp", panel_Vmp*opts.string_length, a.V(panel_Imp*2))
    check("Imp", panel_Imp*2, a.I(panel_Vmp*opts.string_length))
    check("bypass", a.I(-1.0*opts.string_length) - a.Isc, 2.0/(opts.num_subpanels*subpanel_Rbp))

    check_IV("array IV", a, tolerance=2.0)

    a2=pvarray([opts.string_length, opts.string_length])
    a2.set_irradiance(500)

    check_IV("array IV 500W/m2", a2, tolerance=1.0)

    print "All OK"


if __name__ == "__main__":
    parser = OptionParser("pv_model [options] <function>")

    functions = {
        "cell_IV"      : cell_IV,
        "cell_VI"      : cell_VI,
        "cell_P"       : cell_P,
        "cell_T"       : cell_T,
        "cell_L"       : cell_L,
        "subpanel_IV"  : subpanel_IV,
        "subpanel_VI"  : subpanel_VI,
        "subpanel_P"   : subpanel_P,
        "panel_IV"     : panel_IV,
        "panel_VI"     : panel_VI,
        "panel_P"      : panel_P,
        "panel_T"      : panel_T,
        "panel_IV_T"   : panel_IV_T,
        "panel_L"      : panel_L,
        "panel_IV_L"   : panel_IV_L,
        "panel-shading": panel_shading,
        "panel-shadingIV": panel_shading_IV,
        "panel-shadingVI": panel_shading_VI,
        "string_P"     : string_P,
        "string_IV"    : string_IV,
        "string_VI"    : string_VI,
        "array_IV"     : array_IV,
        "array_VI"     : array_VI,
        "2MPP-loss"    : mpp_loss,
        "2MPP-shading" : shaded_2mpp,
        "1MPP-1string" : shaded_1mpp_1string,
        "1MPP-2string" : shaded_1mpp_2string,
        "1MPP-3string" : shaded_1mpp_3string,
        "1MPP-4string" : shaded_1mpp_4string,
        "1MPP-halfshaded" : shaded_1mpp_half,
        "MPPT_sfr"     : mppt_sfr,
        "stringing_tridge"     : stringing_tridge,
        "pv_82"        : pv_82,
        "selftest"     : selftest
        }

    parser.add_option("", "--panel", dest="panel_type", help="panel type")
    parser.add_option("", "--inverter", dest="inverter_type", help="inverter type")
    parser.add_option("", "--string-length", dest="string_length", help="modules per string",
                      type='int', default=11)
    parser.add_option("", "--panel-Isc", dest="panel_Isc", help="panel short circuit current", type='float')
    parser.add_option("", "--panel-Voc", dest="panel_Voc", help="panel open circuit voltage", type='float')
    parser.add_option("", "--panel-Imp", dest="panel_Imp", help="max power panel current", type='float')
    parser.add_option("", "--panel-Vmp", dest="panel_Vmp", help="max power panel voltage", type='float')
    parser.add_option("", "--panel-Roc", dest="panel_Roc", help="panel open circuit resistance (Ohm)",
                      type='float', default=50)
    parser.add_option("", "--panel-Rsc", dest="panel_Rsc", help="panel short circuit resistance (Ohm)",
                      type='float', default=100)
    parser.add_option("", "--subpanel-Rbp", dest="subpanel_Rbp", help="subpanel bypass diode resistance (Ohm)",
                      type='float', default=0.1)
    parser.add_option("", "--inverter-MPP-min", dest="inverter_MPP_min", help="inverter MPP minimum", type='int')
    parser.add_option("", "--inverter-MPP-max", dest="inverter_MPP_max", help="inverter MPP maximum", type='int')
    parser.add_option("", "--temperature", dest="temperature", help="panel temperature (C)",type='float', default=25)
    parser.add_option("", "--irradiance", dest="irradiance", help="light level (W/m2)",type='float', default=1000)
    parser.add_option("", "--shade-factor", dest="shade_factor", help="shading factor",type='float', default=0.05)
    parser.add_option("", "--subpanels", dest="num_subpanels", help="subpanels per panel",
                      type='int', default=3)
    parser.add_option("", "--cells-per-subpanel", dest="cells_per_subpanel", help="cells per subpanel",
                      type='int', default=1)
    parser.add_option("", "--blocking-diodes", dest="blocking_diodes",
                      action="store_true", help="enable blocking diodes", default=False)
    parser.add_option("", "--plot-steps", dest="plot_steps", help="steps on plots", type='int', default=100)

    (opts, args) = parser.parse_args()
    if len(args) != 1:
        print("Available functions are:")
        k = functions.keys()
        k.sort()
        for f in k:
            print("\t%-14s  %s" % (f, functions[f].__doc__))
        parser.error("You must supply a function")
    f = args[0]
    if not f in functions:
        parser.error("Unrecognised function '%s'" % f)

    inverter_range = [x for x in inverter_range]

    if opts.panel_type:
        (panel_Isc, panel_Voc, panel_Imp, panel_Vmp, panel_VT, panel_IT) = panels[opts.panel_type]
    if opts.inverter_type:
        inverter_range = inverters[opts.inverter_type]
        inverter_range = [x for x in inverter_range]
    if opts.panel_Isc:
        panel_Isc = float(opts.panel_Isc)
    if opts.panel_Voc:
        panel_Voc = float(opts.panel_Voc)
    if opts.panel_Imp:
        panel_Imp = float(opts.panel_Imp)
    if opts.panel_Vmp:
        panel_Vmp = float(opts.panel_Vmp)
    if opts.inverter_MPP_min:
        inverter_range[0] = float(opts.inverter_MPP_min)
    if opts.inverter_MPP_max:
        inverter_range[1] = float(opts.inverter_MPP_max)

    cell_scale = opts.num_subpanels * opts.cells_per_subpanel
    cell_Voc = panel_Voc / cell_scale
    cell_Vmp = panel_Vmp / cell_scale
    cell_VT  = panel_VT  / cell_scale
    subpanel_Rbp = (opts.subpanel_Rbp*3) / opts.num_subpanels

    try:
        functions[f]()
    except KeyboardInterrupt:
        sys.exit(1)