/******************************************************************************
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.
******************************************************************************/

 /* SPHereBuiLDer SiMmetric */
 /* divisione aree proporzionate */
 /* costruzione sfere simmetriche rispetto a z */

#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))

void main (void)
{
   int n_par,n_mer,n_elem;
   int m,n;
   TIPOD *punto;
   FILE *fp;
   int j;
   TIPOD y,alfa;
   char nomein[20];

   printf("sphere builder as TIPOD\n\n");
   printf("Nø paralleli (min 1) :");
   scanf("%d",&n_par);
   printf("Nø meridiani (min 2) :");
   scanf("%d",&n_mer);

   n_par++ ;

   n_elem = n_par * n_mer ;

   printf("%d elementi\n",n_elem);

   punto = (TIPOD *)calloc( 12*n_elem , sizeof(TIPOD) ) ;

   for(n=0;n<n_par;n++)
      for(m=0;m<n_mer;m++)
      {
         y = cos( n * PI / (n_par*2) );
         /* y=y*y*y;  */
         y = sqrt(fabs(y))*y/fabs(y);
         /* y = -(2.0*n/n_par) + 1.0 ;   */
         alfa = m * 2*PI / n_mer ;

         punto(n*n_mer+m,0,2)=y;
         punto(n*n_mer+m,0,0)=sqrt(1-y*y)*sin(alfa);
         punto(n*n_mer+m,0,1)=-sqrt(1-y*y)*cos(alfa);

         alfa = (m+1) * 2*PI / n_mer ;

         punto(n*n_mer+m,1,2)=y;
         punto(n*n_mer+m,1,0)=sqrt(1-y*y)*sin(alfa);
         punto(n*n_mer+m,1,1)=-sqrt(1-y*y)*cos(alfa);

         y= cos( (n+1) * PI / (n_par*2) );
         /*y=y*y*y; */
         y = sqrt(fabs(y))*y/fabs(y);
         /*y = -(2.0*(n+1.0)/n_par) + 1.0 ;       */

         punto(n*n_mer+m,2,2)=y;
         punto(n*n_mer+m,2,0)=sqrt(1-y*y)*sin(alfa);
         punto(n*n_mer+m,2,1)=-sqrt(1-y*y)*cos(alfa);

         alfa = m * 2*PI / n_mer ;

         punto(n*n_mer+m,3,2)=y;
         punto(n*n_mer+m,3,0)=sqrt(1-y*y)*sin(alfa);
         punto(n*n_mer+m,3,1)=-sqrt(1-y*y)*cos(alfa);
      }

   printf("nome del file (max 8 char) :");
   scanf("%s",nomein);
   nomein[8]=0;


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

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