home *** CD-ROM | disk | FTP | other *** search
/ Frozen Fish 1: Amiga / FrozenFish-Apr94.iso / bbs / alib / d1xx / d150 / airfoil.lha / Airfoil / airfoil.c < prev    next >
C/C++ Source or Header  |  1987-06-15  |  13KB  |  489 lines
  1. /********************************************************************
  2. *                                                                   *
  3. * Joukowski Airfoil Generator with Streamline and Pressure          *
  4. * Distribution Algorithms                                           *
  5. *                                                                   *
  6. * Written by:   Russell Leighton      15 March 1987                 *
  7. *               Lancaster, CA                                       *
  8. * Modified by:  David Foster          19 June 1988                  *
  9. *               Rochester, MI                                       *
  10. *               To include effect of induced circulation            *
  11. *               by a first order approximation.                     *
  12. ********************************************************************/
  13.  
  14. #include "airfoil.h"
  15.  
  16. main()
  17. {
  18.    float rs,theta,h,vel,eta,S;
  19.    int psi;
  20.    ULONG MessageClass;
  21.  
  22.    open_things();
  23.    do_about();
  24.  
  25.    /****************/
  26.    /* Loop forever */
  27.    /****************/
  28.  
  29.    for(;;)
  30.    {
  31.       while (Continue)
  32.       {
  33.          /****************************************************/
  34.          /* Wait, initially and after each plot, for user to */
  35.          /* bring up the double menu requester               */
  36.          /****************************************************/
  37.  
  38.          Wait(1L<<w->UserPort->mp_SigBit);
  39.          if((message = (struct IntuiMessage *)GetMsg(w->UserPort)) != NULL)
  40.          {
  41.             MessageClass = message->Class;
  42.             ReplyMsg(message);
  43.             if (MessageClass == REQVERIFY)
  44.             {
  45.                do_request();
  46.                break;
  47.             }
  48.          } /* end if */
  49.       } /* end while */
  50.  
  51.       do_init();
  52.  
  53.       if(mode)
  54.       {
  55.  
  56.          /********************/
  57.          /* Plot Streamlines */
  58.          /********************/
  59.  
  60.          FILL = TRUE;
  61.          for (psi = 12;psi > 0;--psi)
  62.          {
  63.             do_mess();
  64.             if(!Continue) break;
  65.  
  66.             PENUP;
  67.             SetAPen(rp,(long)(psi+1));
  68.             SetBPen(rp,(long)(psi+1));
  69.  
  70.             for (theta = 0.015;theta < TWO_PI;theta += PI/100)
  71.             {
  72.                vel = psi/(4.*velocity*r*sin(theta));
  73.                S = (1.+sin(alpha)/sin(theta));
  74.                vel = abs(S/(2.*vel));
  75.                eta = sqrt(vel*vel+1.0)-vel;
  76.                rs = r*(1.+eta)/(1.-eta);
  77.                transform(rs,theta);
  78.             } /* end for */
  79.             AreaEnd(rp);
  80.          } /* end for */
  81.  
  82.          /*******************************/
  83.          /* Plot Stagnation Streamlines */
  84.          /*******************************/
  85.  
  86.          do_mess();
  87.          if(Continue)
  88.          {
  89.             FILL = FALSE;
  90.             PENUP;
  91.             SetAPen(rp,1L);
  92.             h = (r-4.0)/40.0;
  93.             theta = -alpha;
  94.  
  95.             for (rs = 4;rs >= r;rs += h)
  96.             {
  97.                transform(rs,theta);
  98.             }
  99.  
  100.             PENUP;
  101.             theta = PI + alpha;
  102.  
  103.             for (rs = 4;rs > r;rs += h)
  104.             {
  105.                transform(rs,theta);
  106.             }
  107.          } /* end if */
  108.       } /* end if */
  109.  
  110.       else
  111.       {
  112.  
  113.          /******************************/
  114.          /* Plot Pressure Distribution */
  115.          /******************************/
  116.  
  117.          do_mess();
  118.          if(Continue)
  119.          {
  120.             FILL = TRUE;
  121.             PENUP;
  122.             SetAPen(rp,2L);
  123.             SetBPen(rp,2L);
  124.  
  125.             for (theta = 0.0;theta <= TWO_PI;theta += PI/100)
  126.             {
  127.                rs = r+(sin(theta)+sin(alpha))*(sin(theta)+sin(alpha));
  128.                transform(rs,theta);
  129.             } /* end for */
  130.             AreaEnd(rp);
  131.          } /* end if */
  132.       } /* end else */
  133.  
  134.       /****************/
  135.       /* Plot Airfoil */
  136.       /****************/
  137.  
  138.       do_mess();
  139.       if(Continue)
  140.       {
  141.          FILL = TRUE;
  142.          PENUP;
  143.          rs = r;
  144.          SetAPen(rp,1L);
  145.          SetBPen(rp,1L);
  146.  
  147.          for (theta = 0.0;theta <= TWO_PI;theta += PI/100)
  148.          {
  149.             transform(rs,theta);
  150.          }
  151.          AreaEnd(rp);
  152.       } /* end if */
  153.    } /* end forever */
  154. } /* end main */
  155.  
  156. do_init()
  157. {
  158.    float a0;
  159.  
  160.    SetAPen(rp,0L);
  161.    RectFill(rp,(long)(XMIN+1),(long)(YMIN+1),(long)(XMAX-1),(long)(YMAX-1));
  162.    SetOPen(rp,1L);
  163.  
  164.    /***********************************************************/
  165.    /* Calculate circle constants (circle center and radius)   */
  166.    /* from airfoil constants through a reverse transformation */
  167.    /***********************************************************/
  168.  
  169.    alpha = (float)angle*PI/180;
  170.    c = (float)camber/25;
  171.    t = (float)thickness/25;
  172.    a = 0;
  173.    b = c/2;
  174.  
  175.    do
  176.    {
  177.       a0 = a;
  178.       a = t*(2*a0+1)/4/sqrt(b*b+2*a0+1);
  179.       b = c*(1+2*a)/(2+2*a);
  180.    }
  181.    while(abs(a-a0) > TOL);
  182.    
  183.    r = sqrt(b*b+(a+1)*(a+1));
  184.    Continue = TRUE;
  185. }
  186.  
  187. do_mess()
  188. {
  189.    ULONG MessageClass;
  190.  
  191.    /******************************************************************/
  192.    /* Check for double menu requester. This can be brought up at any */
  193.    /* time but can not be displayed unless the message is replied    */
  194.    /* to, therefore this routine must be called periodically.        */
  195.    /******************************************************************/
  196.  
  197.    if((message = (struct IntuiMessage *)GetMsg(w->UserPort)) != NULL)
  198.    {
  199.       MessageClass = message->Class;
  200.       ReplyMsg(message);
  201.       if(MessageClass == REQVERIFY) do_request();
  202.    }
  203.    ixo = XCEN;
  204.    iyo = YCEN;
  205. }
  206.  
  207. transform(rs,theta)
  208.  
  209. float rs,theta;
  210. {
  211.    float x,y,z,u,v;
  212.    int PLOT;
  213.    long ix,iy,cx,cy;
  214.  
  215.    /********************************************************************/
  216.    /* This is the Joukowski transformation routine. This is also where */
  217.    /* the plotting is done. Plotting is usually done using the area    */
  218.    /* fill routines (AreaDraw & AreaMove). If FILL is true then the    */
  219.    /* points are used to build up the area shape to be filled. See the */
  220.    /* Rom Kernal manual containing the graphics primatives for more    */
  221.    /* details. If FILL is false then normal plotting is done.          */
  222.    /********************************************************************/
  223.  
  224.    x = rs*cos(theta-alpha)+a;
  225.    y = rs*sin(theta-alpha)+b;
  226.    z = 1/(x*x+y*y);
  227.    u = x*(1+z)*cos(alpha)-y*(1-z)*sin(alpha);
  228.    v = y*(1-z)*cos(alpha)+x*(1+z)*sin(alpha);
  229.  
  230.    ix = (long)(SFAC*u+XCEN);
  231.    iy = (long)(YCEN-SFAC*v);
  232.  
  233.    if(FILL)
  234.    {
  235.       PLOT = FALSE;
  236.       cx = ix;
  237.       cy = iy;
  238.  
  239.       /*************************************************************/
  240.       /* This subroutine also clips the display area, hence the    */
  241.       /* following statements. Contrary to the what the Rom Kernal */
  242.       /* manual states, all clipping must be done by the program   */
  243.       /* if the area fill routines are used otherwise a nasty      */
  244.       /* crash will result if plotting takes place out of bounds.  */
  245.       /* Only the x values are clipped accurately since, for this  */
  246.       /* routine only the x values at the boundary need to be      */
  247.       /* accurate. The PEN parameter indicates the pen status for  */
  248.       /* moves and draws as set by either PENUP or PENDOWN.        */
  249.       /*************************************************************/
  250.  
  251.       if((ix <= XMIN) && (ixo > XMIN))
  252.       {
  253.          cy += (float)(iyo-iy)*(XMIN-ix)/(ixo-ix);
  254.          cx = XMIN;
  255.       }
  256.       else if((ixo <= XMIN) && (ix > XMIN))
  257.       {
  258.          iyo += (float)(iy-iyo)*(XMIN-ixo)/(ix-ixo);
  259.          ixo = XMIN;
  260.          AreaDraw(rp,ixo,iyo);
  261.          PLOT = TRUE;
  262.       }
  263.       else if((ix >= XMAX) && (ixo < XMAX))
  264.       {
  265.          cy += (float)(iyo-iy)*(XMAX-ix)/(ixo-ix);
  266.          cx = XMAX;
  267.       }
  268.       else if((ixo >= XMAX) && (ix < XMAX))
  269.       {
  270.          iyo += (float)(iy-iyo)*(XMAX-ixo)/(ix-ixo);
  271.          ixo = XMAX;
  272.          AreaDraw(rp,ixo,iyo);
  273.          PLOT = TRUE;
  274.       }
  275.  
  276.       if((ixo > XMIN) && (ixo < XMAX)) PLOT = TRUE;
  277.       if(cy < YMIN) cy = YMIN;
  278.       if(cy > YMAX) cy = YMAX;
  279.  
  280.       if(PLOT)
  281.       {
  282.          if(PEN) AreaDraw(rp,cx,cy);
  283.          else { AreaMove(rp,cx,cy); PENDOWN; }
  284.       }
  285.       ixo = ix;
  286.       iyo = iy;
  287.    }
  288.    else
  289.    {
  290.       if(PEN) Draw(rp,ix,iy);
  291.       else { Move(rp,ix,iy); PENDOWN; }
  292.    }
  293. }
  294.  
  295. do_request()
  296. {
  297.    ULONG MessageClass;
  298.  
  299.    /***************************************************************/
  300.    /* This subroutine is activated when the double menu requester */
  301.    /* is brought up. It handles all input from this requester.    */
  302.    /***************************************************************/
  303.  
  304.    if(title) ShowTitle(s,FALSE);
  305.  
  306.    for (;;)
  307.    {
  308.       Wait(1L<<w->UserPort->mp_SigBit);
  309.       if((message = (struct IntuiMessage *)GetMsg(w->UserPort)) != NULL)
  310.       {
  311.          MessageClass = message->Class;
  312.          ReplyMsg(message);
  313.          switch (MessageClass)
  314.          {
  315.             case GADGETUP    :
  316.  
  317.             case GADGETDOWN  : if (do_gadgets(message) == CLOSE_GAD)
  318.                                {
  319.                                   if(title) ShowTitle(s,TRUE);
  320.                                   return;
  321.                                }
  322.                                break;
  323.          } /* end switch */
  324.       } /* end if */
  325.    } /* end forever */
  326. }
  327.  
  328. int do_gadgets(mes)
  329.  
  330. struct IntuiMessage *mes;
  331. {
  332.    struct Gadget *igad;
  333.    int gadgid;
  334.  
  335.    /*********************************************/
  336.    /* This subroutine handles all gadget input. */
  337.    /*********************************************/
  338.  
  339.    igad = (struct Gadget *)mes->IAddress;
  340.    gadgid = igad->GadgetID;
  341.  
  342.    switch(gadgid)
  343.    {
  344.       case CAMBER_GAD   : camber = (ULONG)camber_string.LongInt;
  345.                           Continue = FALSE;
  346.                           break;
  347.  
  348.       case THICK_GAD    : thickness = (ULONG)thick_string.LongInt;
  349.                           Continue = FALSE;
  350.                           break;
  351.  
  352.       case ANGLE_GAD    : angle = (ULONG)angle_string.LongInt;
  353.                           Continue = FALSE;
  354.                           break;
  355.  
  356.       case VELOCITY_GAD : velocity = (ULONG)(16-(velocity_prop.HorizPot >> 12));
  357.                           Continue = FALSE;
  358.                           break;
  359.  
  360.       case AIRFOIL_GAD  : if(mode) mode = FALSE;
  361.                           else mode = TRUE;
  362.                           Continue = FALSE;
  363.                           break;
  364.  
  365.       case TITLE_GAD    : if(title) title = FALSE;
  366.                           else title = TRUE;
  367.                           break;
  368.  
  369.       case CLOSE_GAD    : break;
  370.  
  371.       case QUIT_GAD     : close_things();
  372.                           exit(0);
  373.                           break;
  374.    }
  375.    return gadgid;
  376. }
  377.  
  378. do_about()
  379. {
  380.    int i;
  381.  
  382.    /*****************************************************/
  383.    /* This subroutine displays the initial information. */
  384.    /*****************************************************/
  385.  
  386.    SetAPen(rp,0L);
  387.    RectFill(rp,(long)(XMIN+1),(long)(YMIN+1),(long)(XMAX-1),(long)(YMAX-1));
  388.    SetAPen(rp,1L);
  389.    for(i = 0;i < 24;i++)
  390.    {
  391.       Move(rp,2L,(long)(8*i+8));
  392.       Text(rp,about[i],(long)strlen(about[i]));
  393.    }
  394. }
  395.  
  396. open_things()
  397. {
  398.    struct Library *OpenLibrary();
  399.    struct Screen *OpenScreen();
  400.    struct Window *OpenWindow();
  401.    struct ViewPort *ViewPortAddress();
  402.    BYTE *AllocRaster();
  403.  
  404.    if(!(GfxBase = (struct GfxBase *)OpenLibrary("graphics.library",0L)))
  405.    {
  406.       printf("no graphics library!!!\n");
  407.       close_things();
  408.       exit(1);
  409.    }
  410.    mask |= GRAPHICS;
  411.  
  412.    if(!(IntuitionBase = (struct IntuitionBase *)
  413.         OpenLibrary("intuition.library",0L)))
  414.    {
  415.       printf("no intuition library!!!\n");
  416.       close_things();
  417.       exit(2);
  418.    }
  419.    mask |= INTUITION;
  420.  
  421.    if (!(s = OpenScreen(&ns)))
  422.    {
  423.       printf("could not open the screen\n");
  424.       close_things();
  425.       exit(3);
  426.    }
  427.    mask |= SCREEN;
  428.  
  429.    nw.Screen = s;
  430.  
  431.    if (!(w = OpenWindow(&nw)))
  432.    {
  433.       printf("could not open the window\n");
  434.       close_things();
  435.       exit(4);
  436.    }
  437.    mask |= WINDOW;
  438.  
  439.    rp = w->RPort;
  440.    vp = ViewPortAddress(w);
  441.  
  442.    InitArea(&area,buffer,250L);
  443.    rp->AreaInfo = &area;
  444.    trp.Size = (long)RASSIZE(640L,400L);
  445.  
  446.    if (!(trp.RasPtr = (BYTE *)AllocRaster(640L,400L)))
  447.    {
  448.       printf("could not allocate memory for area fills\n");
  449.       close_things();
  450.       exit(5);
  451.    }
  452.    mask |= AREAFILL;
  453.  
  454.    rp->TmpRas = &trp;
  455.  
  456.    LoadRGB4(vp,colors,16L);
  457.  
  458.    InitRequester(&AirRequest);
  459.    AirRequest.LeftEdge = 0L;
  460.    AirRequest.TopEdge = 0L;
  461.    AirRequest.Width = 200L;
  462.    AirRequest.Height = 150L;
  463.    AirRequest.RelLeft = -100L;
  464.    AirRequest.RelTop = -75L;
  465.    AirRequest.ReqGadget = &quit_gad;
  466.    AirRequest.Flags = POINTREL;
  467.    AirRequest.BackFill = 0;
  468.    
  469.    if (!(SetDMRequest(w,&AirRequest)))
  470.    {
  471.       printf("could not set Requester\n");
  472.       close_things();
  473.       exit(6);
  474.    }
  475.    mask |= DMREQUEST;
  476.  
  477.    ShowTitle(s,FALSE);
  478. }
  479.  
  480. close_things()
  481. {
  482.    if(mask & DMREQUEST) ClearDMRequest(w,&AirRequest);
  483.    if(mask & AREAFILL)  FreeRaster(trp.RasPtr,640L,400L);
  484.    if(mask & WINDOW)    CloseWindow(w);
  485.    if(mask & SCREEN)    CloseScreen(s);
  486.    if(mask & GRAPHICS)  CloseLibrary(GfxBase);
  487.    if(mask & INTUITION) CloseLibrary(IntuitionBase);
  488. }
  489.