python, plot planck curves looping through arrays
Question: I try to get myself familiar with programming in python but have just started
and struggling with the following problem. Maybe someone can give me a hint
how to proceed or where I can look for a nice solution.

I'd like to plot planck curves for 132 wavelength in 6 different temperatures
via a loop in a loop. The function planckwavel receives two parameters,
wavelength and temperature, which I separated in two loops.

I so far managed to use lists, which worked, however probably not solved in an
elegant way:

    
    
    plancks = []
    temp = [280, 300, 320, 340, 360, 380]
    temp_len = len(temp)
    
    ### via fun planckwavel
    for i in range(temp_len):
        t_list = [] # list nach jeder j schleife wieder leeren
        for j in range(wl_centers_ar.shape[0]):
            t = planckwavel(wl_centers_ar[j],temp[i])
            t_list.append(t)      
        plancks.append(t_list)
    
    
    ### PLOT Planck curves
    
    plancks = np.array(plancks).T # convert list to array and transpose
    
    view_7 = plt.figure(figsize=(8.5, 4.5))
    plt.plot(wl_centers_ar,plancks)
    plt.xticks(rotation='vertical')
    

But I would like to use arrays insted of lists, as I like to continue
afterwards with huge more dimensional images. So I tried the same with arrays
but unfortunately failed with this code:

    
    
    plancks_ar = zeros([132,6], dtype=float ) # create array and fill with zeros
    temp_ar = array([273, 300, 310, 320, 350, 373])
    
    for i in range(temp_ar.shape[0]):
        t_ar = np.zeros(plancks_ar.shape[0])
        for j in range(plancks_ar.shape[0]):
            t = planck(wl_centers_ar[j]*1e-6,temp[1])/10**6
            np.append(t_ar,t)        
        np.append(plancks_ar, t_ar)
    
    plt.plot(wl_centers_ar,plancks)
    

I would be very thankful, if someone can give me some advice.

Thanx, best regards,

peter


Answer: I think you're asking about how to use NumPy's [broadcasting and
vectorization](http://docs.scipy.org/doc/numpy/user/basics.broadcasting.html).
Here's a way to remove the explicit Python loops:

    
    
    import numpy as np
    
    # Some physical constants we'll need
    h, kB, c = 6.626e-34, 1.381e-23, 2.998e8
    
    def planck(lam, T):
        # The Planck function, using NumPy vectorization
        return 2*h*c**2/lam**5 / (np.exp(h*c/lam/kB/T) - 1)
    
    # wavelength array, 3 - 300 um
    lam = np.linspace(3, 75, 132)
    # temperature array
    T = np.array([280, 300, 320, 340, 360, 380])
    # Remember to convert wavelength from um to m
    pfuncs = planck(lam * 1.e-6, T[:,None])
    
    import pylab
    for pfunc in pfuncs:
        pylab.plot(lam, pfunc)
    pylab.show()
    

[![enter image description
here](http://i.stack.imgur.com/tfWPV.png)](http://i.stack.imgur.com/tfWPV.png)

We want to calculate `planck` for each wavelength and for each T, so we need
to broadcast the calculation over the two arrays. Following the rules laid out
in the documentation linked to above, we can do that by adding a new axis to
the temperature array (with `T[:, None]`):

    
    
    lam:       132
    T      6 x   1
    --------------
           6 x 132
    

The final dimension of `T[:, None]` is 1, so the 132 values of `lam` can be
broadcast across it to produce a `6 x 132` array: 6 rows (one for each T) of
132 values (the wavelengths).

