#include 

/* plot a continuous frequency response in two Windows */
freqsplot(h, w, winmag, winphase)
{
        local m, p, wmag, wphase, wmnum, wpnum, win1, win2;

        /* parse args */
        if (argc < 4)
        {
                if (argc < 3)
                {
                        if (argc < 2)
                        {
                                if (argc < 1) error("freqsplot - complex frequency response required");
                                
                                w = xvals(h);
                        }
                        
                        winmag = refwin(W0);
                }

                winphase = {};

                if (iswindow(winmag))
                {
                        wpnum = getwnum(winmag) + 1;

                        if (wpnum > 0 && wpnum <= numwindows())
                        {
                                winphase = castwindow(getwnum(winmag) + 1);
                        }
                }
        }

        if (iswindow(winmag))
        {
                /* window numbers for casting */
                wmnum = getwnum(winmag);

                /* window for mag response */
                wmag = castwindow(wmnum);

                /* magnitude */
                m = mag(h);

                /* clear old formula */
                protect(wmag, "Magnitude Response");
                clear(wmag);

                /* plot as lines */
                setplotstyle(wmag, 0);

                /* log-log axes */
                setylog(wmag, 1, 0, 1);
                setxlog(wmag, 1, 0, 1);

                /* label axes */
                scales(wmag, 2);

                /* grids */
                setgridstyle(wmag, 1, 3);
                setgridstyle(wmag, 2, 3);

                /* assign to window */
                win1 = castwindow(wmnum);
                win1 = xy(w, m);

                /* xy labels */
                label(wmag, "Magnitude Response");
                setylabel(wmag, "Magnitude");
                setxlabel(wmag, sprintf("Frequency (%s)", getvunits(w)));
        }

        if (iswindow(winphase))
        {
                /* window numbers for casting */
                wpnum = getwnum(winphase);

                /* window for phase response */
                wphase = castwindow(wpnum);

                /* phase in degrees */
                p = unwrap(phase(h)) * 180 / pi;
                setvunits(p, "Degrees");

                /* clear old formula */
                protect(wphase, "Phase Response");
                clear(wphase);

                /* plot as lines */
                setplotstyle(wphase, 0);

                /* linear-log axes */
                setylog(wphase, 0);
                setxlog(wphase, 1, 0, 1);

                /* label axes */
                scales(wphase, 2);

                /* grids */
                setgridstyle(wphase, 1, 3);
                setgridstyle(wphase, 2, 3);

                /* assign to window */ 
                win2 = castwindow(wpnum);
                win2 =  xy(w, p);

                /* xy labels */
                label(wphase, "Phase Response");
                setylabel(wphase, "Phase (degrees)");
                setxlabel(wphase, sprintf("Frequency (%s)", getvunits(w)));
        }
}