Plot the geometry of the column model, optionally with somatic locations
    and optionally with reconstructed neurons
    
    kwargs:
    ::
        ax : matplotlib.axes.AxesSubplot
        aspect : str
            matplotlib.axis argument
        isometricangle : float
            pseudo-3d view angle
        plot_somas : bool
            plot soma locations
        plot_morphos : bool
            plot full morphologies
        num_unitsE : int
            number of excitatory morphos plotted per population
        num_unitsI : int
            number of inhibitory morphos plotted per population
        clip_dendrites : bool
            draw dendrites outside of axis
        mainpops : bool
            if True, plot only main pops, e.g. b23 and nb23 as L23I
    
    return:
    ::
        axis : list
            the plt.axis() corresponding to input aspect

def plot_population(ax,
                    populationParams,
                    electrodeParams,
                    layerBoundaries,
                    aspect='equal',
                    isometricangle=np.pi/12,
                    X=['EX', 'IN'],
                    markers=['^', 'o'],
                    colors=['r', 'b'],
                    layers = ['upper', 'lower'],
                    title='positions'):
    '''
    Plot the geometry of the column model, optionally with somatic locations
    and optionally with reconstructed neurons
    
    kwargs:
    ::
        ax : matplotlib.axes.AxesSubplot
        aspect : str
            matplotlib.axis argument
        isometricangle : float
            pseudo-3d view angle
        plot_somas : bool
            plot soma locations
        plot_morphos : bool
            plot full morphologies
        num_unitsE : int
            number of excitatory morphos plotted per population
        num_unitsI : int
            number of inhibitory morphos plotted per population
        clip_dendrites : bool
            draw dendrites outside of axis
        mainpops : bool
            if True, plot only main pops, e.g. b23 and nb23 as L23I
    
    return:
    ::
        axis : list
            the plt.axis() corresponding to input aspect
    '''

    remove_axis_junk(ax, ['right', 'bottom', 'left', 'top'])

    
    # DRAW OUTLINE OF POPULATIONS 
    ax.xaxis.set_ticks([])
    ax.yaxis.set_ticks([])

    #contact points
    ax.plot(electrodeParams['x'],
            electrodeParams['z'],
            '.', marker='o', markersize=5, color='k', zorder=0)

    #outline of electrode       
    x_0 = np.array(electrodeParams['r_z'])[1, 1:-1]
    z_0 = np.array(electrodeParams['r_z'])[0, 1:-1]
    x = np.r_[x_0[-1], x_0[::-1], -x_0[1:], -x_0[-1]]
    z = np.r_[100, z_0[::-1], z_0[1:], 100]
    ax.fill(x, z, color=(0.5, 0.5, 0.5), lw=None, zorder=-0.1)

    #outline of populations:
    #fetch the population radius from some population
    r = populationParams[populationParams.keys()[0]]['radius']

    theta0 = np.linspace(0, np.pi, 20)
    theta1 = np.linspace(np.pi, 2*np.pi, 20)
    
    zpos = np.r_[np.array(layerBoundaries)[:, 0],
                 np.array(layerBoundaries)[-1, 1]]
    
    for i, z in enumerate(np.mean(layerBoundaries, axis=1)):
        ax.text(r, z, ' %s' % layers[i],
                va='center', ha='left', rotation='vertical')

    for i, zval in enumerate(zpos):
        if i == 0:
            ax.plot(r*np.cos(theta0),
                    r*np.sin(theta0)*np.sin(isometricangle)+zval,
                    color='k', zorder=-r, clip_on=False)
            ax.plot(r*np.cos(theta1),
                    r*np.sin(theta1)*np.sin(isometricangle)+zval,
                    color='k', zorder=r, clip_on=False)
        else:
            ax.plot(r*np.cos(theta0),
                    r*np.sin(theta0)*np.sin(isometricangle)+zval,
                    color='gray', zorder=-r, clip_on=False)
            ax.plot(r*np.cos(theta1),
                    r*np.sin(theta1)*np.sin(isometricangle)+zval,
                    color='k', zorder=r, clip_on=False)
    
    ax.plot([-r, -r], [zpos[0], zpos[-1]], 'k', zorder=0, clip_on=False)
    ax.plot([r, r], [zpos[0], zpos[-1]], 'k', zorder=0, clip_on=False)
    
    #plot a horizontal radius scalebar
    ax.plot([0, r], [z_0.min()]*2, 'k', lw=2, zorder=0, clip_on=False)
    ax.text(r / 2., z_0.min()-100, 'r = %i $\mu$m' % int(r), ha='center')
    
    #plot a vertical depth scalebar
    ax.plot([-r]*2, [z_0.min()+50, z_0.min()-50],
        'k', lw=2, zorder=0, clip_on=False)
    ax.text(-r, z_0.min(), r'100 $\mu$m', va='center', ha='right')
    
    ax.set_yticks([])
    ax.set_yticklabels([])
    
    #fake ticks:
    for pos in zpos:
        ax.text(-r, pos, 'z=%i-' % int(pos), ha='right', va='center')
 
    ax.set_title(title)
   
    axis = ax.axis(ax.axis(aspect))