/* 0015grid.c  wing in wind tunnel, plot, get left and right grid         */
/*             wing defined by y=f(x), sequence of points around boundary */
/*             pitch about maximum chord, compute x,y,dx,dy               */
/*             for grid position, slope and normal                        */
/*             write output file  0015a**.dat  ** is angle of attack      */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <GL/glut.h>
#include <math.h>

static double xmin = 0.0;
static double xmax = 1.0;
static double xoff = 0.0; /* 2.0; works with/without offset*/
static double yoff = 0.0; /* 2.0; */
static double ymin; /* computed, depends on angle of attack */
static double ymax;
static double xmins;
static double xmaxs;
static double xtop, ytop, xbot, ybot; /* special points */
static double grid = 0.02; /* grid spacing, may be computed */
static int written = 0;
static int maxwrite = 2;
static int width = 600;
static int height = 600; /* upper part for title */
static double Pi = M_PI;
#define npts 200
static double xp[npts];
static double yp[npts];
#define outpts 100
static double leftp[outpts]; /* xinter */
static double rightp[outpts];
static double ygrid[outpts];
static int npts_lr;
static double angle_deg = 22.0; /* may be read in */
static char angle_ch[] = "22";
static char fname[] = "0015a22.dat";
static FILE* fout;

#undef min
#define min(a,b) ((a)<(b)?(a):(b))
#undef max
#define max(a,b) ((a)>(b)?(a):(b))

static double f(double x) /* basic symmetrical wing, x=0.0,1.0, need top/bot */
{
  return (15./20.)*(0.29690*sqrt(x)-0.12600*x-0.35160*x*x+
		    0.28430*x*x*x-0.10150*x*x*x*x);
} /* end f */

/* rotate wing coordinates based on angle of attack */
static void pitch(double angle, double x, double y, double *ax, double *ay)
{
  double xc;
  xc = x-0.3; /* rotate about maximum chord */
  *ax = 0.3+xc*cos(angle)+y*sin(angle);
  *ay = -xc*sin(angle)+y*cos(angle);
} /* end pitch */

static void yminmax() /* uses and sets globals for min/max */
{
  int j;
  ymax = yp[0];
  ymin = yp[0];
  xmaxs = xp[0];
  xmins = xp[0];
  ytop = yp[0]; /* for top and bottom grid point */
  xtop = xp[0];
  for(j=1; j<npts; j++)
  {
    ymax = max(ymax,yp[j]);
    ymin = min(ymin,yp[j]);
    xmaxs = max(xmaxs,xp[j]);
    xmins = min(xmins,xp[j]);
    if(yp[j]>ytop) { ytop=yp[j]; xtop = xp[j];}
  }
  grid = (ymax-ymin)/11.0; /* just first time */  
  printf("ymin=%f,  ymax=%f  \n", ymin,  ymax);
  printf("xmins=%f, xmaxs=%f \n", xmins, xmaxs);
  printf("xtop=%f,  ytop=%f  \n", xtop,  ytop);
  printf("grid step=%f \n", grid);
} /* end yminmax */

static void build_bound(double angle_deg) /* build wing boundary */
{
  int j;
  double x1, y1, x2, y2;
  double ax1, ay1, ax2, ay2;     /* top, positive surface */
  double nax1, nay1, nax2, nay2; /* negative surface */
  int n;
  double dx;
  double angle;

  angle = Pi*angle_deg/180.0;
  n = npts/2;
  dx = (xmax-xmin)/(double)(n-1);
  for(j=0; j<n; j++)
  {
    x1=xmin+j*dx;
    y1=f(x1);
    pitch(angle, x1, y1, &ax1, &ay1);
    pitch(angle, x1, -y1, &nax1, &nay1);
    xp[j] = nax1;
    yp[j] = nay1;
    xp[npts-1-j] = ax1; /* run bottom backwards */
    yp[npts-1-j] = ay1; 
  }
  pitch(angle, 1.0, 0.0, &xbot, &ybot);
  printf("build bound xbot=%f, ybot=%f \n", xbot, ybot);
} /* end build_bound */

void draw() /* wing side view */
{
  int j;

  glColor3f(1.0, 0.0, 0.0);
  for(j=1; j<npts; j++)
  {
    glBegin(GL_LINES);
      glVertex2f(xoff+xp[j-1],yoff+yp[j-1]);
      glVertex2f(xoff+xp[j],yoff+yp[j]);
    glEnd();
  }
} /* end draw */

void draw_lr() /* left/right grid */
{
  int j;

  glColor3f(0.0, 0.0, 0.0);
  for(j=0; j<npts_lr; j++)
  {
    glBegin(GL_POINTS);
      glVertex2f(xoff+leftp[j],yoff+ygrid[j]);
      glVertex2f(xoff+rightp[j],yoff+ygrid[j]);
    glEnd();
  }
} /* end draw_lr */

int inside(double x1, double y1)
                     /* an even number of crossings means outside */
{
  int j;
  int c=0; /* out, setting this one, plots points outside */
  double xinter;
  double x=x1-xoff;
  double y=y1-yoff;

  if(x<xmins) return 0; /* obvious outside */
  if(x>xmaxs) return 0;
  if(y<ymin)  return 0;
  if(y>ymax)  return 0;
  for(j=0; j<npts-1; j++) /* scan for x crossings */
  {
    if(y>max(yp[j],yp[j+1])) continue;  /* y high   */
    if(y<min(yp[j],yp[j+1])) continue;  /* y low    */
    if(x>max(xp[j],xp[j+1])) continue;  /* x beyond */
    xinter = xp[j]+(y-yp[j])*(xp[j+1]-xp[j])/(yp[j+1]-yp[j]);
    if(x>xinter) continue;              /* x beyond */
    c = 1-c; /* crossing */
  }
  return c;
} /* end inside */

void display(void)
{
  double x, y, xlow, xhigh;
  char title[] = "grid for ** degree angle of attack";
  char *p;
  int i, j;

  title[9] = angle_ch[0];
  title[10] = angle_ch[1];

  /* clear window */
  glClear(GL_COLOR_BUFFER_BIT); 
  glLoadIdentity ();

  /* Draw wings */
  draw();

  /* fill inside wing, PDE will use points outside */
  glColor3f(0.0, 0.0, 1.0);
  glPointSize(2.0);
  for(i=1; i<npts_lr-1; i++) /* not end points for plot */
  {
    y = ygrid[i];
    xlow = leftp[i];
    xhigh = rightp[i];
    for(j=0; j<npts_lr; j++)
    {
      x = leftp[j];
      if(x>xlow && x<xhigh)
      {
        glBegin(GL_POINTS);
          glVertex2f(x,y);
        glEnd();
      }
      x = rightp[j];
      if(x>xlow && x<xhigh)
      {
        glBegin(GL_POINTS);
          glVertex2f(x,y);
        glEnd();
      }
    }
  }
  glPointSize(4.0);
  draw_lr();
  glColor3f(0.0, 1.0, 0.0);
  glPointSize(6.0);
  glBegin(GL_POINTS); /* center at maximum chord */
    glVertex2f(xoff+0.3,yoff+0.0);
  glEnd();

  /* draw text, in its own context */
  glPushMatrix();
  glLoadIdentity ();
  glColor3f(0.0, 0.0, 0.0);
  glEnable(GL_LINE_SMOOTH);
  glTranslatef(0.05+xoff, 0.4+yoff, 0.0);
  glScalef(0.0004, 0.0004, 0.0);
  for(p=title; *p; p++)
    glutStrokeCharacter(GLUT_STROKE_ROMAN, *p);
  glPopMatrix();
  
  glFlush(); 

} /* end display */

/* This routine handles mouse events */
static void mouse(int button, int state, int x, int y)
{
  float wx, wy;
  
  if (button == GLUT_RIGHT_BUTTON && state == GLUT_DOWN)
  {
    exit(0);
  }
  glutPostRedisplay();
} /* end mouse */

/* This routine handles window resizes */
void reshape(int w, int h)
{
  width = w;
  height = h;
  /* Set the transformations */
  glMatrixMode(GL_PROJECTION);
  glLoadIdentity();
  glOrtho(xoff, xoff+1.00, yoff-0.5, yoff+0.5, -1.0, 1.0);
  glMatrixMode(GL_MODELVIEW);
  glViewport(0, 0, w, h);
} /* end reshape */

static void compute_lr()
{
  double x, y, xinter;
  int j, k, i=0;

  printf("   Y          XLEFT     XRIGHT \n");
  ygrid[i] = ybot+yoff; /* top point */
  leftp[i] = xbot;
  rightp[i] = xbot;
  printf(" %f  %f  %f \n", ygrid[i], leftp[i], rightp[i]);
  i = 1;
  for(y=ymin+yoff+grid/2.0; y<ymax+yoff; y=y+grid)
  {
    leftp[i] = 1000.0;
    for(j=0; j<npts-1; j++)
    {
      if(y>max(yp[j],yp[j+1])) continue;  /* y high   */
      if(y<min(yp[j],yp[j+1])) continue;  /* y low    */
      xinter = xp[j]+(y-yp[j])*(xp[j+1]-xp[j])/(yp[j+1]-yp[j]);
      leftp[i] = min(leftp[i],xinter);
    }
    rightp[i] = -1000.0;
    for(j=0; j<npts-1; j++)
    {
      if(y>max(yp[j],yp[j+1])) continue;  /* y high   */
      if(y<min(yp[j],yp[j+1])) continue;  /* y low    */
      xinter = xp[j]+(y-yp[j])*(xp[j+1]-xp[j])/(yp[j+1]-yp[j]);
      rightp[i] = max(rightp[i],xinter);
    }
    ygrid[i] = y;
    if(leftp[i]>=rightp[i]) printf("error on next line \n");
    printf(" %f  %f  %f \n", ygrid[i], leftp[i], rightp[i]);
    i++;
  }
  ygrid[i] = ytop+yoff; /* bottom point */
  leftp[i] = xtop;
  rightp[i] = xtop;
  printf(" %f  %f  %f \n", ygrid[i], leftp[i], rightp[i]);
  i++;
  npts_lr = i;
  printf("\nraw points i, x, y \n");
  for(i=0; i<npts; i++) printf("%3d  %f  %f \n", i, xp[i], yp[i]);
} /* end compute_lr */

static void write_lr()
{
  int i;

  fname[5] = angle_ch[0];
  fname[6] = angle_ch[1];
  fname[4]++;
  fout = fopen(fname, "w");
  for(i=0; i<npts_lr; i++)
  {
    fprintf(fout, " %f  %f  %f \n", ygrid[i], leftp[i], rightp[i]);
  }
  fclose(fout);
  written++;
} /* end write_lr */

static void keyboard(unsigned char key, int x, int y)
{
  if(key=='s' || key=='S') /* smaller grid */
  { 
    grid = grid/2.0;
    compute_lr();
    if(written<maxwrite) write_lr();
  }
  if(key=='l' || key=='L') /* larger grid */
  {
    grid = grid*2.0;
    written = maxwrite+1;
  }
    glutPostRedisplay();
} /* end keyboard */

static void init()
{
  int ip;
  double x,y,x1,y1,x2,y2,x3,y3;

  printf("0015grid.c running Pi=%f \n", Pi);
  /* set clear color to white */
  glClearColor (1.0, 1.0, 1.0, 0.0);
  /* set fill  color to black */
  glColor3f(0.0, 0.0, 0.0);

  glMatrixMode (GL_PROJECTION);
  glLoadIdentity ();
  glOrtho(xoff, xoff+1.00, yoff-0.5, yoff+0.5, -1.0, 1.0);
  glMatrixMode (GL_MODELVIEW);
  glViewport(0, 0, width, height);
  build_bound(angle_deg);
  yminmax();
  compute_lr();
  write_lr();
} /* end init */

int main(int argc, char* argv[])
{
  glutInit(&argc,argv);
  glutInitDisplayMode (GLUT_SINGLE | GLUT_RGB);  
  glutInitWindowSize(width, height);
  glutInitWindowPosition(20,20); 
  
  glutCreateWindow(argv[0]); 
  glutDisplayFunc(display);
  glutReshapeFunc(reshape);
  glutMouseFunc(mouse);
  glutKeyboardFunc(keyboard);
  if(argc>1)
  {
    angle_ch[0] = argv[1][0]; /* use exactly 2 digits. e.g. 00 */
    angle_ch[1] = argv[1][1];
    angle_deg = (double)atoi(angle_ch);
  }  
  init();
  glutMainLoop();
  return 0;
} /* end main of 0015grid.c */