/************************************************************************/
/*                             PolePhase.c                              */
/*                          Alan Swithenbank                            */
/*                            October 2002                              */
/************************************************************************/
/* This program calculates and displays relative waveforms generated by */
/* pole magnets passing the first armature coils in an axial type three */
/* phase machine having 12 armature coils spaced at 30 degrees around   */
/* its armature frame and 16 pole magnets spaced at 22.5 degrees around */
/* its rotor face. (Oddly familiar, that...)                            */
/*                                                                      */
/* In general, 3-phase waveforms, represented by square waves, look as  */
/* depicted in the diagram below (Figure 1):                            */
/*                                                                      */
/*                                                                      */
/*                         electrical degrees                           */
/*                                                                      */
/*                0    60    120   180   240   300   360                */
/*                |     |     |     |     |     |     |                 */
/*                                                                      */
/*                -------------------                 -                 */
/*    Phase A:    |                 |                 |                 */
/*                |                 |                 |                 */
/*                -                 -------------------                 */
/*                                                                      */
/*                            -------------------                       */
/*    Phase B:                |                 |                       */
/*                            |                 |                       */
/*                -------------                 -------                 */
/*                                                                      */
/*                -------                 -------------                 */
/*    Phase C:          |                 |                             */
/*                      |                 |                             */
/*                      -------------------                             */
/*                                                                      */
/*                                                                      */
/*                 Figure 1: 3-phase square waves                       */
/*                                                                      */
/*                                                                      */
/* That is, 3 waveforms, A,B,C, with phase A leading phase B by 120     */
/* degrees, and phase B leading phase C by 120 degrees. Usually these   */
/* waveforms are represented by cosine functions with their phase angle */
/* differences set relative to phase A, where phase A has no relative   */
/* phase angle, and hence, phase A has maximum amplitude at the phi = 0 */
/* electrical degree position, i.e.,                                    */
/*                                                                      */
/*                                                                      */
/*                Phase A:     A = Ma*cos(phi)                          */
/*                Phase B:     B = Mb*cos(phi + b)                      */
/*                Phase C:     C = Mb*cos(phi + c)                      */
/*                                                                      */
/*                                                                      */
/* Where A,B,C are magnitudes calculated for the waveform at a phase    */
/* angle of phi electrical degrees, Ma,Mb,Mc are the maximum positive   */
/* magnitudes for each waveform, and b,c are phase angles of waveforms  */
/* B and C relative to waveform A, respectively.                        */
/*                                                                      */
/* To invoke the condition of no included phase A phase angle, note in  */
/* Figure 1 the points in the A and B phase waveforms where both are    */
/* heading positive with the same magnitude are at phi=360 (=0) degrees */
/* for phase A, with phase B leading 240 electrical degrees at phi=120  */
/* degrees. Similarly for phase C, this occurs at phi=360 (=0) degrees  */
/* for phase A with phase C leading 120 electrical degrees at phi=240.  */
/* The lead angles for phase are the ones required for b and c, so, the */
/* final three 3-phase cosine expressions are:                          */      
/*                                                                      */
/*                                                                      */
/*                Phase A:     A = Ma*cos(phi)                          */
/*                Phase B:     B = Mb*cos(phi + 240)                    */
/*                Phase C:     C = Mc*cos(phi + 120)                    */
/*                                                                      */
/*                                                                      */
/* These are the equations used in this program. However a further look */
/* at phi is required. The phrase "electrical degrees" used frequently  */
/* above is an important distinction for electical machines. 360 degree */
/* mechanical rotation of an electrical machine in general does not     */
/* equate to a single 360 phase rotation for the machine waveforms. For */
/* a motor/generator type machine 360 degrees phase change in waveforms */
/* occurs each time the armature coil to pole magnet alignment pattern  */
/* repeats itself with respect to N and S poles. (This is true whether  */
/* the machine is being used as a motor or as a generator.) Since there */
/* are multiple pole magnet to armature coil alignments with a single   */
/* 360 degree mechanical rotation of the machine there are multiple 360 */
/* electrical degree phase rotations of the waveforms or each machine   */
/* rotation. The number of 360 elecrical degree cycles per 360 degree   */
/* mechanical machine rotations depends on the number and arrangement   */
/* of pole magnets and armature coils in the machine.                   */
/*                                                                      */
/* For the machine under discussion the zero mechanical degree position */
/* is where pole magnet 0 (a North pole) aligns with the first phase A  */
/* armature coil. This is also the first electrical 0 position, and the */
/* position where the magnitude of the phase A waveform is a maximum    */
/* when the machine is in operation. The relative degrees displacement  */
/* of armature coils and pole magnets for the machine passing though    */
/* its mechanical zero degree position are depicted as a linear spread  */
/* in the diagram below (Figure 2):                                     */
/*                                                                      */
/*                                                                      */
/*                      mechanical                                      */
/*                          0                                           */
/*                                                                      */
/*                          |    30 deg    |    30 deg   |              */
/*   coil:              ... A              B             C ...          */
/*                          |              |             |              */
/*                                                                      */
/*               |          |          |          |          |          */
/*               S          N          S          N          S          */
/*               | 22.5 deg | 22.5 deg | 22.5 deg | 22.5 deg |          */
/*                                                                      */
/*   pole:   ... 15         0          1          2          3 ...      */
/*                                                                      */
/*                                     | n |  o   |   p  | q |          */
/*                                                                      */
/*        -----> rotation               n =  7.5 degrees                */
/*                                      o = 15.0 degrees                */
/*                                      p = 15.0 degrees                */
/*                                      q =  7.5 degress                */ 
/*                                                                      */
/*                                                                      */
/*   Figure 2: Relative degree offset positions of rotor pole magnets   */
/*             to the first three armature phase coils in a 16 pole     */
/*             magnet, 12 armature coil electrical machine passing      */
/*             through its mechanical zero degree position.             */
/*                                                                      */
/*                                                                      */
/* From the diagram in Figure 2, it can be seen that for each rotation  */
/* of 45 mechanical degrees rotation the N-S pole to coil alignment     */
/* is the same as the alignment in the mechanical 0 degrees position.   */
/* So, 45 degrees of mechanical rotation equals 360 electrical degrees  */
/* rotation for this machine. Thus, there are 360/45 = 8 electrical     */
/* phase rotations for each machine rotation. So, if machine rotation   */
/* degrees are used for angle phi, that angle must be multiplied by 8   */
/* to obtain the value for use in the phase waveform equations and thus */
/* generate the correct 3-phase waveforms for the electrical machine    */
/* under discussion, i.e.,                                              */
/*                                                                      */
/*                                                                      */
/*                Phase A:     A = Ma*cos(8*phi)                        */
/*                Phase B:     B = Mb*cos(8*phi + 240)                  */
/*                Phase C:     C = Mc*cos(8*phi + 120)                  */
/*                                                                      */
/*                                                                      */ 
/* The phase angle factor will, of course, change for an electrical     */
/* machine with a different pole magnet and armature coil arrangment    */
/* than used here.                                                      */
/*                                                                      */
/* Another thing to notice are that the phase waveform zero magnitude   */
/* position is when a N pole magnet centerline and a S pole magnet      */
/* centerline are equally distant from opposite sides of an armature    */
/* coil centerline; which, in the case of Figure 2 is 11.25 degrees.    */
/* Also, at the mechanical 0 position the closest displacement of any   */
/* nonaligned pole magnet's centerline from a phase coil centerline is  */
/* 7.5 degrees. So, in a program uses discrete degree steps to simulate */
/* rotation, the step size should be chosen as a factor of 7.5 degrees  */
/* to guarantee N and S pole magnet centerlines will align with all     */
/* armature coil centerlines at some point during the rotation process; */
/* guaranteeing full amplitude phase value calculations both positive   */
/* and negative for all phases. Here a step value of 1.5 mechanical     */
/* degrees was chosen.                                                  */
/*                                                                      */
/* Each coil in a phase coil set adds to the phase waveform magnitude,  */
/* but does not affect waveform phase offsets. That is, the relative    */
/* alignment of phases does not change with magnitude. So, in the code  */
/* below, the Ma,Mb,Mc amplitde factors were ignored, and calculations  */
/* made just for the [-1,+1] range of the cosine function.              */
/*                                                                      */
/*                                                                      */
/* NOTE: Requires Workbench 2.0 or greater.                             */
/*                                                                      */
/* NOTE: Window pixels count increasing right and increasing down.      */
/*                                                                      */
/* NOTE: While working on finalizing the drawplots() function, found a  */
/*       reason for the program hanging after opening the plot display  */
/*       window when it is started from Workbench. This was due to the  */
/*       stack size set too small in the Workbench icon information.    */
/*       Setting stack size from the default 4096 to 20000 solved this. */
/*       (Actually, stack issues should have come to mind from the old  */
/*       FORTRAN programming experiences...DOH!...)                     */
/*                                                                      */
/* NOTE: Updating the drawplots() function included a change in the     */
/*       PlotInfo structure and moving its definition and typedef to    */
/*       my_headers/GlobalStruct.h. So, there is a conflict with the    */
/*       local definition of PlotInfo structures here. To recompile     */
/*       this code the PlotInfo structure definition and its typedef    */
/*       need to be temporarily commented out of GlobalStruct.h.        */
/*                                                                      */
/* NOTE: C programming always sucks for some obscure reason or another. */
/*                                                                      */
/************************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

/*=======================================================================*/
/*================== Amiga System Includes And Typdefs ==================*/
/*=======================================================================*/

#include "thin_man_2:my_headers/SystemIncludes.h"

/*=======================================================================*/
/*============= Global Visibility Structures And Variables ==============*/
/*=======================================================================*/

#include "Thin_Man_2:my_headers/GlobalStruct.h"

struct PlotInfo
{
  BOOL          Xforce;               /* force X axis scale to minX maxX */
  BOOL          Yforce;               /* force Y axis scale to minY maxY */ 
  char         *xtitle;               /* pointer to X axis title string  */
  char         *ytitle;               /* pointer to Y axis title string  */
  float         datstp;               /* X axis step value               */
  float         maxY;                 /* maximum Y data value            */
  float         minY;                 /* minimum Y data value            */
  float         maxX;                 /* maximum X data value            */
  float         minX;                 /* minimum X data value            */
  int           numcol;               /* number of data values per row   */
  int           numrow;               /* number of data rows             */
  int           dsphgt;               /* display window height           */
  int           dspwth;               /* display window width            */
  int          *Clrs;                 /* plot line colors                */
  RastPort     *Rst;                  /* program window RastPort pointer */
  ULONG         CLI;                  /* program I/O CLI handle pointer  */
  void         *datary;               /* void pointer for Y data array   */
  Window       *Wnd;                  /* program window struct pointer   */
};

typedef struct PlotInfo PlotInfo;

enum pnum {axis=0,text,line0,line1,line2,line3};  /* pen abstractions   */
enum phse {A=0,B,C};                              /* phase abstractions */

/*=======================================================================*/
/*========================= Function Prototypes =========================*/
/*=======================================================================*/

/*=== System Specific:                                                   */

#include "Thin_Man_2:my_headers/AmigaIO.h"
 
BOOL drawplots(PlotInfo *);  /* plot drawing function                    */
BOOL my_main(BOOL);          /* main code function prototype             */

/*============== ANSI:                                                   */

#include "Thin_Man_2:my_headers/Extensions.h"

/*=======================================================================*/
/*=================== Miscellaneous Preprocessor Stuff ==================*/
/*=======================================================================*/

#define _debug FALSE

/*=======================================================================*/
/*====================== Amiga CLI Startup Routine ======================*/
/*=======================================================================*/

void main(int ac, char *av[])       /* main() is for CLI startup         */
{                                   /* known to StormC so no prototype   */

/*=== open libraries: ===*/

  IntuitionBase = (Library *)OpenLibrary("intuition.library",0);
  if (IntuitionBase == NULL)        /* Intuition library open failure    */
    { exit(10); }    
 
  GfxBase = (struct GfxBase *)OpenLibrary("graphics.library",0);
  if (GfxBase == NULL)              /* Graphics library open failure     */
    {     
      CloseLibrary(IntuitionBase);
      exit(20);
    }

/*=== call main routine: ===*/

  if (!my_main(FALSE))
    {
      DisplayBeep(0);               /* beep if main code function failed */
    } 
   
/*=== close libraries: ===*/

  CloseLibrary((Library *)GfxBase);
  CloseLibrary(IntuitionBase);

  exit(0);
}

/*=======================================================================*/
/*==================== Amiga Workbench Startup Routine ==================*/
/*=======================================================================*/

void wbmain(struct WBStartup *wbs)  /* wbmain() is for Workbench startup */
{                                   /* known to StormC so no prototype   */
    
/*=== open libraries: ===*/

  IntuitionBase = (Library *)OpenLibrary("intuition.library",0);
  if (IntuitionBase == NULL)        /* Intuition library open failure    */
    { exit(10); }    
 
  GfxBase = (struct GfxBase *)OpenLibrary("graphics.library",0);
  if (GfxBase == NULL)              /* Graphics library open failure     */
    {     
      CloseLibrary(IntuitionBase);
      exit(20);
    }

/*=== call main routine: ===*/

  if (!my_main(TRUE))
    {
      DisplayBeep(0);               /* beep if main code function failed */
    }

/*=== close libraries: ===*/

  CloseLibrary((Library *)GfxBase);
  CloseLibrary(IntuitionBase);

  exit(0);
}

/*=======================================================================*/
/*========================= Main Code Function ==========================*/
/*=======================================================================*/

BOOL my_main(BOOL WBStart) 
{
  char      TmpStr[100];              /* temporary use string            */
  char      XTitle[100];              /* string for X axis title         */
  char      YTitle[100];              /* string for Y axis title         */
  char      FilStr[32];               /* output data file name string    */

  FILE     *phasedata;                /* file handle for data output     */

  float     deg2rad = PI/180;         /* degrees to radians conversion   */

const float step = 1.5;               /* pole rotation step              */
const int   numstep=1+(int)(360/step);/* number of rotation steps        */
const int   numphse = 3;              /* number of phases                */
  float     phases[numstep][numphse]; /* phase data storage (pseudo 2D)  */
  float     tmpval;                   /* multiuse float                  */

  int       PltClr[16];               /* pen colors for plot lines       */
  int       numpole = 16;             /* number of pole magnets          */
  int       k;                        /* multiuse index                  */
  int       MonWth;                   /* monitor screen pixel width      */
  int       MonHgt;                   /* monitor screen pixel height     */
  int       WndHgt;                   /* program window height           */
  int       WndWth;                   /* program window width            */

  PlotInfo  PlotData;                 /* plot axes information struct    */

  RastPort *PgmRst;                   /* program window RastPort pointer */
  Window   *PgmWnd;                   /* program window struct pointer   */

  BOOL          closeMe;              /* IDCMP while loop variable       */
  ULONG         msgClass;             /* IDCMP flag (Amiga LONG)         */  
  IntuiMessage *WndoMsg;              /* pointer to IntuiMessage struct  */

/*=========================== INITIALIZATION ============================*/

/*=== set monitor screen width and height: ===*/

  MonWth = 720;                          /* use overscan width to center */
  MonHgt = 440;

  WndWth = MonWth*3/4;
  WndHgt = MonHgt*3/4;

/*=== open program display window: ===*/

  PgmWnd = OpenWindowTags(0,
    WA_Title,         (ULONG)"PolePhase Display Window: quit with gadget",
    WA_Top,           WndHgt/5,
    WA_Left,          (MonWth-WndWth)/2,
    WA_Width,         WndWth,
    WA_Height,        WndHgt,
    WA_GimmeZeroZero, (ULONG)TRUE,
    WA_CloseGadget,   (ULONG)TRUE,
    WA_IDCMP,         IDCMP_CLOSEWINDOW|IDCMP_REFRESHWINDOW,
    TAG_END );

  PgmRst = PgmWnd->RPort;

/*============================= MAIN CODE ===============================*/

/*=== initialize: ===*/

/* set plot line colors: */

  PltClr[axis]  =  1;   
  PltClr[text]  =  7;
  PltClr[line0] =  3;
  PltClr[line1] =  4;
  PltClr[line2] =  5;
  PltClr[line3] =  6;

/* open a data output file: */
  
  sprintf(FilStr,"RAM:polephase.data");

  if (!(phasedata = fopen(FilStr,"w"))) {
    return FALSE; 
  }

  sprintf(TmpStr,"\nOutput data file: %s\n\n",FilStr);
  fprintf(phasedata,TmpStr);

  sprintf(TmpStr,"pole rotation step: %.3f degrees\n\n",step);
  fprintf(phasedata,TmpStr);

  sprintf(TmpStr,"            deg-m    deg-e   Phase A  Phase B  Phase C\n\n");
  fprintf(phasedata,TmpStr);

/*=== do calculations: ===*/

  tmpval = 0;

  for (k=0;k<numstep;k++) {

    sprintf(TmpStr,"step %3i: ",k);
    fprintf(phasedata,TmpStr);
    
    phases[k][A] = cos(deg2rad*8*(float)k*step);
    phases[k][B] = cos(deg2rad*(8*(float)k*step + 240)); 
    phases[k][C] = cos(deg2rad*(8*(float)k*step + 120)); 

    sprintf(TmpStr,"%7.3f%9.3f%9.3f%9.3f%9.3f\n",
            (float)k*step,tmpval,phases[k][A],phases[k][B],phases[k][C]);   
    fprintf(phasedata,TmpStr);

    tmpval+=(8*step);
    if (tmpval>360) { tmpval-=360; }
  }

/*=== save data ===*/

  fclose(phasedata);

/*=== plot all phase data: ===*/

  strcpy(XTitle,"Mechanical Degrees");
  PlotData.xtitle  =   XTitle;        /* pointer to X axis title string  */
  strcpy(YTitle,"Phase Magnitude");
  PlotData.ytitle  =   YTitle;        /* pointer to Y axis title string  */
  PlotData.datary  =   &phases;       /* Y data array pointer            */
  PlotData.datstp  =   step;          /* X axix step value               */
  PlotData.maxY    =   1.0;           /* maximum Y data value            */
  PlotData.minY    =  -1.0;           /* minimum Y data value            */
  PlotData.maxX    = 360.0;           /* maximum X data value            */
  PlotData.minX    =   0.0;           /* minimum X data value            */
  PlotData.numcol  =   numphse;       /* number of data values per row   */
  PlotData.numrow  =   numstep;       /* number of data rows             */
  PlotData.dsphgt  =   WndHgt;        /* display window height           */
  PlotData.dspwth  =   WndWth;        /* display window width            */
  PlotData.Rst     =   PgmRst;        /* program window RastPort pointer */
  PlotData.Wnd     =   PgmWnd;        /* program window struct pointer   */
  PlotData.Clrs    =   PltClr;        /* plot pen color numbers          */
  PlotData.Xforce  =   TRUE;          /* force X axis to be minX maxX    */
  PlotData.Yforce  =   TRUE;          /* force Y axis to be minY maxY    */

  drawplots(&PlotData);

/*=== wait for user to examine plot: ===*/

  closeMe = FALSE;

  while(!closeMe) { 
    Wait( 1 << PgmWnd->UserPort->mp_SigBit );
    WndoMsg = (IntuiMessage *)GetMsg(PgmWnd->UserPort);
    if (WndoMsg) {
      msgClass = WndoMsg->Class; 
      ReplyMsg((Message *)WndoMsg); 
      if (msgClass==CLOSEWINDOW)
        { closeMe = TRUE; }
    }
  }                    
                 
/*============================== BAIL OUT ===============================*/

  CloseWindow(PgmWnd);

  return TRUE;

}

/*=======================================================================*/
/*===================== Local Function Definitions ======================*/
/*=======================================================================*/

BOOL drawplots(PlotInfo *PlotData)
{ 

/* This is a more generic version of drawplots() than from TurbTorque.c. */
/* It uses a data array from a PlotInfo structure and plots information  */
/* based on axes and scale data also taken from the PlotInfo structure.  */    
/* The data array pointed to by the PlotInfo structure datary parameter  */
/* datary is assumed to be two dimensional, with the assignment made to  */
/* the address of the array pointer (e.g., a pointer to Ydata[193][3]    */
/* would be placed in the PlotInfo structure PlotData via the statement: */
/* PlotData.datary = &Ydata, with datary given the type pointer to void  */
/* in the PlotInfo structure definition). The numcol and numrow values   */
/* define where to parse out proper row and column data from the data    */
/* array. It is assumed the data are stored so each row represents one   */
/* time step across all data, i.e., each column gives one parameter's    */
/* values in the array. So, e.g., 3 parameters, A,B,C, have values saved */
/* in a two dim array[numrow][numcol] as:                                */
/*                                                                       */
/*                                                                       */
/*       A0  B0 C0  <-->  array[0][0] array[0][1] array[0][2]            */
/*       A1  B1 C1  <-->  array[1][0] array[1][1] array[1][2]            */
/*       A2  B2 C2  <-->  array[2][0] array[2][1] array[2][2]            */
/*       A3  B3 C3  <-->  array[3][0] array[3][1] array[3][2]            */ 
/*       A4  B4 C4  <-->  array[4][0] array[4][1] array[4][2]            */
/*                                                                       */
/*                               .                                       */
/*                               .                                       */
/*                               .                                       */
/*                                                                       */
/*                                                                       */
/* So, here, numcol is 3 for the three different parameters, and numrow  */
/* gives the number of data values in each data column. (This probably   */
/* will make more sense to FORTRAN programmers than to C programmers.)   */
/* Thus, sequential data values for a particular parameter are accessed  */
/* in the data array by indexing from from the parameter's first data    */
/* position by the numcol value. The first values for each parameter are */
/* respectively, the first three values in the data array, i.e., e.g.,   */
/*                                                                       */
/*                                                                       */
/*                parameter A0 is array[0 + 0*numcol],                   */
/*                parameter A1 is array[0 + 1*numcol],                   */
/*                parameter A2 is array[0 + 2*numcol], ...               */
/*                                                                       */
/*                parameter B0 is array[1 + 0*numcol],                   */
/*                parameter B1 is array[1 + 1*numcol],                   */
/*                parameter B2 is array[1 + 2*numcol], ...               */
/*                                                                       */
/*                parameter C0 is array[2 + 0*numcol],                   */
/*                parameter C1 is array[2 + 1*numcol],                   */
/*                parameter C2 is array[2 + 2*numcol], ...               */
/*                                                                       */
/*                                                                       */
/* That is, to access an element from a given column, index it by its    */
/* column value plus its row value multiplied by the number of columns,  */
/*                                                                       */
/* For this access to work the data must be saved by the calling program */
/* in a 2D array that is dimensioned to be exactly the numrow by numcol  */
/* values used here. This can be done by using const int parameters to   */
/* declare the array, e.g., as in PolePhase.c here:                      */
/*                                                                       */
/*                                                                       */
/*         const float  step = 1.875;                                    */
/*         const int    numstep=1+(int)(360/step);   (the row count)     */
/*         const int    numphse = 3;                 (the column count)  */
/*               float  phases[numstep][numphse];                        */
/*                                                                       */
/*                                                                       */
/* Without the const modifiers, the array declaration will fail.         */
/*                                                                       */
/*                                                                       */
/* NOTE: there is a half-assed attempt at compensation for min values    */
/*       less than zero in the plot factor calculations, but assumes     */
/*       maximum values are greater than zero.                           */
/*                                                                       */
/* NOTE: currently this function uses only implict (defined step) X axes */
/*       values. Could add an array for explict (declared value) X axes. */
/*                                                                       */
/* NOTE: screen pixel Y values increase down.                            */
/*                                                                       */
/*                                                                       */
/*                                                                       */

/*=== declare and set parameters: ===*/

  char    TmpStr[100];         /* temporary use text string              */

  float  *ydata;               /* pointer to two dimensional array data  */
  float   maxY,minY,maxX,minX; /* max/min X and Y values from data set   */
  float   diffX,diffY;         /* used in axes pixel factor calculations */
  float   Xlen,Ylen;           /* x and y axes length parameters         */
  float   Xzero,Yzero;         /* plot axes (0,0) pixel position         */
  float   Xnext,Ynext;         /* plotting position parameters           */
  float   Xshift,Yshift;       /* use for shifted axes when minimums < 0 */
  float   Xplcs;               /* decimal places for X axis labeling     */
  float   Xaxis;               /* rouned up maximum X axis value         */
  float   Xfac;                /* X value to pixels conversion factor    */
  float   Yplcs;               /* decimal places for Y axis labeling     */
  float   Yaxis;               /* rouned up maximum Y axis value         */
  float   Yfac;                /* Y value to pixels conversion factor    */
  float   step;                /* X axis step value                      */
  float   tmpX,tmpY1,tmpY2;    /* multiuse X and Y parameters            */

  int     l,m;                 /* multiuse indices                       */
  int     numcol;              /* number of data array cols (X extent)   */
  int     numrow;              /* number of data array rows (Y extent)   */
  int     WndWth,WndHgt;       /* plotting window pixel width and height */
  int    *PltClr;              /* pointer to plot pen color value array  */

  RastPort *PgmRst;            /* plotting window rastport pointer       */

  ydata  = PlotData->datary;
  step   = PlotData->datstp;
  numrow = PlotData->numrow;
  numcol = PlotData->numcol;
  WndWth = PlotData->dspwth;
  WndHgt = PlotData->dsphgt;
  maxX   = PlotData->maxX;
  minX   = PlotData->minX;
  maxY   = PlotData->maxY;
  minY   = PlotData->minY;
  PgmRst = PlotData->Rst;
  PltClr = PlotData->Clrs;

  Xlen = 3*(float)WndWth/4;       /* set x and y axes to 3/4 window size */
  Ylen = 3*(float)WndHgt/4;
  
/*=== calculate axes pixel conversion factors: ===*/

/* set pixel conversion factors so can label axes rationally:            */

  if (minY<0) {                   /* get maximum Y axis extent           */
    diffY = maxY - minY;           
  } else {
    diffY = maxY;
  }
  if (!(PlotData->Yforce)) {      /* generate reasonable Y scale values  */
    Yplcs = ceil(log10(diffY));   /* round up pwr of 10 for Y axis value */
    Yaxis = pow(10,Yplcs);        /* calculate rounded up max axis value */
    if ((Yaxis/diffY>4)) {        /* reduce max axis where reasonable    */
      Yaxis=Yaxis/4;              /* (rounded up power could be way big) */
    } else if(Yaxis/diffY>3) {
      Yaxis=Yaxis/3;
    } else if(Yaxis/diffY>2) {
      Yaxis=Yaxis/2;
    }
  } else {                        /* force Y axis scale to be minY maxY  */
    Yaxis = diffY;
  }
  Yfac  = Ylen/Yaxis;             /* Y value to pixels conversion factor */
         
  if (minX<0) {                   /* get maximum Y axis extent           */
    diffX= maxX - minX;           
  } else {
    diffX = maxX;
  }
  if (!(PlotData->Xforce)) {      /* generate reasonalbe X scale values  */
    Xplcs = ceil(log10(diffX));   /* round up pwr of 10 for X axis value */
    Xaxis = pow(10,Xplcs);        /* calculate rounded up max axis value */
    if ((Xaxis/diffX>4)) {        /* reduce max axis where reasonable    */
      Xaxis=Xaxis/4;              /* (rounded up power could be way big) */
    } else if(Xaxis/diffX>3) {
      Xaxis=Xaxis/3;
    } else if(Xaxis/diffX>2) {
      Xaxis=Xaxis/2;
    }
  } else {                        /* force X axis scale to be minX minY  */  
    Xaxis = diffX;
  }
  Xfac = Xlen/Xaxis;              /* X value to pixels conversion factor */

/*=== calculate axes (0,0) pixel location: ===*/

/* note: X/Y axes lengths are 3/4 of corresponding window sizes          */
/* note: Y pixels incease down                                           */

  if (minX<0) {
    Xzero = (float)WndWth/8 + fabs(minX/(maxX-minX))*Xlen;
  } else {
    Xzero = (float)WndWth/8;  
  }

  if (minY<0) {
    Yzero = 7*(float)WndHgt/8 - fabs(minY/(maxY-minY))*Ylen;
  } else {
    Yzero= 7*(float)WndHgt/8;
  }

/*=== draw axes: ===*/
  
/* if minimums less than zero, adjust axis endpoint shift accordingly:   */

  if (minX>=0) {
    Xshift = Xzero;
  } else {
    Xshift = Xzero-fabs(minX/(maxX-minX))*Xlen;
  }

  if (minY>=0) {
    Yshift = Yzero;
  } else {
    Yshift = Yzero+fabs(minY/(maxY-minY))*Ylen;
  }
  
/* set for solid axes and tic lines:                                     */

  SetDrPt(PgmRst,0xFFFF);          
  SetAPen(PgmRst,PltClr[axis]);

/* draw X axes:                                                          */                                    

  Move(PgmRst,(int)Xshift,(int)Yzero);              
  Draw(PgmRst,(int)(Xshift+Xlen),(int)Yzero);

/* draw Y axes:                                                          */                                    

  Move(PgmRst,(int)Xzero,(int)Yshift);                   
  Draw(PgmRst,(int)Xzero,(int)Yshift-(int)Ylen);
  
/*=== put tic marks on axes ===*/

/* X axis tics:                                                          */

  Move(PgmRst,(int)Xshift,(int)Yzero); /* screen Y is pos down */ 
  Draw(PgmRst,(int)Xshift,(int)Yzero+8);        
  Move(PgmRst,(int)(Xshift+Xlen/2.0),(int)Yzero);
  Draw(PgmRst,(int)(Xshift+Xlen/2.0),(int)Yzero+8);
  Move(PgmRst,(int)(Xshift+Xlen),(int)Yzero);
  Draw(PgmRst,(int)(Xshift+Xlen),(int)Yzero+8);

/* Y axis tics:                                                          */

  Move(PgmRst,(int)Xzero,(int)Yshift);
  Draw(PgmRst,(int)Xzero-8,(int)Yshift);
  Move(PgmRst,(int)Xzero,(int)(Yshift-Ylen/2.0));
  Draw(PgmRst,(int)Xzero-8,(int)(Yshift-Ylen/2.0));
  Move(PgmRst,(int)Xzero,(int)(Yshift-Ylen));
  Draw(PgmRst,(int)Xzero-8,(int)(Yshift-Ylen));

/*=== label plot: ===*/

/* note: all labeling offsets assume topaz 8 font                        */

/* note: tic labels relative to calculated axes length, not min/max data */

  SetAPen(PgmRst,PltClr[text]);                   /* set text pen color  */

/* X axis tic labels:                                                    */

  if (Xshift==Xzero) {                       /* minX >= 0           */

    tmpX = minX + Xlen/Xfac;                      /* max X on axes value */

    if (minY>=0) {                                /* don't put on Y axis */
      sprintf(TmpStr,"%.1f",minX);
      Move(PgmRst,(int)(Xshift-4*strlen(TmpStr)),(int)(Yzero+18));
      Text(PgmRst,TmpStr,strlen(TmpStr));
    }

    sprintf(TmpStr,"%.1f",minX+(tmpX-minX)/2);
    Move(PgmRst,(int)(Xshift+Xlen/2.0-4*strlen(TmpStr)),(int)(Yzero+18));
    Text(PgmRst,TmpStr,strlen(TmpStr));

    sprintf(TmpStr,"%.1f",tmpX);
    Move(PgmRst,(int)(Xshift+Xlen-4*strlen(TmpStr)),(int)(Yzero+18));
    Text(PgmRst,TmpStr,strlen(TmpStr));


  } else {                                        /* minX < 0            */


// don't have to worry about this case here...for now...


  }


/* Y axis tic labels:                                                    */

  if (Yshift==Yzero) {                            /* minY >= 0           */


// don't have to worry about this case here...for now...


  } else {                                        /* minY < 0            */

    tmpY1 = (Yzero-Yshift)/Yfac;       /* min Y axis value (is negative) */
    tmpY2 = Ylen/Yfac+tmpY1;           /* max Y axis value (is positive) */

    sprintf(TmpStr,"%.2f",tmpY1);
    Move(PgmRst,(int)(Xzero-10-8*strlen(TmpStr)),(int)(Yshift+4));
    Text(PgmRst,TmpStr,strlen(TmpStr));

    sprintf(TmpStr,"%.2f",tmpY1+(tmpY2-tmpY1)/2);
    Move(PgmRst,(int)(Xzero-10-8*strlen(TmpStr)),(int)(Yshift+4-Ylen/2.0));
    Text(PgmRst,TmpStr,strlen(TmpStr));

    sprintf(TmpStr,"%.2f",tmpY2);
    Move(PgmRst,(int)(Xzero-10-8*strlen(TmpStr)),(int)(Yshift+4-Ylen));
    Text(PgmRst,TmpStr,strlen(TmpStr));

  }

/* X axis label:                                                         */

  Move(PgmRst,(int)Xshift+((int)Xlen-strlen(TmpStr)*8)/2,WndHgt-25);
  sprintf(TmpStr,PlotData->xtitle);
  Text(PgmRst,TmpStr,strlen(TmpStr));

/* Y axis label:                                                         */

  Move(PgmRst,35,25);
  sprintf(TmpStr,PlotData->ytitle);
  Text(PgmRst,TmpStr,strlen(TmpStr));

/*=== plot data: ===*/

  for (m=0;m<numcol;m++) {               /* plot data for each parameter */
                                                                                
    SetAPen(PgmRst,PltClr[line0+m]);     /* change each plot line color  */    

    Xnext = Xzero+minX*Xfac;             /* set position to first point  */
    Ynext = Yzero-ydata[m]*Yfac;         /* screen Y positive down       */
   
    Move(PgmRst,(int)Xnext,(int)Ynext); 

    for (l=1;l<numrow;l++) {             /* step thru data from point 2  */         
    
      Xnext = Xnext+step*Xfac;
      Ynext = Yzero - ydata[m+l*numcol]*Yfac;  /* screen Y positive down */

      Draw(PgmRst,(int)Xnext,(int)Ynext);

    }
  }

  return TRUE;

}
