#**
#**   Intermodulation Products
#**                                         In this example:    f1 = 5 MHz
#**   |                                                         f2 = 6 MHz
#**   |
#**   |
#**   |                                              [2]  Second Order IMD Products
#**   |                  f1   f2                     [3]  Third  Order IMD Products
#**   |                   :   :
#**   |                   :   :
#**   |                   :   :
#**   |  [2]              :   :                 [2]
#**   |                   :   :                                     [3] [3]
#**   | f2-f1             :   :                f1+f2
#**   |   :               :   :                  :               2f1+f2 2f2+f1
#**   |   :          [3]  :   :  [3]        2f1  :  2f2              :   :
#**   |   :               :   :              :   :   :          3f1  :   :  3f2
#**   |   :        2f1-f2 :   :2f2-f1        :   :   :           :   :   :   :
#**  dc   :           :   :   :   :          :   :   :           :   :   :   :
#**   :   :           :   :   :   :          :   :   :           :   :   :   :
#**   :___:___________:___:___:___:__________:___:___:___________:___:___:___:_____
#**
#**   0   1           4   5   6   7         10  11  12          15  16  17  18  MHz
#**
#**   !                   !   !              !       !           !           !
#** offset             FUNDAMENTALS       second harmonics      third harmonics
#**
#**
#**
#** The Harmonics:
#**
#**    y[t] = k0 + k1*x[t] + k2*x[t]^2 + k3*x[t]^3       where x[t] = A*Cos[wt]
#**
#**    y[t] = k0 + 0.50*k2*A^2
#**              + (A*k1 + 0.75*A^3*k3) * Cos[  wt]
#**              + 0.50*k2*A^2          * Cos[2*wt]
#**              + 0.25*k3*A^3          * Cos[3*wt]
#**
#**
#** The Intermodulation Products:
#**
#**    y[t] = k0 + k1*x[t] + k2*x[t]^2 + k3*x[t]^3       where x[t] = A0 + A1*Cos[w1t] + A2*Cos[w2t]
#**
#**    y[t] = k0 + A0*k1 + (A0^2 + 0.50*A1^2 + 0.50*A2^2)*k2 + (A0^3 + 1.50*A0*A1^2 + 1.50*A0*A2^2)*k3
#**              + A1*k1 + 2*A0*A1*k2 + (0.75*A1^3 + 1.50*A1*A2^2 + 3*A0^2*A1)*k3          * Cos[  w1t]
#**              + A2*k1 + 2*A0*A2*k2 + (0.75*A2^3 + 1.50*A2*A1^2 + 3*A0^2*A2)*k3          * Cos[  w2t]
#**              + 0.50*A1^2*k2 + 1.5*A0*A1^2*k3                                           * Cos[2*w1t]
#**              + 0.50*A2^2*k2 + 1.5*A0*A2^2*k3                                           * Cos[2*w2t]
#**              + 0.25*A1^3*k3                                                            * Cos[3*w1t]
#**              + 0.25*A2^3*k3                                                            * Cos[3*w2t]
#**              + A1*A2*k2 + 3*A0*A1*A2*k3                                                * Cos[  w1t-  w2t]
#**              + A1*A2*k2 + 3*A0*A1*A2*k3                                                * Cos[  w1t+  w2t]
#**              + 0.75*A1*A2^2*k3                                                         * Cos[  w1t-2*w2t]
#**              + 0.75*A2*A1^2*k3                                                         * Cos[2*w1t-  w2t]
#**              + 0.75*A2*A1^2*k3                                                         * Cos[2*w1t+  w2t]
#**              + 0.75*A1*A2^2*k3                                                         * Cos[  w1t+2*w2t]
#**
#**
#** NB:  As distorsion is coming from signal compression at high amplitudes or smooth clipping
#**      parameters, at least one of the parameters k2 or k3 is negative.