/******************************************************************************
Eulero codice per il calcolo del flusso potenziale su un corpo generico 
non portante tridimensionale - dal metodo di Hess Smith 

Copyright (C) 1999 Matteo Lucarelli - matteo@matteolucarelli.net

This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License
as published by the Free Software Foundation; either version 2
of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
******************************************************************************/

 /* NACA 0021 builder */
 /* solido di rotazione con sezione NACA 0021 */
 /* solo parte con z negative */
 /* con aggiunta ellittica */

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

#define TIPOD float

#define PI 3.1415926536
#define punto(a,b,c) *(punto+(a)*12+(b)*3+(c))

#define ns 21    /* numero settori */
#define np 54    /* numero spessori */

TIPOD r,a;
int i,j ;
TIPOD valx[np]={ 0,0.1,0.2,0.4,0.8,1.25,1.6,2,2.5,3,4,5,5.7,6.5,7.5,8.2,9,10,
                 11.5,13.5,15,
                 16.5,18,20,21.5,23,25,26.5,28,30,33,36,40,43,46,50,53,56,60,63,
                 66,70,73,76,80,83,86,90,92,93.5,95,96.5,98,100 };
TIPOD valr[np]={ 0,0.341,0.668,1.284,2.358,3.315,3.697,4.106,4.576,4.946,5.624,
                 6.221,6.565,
                 6.932,7.35,7.604,7.87,8.195,8.592,9.055,9.354,9.594,9.805,
                 10.04,10.17,10.28,10.39,10.44,10.48,10.50,10.45,10.36,
                 10.15,9.929,9.667,9.265,8.912,8.533,7.986,7.539,7.071,
                 6.412,5.890,5.347,4.591,4.000,3.386,2.534,2.093,1.755,
                 1.412,1.061,0.705,0.221	 };

void main(void)
{
   TIPOD inserto(TIPOD x,TIPOD y);
   TIPOD *punto;
   int nelem ;
   FILE *fp;

   nelem = (np-1)*ns ;
   punto = (TIPOD *)calloc((size_t)(12*nelem) , sizeof(TIPOD) ) ;


   for (i=0;i<(np-1);i++)
      for (j=0;j<ns;j++)
      {
         r=valr[i];
         a=(PI/ns*j);

         punto(i*ns+j,0,0)=valx[i];
         punto(i*ns+j,0,2)=-r*sin(a);
         if (cos(a)>0)
         punto(i*ns+j,0,1)=r*cos(a)+inserto(punto(i*ns+j,0,0),punto(i*ns+j,0,2));
         else  punto(i*ns+j,0,1)=r*cos(a);

         a=(PI/ns*(j+1));

         punto(i*ns+j,1,0)=valx[i];
         punto(i*ns+j,1,2)=-r*sin(a);
         if (cos(a)>0)
         punto(i*ns+j,1,1)=r*cos(a)+inserto(punto(i*ns+j,1,0),punto(i*ns+j,1,2));
         else  punto(i*ns+j,1,1)=r*cos(a);
         r=valr[i+1];

         punto(i*ns+j,2,0)=valx[i+1];
         punto(i*ns+j,2,2)=-r*sin(a);
         if (cos(a)>0)
         punto(i*ns+j,2,1)=r*cos(a)+inserto(punto(i*ns+j,2,0),punto(i*ns+j,2,2));
         else  punto(i*ns+j,2,1)=r*cos(a);
         a=(PI/ns*j);

         punto(i*ns+j,3,0)=valx[i+1];
         punto(i*ns+j,3,2)=-r*sin(a);
         if (cos(a)>0)
         punto(i*ns+j,3,1)=r*cos(a)+inserto(punto(i*ns+j,3,0),punto(i*ns+j,3,2));
         else  punto(i*ns+j,3,1)=r*cos(a);
      }

   if((fp=fopen("NAC2154i.din" , "wb" ))==NULL)
   {   printf("errore apertura file");
       exit(1);
   }
   fwrite(&nelem,sizeof(int),1,fp);     printf("eseguito\n");

   fwrite(punto,sizeof(TIPOD),(size_t)(nelem*12),fp);
   fclose(fp);
}

TIPOD inserto(TIPOD x,TIPOD z)
{
   TIPOD y=0;

   if ((x>25)&&(x<50)&&(z>-4.17)&&(z<4.17))
      y=4*sqrt((1.0-z*z/17.389)*(1.0-(x-37.5)*(x-37.5)/156.25));

   return  y;
}