/********************************************************************
*                                                                   *
* Joukowski Airfoil Generator with Streamline and Pressure          *
* Distribution Algorithms                                           *
*                                                                   *
* Written by:   Russell Leighton      15 March 1987                 *
*               Lancaster, CA                                       *
* Modified by:  David Foster          19 June 1988                  *
*               Rochester, MI                                       *
*               To include effect of induced circulation            *
*               by a first order approximation.                     *
********************************************************************/

#include "airfoil.h"

main()
{
   float rs,theta,h,vel,eta,S;
   int psi;
   ULONG MessageClass;

   open_things();
   do_about();

   /****************/
   /* Loop forever */
   /****************/

   for(;;)
   {
      while (Continue)
      {
         /****************************************************/
         /* Wait, initially and after each plot, for user to */
         /* bring up the double menu requester               */
         /****************************************************/

         Wait(1L<<w->UserPort->mp_SigBit);
         if((message = (struct IntuiMessage *)GetMsg(w->UserPort)) != NULL)
         {
            MessageClass = message->Class;
            ReplyMsg(message);
            if (MessageClass == REQVERIFY)
            {
               do_request();
               break;
            }
         } /* end if */
      } /* end while */

      do_init();

      if(mode)
      {

         /********************/
         /* Plot Streamlines */
         /********************/

         FILL = TRUE;
         for (psi = 12;psi > 0;--psi)
         {
            do_mess();
            if(!Continue) break;

            PENUP;
            SetAPen(rp,(long)(psi+1));
            SetBPen(rp,(long)(psi+1));

            for (theta = 0.015;theta < TWO_PI;theta += PI/100)
            {
               vel = psi/(4.*velocity*r*sin(theta));
			   S = (1.+sin(alpha)/sin(theta));
			   vel = abs(S/(2.*vel));
               eta = sqrt(vel*vel+1.0)-vel;
               rs = r*(1.+eta)/(1.-eta);
               transform(rs,theta);
            } /* end for */
            AreaEnd(rp);
         } /* end for */

         /*******************************/
         /* Plot Stagnation Streamlines */
         /*******************************/

         do_mess();
         if(Continue)
         {
            FILL = FALSE;
            PENUP;
            SetAPen(rp,1L);
            h = (r-4.0)/40.0;
            theta = -alpha;

            for (rs = 4;rs >= r;rs += h)
            {
               transform(rs,theta);
            }

            PENUP;
            theta = PI + alpha;

            for (rs = 4;rs > r;rs += h)
            {
               transform(rs,theta);
            }
         } /* end if */
      } /* end if */

      else
      {

         /******************************/
         /* Plot Pressure Distribution */
         /******************************/

         do_mess();
         if(Continue)
         {
            FILL = TRUE;
            PENUP;
            SetAPen(rp,2L);
            SetBPen(rp,2L);

            for (theta = 0.0;theta <= TWO_PI;theta += PI/100)
            {
               rs = r+(sin(theta)+sin(alpha))*(sin(theta)+sin(alpha));
               transform(rs,theta);
            } /* end for */
            AreaEnd(rp);
         } /* end if */
      } /* end else */

      /****************/
      /* Plot Airfoil */
      /****************/

      do_mess();
      if(Continue)
      {
         FILL = TRUE;
         PENUP;
         rs = r;
         SetAPen(rp,1L);
         SetBPen(rp,1L);

         for (theta = 0.0;theta <= TWO_PI;theta += PI/100)
         {
            transform(rs,theta);
         }
         AreaEnd(rp);
      } /* end if */
   } /* end forever */
} /* end main */

do_init()
{
   float a0;

   SetAPen(rp,0L);
   RectFill(rp,(long)(XMIN+1),(long)(YMIN+1),(long)(XMAX-1),(long)(YMAX-1));
   SetOPen(rp,1L);

   /***********************************************************/
   /* Calculate circle constants (circle center and radius)   */
   /* from airfoil constants through a reverse transformation */
   /***********************************************************/

   alpha = (float)angle*PI/180;
   c = (float)camber/25;
   t = (float)thickness/25;
   a = 0;
   b = c/2;

   do
   {
      a0 = a;
      a = t*(2*a0+1)/4/sqrt(b*b+2*a0+1);
      b = c*(1+2*a)/(2+2*a);
   }
   while(abs(a-a0) > TOL);
   
   r = sqrt(b*b+(a+1)*(a+1));
   Continue = TRUE;
}

do_mess()
{
   ULONG MessageClass;

   /******************************************************************/
   /* Check for double menu requester. This can be brought up at any */
   /* time but can not be displayed unless the message is replied    */
   /* to, therefore this routine must be called periodically.        */
   /******************************************************************/

   if((message = (struct IntuiMessage *)GetMsg(w->UserPort)) != NULL)
   {
      MessageClass = message->Class;
      ReplyMsg(message);
      if(MessageClass == REQVERIFY) do_request();
   }
   ixo = XCEN;
   iyo = YCEN;
}

transform(rs,theta)

float rs,theta;
{
   float x,y,z,u,v;
   int PLOT;
   long ix,iy,cx,cy;

   /********************************************************************/
   /* This is the Joukowski transformation routine. This is also where */
   /* the plotting is done. Plotting is usually done using the area    */
   /* fill routines (AreaDraw & AreaMove). If FILL is true then the    */
   /* points are used to build up the area shape to be filled. See the */
   /* Rom Kernal manual containing the graphics primatives for more    */
   /* details. If FILL is false then normal plotting is done.          */
   /********************************************************************/

   x = rs*cos(theta-alpha)+a;
   y = rs*sin(theta-alpha)+b;
   z = 1/(x*x+y*y);
   u = x*(1+z)*cos(alpha)-y*(1-z)*sin(alpha);
   v = y*(1-z)*cos(alpha)+x*(1+z)*sin(alpha);

   ix = (long)(SFAC*u+XCEN);
   iy = (long)(YCEN-SFAC*v);

   if(FILL)
   {
      PLOT = FALSE;
      cx = ix;
      cy = iy;

      /*************************************************************/
      /* This subroutine also clips the display area, hence the    */
      /* following statements. Contrary to the what the Rom Kernal */
      /* manual states, all clipping must be done by the program   */
      /* if the area fill routines are used otherwise a nasty      */
      /* crash will result if plotting takes place out of bounds.  */
      /* Only the x values are clipped accurately since, for this  */
      /* routine only the x values at the boundary need to be      */
      /* accurate. The PEN parameter indicates the pen status for  */
      /* moves and draws as set by either PENUP or PENDOWN.        */
      /*************************************************************/

      if((ix <= XMIN) && (ixo > XMIN))
      {
         cy += (float)(iyo-iy)*(XMIN-ix)/(ixo-ix);
         cx = XMIN;
      }
      else if((ixo <= XMIN) && (ix > XMIN))
      {
         iyo += (float)(iy-iyo)*(XMIN-ixo)/(ix-ixo);
         ixo = XMIN;
         AreaDraw(rp,ixo,iyo);
         PLOT = TRUE;
      }
      else if((ix >= XMAX) && (ixo < XMAX))
      {
         cy += (float)(iyo-iy)*(XMAX-ix)/(ixo-ix);
         cx = XMAX;
      }
      else if((ixo >= XMAX) && (ix < XMAX))
      {
         iyo += (float)(iy-iyo)*(XMAX-ixo)/(ix-ixo);
         ixo = XMAX;
         AreaDraw(rp,ixo,iyo);
         PLOT = TRUE;
      }

      if((ixo > XMIN) && (ixo < XMAX)) PLOT = TRUE;
      if(cy < YMIN) cy = YMIN;
      if(cy > YMAX) cy = YMAX;

      if(PLOT)
      {
         if(PEN) AreaDraw(rp,cx,cy);
         else { AreaMove(rp,cx,cy); PENDOWN; }
      }
      ixo = ix;
      iyo = iy;
   }
   else
   {
      if(PEN) Draw(rp,ix,iy);
      else { Move(rp,ix,iy); PENDOWN; }
   }
}

do_request()
{
   ULONG MessageClass;

   /***************************************************************/
   /* This subroutine is activated when the double menu requester */
   /* is brought up. It handles all input from this requester.    */
   /***************************************************************/

   if(title) ShowTitle(s,FALSE);

   for (;;)
   {
      Wait(1L<<w->UserPort->mp_SigBit);
      if((message = (struct IntuiMessage *)GetMsg(w->UserPort)) != NULL)
      {
         MessageClass = message->Class;
         ReplyMsg(message);
         switch (MessageClass)
         {
            case GADGETUP    :

            case GADGETDOWN  : if (do_gadgets(message) == CLOSE_GAD)
                               {
                                  if(title) ShowTitle(s,TRUE);
                                  return;
                               }
                               break;
         } /* end switch */
      } /* end if */
   } /* end forever */
}

int do_gadgets(mes)

struct IntuiMessage *mes;
{
   struct Gadget *igad;
   int gadgid;

   /*********************************************/
   /* This subroutine handles all gadget input. */
   /*********************************************/

   igad = (struct Gadget *)mes->IAddress;
   gadgid = igad->GadgetID;

   switch(gadgid)
   {
      case CAMBER_GAD   : camber = (ULONG)camber_string.LongInt;
                          Continue = FALSE;
                          break;

      case THICK_GAD    : thickness = (ULONG)thick_string.LongInt;
                          Continue = FALSE;
                          break;

      case ANGLE_GAD    : angle = (ULONG)angle_string.LongInt;
                          Continue = FALSE;
                          break;

      case VELOCITY_GAD : velocity = (ULONG)(16-(velocity_prop.HorizPot >> 12));
                          Continue = FALSE;
                          break;

      case AIRFOIL_GAD  : if(mode) mode = FALSE;
                          else mode = TRUE;
                          Continue = FALSE;
                          break;

      case TITLE_GAD    : if(title) title = FALSE;
                          else title = TRUE;
                          break;

      case CLOSE_GAD    : break;

      case QUIT_GAD     : close_things();
                          exit(0);
                          break;
   }
   return gadgid;
}

do_about()
{
   int i;

   /*****************************************************/
   /* This subroutine displays the initial information. */
   /*****************************************************/

   SetAPen(rp,0L);
   RectFill(rp,(long)(XMIN+1),(long)(YMIN+1),(long)(XMAX-1),(long)(YMAX-1));
   SetAPen(rp,1L);
   for(i = 0;i < 24;i++)
   {
      Move(rp,2L,(long)(8*i+8));
      Text(rp,about[i],(long)strlen(about[i]));
   }
}

open_things()
{
   struct Library *OpenLibrary();
   struct Screen *OpenScreen();
   struct Window *OpenWindow();
   struct ViewPort *ViewPortAddress();
   BYTE *AllocRaster();

   if(!(GfxBase = (struct GfxBase *)OpenLibrary("graphics.library",0L)))
   {
      printf("no graphics library!!!\n");
      close_things();
      exit(1);
   }
   mask |= GRAPHICS;

   if(!(IntuitionBase = (struct IntuitionBase *)
        OpenLibrary("intuition.library",0L)))
   {
      printf("no intuition library!!!\n");
      close_things();
      exit(2);
   }
   mask |= INTUITION;

   if (!(s = OpenScreen(&ns)))
   {
      printf("could not open the screen\n");
      close_things();
      exit(3);
   }
   mask |= SCREEN;

   nw.Screen = s;

   if (!(w = OpenWindow(&nw)))
   {
      printf("could not open the window\n");
      close_things();
      exit(4);
   }
   mask |= WINDOW;

   rp = w->RPort;
   vp = ViewPortAddress(w);

   InitArea(&area,buffer,250L);
   rp->AreaInfo = &area;
   trp.Size = (long)RASSIZE(640L,400L);

   if (!(trp.RasPtr = (BYTE *)AllocRaster(640L,400L)))
   {
      printf("could not allocate memory for area fills\n");
      close_things();
      exit(5);
   }
   mask |= AREAFILL;

   rp->TmpRas = &trp;

   LoadRGB4(vp,colors,16L);

   InitRequester(&AirRequest);
   AirRequest.LeftEdge = 0L;
   AirRequest.TopEdge = 0L;
   AirRequest.Width = 200L;
   AirRequest.Height = 150L;
   AirRequest.RelLeft = -100L;
   AirRequest.RelTop = -75L;
   AirRequest.ReqGadget = &quit_gad;
   AirRequest.Flags = POINTREL;
   AirRequest.BackFill = 0;
   
   if (!(SetDMRequest(w,&AirRequest)))
   {
      printf("could not set Requester\n");
      close_things();
      exit(6);
   }
   mask |= DMREQUEST;

   ShowTitle(s,FALSE);
}

close_things()
{
   if(mask & DMREQUEST) ClearDMRequest(w,&AirRequest);
   if(mask & AREAFILL)  FreeRaster(trp.RasPtr,640L,400L);
   if(mask & WINDOW)    CloseWindow(w);
   if(mask & SCREEN)    CloseScreen(s);
   if(mask & GRAPHICS)  CloseLibrary(GfxBase);
   if(mask & INTUITION) CloseLibrary(IntuitionBase);
}