def seis_sim(sp, amp_ratio=1.5, flength=False, phaseout='all'):
    """
    Generate a simulated seismogram from a given S-P time.

    Will generate spikes separated by a given S-P time, which are then
    convolved with a decaying sine function.  The P-phase is simulated by a
    positive spike of value 1, the S-arrival is simulated by a decaying
    boxcar of maximum amplitude 1.5.  These amplitude ratios can be altered by
    changing the amp_ratio, which is the ratio S amplitude:P amplitude.

    .. note::
        In testing this can achieve 0.3 or greater cross-correlations with
        data.

    :type sp: int
    :param sp: S-P time in samples
    :type amp_ratio: float
    :param amp_ratio: S:P amplitude ratio
    :type flength: int
    :param flength: Fixed length in samples, defaults to False
    :type phaseout: str
    :param phaseout:
        Either 'P', 'S' or 'all', controls which phases to cut around, defaults
        to 'all'. Can only be used with 'P' or 'S' options if flength
        is set.

    :returns: Simulated data.
    :rtype: :class:`numpy.ndarray`
    """
    if flength and 2.5 * sp < flength and 100 < flength:
        additional_length = flength
    elif 2.5 * sp < 100.0:
        additional_length = 100
    else:
        additional_length = 2.5 * sp
    synth = np.zeros(int(sp + 10 + additional_length))
    # Make the array begin 10 samples before the P
    # and at least 2.5 times the S-P samples after the S arrival
    synth[10] = 1.0  # P-spike fixed at 10 samples from start of window
    # The length of the decaying S-phase should depend on the SP time,\
    # Some basic estimations suggest this should be atleast 10 samples\
    # and that the coda should be about 1/10 of the SP time
    S_length = 10 + int(sp // 3)
    S_spikes = np.arange(amp_ratio, 0, -(amp_ratio / S_length))
    # What we actually want, or what appears better is to have a series of\
    # individual spikes, of alternating polarity...
    for i in range(len(S_spikes)):
        if i in np.arange(1, len(S_spikes), 2):
            S_spikes[i] = 0
        if i in np.arange(2, len(S_spikes), 4):
            S_spikes[i] *= -1
    # Put these spikes into the synthetic
    synth[10 + sp:10 + sp + len(S_spikes)] = S_spikes
    # Generate a rough damped sine wave to convolve with the model spikes
    sine_x = np.arange(0, 10.0, 0.5)
    damped_sine = np.exp(-sine_x) * np.sin(2 * np.pi * sine_x)
    # Convolve the spike model with the damped sine!
    synth = np.convolve(synth, damped_sine)
    # Normalize snyth
    synth = synth / np.max(np.abs(synth))
    if not flength:
        return synth
    else:
        if phaseout in ['all', 'P']:
            synth = synth[0:flength]
        elif phaseout == 'S':
            synth = synth[sp:]
            if len(synth) < flength:
                # If this is too short, pad
                synth = np.append(synth, np.zeros(flength - len(synth)))
            else:
                synth = synth[0:flength]
        return synth
Generate a simulated seismogram from a given S-P time.

    Will generate spikes separated by a given S-P time, which are then
    convolved with a decaying sine function.  The P-phase is simulated by a
    positive spike of value 1, the S-arrival is simulated by a decaying
    boxcar of maximum amplitude 1.5.  These amplitude ratios can be altered by
    changing the amp_ratio, which is the ratio S amplitude:P amplitude.

    .. note::
        In testing this can achieve 0.3 or greater cross-correlations with
        data.

    :type sp: int
    :param sp: S-P time in samples
    :type amp_ratio: float
    :param amp_ratio: S:P amplitude ratio
    :type flength: int
    :param flength: Fixed length in samples, defaults to False
    :type phaseout: str
    :param phaseout:
        Either 'P', 'S' or 'all', controls which phases to cut around, defaults
        to 'all'. Can only be used with 'P' or 'S' options if flength
        is set.

    :returns: Simulated data.
    :rtype: :class:`numpy.ndarray`