Zero Crossings

S-functions model zero crossings using the mode work vector (or a DWork vector configured as a mode vector) and the continuous zero-crossing vector. Whether the S-function uses mode or DWork vectors, the concept and implementation are the same. For an example using DWork vectors to model zero crossings, see DWork Mode Vector in the “Using Work Vectors” section. The remainder of this section uses mode vectors to model zero crossings.

Elements of the mode vector are integer values. You specify the number of mode vector elements in mdlInitializeSizes, using ssSetNumModes(S,num). You can then access the mode vector using ssGetModeVector. The mode vector values determine how the mdlOutputs routine operates when the solvers are homing in on zero crossings. The Simulink^® solvers track the zero crossings or state events (i.e., discontinuities in the first derivatives) of some signal, usually a function of an input to your S-function, by looking at the continuous zero crossings. Register the number of continuous zero crossings in mdlInitializeSizes, using ssSetNumNonsampledZCs(S, num), then include an mdlZeroCrossings routine to calculate the continuous zero crossings. The S-function sfun_zc_sat.c contains a zero-crossing example. The remainder of this section describes the portions of this S-function that pertain to zero-crossing detection. For a full description of this example, see Zero-Crossing Detection.

First, mdlInitializeSizes specifies the sizes for the mode and continuous zero-crossing vectors using the following lines of code.

ssSetNumModes(S, DYNAMICALLY_SIZED);
ssSetNumNonsampledZCs(S, DYNAMICALLY_SIZED);

Since the number of modes and continuous zero crossings is dynamically sized, mdlSetWorkWidths must initialize the actual size of these vectors. In this example, shown below, there is one mode vector for each output element and two continuous zero crossings for each mode. In general, the number of continuous zero crossings needed for each mode depends on the number of events that need to be detected. In this case, each output (mode) needs to detect when it hits the upper or the lower bound, hence two continuous zero crossings per mode.

static void mdlSetWorkWidths(SimStruct *S)
{
    int nModes;
    int nNonsampledZCs;

    nModes         = numOutput;
    nNonsampledZCs = 2 * numOutput;
    
    ssSetNumModes(S,nModes);
    ssSetNumNonsampledZCs(S,nNonsampledZCs);
}

Next, mdlOutputs determines which mode the simulation is running in at the beginning of each major time step. The method stores this information in the mode vector so it is available when calculating outputs at both major and minor time steps.

/* Get the mode vector */
int_T *mode = ssGetModeVector(S);

    /* Specify three possible mode values.*/
    enum { UpperLimitEquation, NonLimitEquation, LowerLimitEquation };

    /* Update the mode vector at the beginning of a major time step */
    if ( ssIsMajorTimeStep(S) ) {
       for ( iOutput = 0; iOutput < numOutput; iOutput++ ) {
            if ( *uPtrs[uIdx] > *upperLimit ) {
               /* Upper limit is reached. */
               mode[iOutput] = UpperLimitEquation;

            } else if ( *uPtrs[uIdx] < *lowerLimit ) {
               /* Lower limit is reached. */
               mode[iOutput] = LowerLimitEquation;

            } else {
               /* Output is not limited. */
                mode[iOutput] = NonLimitEquation;
            }

            /* Adjust indices to give scalar expansion. */
            uIdx       += uInc;
            upperLimit += upperLimitInc;
            lowerLimit += lowerLimitInc;
        }
        /* Reset index to input and limits. */
        uIdx       = 0;
        upperLimit = mxGetPr( P_PAR_UPPER_LIMIT );
        lowerLimit = mxGetPr( P_PAR_LOWER_LIMIT );

    } /* end IsMajorTimeStep */

Output calculations in mdlOutputs are done based on the values stored in the mode vector.

for ( iOutput = 0; iOutput < numOutput; iOutput++ ) {
    if ( mode[iOutput] == UpperLimitEquation ) {
       /* Output upper limit. */
       *y++ = *upperLimit;

    } else if ( mode[iOutput] == LowerLimitEquation ) {
       /* Output lower limit. */
       *y++ = *lowerLimit;

    } else {
       /* Output is equal to input */
       *y++ = *uPtrs[uIdx];
    }

After outputs are calculated, the Simulink engine calls mdlZeroCrossings to determine if a zero crossing has occurred. A zero crossing is detected if any element of the continuous zero-crossing vector switches from negative to positive, or positive to negative. If this occurs, the simulation modifies the step size and recalculates the outputs to try to locate the exact zero crossing. For this example, the values for the continuous zero-crossing vectors are calculated as shown below.

static void mdlZeroCrossings(SimStruct *S)
{
    int_T             iOutput;
    int_T             numOutput = ssGetOutputPortWidth(S,0);
    real_T            *zcSignals = ssGetNonsampledZCs(S);
    InputRealPtrsType uPtrs      = ssGetInputPortRealSignalPtrs(S,0);

    /* Set index and increment for the input signal, upper limit, and lower 
     * limit parameters so that each gives scalar expansion if needed. */
    int_T  uIdx          = 0;
    int_T  uInc          = ( ssGetInputPortWidth(S,0) > 1 );
    const real_T *upperLimit   = mxGetPr( P_PAR_UPPER_LIMIT );
    int_T  upperLimitInc = ( mxGetNumberOfElements( P_PAR_UPPER_LIMIT ) > 1 );
    const real_T *lowerLimit   = mxGetPr( P_PAR_LOWER_LIMIT );
    int_T  lowerLimitInc = ( mxGetNumberOfElements( P_PAR_LOWER_LIMIT ) > 1 );

    /*Check if the input has crossed an upper or lower limit */
    for ( iOutput = 0; iOutput < numOutput; iOutput++ ) {
        zcSignals[2*iOutput] = *uPtrs[uIdx] - *upperLimit;
        zcSignals[2*iOutput+1] = *uPtrs[uIdx] - *lowerLimit;

        /* Adjust indices to give scalar expansion if needed */
        uIdx       += uInc;
        upperLimit += upperLimitInc;
        lowerLimit += lowerLimitInc;
    }
}