homopuls.c
Upload User: loeagle
Upload Date: 2013-03-02
Package Size: 1236k
Code Size: 14k
Development Platform:

Matlab

  1. /* $Revision: 1.23 $ */
  2. /*
  3.  * HOMOPULS   A Simulink homonic pulse generator.
  4.  *
  5.  *           Syntax:  [sys, x0] = homopuls(t,x,u,flag,sample,divider,offset)
  6.  *
  7.  * Wes Wang  8/18/1994 revised 1/24/1996.
  8.  * Copyright 1996-2001 The MathWorks, Inc.
  9.  */
  10. #define S_FUNCTION_NAME homopuls
  11. #ifdef MATLAB_MEX_FILE
  12. #include "mex.h"      /* needed for declaration of mexErrMsgTxt */
  13. #endif
  14. /*
  15.  * need to include simstruc.h for the definition of the SimStruct and
  16.  * its associated macro definitions.
  17.  */
  18. #include "simstruc.h"
  19. #include "tmwtypes.h"
  20. #include <math.h>
  21. /* For RTW */
  22. #if defined(RT) || defined(NRT)  
  23. #undef  mexPrintf
  24. #define mexPrintf printf
  25. #endif
  26. /*
  27.  * Defines for easy access of the input parameters
  28.  */
  29. #define NUM_ARGS     3
  30. #define SAMPLE_TIME       ssGetArg(S,0)
  31. #define DIVIDER_EACH      ssGetArg(S,1)
  32. #define OFFSET_EACH       ssGetArg(S,2)
  33. /*
  34.  * mdlInitializeSizes - called to initialize the sizes array stored in
  35.  *                      the SimStruct.  The sizes array defines the
  36.  *                      characteristics (number of inputs, outputs,
  37.  *                      states, etc.) of the S-Function.
  38.  */
  39. static int_T isclose(real_T a, real_T b)
  40. {
  41. #define EPS       4.656612875245797e-10    /*  reciprocal of 2^31-1 */
  42.      real_T diff;
  43. #ifdef MATLAB_MEX_FILE
  44.     if (0) {
  45.       char_T err_msg[255];
  46.       sprintf(err_msg, "a %f, b %f n", a, b);
  47.       mexPrintf(err_msg);
  48.     } 
  49. #endif
  50.     
  51.      diff = fabs(a - b);
  52.      if (diff * 1.e-5 /(fabs(a) + EPS) < EPS)
  53.        return(1);
  54.      else
  55.        return(0);
  56. }
  57. static void mdlInitializeSizes(SimStruct *S)
  58. {
  59.   /*
  60.    * Set-up size information.
  61.    */ 
  62.   int_T NumSampleTime;
  63.   if (ssGetNumArgs(S) == NUM_ARGS) {
  64.     int_T dividerSize;
  65.     if ((mxGetN(SAMPLE_TIME) * mxGetM(SAMPLE_TIME)) > 1) {
  66. #ifdef MATLAB_MEX_FILE
  67.       mexErrMsgTxt("The sample time is a scalar.");
  68. #endif
  69.     }
  70.     
  71.     if ((mxGetN(DIVIDER_EACH) != 1) && (mxGetM(DIVIDER_EACH) != 1)) {
  72. #ifdef MATLAB_MEX_FILE
  73.       mexErrMsgTxt("The divider must be a vector");
  74. #endif
  75.     }
  76.                
  77.     if ((mxGetN(OFFSET_EACH) != 1) && (mxGetM(OFFSET_EACH) != 1)) {
  78. #ifdef MATLAB_MEX_FILE
  79.       mexErrMsgTxt("The offset must be a vector");
  80. #endif
  81.     }
  82.                
  83.     if ((mxGetN(OFFSET_EACH) * (mxGetM(OFFSET_EACH))) != ((mxGetN(DIVIDER_EACH) * (mxGetM(DIVIDER_EACH))))) {
  84. #ifdef MATLAB_MEX_FILE
  85.       mexErrMsgTxt("Divider and Offset must have the same length");
  86. #endif
  87.     }
  88.     
  89.     dividerSize = mxGetN(DIVIDER_EACH) * mxGetM(DIVIDER_EACH);
  90.     NumSampleTime = dividerSize;
  91.     if (dividerSize > 0) {
  92.        real_T sampleTime, sampleTimeI, offsetTimeI;
  93.        int_T    adj, i;
  94.        real_T *past_sample, *past_offset;
  95.        
  96.        past_sample = (real_T *)calloc(dividerSize, sizeof(real_T));
  97.        past_offset = (real_T *)calloc(dividerSize, sizeof(real_T));
  98.        sampleTime = mxGetPr(SAMPLE_TIME)[0];
  99.        NumSampleTime = 0;
  100.        for (i=0; i < dividerSize; i++) {
  101.          offsetTimeI = mxGetPr(OFFSET_EACH)[i];
  102.          sampleTimeI = sampleTime/2./mxGetPr(DIVIDER_EACH)[i];
  103.          if (sampleTimeI <= 0) {
  104. #ifdef MATLAB_MEX_FILE
  105.      mexErrMsgTxt("Sample time must be positive number.");
  106. #endif
  107.          }
  108.          while  (offsetTimeI < 0)
  109.             offsetTimeI += sampleTimeI;
  110.          adj = (int_T)(offsetTimeI / sampleTimeI);
  111.          offsetTimeI = offsetTimeI - ((real_T)adj) * sampleTimeI;
  112. #ifdef MATLAB_MEX_FILE
  113.          if (0) {
  114.            char_T err_msg[255];
  115.            sprintf(err_msg, "Iteration %d, offsetTimeI %f, sampleTimeI %f n", i, offsetTimeI, sampleTimeI);
  116.            mexPrintf(err_msg);
  117.          }     
  118. #endif
  119.  if ((isclose(offsetTimeI, 0.0)) || (isclose(offsetTimeI, sampleTimeI)))
  120.    offsetTimeI = 0;
  121.  if (i > 0) {
  122.    int_T test_flag, ii;
  123.    test_flag = 1;
  124.    for (ii = 0; ii<NumSampleTime; ii++) {
  125. #ifdef MATLAB_MEX_FILE
  126.           if (0) {
  127.                char_T err_msg[255];
  128.                sprintf(err_msg, "sampleTimeI %f, past_sample[ii] %f, offsetTimeI %f, past_offset[ii] %f n", sampleTimeI, past_sample[ii], offsetTimeI, past_offset[ii]);
  129.                mexPrintf(err_msg);
  130.   }
  131. #endif
  132.      if ((isclose(sampleTimeI, past_sample[ii])) && (isclose(offsetTimeI, past_offset[ii])))
  133.         test_flag = 0;
  134.    }
  135.    if (test_flag) {
  136.      past_sample[NumSampleTime] = sampleTimeI;
  137.              past_offset[NumSampleTime] = offsetTimeI;
  138.      NumSampleTime++;
  139.    }
  140.  } else {
  141.    past_sample[NumSampleTime] = sampleTimeI;
  142.    past_offset[NumSampleTime] = offsetTimeI;
  143.            NumSampleTime++;
  144.  }
  145.        }
  146.        free(past_sample);
  147.        free(past_offset);
  148.     }
  149. #ifdef MATLAB_MEX_FILE     
  150.     if (0) {
  151.       char_T err_msg[255];
  152.       sprintf(err_msg, "NumSampleTime %d: n", NumSampleTime);
  153.       mexPrintf(err_msg);
  154.     }
  155. #endif
  156.     ssSetNumContStates(    S, 0);
  157.     ssSetNumDiscStates(    S, 0);
  158.     ssSetNumInputs(        S, 0);
  159.     ssSetNumOutputs(       S, dividerSize);
  160.     ssSetDirectFeedThrough(S, 0);
  161.     ssSetNumInputArgs(     S, NUM_ARGS);
  162.     ssSetNumSampleTimes(   S, NumSampleTime);
  163.     ssSetNumRWork(         S, 4*dividerSize+2); /* store the timing */
  164.     ssSetNumIWork(         S, dividerSize); /* store the last time access. */
  165.     ssSetNumPWork(         S, 0);
  166.   } else {
  167. #ifdef MATLAB_MEX_FILE
  168.     char_T err_msg[256];
  169.     sprintf(err_msg, "Wrong number of input arguments passed to S-function MEX-file.n %d input arguments were passed in when expecting %d input arguments.n", ssGetNumArgs(S) + 4, NUM_ARGS + 4);
  170.     mexErrMsgTxt(err_msg);
  171. #endif
  172.   }
  173. }
  174. /*
  175.  * mdlInitializeSampleTimes - initializes the array of sample times stored in
  176.  *                            the SimStruct associated with this S-Function.
  177.  */
  178. static void mdlInitializeSampleTimes(SimStruct *S)
  179. {
  180.   int_T dividerSize, NumSampleTime;
  181.   int_T NumSmpTm = ssGetNumSampleTimes(S);
  182.   dividerSize = mxGetN(DIVIDER_EACH) * mxGetM(DIVIDER_EACH);
  183.   NumSampleTime = dividerSize;
  184.   if (dividerSize > 0) {
  185.     real_T sampleTime, sampleTimeI, offsetTimeI;
  186.     int_T    adj, i;
  187.     real_T *past_sample, *past_offset;
  188.        
  189.     past_sample = (real_T *)calloc(dividerSize, sizeof(real_T));
  190.     past_offset = (real_T *)calloc(dividerSize, sizeof(real_T));
  191.     sampleTime = mxGetPr(SAMPLE_TIME)[0];
  192.     NumSampleTime = 0;
  193.     for (i=0; i < dividerSize; i++) {
  194.       offsetTimeI = mxGetPr(OFFSET_EACH)[i];
  195.       sampleTimeI = sampleTime/2./mxGetPr(DIVIDER_EACH)[i];
  196.       while  (offsetTimeI < 0)
  197.          offsetTimeI += sampleTimeI;
  198.       adj = (int_T)(offsetTimeI / sampleTimeI);
  199.       offsetTimeI = offsetTimeI - ((real_T)adj) * sampleTimeI;
  200.       if ((isclose(offsetTimeI, 0.0)) || (isclose(offsetTimeI, sampleTimeI)))
  201. offsetTimeI = 0;
  202.       while  (offsetTimeI > sampleTimeI)
  203.         offsetTimeI -= sampleTimeI;
  204.       if (i > 0) {
  205. int_T test_flag, ii;
  206. test_flag = 1;
  207. for (ii = 0; ii<NumSampleTime; ii++) {
  208.   if ((isclose(sampleTimeI, past_sample[ii])) && (isclose(offsetTimeI, past_offset[ii])))
  209.     test_flag = 0;
  210. }
  211. if (test_flag) {
  212.   past_sample[NumSampleTime] = sampleTimeI;
  213.           past_offset[NumSampleTime] = offsetTimeI;
  214. #ifdef MATLAB_MEX_FILE
  215.   if (0) {
  216.       char_T err_msg[255];
  217.       sprintf(err_msg, "Sample Time: %f, Offset Time: %fn", sampleTimeI, offsetTimeI);
  218.       mexPrintf(err_msg);
  219.   } 
  220. #endif    
  221.           ssSetSampleTimeEvent(S, NumSampleTime, sampleTimeI);
  222.           ssSetOffsetTimeEvent(S, NumSampleTime, offsetTimeI);
  223.   NumSampleTime++;
  224. }
  225.       } else {
  226. past_sample[NumSampleTime] = sampleTimeI;
  227. past_offset[NumSampleTime] = offsetTimeI;
  228. #ifdef MATLAB_MEX_FILE
  229. if (0) {
  230.     char_T err_msg[255];
  231.     sprintf(err_msg, "Sample Time: %f, Offset Time: %f", sampleTimeI, offsetTimeI);
  232.     mexPrintf(err_msg);
  233.         }     
  234. #endif
  235.         ssSetSampleTimeEvent(S, NumSampleTime, sampleTimeI);
  236.         ssSetOffsetTimeEvent(S, NumSampleTime, offsetTimeI);
  237.         NumSampleTime++;
  238.       }
  239.       if (NumSampleTime > NumSmpTm) {
  240. #ifdef MATLAB_MEX_FILE
  241.         mexErrMsgTxt("Check your sample time and offset time. It is not consistant.");
  242. #endif
  243.       }
  244.     }
  245.     free(past_sample);
  246.     free(past_offset);
  247.   }
  248. #ifdef MATLAB_MEX_FILE
  249.   if (0) {
  250.       char_T err_msg[255];
  251.       sprintf(err_msg, "NumSampleTime %d: n", NumSampleTime);
  252.       mexPrintf(err_msg);
  253.   }
  254. #endif
  255.     ssSetNumContStates(    S, 0);
  256.     ssSetNumDiscStates(    S, 0);
  257. }
  258. /*
  259.  * mdlInitializeConditions - initializes the states for the S-Function
  260.  */
  261. static void mdlInitializeConditions(real_T *x0, SimStruct *S)
  262. {
  263.   int_T dividerSize = mxGetN(DIVIDER_EACH) * mxGetM(DIVIDER_EACH);
  264.   int_T offsetSize  = mxGetN(OFFSET_EACH) * mxGetM(OFFSET_EACH);
  265.   real_T *hit_base   = ssGetRWork(S);
  266.   real_T *next_hit   = ssGetRWork(S) + 1;    
  267.   real_T *last_value = ssGetRWork(S) + dividerSize + 1;
  268.   real_T *increment  = ssGetRWork(S) + dividerSize * 2 + 1;
  269.   real_T *adj_offset = ssGetRWork(S) + dividerSize * 3 + 1;
  270.   real_T *tolerance  = ssGetRWork(S) + dividerSize * 4 + 1;
  271.   int_T  *reverse    = ssGetIWork(S);        
  272.   real_T  sampleTime, offsetTime, offsetMin, tmp, tol;
  273.   int_T   i, adj;
  274.   
  275.   dividerSize = (dividerSize < offsetSize) ? dividerSize : offsetSize;
  276.   offsetMin = mxGetPr(OFFSET_EACH)[0];
  277.   sampleTime = mxGetPr(SAMPLE_TIME)[0];
  278.   /* 
  279.    * Initialize the last_accs to all zeros.
  280.    */
  281.   tol = sampleTime;
  282.   for (i = 0; i < dividerSize; i++) {
  283.     *last_value++ = 0.;
  284.     
  285.     offsetTime = mxGetPr(OFFSET_EACH)[i];        
  286.     offsetMin = (offsetMin < offsetTime) ? offsetMin : offsetTime;
  287.     
  288.     next_hit[i] = offsetTime;
  289.     
  290.     tmp = sampleTime/2./mxGetPr(DIVIDER_EACH)[i];
  291.     increment[i] = tmp;
  292.     
  293.     tol = (tol < tmp) ? tol : tmp;
  294.     
  295.     adj = (int_T)(offsetTime / tmp);
  296.     tmp = offsetTime - ((real_T)adj) * tmp;
  297.     /*
  298.     if (isclose(tmp, 0.0)) {
  299.       reverse[i] = 0;
  300.     } else if (isclose(tmp, offsetTime - ((real_T)(((int_T)(offsetTime/tmp/2)))) * tmp * 2)) {
  301.     mexPrintf("offsetTime=%f, increment[%d]=%f.",offsetTime,i,increment[i]);
  302.     */
  303.     if (increment[i] > offsetTime - 2 * increment[i] * floor(offsetTime/2./increment[i])) {
  304.       reverse[i] = 0;
  305.     } else {
  306.       reverse[i] = 1;
  307.     }
  308. #ifdef MATLAB_MEX_FILE
  309.     if (0) {
  310.       char_T err_msg[255];
  311.       sprintf(err_msg, "reverse[%d]=%d; n", i, reverse[i]);
  312.       mexPrintf(err_msg);
  313.     }
  314. #endif
  315.     adj_offset[i] = tmp;
  316.     if (tmp > 0)
  317.       tol = (tol < tmp) ? tol : tmp;        
  318.   }
  319. #ifdef MATLAB_MEX_FILE
  320.   if (0) {
  321.     char_T err_msg[255];
  322.     sprintf(err_msg, "n");
  323.     mexPrintf(err_msg);
  324.   }
  325. #endif
  326.   *hit_base = offsetMin + sampleTime;
  327.   *hit_base = sampleTime;
  328.   *tolerance = tol / 100;
  329. }
  330. /*
  331.  * mdlOutputs - computes the outputs of the S-Function
  332.  */
  333. static void mdlOutputs(real_T *y, const real_T *x, const real_T *u, SimStruct *S, int_T tid)
  334. {
  335.   int_T dividerSize = mxGetN(DIVIDER_EACH) * mxGetM(DIVIDER_EACH);
  336.   real_T *hit_base      = ssGetRWork(S);
  337.   real_T *next_hit      = ssGetRWork(S) + 1;    
  338.   real_T *last_value    = ssGetRWork(S) + dividerSize + 1;
  339.   real_T *increment     = ssGetRWork(S) + dividerSize * 2 + 1;
  340.   real_T *adj_offset    = ssGetRWork(S) + dividerSize * 3 + 1;
  341.   real_T *tolerance     = ssGetRWork(S) + dividerSize * 4 + 1;
  342.   int_T    *reverse = ssGetIWork(S);
  343.   
  344.   real_T time_clock, tol;
  345.   int_T i;
  346.   
  347.   tol = *tolerance;
  348.   
  349.   time_clock = ssGetT(S) + tol;
  350.   for (i = 0; i < dividerSize; i++) {
  351.     if (time_clock >= next_hit[i]) {
  352.       int_T    num;
  353.       real_T sampleTime;
  354.       
  355.       sampleTime = mxGetPr(SAMPLE_TIME)[0];        
  356.       num = (int_T)((time_clock - *hit_base - adj_offset[i] - sampleTime) / increment[i]);
  357.       num = num%2;
  358.       if (num)
  359. last_value[i] = 1.;
  360.       else
  361. last_value[i] = 0.;
  362. #ifdef MATLAB_MEX_FILE
  363.       if (0) {
  364. char_T err_msg[255];
  365. sprintf(err_msg, "reverse[%d]=%d; last_value=%d; n", i, reverse[i], (int_T)last_value[i]);
  366. mexPrintf(err_msg);
  367.       }
  368. #endif
  369.       if (reverse[i])
  370. last_value[i] = ((int_T)last_value[i] + 1) % 2;
  371. #ifdef MATLAB_MEX_FILE
  372.       if (0) {
  373. char_T err_msg[255];
  374. sprintf(err_msg, "new_last_value=%d;n ", (int_T)last_value[i]);
  375. mexPrintf(err_msg);
  376.       }
  377. #endif
  378.       next_hit[i] = next_hit[i] + increment[i];
  379.     }
  380.     y[i] = last_value[i];
  381.   }
  382.   
  383.   /*    sprintf(err_msg, "n ");
  384. mexPrintf(err_msg);
  385. for (i = 0; i < dividerSize; i++) {
  386. sprintf(err_msg, "output_y[%d]= %f, ", i, y[i]);
  387. mexPrintf(err_msg);
  388. }
  389. */
  390.   /* adjust the "current hit" every sample cycle.    */
  391.   if (time_clock > *hit_base){
  392.     real_T sampleTime;
  393.     sampleTime = mxGetPr(SAMPLE_TIME)[0];        
  394.     for (i = 0; i < dividerSize; i++) {
  395.       if (time_clock > mxGetPr(OFFSET_EACH)[i]) {
  396. if (adj_offset[i] <= 0) {
  397.   next_hit[i] = *hit_base + increment[i];
  398. } else {
  399.   next_hit[i] = *hit_base + adj_offset[i];
  400. }
  401.       }
  402.     }
  403.     *hit_base = *hit_base +sampleTime;        
  404.   }    
  405. }
  406. /*
  407.  * mdlUpdate - computes the discrete states of the S-Function
  408.  */
  409. static void mdlUpdate(real_T *x, const real_T *u, SimStruct *S, int_T tid)
  410. {
  411. }
  412. /*
  413.  * mdlDerivatives - computes the derivatives of the S-Function
  414.  */
  415. static void mdlDerivatives(real_T *dx, const real_T *x, const real_T *u, SimStruct *S, int_T tid)
  416. {
  417. }
  418. /*
  419.  * mdlTerminate - called at termination of model execution.
  420.  */
  421. static void mdlTerminate(SimStruct *S)
  422. {
  423. }
  424. #ifdef     MATLAB_MEX_FILE    /* Is this file being compiled as a MEX-file? */
  425. #include "simulink.c"      /* MEX-File interface mechanism */
  426. #else
  427. #include "cg_sfun.h"       /* Code generation registration function */
  428. #endif