def snfpsf(wave, psfparams, header, psftype):
    """Create a 3-d PSF based on SNFactory-specific parameterization of
    Gaussian + Moffat PSF parameters and ADR."""

    # Get Gaussian+Moffat parameters at each wavelength.
    relwave = wave / REFWAVE - 1.0
    ellipticity = abs(psfparams[0]) * np.ones_like(wave)
    alpha = np.abs(psfparams[1] +
                   psfparams[2] * relwave +
                   psfparams[3] * relwave**2)

    # correlated parameters (coefficients determined externally)
    sigma = 0.545 + 0.215 * alpha  # Gaussian parameter
    beta  = 1.685 + 0.345 * alpha  # Moffat parameter
    eta   = 1.040 + 0.0   * alpha  # gaussian ampl. / moffat ampl.

    # Atmospheric differential refraction (ADR): Because of ADR,
    # the center of the PSF will be different at each wavelength,
    # by an amount that we can determine (pretty well) from the
    # atmospheric conditions and the pointing and angle of the
    # instrument. We calculate the offsets here as a function of
    # observation and wavelength and input these to the model.

    # Correction to parallactic angle and airmass for 2nd-order effects
    # such as MLA rotation, mechanical flexures or finite-exposure
    # corrections. These values have been trained on faint-std star
    # exposures.
    #
    # `predict_adr_params` uses 'AIRMASS', 'PARANG' and 'CHANNEL' keys
    # in input dictionary.
    delta, theta = Hyper_PSF3D_PL.predict_adr_params(header)

    # check for crazy values of pressure and temperature, and assign default
    # values.
    pressure = header.get('PRESSURE', 617.)
    if not 550. < pressure < 650.:
        pressure = 617.
    temp = header.get('TEMP', 2.)
    if not -20. < temp < 20.:
        temp = 2.

    adr = ADR(pressure, temp, lref=REFWAVE, delta=delta, theta=theta)
    adr_refract = adr.refract(0, 0, wave, unit=SPAXEL_SIZE)