def step_generator(intervals, parameter):
    """




    :param intervals: the 1- or 2-D intervals to be used
    :param parameter:
    """

    intervals = np.atleast_2d(intervals)

    # need to make sure the shape is right
    # assert self._intervals.shape


    if intervals.shape[0] > 1 and intervals.shape[1] == 2:

        n_intervals = intervals.shape[0]




        is_2d = True

    elif intervals.shape[0] == 1:

        n_intervals = intervals.shape[1]
        intervals = intervals[0]

        is_2d = False

    else:

        raise RuntimeError("These intervals are not yet supported")

    parameter_min = parameter.min_value
    parameter_max = parameter.max_value
    initial_value = parameter.value

    if is_2d:

        func = StepFunctionUpper()

        for i in range(n_intervals - 1):

            func += StepFunctionUpper()

        for i, interval in enumerate(intervals):

            i=i+1

            func.free_parameters['value_%d' % i].value = initial_value
            func.free_parameters['value_%d' % i].min_value = parameter_min
            func.free_parameters['value_%d' % i].max_value = parameter_max

            func.parameters['upper_bound_%d' % i].value = interval[1]

            func.parameters['lower_bound_%d' % i].value = interval[0]


    else:

        func = DiracDelta()

        for i in range(n_intervals - 1):

            func += DiracDelta()

        for i, interval in enumerate(intervals):

            i=i+1

            func.free_parameters['k_%d' % i].value = initial_value
            func.free_parameters['k_%d' % i].min_value = parameter_min
            func.free_parameters['k_%d' % i].max_value = parameter_max

            func.parameters['zero_point_%d' % i].value = interval



    return func



# how it would work:


triggerName = 'bn080916009'
ra = 121.8
dec = -61.3




data_list = DataList(nai6,nai7,nai9,bgo1,nai6_2,nai7_2,nai9_2,bgo1_2 )







cpl = Cutoff_powerlaw()

syn = SynchrotronNumerical()





GRB = PointSource( triggerName, ra, dec, spectral_shape=cpl )

model = Model( GRB )


#syn.B=1E7
#syn.K = 1E25
#syn.gamma_cool=7E4
#syn.gamma_max = 1E6
time = IndependentVariable("time",13.0, unit='s')
model

law = Powerlaw(K=500,index=-1)

#law = Line(a=-0.02,b=-2.0)
model.add_independent_variable(time)
model.link(model.bn080916009.spectrum.main.Cutoff_powerlaw.xc,time,law)

ff = step_generator([12.5,17.5],cpl.index)



model.link(model.bn080916009.spectrum.main.Cutoff_powerlaw.index,time,ff)


ff = step_generator([12.5,17.5],cpl.K)



model.link(model.bn080916009.spectrum.main.Cutoff_powerlaw.K,time,ff)

nai6.external_property(time,12.5)
nai7.external_property(time,12.5)
nai9.external_property(time,12.5)
bgo1.external_property(time,12.5)


nai6_2.external_property(time,17.5)
nai7_2.external_property(time,17.5)
nai9_2.external_property(time,17.5)
bgo1_2.external_property(time,17.5)




class StepFunctionUpper(Function1D):
    r"""
        description :

            A function which is constant on the interval lower_bound - upper_bound and 0 outside the interval. The
            upper interval is open.

        latex : $ f(x)=\begin{cases}0 & x < \text{lower_bound} \\\text{value} & \text{lower_bound} \le x \le \text{upper_bound} \\ 0 & x > \text{upper_bound} \end{cases}$

        parameters :

            lower_bound :

                desc : Lower bound for the interval
                initial value : 0
                min : -np.inf
                max : np.inf
                fix : yes

            upper_bound :

                desc : Upper bound for the interval
                initial value : 1
                min : -np.inf
                max : np.inf
                fix : yes

            value :

                desc : Value in the interval
                initial value : 1.0

        tests :
            - { x : 0.5, function value: 1.0, tolerance: 1e-20}
            - { x : -0.5, function value: 0, tolerance: 1e-20}

        """

    __metaclass__ = FunctionMeta

    def _set_units(self, x_unit, y_unit):
        # Lower and upper bound has the same unit as x
        self.lower_bound.unit = x_unit
        self.upper_bound.unit = x_unit

        # value has the same unit as y
        self.value.unit = y_unit

    def evaluate(self, x, lower_bound, upper_bound, value):
        # The value * 0 is to keep the units right

        result = np.zeros(x.shape) * value * 0

        idx = (x >= lower_bound) & (x < upper_bound)
        result[idx] = value

        return result