@@ -15,6 +15,8 @@ import {
1515 functionFromPolynomialTermsArray ,
1616 areAllTfTermsNumbers ,
1717 roundDecimal ,
18+ toleranceNumericalAnalysisSmall ,
19+ tolerancePhaseAdjustmentLarge ,
1820} from "../../util/commons.js" ;
1921import { logMessages } from "../../util/loggingService.js" ;
2022import {
@@ -159,8 +161,16 @@ export default class BodePlot {
159161
160162 //compute expected steep phase shifts due to polynomial zeros & poles at y-axis
161163 const expectedSteepPhaseShiftsMap = new Map ( ) ;
162- const zerosAtYAxis = this . #zeros. filter ( ( x ) => x [ 0 ] === 0 && x [ 1 ] !== 0 ) ;
163- const polesAtYAxis = this . #poles. filter ( ( x ) => x [ 0 ] === 0 && x [ 1 ] !== 0 ) ;
164+ const zerosAtYAxis = this . #zeros. filter (
165+ ( x ) =>
166+ ( x [ 0 ] === 0 && x [ 1 ] !== 0 ) ||
167+ ( Math . abs ( x [ 0 ] ) < 1 && Math . abs ( x [ 1 ] / x [ 0 ] ) > 15 )
168+ ) ;
169+ const polesAtYAxis = this . #poles. filter (
170+ ( x ) =>
171+ ( x [ 0 ] === 0 && x [ 1 ] !== 0 ) ||
172+ ( Math . abs ( x [ 0 ] ) < 1 && Math . abs ( x [ 1 ] / x [ 0 ] ) > 15 )
173+ ) ;
164174
165175 zerosAtYAxis . forEach ( ( x ) => {
166176 const value = expectedSteepPhaseShiftsMap . get ( Math . abs ( x [ 1 ] ) ) ;
@@ -186,7 +196,8 @@ export default class BodePlot {
186196 let newMagnitudeValue ;
187197 let newPhaseValue ;
188198 let lastMagnitudeValue = magnitude ( this . #wMin) ;
189- let lastPhaseValue = ( 180 / Math . PI ) * phase ( this . #wMin) ;
199+ const firstPhaseValue = ( 180 / Math . PI ) * phase ( this . #wMin) ;
200+ let lastPhaseValue = firstPhaseValue ;
190201
191202 const variableStep = new VariableStep ( ) ;
192203 const phaseUnwrapper = new PhaseUnwrapper ( expectedSteepPhaseShiftsMap ) ;
@@ -256,22 +267,49 @@ export default class BodePlot {
256267 ? Math . max ( ...magnitudeCurvePoints )
257268 : undefined ;
258269
259- //adjust phase based on expected phase value at wMax according to the transfer function polynomials
260- const zerosAtPositiveHalfplane = this . #zeros. filter ( ( x ) => x [ 0 ] > 0 ) . length ;
261- const zerosAtNegativeHalfplaneOrYAxis = this . #zeros. filter (
262- ( x ) => x [ 0 ] <= 0
263- ) . length ;
270+ //calculate pre-adjusted phase min & max values
271+ const minPhaseTemp =
272+ this . #phaseCurvePoints. length > 0
273+ ? Math . min ( ...this . #phaseCurvePoints. map ( ( x ) => x [ 1 ] ) )
274+ : undefined ;
275+ const maxPhaseTemp =
276+ this . #phaseCurvePoints. length > 0
277+ ? Math . max ( ...this . #phaseCurvePoints. map ( ( x ) => x [ 1 ] ) )
278+ : undefined ;
279+
280+ let expectedPhaseValueAtWmax = lastPhaseValue ;
281+ const tol = tolerancePhaseAdjustmentLarge ;
264282
265- const expectedPhaseValueAtWmax =
266- - 90 *
267- ( this . #denominatorTermsArray. length -
268- 1 +
269- zerosAtPositiveHalfplane -
270- zerosAtNegativeHalfplaneOrYAxis ) ;
283+ //adjust expected phase value at wMax - case 1
284+ if ( minPhaseTemp > 180 - tol && maxPhaseTemp > minPhaseTemp + 2 * tol ) {
285+ expectedPhaseValueAtWmax = lastPhaseValue - 360 ;
286+ } else if (
287+ maxPhaseTemp < - 180 + tol &&
288+ maxPhaseTemp > minPhaseTemp + 2 * tol
289+ ) {
290+ expectedPhaseValueAtWmax = lastPhaseValue + 360 ;
291+ }
292+
293+ //adjust expected phase value at wMax - case 2
294+ //(display absolute phase values beyond 180, as an exception in this case)
295+ const polesAtYAxisOrOrigin = this . #poles. filter ( ( x ) => x [ 0 ] === 0 ) . length ;
296+ const zerosAtYAxisOrOrigin = this . #zeros. filter ( ( x ) => x [ 0 ] === 0 ) . length ;
297+ if (
298+ this . #poles. length === polesAtYAxisOrOrigin &&
299+ this . #zeros. length === 0
300+ ) {
301+ expectedPhaseValueAtWmax = - 90 * polesAtYAxisOrOrigin ;
302+ } else if (
303+ this . #zeros. length === zerosAtYAxisOrOrigin &&
304+ this . #poles. length === 0
305+ ) {
306+ expectedPhaseValueAtWmax = 90 * zerosAtYAxisOrOrigin ;
307+ }
271308
272- if ( lastPhaseValue > expectedPhaseValueAtWmax + 10 ) {
309+ //apply adjustment
310+ if ( lastPhaseValue > expectedPhaseValueAtWmax + tol ) {
273311 const factor = Math . ceil (
274- ( Math . abs ( lastPhaseValue - expectedPhaseValueAtWmax ) - 10 ) / 180
312+ ( Math . abs ( lastPhaseValue - expectedPhaseValueAtWmax ) - tol ) / 180
275313 ) ;
276314 logMessages (
277315 [
@@ -287,9 +325,9 @@ export default class BodePlot {
287325 this . #phaseCurvePoints. forEach (
288326 ( _ , i ) => ( this . #phaseCurvePoints[ i ] [ 1 ] -= factor * 180 )
289327 ) ;
290- } else if ( lastPhaseValue + 10 < expectedPhaseValueAtWmax ) {
328+ } else if ( lastPhaseValue + tol < expectedPhaseValueAtWmax ) {
291329 const factor = Math . ceil (
292- ( Math . abs ( expectedPhaseValueAtWmax - lastPhaseValue ) - 10 ) / 180
330+ ( Math . abs ( expectedPhaseValueAtWmax - lastPhaseValue ) - tol ) / 180
293331 ) ;
294332 logMessages (
295333 [
0 commit comments