While working on Exercise 6.5 of Ch06 in Dr. Middlebrook's D-OA method, I tried to make bode plot of the transfer function:
bodeplot[s/100+100/s*(1+10/s)] (input to wolframalpha)
in J
Somehow the J code phase plot doesn't agree with Mathematica's result, though the magnitude plot matches fine.
Anything wrong with my J code?
Af =: 4 : 0"_ 0
s=.0j1*y
'w q'=.x
f=.(s%w) + (w%s)*(1+w%q*s)
20*10^. | f
)
Pf =: 4 : 0"_ 0
s=.0j1*y
'w q'=.x
f=.(s%w) + (w%s)*(1+w%q*s)
(180%o.1)* 1{ *. f
)
load 'plot'
plot (; (100 10 Af (10 ^ ]))) 0.02*i.200
plot (; (100 10 Pf (10 ^ ]))) 0.02*i.200
To be more general, say a complex variable on the unit circle in the complex plane z = cos x + I sin x
If we plot its phase angle, there will be a jump at 180 degree (from 180 to -180)
z_unit_circle =. ((2 o. ]) + (0j1 * (1 o.]))) @ (180 %~ o.)
plot (180%o.1)*1{"1 *. z_unit_circle i.360
I think that's what happens when phase angle goes around 180 or -180 in the earlier J bode plot.
To avoid this jump, we can make use of the relationship Tan(Im(z)/Re(z)) = Tan(-180 + Im(z)/Re(z)), i.e. to turn -180 before hand.
phase_angle =. _180 + (180 % o.1) * (_3 o. %~/) @ +.
Pf =: 4 : 0"_ 0
s=.0j1*y
'w q'=.x
f=.(s%w) + (w%s)*(1+w%q*s)
phase_angle f
)
plot (; (100 10 Pf (10 ^ ]))) 0.02*i.200
This is essentially the same as the answer provided by Eelvex.
However this phase_angle[z] has more jumps than Arg[z]
plot phase_angle"1 z_unit_circle i.360
So my question is how to make the correct bode plot in J. In other words, knowing the phase angle goes from 3rd quadrant into 2nd quadrant, thus -180 before hand
