1 
2 /*
3  * TRANSLINE.cpp - base for a transmission line implementation
4  *
5  * Copyright (C) 2005 Stefan Jahn <stefan@lkcc.org>
6  * Modified for Kicad: 2018 jean-pierre.charras
7  *
8  * This program is free software; you can redistribute it and/or
9  * modify it under the terms of the GNU General Public License
10  * as published by the Free Software Foundation; either version 2
11  * of the License, or (at your option) any later version.
12  *
13  * This program is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
16  * GNU General Public License for more details.
17  *
18  * You should have received a copy of the GNU General Public License
19  * along with this package; see the file COPYING.  If not, write to
20  * the Free Software Foundation, Inc., 51 Franklin Street - Fifth Floor,
21  * Boston, MA 02110-1301, USA.
22  *
23  */
24 
25 #include <cmath>
26 #include <limits>
27 
28 #include <wx/colour.h>
29 #include <wx/settings.h>
30 
31 #include "transline.h"
32 #include "units.h"
33 
34 #ifndef INFINITY
35 #define INFINITY std::numeric_limits<double>::infinity()
36 #endif
37 
38 
39 #ifndef M_PI_2
40 #define M_PI_2 ( M_PI / 2 )
41 #endif
42 
43 
44 // Functions to Read/Write parameters in pcb_calculator main frame:
45 // They are wrapper to actual functions, so all transline functions do not
46 // depend on Graphic User Interface
47 void SetPropertyInDialog( enum PRMS_ID aPrmId, double value );
48 
49 /* Puts the text into the given result line.
50 */
51 void SetResultInDialog( int line, const char* text );
52 
53 /* print aValue into the given result line.
54 */
55 void SetResultInDialog( int aLineNumber, double aValue, const char* aText );
56 
57 /* Returns a named property value. */
58 double GetPropertyInDialog( enum PRMS_ID aPrmId );
59 
60 // Returns true if the param aPrmId is selected
61 // Has meaning only for params that have a radio button
62 bool IsSelectedInDialog( enum PRMS_ID aPrmId );
63 
64 /** Function SetPropertyBgColorInDialog
65  *  Set the background color of a parameter
66  *  @param aPrmId = param id to set
67  *  @param aCol = new color
68  */
69 void SetPropertyBgColorInDialog( enum PRMS_ID aPrmId, const KIGFX::COLOR4D* aCol );
70 
71 
72 /* Constructor creates a transmission line instance. */
TRANSLINE()73 TRANSLINE::TRANSLINE()
74 {
75     m_parameters[MURC_PRM] = 1.0;
76     m_Name                 = nullptr;
77     ang_l                  = 0.0;       // Electrical length in angle
78     len                    = 0.0;       // length of line
79     er_eff                 = 1.0;       // effective dielectric constant
80     Init();
81 }
82 
83 
84 /* Destructor destroys a transmission line instance. */
~TRANSLINE()85 TRANSLINE::~TRANSLINE()
86 {
87 }
88 
89 
Init()90 void TRANSLINE::Init()
91 {
92     wxColour wxcol = wxSystemSettings::GetColour( wxSYS_COLOUR_WINDOW );
93     okCol          = KIGFX::COLOR4D( wxcol );
94     okCol.r        = wxcol.Red() / 255.0;
95     okCol.g        = wxcol.Green() / 255.0;
96     okCol.b        = wxcol.Blue() / 255.0;
97     int i;
98     // Initialize these variables mainly to avoid warnings from a static analyzer
99     for( i = 0; i < EXTRA_PRMS_COUNT; ++i )
100     {
101         m_parameters[i] = 0;
102     }
103 }
104 
105 
106 /* Sets a named property to the given value, access through the
107  *  application.
108  */
setProperty(enum PRMS_ID aPrmId,double value)109 void TRANSLINE::setProperty( enum PRMS_ID aPrmId, double value )
110 {
111     SetPropertyInDialog( aPrmId, value );
112 }
113 
114 
115 /*
116  *Returns true if the param aPrmId is selected
117  * Has meaning only for params that have a radio button
118  */
isSelected(enum PRMS_ID aPrmId)119 bool TRANSLINE::isSelected( enum PRMS_ID aPrmId )
120 {
121     return IsSelectedInDialog( aPrmId );
122 }
123 
124 
125 /* Puts the text into the given result line.
126 */
setResult(int line,const char * text)127 void TRANSLINE::setResult( int line, const char* text )
128 {
129     SetResultInDialog( line, text );
130 }
131 
132 
setResult(int line,double value,const char * text)133 void TRANSLINE::setResult( int line, double value, const char* text )
134 {
135     SetResultInDialog( line, value, text );
136 }
137 
138 
139 /* Returns a property value. */
getProperty(enum PRMS_ID aPrmId)140 double TRANSLINE::getProperty( enum PRMS_ID aPrmId )
141 {
142     return GetPropertyInDialog( aPrmId );
143 }
144 
145 
146 /** @function getProperties
147  *
148  *  Get all properties from the UI. Computes some extra ones.
149  **/
getProperties()150 void TRANSLINE::getProperties()
151 {
152     for( int i = 0; i < DUMMY_PRM; ++i )
153     {
154         m_parameters[i] = getProperty( (PRMS_ID) i );
155         setErrorLevel( (PRMS_ID) i, TRANSLINE_OK );
156     }
157 
158     m_parameters[SIGMA_PRM]       = 1.0 / getProperty( RHO_PRM );
159     m_parameters[EPSILON_EFF_PRM] = 1.0;
160     m_parameters[SKIN_DEPTH_PRM]  = skin_depth();
161 }
162 
163 
164 /** @function checkProperties
165  *
166  *  Checks the input parameters (ie: negative length).
167  *  Does not check for incompatibility between values as this depends on the line shape.
168  **/
checkProperties()169 void TRANSLINE::checkProperties()
170 {
171     // Do not check for values that are results of analyzing / synthesizing
172     // Do not check for transline specific incompatibilities ( like " conductor height should be lesser than dielectric height")
173     if( !std::isfinite( m_parameters[EPSILONR_PRM] ) || m_parameters[EPSILONR_PRM] <= 0 )
174         setErrorLevel( EPSILONR_PRM, TRANSLINE_WARNING );
175 
176     if( !std::isfinite( m_parameters[TAND_PRM] ) || m_parameters[TAND_PRM] < 0 )
177         setErrorLevel( TAND_PRM, TRANSLINE_WARNING );
178 
179     if( !std::isfinite( m_parameters[RHO_PRM] ) || m_parameters[RHO_PRM] < 0 )
180         setErrorLevel( RHO_PRM, TRANSLINE_WARNING );
181 
182     if( !std::isfinite( m_parameters[H_PRM] ) || m_parameters[H_PRM] < 0 )
183         setErrorLevel( H_PRM, TRANSLINE_WARNING );
184 
185     if( !std::isfinite( m_parameters[TWISTEDPAIR_TWIST_PRM] )
186             || m_parameters[TWISTEDPAIR_TWIST_PRM] < 0 )
187         setErrorLevel( TWISTEDPAIR_TWIST_PRM, TRANSLINE_WARNING );
188 
189     if( !std::isfinite( m_parameters[STRIPLINE_A_PRM] ) || m_parameters[STRIPLINE_A_PRM] <= 0 )
190         setErrorLevel( STRIPLINE_A_PRM, TRANSLINE_WARNING );
191 
192     if( !std::isfinite( m_parameters[H_T_PRM] ) || m_parameters[H_T_PRM] <= 0 )
193         setErrorLevel( H_T_PRM, TRANSLINE_WARNING );
194 
195     // How can we check ROUGH_PRM ?
196 
197     if( !std::isfinite( m_parameters[MUR_PRM] ) || m_parameters[MUR_PRM] < 0 )
198         setErrorLevel( MUR_PRM, TRANSLINE_WARNING );
199 
200     if( !std::isfinite( m_parameters[TWISTEDPAIR_EPSILONR_ENV_PRM] )
201             || m_parameters[TWISTEDPAIR_EPSILONR_ENV_PRM] <= 0 )
202         setErrorLevel( TWISTEDPAIR_EPSILONR_ENV_PRM, TRANSLINE_WARNING );
203 
204     if( !std::isfinite( m_parameters[MURC_PRM] ) || m_parameters[MURC_PRM] < 0 )
205         setErrorLevel( MURC_PRM, TRANSLINE_WARNING );
206 
207     if( !std::isfinite( m_parameters[FREQUENCY_PRM] ) || m_parameters[FREQUENCY_PRM] <= 0 )
208         setErrorLevel( FREQUENCY_PRM, TRANSLINE_WARNING );
209 }
210 
analyze()211 void TRANSLINE::analyze()
212 {
213     getProperties();
214     checkProperties();
215     calcAnalyze();
216     showAnalyze();
217     show_results();
218 }
219 
synthesize()220 void TRANSLINE::synthesize()
221 {
222     getProperties();
223     checkProperties();
224     calcSynthesize();
225     showSynthesize();
226     show_results();
227 }
228 
229 
230 /**
231  * @function skin_depth
232  * calculate skin depth
233  *
234  * \f$ \frac{1}{\sqrt{ \pi \cdot f \cdot \mu \cdot \sigma }} \f$
235  */
236 #include <cstdio>
skin_depth()237 double TRANSLINE::skin_depth()
238 {
239     double depth;
240     depth = 1.0
241             / sqrt( M_PI * m_parameters[FREQUENCY_PRM] * m_parameters[MURC_PRM] * MU0
242                     * m_parameters[SIGMA_PRM] );
243     return depth;
244 }
245 
246 
247 /* *****************************************************************
248 **********                                             **********
249 **********          mathematical functions             **********
250 **********                                             **********
251 ***************************************************************** */
252 
253 #define NR_EPSI 2.2204460492503131e-16
254 
255 /* The function computes the complete elliptic integral of first kind
256  *  K() and the second kind E() using the arithmetic-geometric mean
257  *  algorithm (AGM) by Abramowitz and Stegun.
258  */
ellipke(double arg,double & k,double & e)259 void TRANSLINE::ellipke( double arg, double& k, double& e )
260 {
261     int iMax = 16;
262 
263     if( arg == 1.0 )
264     {
265         k = INFINITY; // infinite
266         e = 0;
267     }
268     else if( std::isinf( arg ) && arg < 0 )
269     {
270         k = 0;
271         e = INFINITY; // infinite
272     }
273     else
274     {
275         double a, b, c, fr, s, fk = 1, fe = 1, t, da = arg;
276         int    i;
277 
278         if( arg < 0 )
279         {
280             fk = 1 / sqrt( 1 - arg );
281             fe = sqrt( 1 - arg );
282             da = -arg / ( 1 - arg );
283         }
284 
285         a  = 1;
286         b  = sqrt( 1 - da );
287         c  = sqrt( da );
288         fr = 0.5;
289         s  = fr * c * c;
290 
291         for( i = 0; i < iMax; i++ )
292         {
293             t = ( a + b ) / 2;
294             c = ( a - b ) / 2;
295             b = sqrt( a * b );
296             a = t;
297             fr *= 2;
298             s += fr * c * c;
299 
300             if( c / a < NR_EPSI )
301                 break;
302         }
303 
304         if( i >= iMax )
305         {
306             k = 0;
307             e = 0;
308         }
309         else
310         {
311             k = M_PI_2 / a;
312             e = M_PI_2 * ( 1 - s ) / a;
313             if( arg < 0 )
314             {
315                 k *= fk;
316                 e *= fe;
317             }
318         }
319     }
320 }
321 
322 
323 /* We need to know only K(k), and if possible KISS. */
ellipk(double k)324 double TRANSLINE::ellipk( double k )
325 {
326     double r, lost;
327 
328     ellipke( k, r, lost );
329     return r;
330 }
331 
332 #define MAX_ERROR 0.000001
333 
334 /**
335  * @function minimizeZ0Error1D
336  *
337  * Tries to find a parameter that minimizes the error ( on Z0 ).
338  * This function only works with a single parameter.
339  * Calls @ref calcAnalyze several times until the error is acceptable.
340  * While the error is unnacceptable, changes slightly the parameter.
341  *
342  * This function does not change Z0 / Angl_L.
343  *
344  * @param avar Parameter to synthesize
345  * @return 'true' if error < MAX_ERROR, else 'false'
346  */
347 
348 
minimizeZ0Error1D(double * aVar)349 bool TRANSLINE::minimizeZ0Error1D( double* aVar )
350 {
351     double Z0_dest, Z0_current, Z0_result, angl_l_dest, increment, slope, error;
352     int    iteration;
353 
354     if( !std::isfinite( m_parameters[Z0_PRM] ) )
355     {
356         *aVar = NAN;
357         return false;
358     }
359 
360     if( ( !std::isfinite( *aVar ) ) || ( *aVar == 0 ) )
361         *aVar = 0.001;
362 
363     /* required value of Z0 */
364     Z0_dest = m_parameters[Z0_PRM];
365 
366     /* required value of angl_l */
367     angl_l_dest = m_parameters[ANG_L_PRM];
368 
369     /* Newton's method */
370     iteration = 0;
371 
372     /* compute parameters */
373     calcAnalyze();
374     Z0_current = m_parameters[Z0_PRM];
375 
376     error = fabs( Z0_dest - Z0_current );
377 
378     while( error > MAX_ERROR )
379     {
380         iteration++;
381         increment = *aVar / 100.0;
382         *aVar += increment;
383         /* compute parameters */
384         calcAnalyze();
385         Z0_result = m_parameters[Z0_PRM];
386         /* f(w(n)) = Z0 - Z0(w(n)) */
387         /* f'(w(n)) = -f'(Z0(w(n))) */
388         /* f'(Z0(w(n))) = (Z0(w(n)) - Z0(w(n+delw))/delw */
389         /* w(n+1) = w(n) - f(w(n))/f'(w(n)) */
390         slope = ( Z0_result - Z0_current ) / increment;
391         slope = ( Z0_dest - Z0_current ) / slope - increment;
392         *aVar += slope;
393 
394         if( *aVar <= 0.0 )
395             *aVar = increment;
396 
397         /* find new error */
398         /* compute parameters */
399         calcAnalyze();
400         Z0_current = m_parameters[Z0_PRM];
401         error      = fabs( Z0_dest - Z0_current );
402 
403         if( iteration > 100 )
404             break;
405     }
406 
407     /* Compute one last time, but with correct length */
408     m_parameters[Z0_PRM]       = Z0_dest;
409     m_parameters[ANG_L_PRM]    = angl_l_dest;
410     m_parameters[PHYS_LEN_PRM] = C0 / m_parameters[FREQUENCY_PRM]
411                                  / sqrt( m_parameters[EPSILON_EFF_PRM] ) * m_parameters[ANG_L_PRM]
412                                  / 2.0 / M_PI; /* in m */
413     calcAnalyze();
414 
415     /* Restore parameters */
416     m_parameters[Z0_PRM]       = Z0_dest;
417     m_parameters[ANG_L_PRM]    = angl_l_dest;
418     m_parameters[PHYS_LEN_PRM] = C0 / m_parameters[FREQUENCY_PRM]
419                                  / sqrt( m_parameters[EPSILON_EFF_PRM] ) * m_parameters[ANG_L_PRM]
420                                  / 2.0 / M_PI; /* in m */
421     return error <= MAX_ERROR;
422 }
423 /**
424  * @function setErrorLevel
425  *
426  * set an error / warning level for a given parameter.
427  *
428  * @see TRANSLINE_OK
429  * @see TRANSLINE_WARNING
430  * @see TRANSLINE_ERROR
431  *
432  * @param aP parameter
433  * @param aErrorLevel Error level
434  */
setErrorLevel(PRMS_ID aP,char aErrorLevel)435 void TRANSLINE::setErrorLevel( PRMS_ID aP, char aErrorLevel )
436 {
437     switch( aErrorLevel )
438     {
439     case TRANSLINE_WARNING:  SetPropertyBgColorInDialog( aP, &warnCol );  break;
440     case TRANSLINE_ERROR:    SetPropertyBgColorInDialog( aP, &errCol );   break;
441     default:                 SetPropertyBgColorInDialog( aP, &okCol );    break;
442     }
443 }
444