def lam2rgb(wav, gamma=0.8, output='norm'):
    """
    Convert a wavelength [nm] to an RGB colour tuple in the visible spectrum.

    This converts a given wavelength of light to an
    approximate RGB colour value with edge attenuation.
    The wavelength must be given in nanometres in the
    range from 380 nm - 750 nm (789 THz - 400 THz).

    Adapted from: http://www.noah.org/wiki/Wavelength_to_RGB_in_Python
    Based on code by Dan Bruton
    http://www.physics.sfasu.edu/astro/color/spectra.html

    Parameters
    ----------
    wav : int or float
        Wavelength of light in nanometres [nm].
    gamma : float, optional
        Gamma correction exponent (default=0.8).
    output : {'norm', 'hex'}, optional
        Specify whether to return an RGB tuple or hexadecimal string and
        RGB tuple (default=(RGB)).

    Returns
    -------
    tuple
        If 'norm' output an RGB tuple with ints in [0, 255] or if 'hex', output
        a hexadecimal string followed by an RGB tuple.
    """

    # == A few notes about colour ==

    # Color   Wavelength(nm) Frequency(THz)
    # Red     620-750        484-400
    # Orange  590-620        508-484
    # Yellow  570-590        526-508
    # Green   495-570        606-526
    # Blue    450-495        668-606
    # Violet  380-450        789-668

    # f is frequency (cycles per second)
    # l (lambda) is wavelength (meters per cycle)
    # e is energy (Joules)
    # h (Plank's constant) = 6.6260695729 x 10^-34 Joule*seconds
    #                      = 6.6260695729 x 10^-34 m^2*kg/seconds
    # c = 299792458 meters per second
    # f = c/l
    # l = c/f
    # e = h*f
    # e = c*h/l

    # List of peak frequency responses for each type of
    # photoreceptor cell in the human eye:
    #     S cone: 437 nm
    #     M cone: 533 nm
    #     L cone: 564 nm
    #     rod:    550 nm in bright daylight, 498 nm when dark adapted.
    #             Rods adapt to low light by becoming more sensitive.
    #             Peak frequency response shifts to 498 nm.

    wav = float(wav)
    if 380 <= wav < 440:
        attenuation = 0.3 + 0.7 * (wav - 380) / (440 - 380)
        R = ((-(wav - 440) / (440 - 380)) * attenuation) ** gamma
        G = 0.0
        B = (1.0 * attenuation) ** gamma
    elif 440 <= wav < 490:
        R = 0.0
        G = ((wav - 440) / (490 - 440)) ** gamma
        B = 1.0
    elif 490 <= wav < 510:
        R = 0.0
        G = 1.0
        B = (-(wav - 510) / (510 - 490)) ** gamma
    elif 510 <= wav < 580:
        R = ((wav - 510) / (580 - 510)) ** gamma
        G = 1.0
        B = 0.0
    elif 580 <= wav < 645:
        R = 1.0
        G = (-(wav - 645) / (645 - 580)) ** gamma
        B = 0.0
    elif 645 <= wav <= 750:
        attenuation = 0.3 + 0.7 * (750 - wav) / (750 - 645)
        R = (1.0 * attenuation) ** gamma
        G = 0.0
        B = 0.0
    else:  # Outside the visible spectrum
        R = 0.0
        G = 0.0
        B = 0.0

    if output == 'norm':
        return (int(R), int(G), int(B))
    elif output == 'hex':
        R *= 255
        R = int(max(0, min(round(R), 255)))
        G *= 255
        G = int(max(0, min(round(G), 255)))
        B *= 255
        B = int(max(0, min(round(B), 255)))
        # return (int(R), int(G), int(B)) # int() truncates towards 0
        return "#{0:02x}{1:02x}{2:02x}".format(R, G, B), (R, G, B)