/************************************************************************/
/*                            TurbTorque.c                              */
/*                          Alan Swithenbank                            */
/*                            October 2002                              */
/************************************************************************/
/* This program calculates power output for a Tesla turbine based on an */
/* analysis by William Tahil in TEBA News Issue #16, 1999. The analysis */
/* determines torque considering the disc surface to represent a spiral */
/* vortex. The shear stress for an infinitesimal area is determined and */
/* integrated across the vortex. From that, torque is determined, and   */
/* power is determined as the product of torque and radial frequency.   */
/* The analysis considers the working fluid to be incompressible and    */
/* inviscid. A better approximation would be to solve the Navier-Stokes */
/* equations for the vortex. Also, there is no accounting for volute    */
/* interactions, rivets, starwashers, etc.                              */
/*                                                                      */
/* This is a semi-interactive routine. Enter a runner diameter and disc */
/* spacing (both in units of inches), and calculations are performed    */
/* and plotted for runners with a vent to disc diameter ratio the same  */
/* as Tesla's 9.75" turbine, with 6, 9, 12, and 15 disc runner stacks,  */
/* (including the inside surfaces of the end plates), over an rpm range */
/* from 2000 to 10000, with air as the working fluid.                   */
/*                                                                      */
/* The apparently pointless use of arrays in the code here and there is */
/* beause this code was produced from a full auto generation program    */
/* (TurbTorque-Auto.c, June 2002, which outputs data for several fixed  */
/* runner diameters and disc spacings), by simply changing the numdia   */
/* and numgap parameter values to 1, and adding the manual input        */
/* features to the main calculation loop.                               */
/*                                                                      */
/* Screen plot data are also output to the file RAM:TurbTorque.data.    */
/*                                                                      */
/* Calculation quantitites:                                             */
/* -----------------------                                              */
/*                                                                      */
/*         R2: disc outer radius                                        */
/*                                                                      */
/*         R1: disc vent radius                                         */
/*                                                                      */
/*        gap: disc separation                                          */
/*                                                                      */
/*          h: gap/2                                                    */
/*                                                                      */
/*        rev: disc revolutions/minute                                  */
/*                                                                      */
/*          w: disc radial velocity                                     */
/*                                                                      */
/*         mu: dynamic viscosity (1.79E-5 Ns/m^2 for air)               */
/*                                                                      */
/*     vtheta: tangential velocity at disc rim (and input air speed)    */
/*                                                                      */
/*      gamma: vortex circulation = V_theta(2Pi)(outer radius)          */
/*                                                                      */
/*          T: surface torque = 3mu(Gamma)(R2^2 - R2^2)/(2h)            */
/*                                                                      */
/*         TD: disc torque = 2T (two surfaces doubles torque)           */
/*                                                                      */
/*         PD: disc power = wTD                                         */
/*                                                                      */
/*         PT: turbine power = PD*(number of discs - 1)                 */
/*                                                                      */
/*             [ -1 disc because only insides of end discs are used]    */
/*                                                                      */
/* Conversion Factors:                                                  */
/* -------------------                                                  */
/*                                                                      */
/*           1 inch = 0.0254 m                                          */
/*                                                                      */
/*            1 rpm = 0.1047 rad/sec                                    */
/*                                                                      */
/*             1 Nm = 1.3558 lbf-ft                                     */
/*                                                                      */
/*    1 Hp-electric = 550.22 lbf-ft/s                                   */
/*                                                                      */
/*                                                                      */
/* NOTE: Requires Workbench 2.0 or greater.                             */
/*                                                                      */
/* NOTE: Window pixels count increasing right and increasing down.      */
/*                                                                      */
/* NOTE: No input error checking! Routine happily crashes on sillyness. */
/*                                                                      */
/* NOTE: While working on finalizing the drawplots() function, found a  */
/*       reason for the apparent need to use a GimmeZeroZero window to  */
/*       prevent the program from hanging up on the data loop when it   */
/*       is started from Workbench. This was due to the stack size set  */
/*       to 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 old FORTRAN programming   */
/*       experiences...DOH!...)                                         */
/*                                                                      */
/* NOTE: The program has been recompiled to use a standard window with  */
/*       proper stack size and runs from Workbench fine. But, updating  */
/*       the drawplots() function included a change in the PlotInfo     */
/*       structure and moving definitions 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: Forgot to include the factor of 2 for inner runner discs which */
/*       have two sides that contribute equally to the torque. Fixed on */
/*       09 October 2005 by adding variable T as defined in comments.   */
/*       Note above about how to recompile this code.                   */
/*                                                                      */
/* 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
{
  float         discdiam;             /* current disc diameter           */
  float         ventratio;            /* vent diam to disc diam ratio    */
  float         MaxHP;                /* maximum horsepower plot Y axis  */
  float        *rpm;                  /* rpm values for plox X axis      */
  float        *HP;                   /* HP values for plot Y axis       */
  float        *space;                /* gap space value array           */
  float        *discs;                /* number of discs value array     */
  float         Xorg;                 /* plot origin pixel X location    */
  float         Yorg;                 /* plot origin pixel X location    */
  float         Xsize;                /* plot X axis pixel length        */
  float         Ysize;                /* plot Y axis pixel length        */
  int          *Clrs;                 /* plot line colors                */
  int           si;                   /* start index for HP array        */
  int           stkcnt;               /* number of stack sizes           */
  int           gapcnt;               /* number of gap sizes             */
  int           diacnt;               /* number of diameter sizes        */
  int           rpmcnt;               /* number of rpm values            */
  int           dsphgt;               /* display window height           */
  int           dspwth;               /* display window width            */
  RastPort     *Rst;                  /* program window RastPort pointer */
  ULONG         CLI;                  /* program I/O CLI handle pointer  */
  Window       *Wnd;                  /* program window struct pointer   */
};

typedef struct PlotInfo PlotInfo;

enum color {gap1=0,gap2,gap3,axis,text};    /* abstractions for plotting */
 
/*=======================================================================*/
/*========================= Function Prototypes =========================*/
/*=======================================================================*/

/*=== System Specific:                                                   */

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

/*============== 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);
    }

  GadToolsBase = (Library *)OpenLibrary("gadtools.library",0);
  if (GadToolsBase == NULL)         /* GadTools library open failure     */
    {     
      CloseLibrary(IntuitionBase);
      CloseLibrary((Library *)GfxBase);
      exit(30);
    }

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

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

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

  CloseLibrary(GadToolsBase);
  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);
    }

  GadToolsBase = (Library *)OpenLibrary("gadtools.library",0);
  if (GadToolsBase == NULL)         /* GadTools library open failure     */
    {     
      CloseLibrary(IntuitionBase);
      CloseLibrary((Library *)GfxBase);
      exit(30);
    }

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

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

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

  CloseLibrary(GadToolsBase);
  CloseLibrary((Library *)GfxBase);
  CloseLibrary(IntuitionBase);
       
  exit(0);
}

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

BOOL my_main(BOOL WBStart) 
{
  PlotInfo      PlotData;             /* plot axes information struct    */

  BOOL          DONE = FALSE;         /* program bail out flag           */

  char          CLIStr[100];          /* CLI control string              */
  char          TmpStr[100];          /* temporary use string            */
  char          FilStr[32];           /* output data file string         */
  FILE         *discdata;             /* file handle for data output     */
  float         D2[]={6,9.75,12,
                      15,20};         /* disc diameter values            */
  float         D1,R2,R1;             /* vent diameter disc/vent radii   */
  float         gap[]={.02,.031,.04}; /* disc separation values          */
  float         h;                    /* (disc separation)/2             */
  float         ndiscs[]={6,9,12,15}; /* number of discs in runner stack */
  float         rev[]={2000,3000,
                       4000,5000,
                       6000,8000,
                       10000};        /* disc rpm                        */
  float         w;                    /* disc rad/sec                    */
  float         vtheta;               /* disc rim tangential velocity    */
  float         gamma;                /* vortex circulation              */
  float         mu;                   /* dynamic viscosity of fluid      */
  float         T;                    /* disc torque (from one side)     */
  float         TD;                   /* disc torque (from both sides)   */
  float         PD,PT[500];           /* disc power, turbine power       */
  float         inches2m = 0.0254;    /* inch to meters conversion       */
  float         rpmtorad = 0.1047;    /* rpm to rad/sec conversion       */
  float         nmtolbft = 1.3558;    /* Nm to lbf-ft conversion         */
  float         hptolfps = 550.22;    /* hp-elec to lbf-ft/s conversion  */
  float         ventfact = 0.487;     /* 9.75" turbine vent to rim ratio */
  float         maxHP;                /* maximum horsepower (plot param) */
  float         Xlen,Ylen;            /* X and Y axes lengths in pixels  */
  float         Xzero,Yzero;          /* pixel position for plot origin  */

  int           i;                    /* multiuse index                  */
  int           numdia = 1;           /* number of different diameters   */
  int           numgap = 1;           /* number of gap sizes             */
  int           numstk = 4;           /* number of disc stack variations */
  int           numrpm = 7;           /* number of used rpm values       */
  int           di,gi,ri,ni,pa;       /* disc paramater array indices    */
  int           PltClr[16];           /* pen colors for plot lines       */
  int           CLIHgt;               /* CLI window height               */
  int           CLITop;               /* CLI window top screen position  */
  int           CLIWth;               /* CLI window width                */
  int           MonWth;               /* monitor screen pixel width      */
  int           MonHgt;               /* monitor screen pixel height     */
  int           WndHgt;               /* program window height           */
  int           WndWth;               /* program window width            */

  RastPort     *PgmRst;               /* program window RastPort pointer */

  ULONG         CLIHnd;               /* program I/O CLI handle pointer  */

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

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

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

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

/*=== open a CLI handle for I/O: ===*/

  CLIHgt = 75;
  CLITop = 25;
  CLIWth = 600;

  strcpy(CLIStr,"CON:");                     /* build CLI control string */
  sprintf(TmpStr,"%i/",(MonWth-CLIWth)/2);   /* from height, top, width  */
  strcat(CLIStr,TmpStr);
  sprintf(TmpStr,"%i/",CLITop);
  strcat(CLIStr,TmpStr);
  sprintf(TmpStr,"%i/",CLIWth);
  strcat(CLIStr,TmpStr);
  sprintf(TmpStr,"%i/",CLIHgt);
  strcat(CLIStr,TmpStr);
  strcat(CLIStr,"TurbTorque CLI");

  CLIHnd = OpnHnd(CLIStr);

  if (CLIHnd==0) { 
    return FALSE; 
  }

/*=== set program window width and height (fit under CLI): ===*/

  WndWth = CLIWth;
  WndHgt = (MonHgt-CLITop-CLIHgt)-50;

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

//NOTE: see bottom of calculation loop for GimmeZeroZero use excuse.
//NOTE: see main comment notes for why no longer needed!

  PgmWnd = OpenWindowTags(0,
     WA_Title,         (ULONG)"TurbTorque Display Window",
     WA_Top,           15+CLIHgt+CLITop,
     WA_Left,          (720-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 X and Y axes lengths to 3/4 of program window sizes: */

  Xlen = 3*WndWth/4;
  Ylen = 3*WndHgt/4;

/* set plot origin location (Y increases down for pixels): */

  Xzero = WndWth/8;
  Yzero= 7*WndHgt/8;

/* set plot line colors: */

  PltClr[axis] =  1;
  PltClr[gap1] =  1;
  PltClr[gap2] =  5;
  PltClr[gap3] =  6;
  PltClr[text] =  2;

/* initialize PlotInfo structure fixed parameters: */

  (&PlotData)->ventratio = ventfact;
  (&PlotData)->rpm       = rev;
  (&PlotData)->HP        = PT;
  (&PlotData)->space     = gap;
  (&PlotData)->discs     = ndiscs;
  (&PlotData)->Xorg      = Xzero;
  (&PlotData)->Yorg      = Yzero;
  (&PlotData)->Xsize     = Xlen;
  (&PlotData)->Ysize     = Ylen;
  (&PlotData)->Clrs      = PltClr;
  (&PlotData)->stkcnt    = numstk;
  (&PlotData)->gapcnt    = numgap;
  (&PlotData)->diacnt    = numdia;
  (&PlotData)->rpmcnt    = numrpm;
  (&PlotData)->dspwth    = WndWth;
  (&PlotData)->dsphgt    = WndHgt;
  (&PlotData)->Rst       = PgmRst;
  (&PlotData)->CLI       = CLIHnd;
  (&PlotData)->Wnd       = PgmWnd;

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

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

  sprintf(TmpStr,"\nOutput data file: %s\n",FilStr);
  WrtCLI(TmpStr,CLIHnd,FALSE,FALSE);
  fprintf(discdata,TmpStr);

/*=== do calculations and plots: ===*/

  mu = .0000179;                          /* dynamic viscosity of air    */
   
  while(!DONE) {

    pa = 0;                               /* set power array start index */

    for (di=0;di<numdia;di++) {

/*** added manual data entry here: ***/

      for (i=0;i<strlen(TmpStr);i++) { TmpStr[i]=' '; }  /* clear string */      

      ReadCL(CLIHnd,TmpStr,"Enter disc diameter (inches): ",TRUE,FALSE);
      D2[di] = atof(TmpStr);

      for (i=0;i<strlen(TmpStr);i++) { TmpStr[i]=' '; }  /* clear string */      
      ReadCL(CLIHnd,TmpStr,"Enter disc gap (inches): ",TRUE,FALSE);
      gap[di] = atof(TmpStr);

/*** end manual data entry ***/

      sprintf(TmpStr,"\n%.3f inch diameter turbine data:",D2[di]);
      fprintf(discdata,TmpStr);
      sprintf(TmpStr,"\nWriting %.3f inch diameter turbine data.",D2[di]);
      WrtCLI(TmpStr,CLIHnd,FALSE,FALSE);
       
      maxHP = 0;                          /* zero maxHP for current disc */
      (&PlotData)->si = pa;               /* data segment start index    */
      (&PlotData)->discdiam = D2[di];     /* current disc diameter       */

      for (gi=0;gi<numgap;gi++) {
        sprintf(TmpStr,"\n  GAP = %.3f\",",gap[gi]);
        fprintf(discdata,TmpStr);

        for (ri=0;ri<7;ri++) {
          sprintf(TmpStr,"\n     RPM = %.3f:",rev[ri]);
          fprintf(discdata,TmpStr);

          for (ni=0;ni<4;ni++) {
            sprintf(TmpStr,"\n        discs = %.3f:",ndiscs[ni]);
            fprintf(discdata,TmpStr);

            D1 = ventfact*D2[di];               /* value in inches       */

            R2 = inches2m*D2[di]/2;             /* values in meters      */
            R1 = inches2m*D1/2;
            h  = inches2m*gap[gi]/2;      
                           
            w = rpmtorad*rev[ri];               /* rpm to rad/sec        */
 
            vtheta = (PI*2*R2)*rev[ri]/60;      /* (m/rev)*(rev/sec=m/s  */
            gamma = vtheta*2*PI*R2;

            T = 3*mu*gamma*(R2*R2 - R1*R1)/h;   /* single suface torque  */
            TD = 2*T;                           /* single disc torque    */
            PD = TD*w;                          /* single disc power     */

            PT[pa] = (ndiscs[ni]-1)*PD;         /* one side of endplates */
            PT[pa] = PT[pa]/nmtolbft/hptolfps;  /* to hp-electric        */

            if (maxHP < PT[pa]) {maxHP=PT[pa];} /* get maximum horsepwr  */

            sprintf(TmpStr," Total HP-electric = %.2f",PT[pa]);
            fprintf(discdata,TmpStr);

            pa++;                               /* update PT array index */

          }
        }
      }

/*=== report max HP for disc size: ===*/

      sprintf(TmpStr,"\n\n  Maximum HP-electric = %.2f\n",maxHP);
      WrtCLI(TmpStr,CLIHnd,FALSE,FALSE);
      fprintf(discdata,TmpStr);

/*=== draw current turbine plots: ===*/

      (&PlotData)->MaxHP = maxHP;       /* max HP for array segment    */

      drawplots(&PlotData);             /* plot data                   */

      pa++;                             /* update power array index    */

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

/* modifed to recycle at user's option rather than just quit */

      for (i=0;i<strlen(TmpStr);i++) { TmpStr[i]=' '; }  /* clear string */
      ReadCL(CLIHnd,TmpStr,"enter quit to quit, hit enter to continue: ",
             TRUE,FALSE);

      if ((strncmp(TmpStr,"quit",4))==0) { 
        DONE = TRUE; 
      } else {
        DONE = FALSE;

/*=== clear previous window image: ===*/

        SetRast(PgmRst,0);          
        RefreshWindowFrame(PgmWnd);

// NOTE: Here, RefreshWindowFrame(PgmWnd); causes hang up on loop 
//       when start with workbench icon outside of StormC environment,
//       So, switched window type to gimmezerozero so frame refresh
//       is not necessary. I don't even want to go there.

// NOTE: See main comment notes for why no longer needed!


      }

    }

  };

/*============================== BAIL OUT ===============================*/

  WrtCLI("\n\n\n\n\n                               BYE BYE!\n\n\n",
         CLIHnd,FALSE,FALSE);

  CloseWindow(PgmWnd);

  fclose(discdata);

  Delay(100);

  Close(CLIHnd);;

  return(DONE);

}

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

BOOL drawplots(PlotInfo *PlotData)
{ 

/* note: pixel count is positive down.                                   */

/* note: draw three lines, one for each gap size, on plot of rpm vs. HP. */
/*       The data contain rpm vs. HP information for the number of discs */
/*       in a turbine. So, do multiple gap related plots, one for each   */
/*       number of discs data set to see everything on the plots. A plot */
/*       of all HP vs. rpm for number of discs data is made for each     */
/*       disc diameter.                                                  */

  char     TmpStr[100];

  float    Xzero,Yzero,Xnext,Ynext;
  float    Xlen,Ylen;
  float   *PT;
  float   *rev;
  float   *gap;
  float    maxHP,HPplcs,HPaxis;
  float    maxRPM,RPMplcs,RPMaxis;
  float    HPfac,RPMfac;
  float    diam,vntfac;
  float   *ndiscs;

  int      pa,i,j,k,l,m,n;
  int      numdia,numgap,numrpm,numstk;
  int      wndwth,wndhgt;

  RastPort *PgmRst;

/*=== pull out fixed data for current power array segment: ===*/

  diam   = PlotData->discdiam;
  vntfac = PlotData->ventratio;
  rev    = PlotData->rpm;
  gap    = PlotData->space;
  ndiscs = PlotData->discs;
  maxHP  = PlotData->MaxHP;
  numdia = PlotData->diacnt;
  numgap = PlotData->gapcnt;
  numrpm = PlotData->rpmcnt;
  numstk = PlotData->stkcnt;
  pa     = PlotData->si;
  PT     = PlotData->HP;
  PgmRst = PlotData->Rst;
  Xzero  = PlotData->Xorg;
  Yzero  = PlotData->Yorg;
  Xlen   = PlotData->Xsize;
  Ylen   = PlotData->Ysize;
  wndwth = PlotData->dspwth;
  wndhgt = PlotData->dsphgt;
  
/*=== calculate axis pixel conversion factors: ===*/

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

  HPplcs = ceil(log10(maxHP));   /* round up pwr of 10 for HP axis value  */ 
  HPaxis = pow(10,HPplcs);       /* calculate rounded up max axis value   */
  if ((HPaxis/maxHP>4)) {        /* reduce max axis where reasonable      */        
    HPaxis=HPaxis/4;
  } else if(HPaxis/maxHP>3) {
    HPaxis=HPaxis/3;
  } else if(HPaxis/maxHP>2) {
    HPaxis=HPaxis/2;
  }
  HPfac  = Ylen/HPaxis;          /* HP value to pixels conversion factor  */

  maxRPM = rev[numrpm-1];          
  RPMplcs = ceil(log10(maxRPM)); /* round up pwr of 10 for RPM axis value */ 
  RPMaxis = pow(10,RPMplcs);     /* calculate rounded up max axis value   */
  if ((RPMaxis/maxRPM>4)) {      /* reduce max axis where reasonable      */        
    RPMaxis=RPMaxis/4;
  } else if(RPMaxis/maxRPM>3) {
    RPMaxis=RPMaxis/3;
  } else if(RPMaxis/maxRPM>2) {
    RPMaxis=RPMaxis/2;
  }
  RPMfac = Xlen/RPMaxis;         /* RPM value to pixels conversion factor */

/*=== draw axes: ===*/

  SetDrPt(PgmRst,0xFFFF);        /* make sure solid line axes and tics    */

  SetAPen(PgmRst,PlotData->Clrs[axis]);

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

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

/* X axis tics:                                                          */

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

/* Y axis tics:                                                          */

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

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

  SetAPen(PgmRst,PlotData->Clrs[gap1]);

  Move(PgmRst,95,50); 
  sprintf(TmpStr,"GAP = %.3f\"",gap[gap1]);
  Text(PgmRst,TmpStr,strlen(TmpStr));

/****  Since just one gap in this code, use only one plot color 

  SetAPen(PgmRst,PlotData->Clrs[gap2]);

  Move(PgmRst,95,62); 
  sprintf(TmpStr,"GAP = %.3f\"",gap[gap2]);
  Text(PgmRst,TmpStr,strlen(TmpStr));

  SetAPen(PgmRst,PlotData->Clrs[gap3]);

  Move(PgmRst,95,74); 
  sprintf(TmpStr,"GAP = %.3f\"",gap[gap3]);
  Text(PgmRst,TmpStr,strlen(TmpStr));

  SetAPen(PgmRst,PlotData->Clrs[text]);

****/
  
  i=150;j=75;
  for (m=0;m<numstk;m++) {           
    switch (m) {                     
      case  0:
        SetDrPt(PgmRst,0xFF00);
        break;
      case  1:                   
        SetDrPt(PgmRst,0xF020);
        break;         
      case  2:                   
        SetDrPt(PgmRst,0xC666);
        break;
      case  3:                   
        SetDrPt(PgmRst,0xFFFF);
        break;
      default:
        SetDrPt(PgmRst,0xFFFF);
        break;
    };
    Move(PgmRst,i,(j+=12));
    Draw(PgmRst,i+55,j);
    Move(PgmRst,i+60,j+2);
    sprintf(TmpStr,"%i discs",(int)ndiscs[m]);
    Text(PgmRst,TmpStr,strlen(TmpStr));
  }

  sprintf(TmpStr,"%.2f\" turbine, %.2f\" vent.",diam,vntfac*diam);
  Move(PgmRst,(wndwth-strlen(TmpStr))/2,25); 
  Text(PgmRst,TmpStr,strlen(TmpStr));

  Move(PgmRst,(int)Xzero-65,(int)(Yzero-Ylen/2.0+4)); 
  strcpy(TmpStr,"HP");
  Text(PgmRst,TmpStr,strlen(TmpStr));

  Move(PgmRst,(int)(Xzero+Xlen/2.0-12),(int)(Yzero+30)); 
  sprintf(TmpStr,"RPM");
  Text(PgmRst,TmpStr,strlen(TmpStr));

  Move(PgmRst,(int)(Xzero-4),(int)(Yzero+18)); 
  sprintf(TmpStr,"0");
  Text(PgmRst,TmpStr,strlen(TmpStr));

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

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

  Move(PgmRst,(int)(Xzero-18),(int)(Yzero+4)); 
  sprintf(TmpStr,"0");
  Text(PgmRst,TmpStr,strlen(TmpStr));

  sprintf(TmpStr,"%.1f",HPaxis/2.0);
  Move(PgmRst,(int)(Xzero-10-8*strlen(TmpStr)),(int)(Yzero+4-Ylen/2.0)); 
  Text(PgmRst,TmpStr,strlen(TmpStr));

  sprintf(TmpStr,"%.1f",HPaxis);
  Move(PgmRst,(int)(Xzero-10-8*strlen(TmpStr)),(int)(Yzero+4-Ylen)); 
  Text(PgmRst,TmpStr,strlen(TmpStr));


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

/* There are HP values in the power array for each of 4 disc stack sizes.*/
/* They are grouped sequentially for each rpm value within gap data sets.*/
/* So, HP values for each gap size must be read from the array mod 4 to  */
/* get the proper data. There are calculations for seven rpm values for  */
/* each gap value. And, there are three gap values. This means there are */
/* 28 HP values for each gap segment in the power array and 84 values in */
/* each disc diameter segment in the power array, arranged as follows:   */

/* NOTE: only using 1 gap value now, but code here is generic.           */               


/*         diameter 1
             gap 1
               rpm 1
                 stack size 1
                 stack size 2
                 stack size 3
                 stack size 4
               rpm 2
                 stack size 1
                 stack size 2
                 stack size 3
                 stack size 4
               rpm 3
                  .
                  .
                  .
               rpm 7
                 stack size 1
                 stack size 2
                 stack size 3
                 stack size 5
             gap 2 
                  .
                  .
                  .
             gap 3
                  .
                  .
                  .
           diameter 2
                  .
                  .
                  .
           diameter 3
                  .
                  .
                  .
           diameter 4
                  .
                  .
                  .               
                                                                         */

  for (m=0;m<numstk;m++) {           /* plot data for each stack size    */

    k = pa + m;                      /* adjust index to stack data start */

    switch (m) {                     /* change line style for stack size */
      case  0:
        SetDrPt(PgmRst,0xFF00);
        break;
      case  1:                   
        SetDrPt(PgmRst,0xF020);
        break;         
      case  2:                   
        SetDrPt(PgmRst,0xC666);
        break;
      case  3:                   
        SetDrPt(PgmRst,0xFFFF);
        break;
      default:
        SetDrPt(PgmRst,0xFFFF);
        break;
    };
                          
    for (j=0;j<numgap;j++) {              /* plot data for each gap      */

      SetAPen(PgmRst,PlotData->Clrs[j]);  /* change color for each gap   */

      Move(PgmRst,(int)Xzero,(int)Yzero); /* set pen position to origin  */
                                        
      n = k + j*numstk*numrpm;            /* set index to gap data start */

      l = 0;                              /* reset rpm index counter     */

      for (i=n;i<(n+numstk*numrpm);i+=numstk) {
    
        Xnext = Xzero + rev[l++]*RPMfac;  /* calculate next point pixels */
        Ynext = Yzero - PT[i]*HPfac;      /* (screen Y is positive down) */

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

      };

    }

  }

  return TRUE;
}