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

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

#define SQ(a) ((a)*(a))
#define MAX(a,b) ( (a) > (b) ? (a) : (b) )
#define MAXNUM 1e+6    /* sostitutivo di infinito */

#define TIPOD float    /******** tipo di dati *********/

#define punto(a,b,c) *(punto+(a)*12+(b)*3+c)
#define prel(a,b,c) *(prel+(a)*8+(b)*2+c)
#define pmedio(a,b) *(pmedio+(a)*3+b)
#define rotazione(a,b,c) *(rotazione+(a)*9+(b)*3+c)
#define np(a,b) *(np+(a)*3+b)
#define minf(a,b) *(minf+(a)*nelem+b)
#define vnex(a) *(vnex+a)
#define intens(a) *(intens+a)

#define d12(a) *(d12+a)
#define d23(a) *(d23+a)
#define d34(a) *(d34+a)
#define d41(a) *(d41+a)
#define m12(a) *(m12+a)
#define m23(a) *(m23+a)
#define m34(a) *(m34+a)
#define m41(a) *(m41+a)
#define etasud12(a) *(etasud12+a)
#define etasud23(a) *(etasud23+a)
#define etasud34(a) *(etasud34+a)
#define etasud41(a) *(etasud41+a)
#define xisud12(a) *(xisud12+a)
#define xisud23(a) *(xisud23+a)
#define xisud34(a) *(xisud34+a)
#define xisud41(a) *(xisud41+a)
#define area(n) *(area+n)

#define assc(a) *(assc+a)
#define relc(a) *(relc+a)

#define npr(a) *(npr+a)

/* variabili globali allocate dinamicamente nella sub loaddati */

char nome[20]=" ";       /* nome del file (.din:dati input .don:dati output) */
int nelem=0;             /* n° di elementi del corpo */
int tiposimm;            /* tipo di simmetria applicabile */
TIPOD vex[3];           /* vettore velocità esterna */
TIPOD *punto=NULL;      /* matrice nelem*4*3 :coordinate vertici degli elementi
                            ordinati in senso orario visti dal lato del flusso*/
TIPOD *pmedio=NULL;     /* matrice nelem*3 : coordinate punti medi */
TIPOD *rotazione=NULL;  /* matrice nelem*3*3 : tensori di rotazione */
TIPOD *np=NULL;         /* matrice nelem*3 : coordinate dei null point */
TIPOD *minf=NULL;       /* matrice nelem*nelem : matrice di influenza */
TIPOD *vnex=NULL;       /* vettore nelem : comp. normale vel. esterna */
TIPOD *intens=NULL;     /* vettore nelem : intensità sorgenti (solution) */
                                                                                                                                                       
/* quantità globali utili */
                                                                                                                                           
TIPOD *prel=NULL;                    /* array nelem*4*2 */
TIPOD *d12=NULL;                     /* vettore (nelem) */
TIPOD *d23=NULL;                     /*    "       "    */
TIPOD *d34=NULL;                     /*    "       "    */
TIPOD *d41=NULL;                     /*    "       "    */
TIPOD *m12=NULL;                     /*    "       "    */
TIPOD *m23=NULL;                     /*    "       "    */
TIPOD *m34=NULL;                     /*    "       "    */
TIPOD *m41=NULL;                     /*    "       "    */
TIPOD *etasud12=NULL;                /*    "       "    */
TIPOD *etasud23=NULL;                /*    "       "    */
TIPOD *etasud34=NULL;                /*    "       "    */
TIPOD *etasud41=NULL;                /*    "       "    */
TIPOD *xisud12=NULL;                 /*    "       "    */
TIPOD *xisud23=NULL;                 /*    "       "    */
TIPOD *xisud34=NULL;                 /*    "       "    */
TIPOD *xisud41=NULL;                 /*    "       "    */
TIPOD *area=NULL;                    /*    "       "    */

/****************************************MAIN**********************************/

void main(void)
{
   int r;
   char c;
   void loaddati(void);
   void calgeo(void);
   void calutili(void);
   void calnp(void);
   void calminf(void);
   void risolvi(void);
   void calvelnp(int r);
   void savedati(void);
   void savealtri(int r);

   printf("EULERO\ncalcolo del flusso potenziale\n");
   
   loaddati();  /* carica  da disco le coordinate dei punti del corpo */
                /* elimina gli elementi di area nulla */
                /* alloca gli array dinamici globali */
   calgeo();    /* calcoli geometrici */
   calutili();  /* calcola costanti utili */
   calnp();     /* calcola i null points */
   calminf();   /* calcola la matrice di influenza */
   risolvi();   /* calcola la vel esterna e risolve il sistema */

   do
   {
      printf("\nInserire una opzione\n");
      printf("  1-uscita\n  2-ricalcolo intesità con altra vel ext\n");
      printf("  3-calcolo velocità nei null point\n");
      printf("  4-salvataggio coordinate proiettate\n");
      printf("  5-salvataggio coordinate null point\n");
      printf("  6-salvataggio matrice di influenza\n");
      printf("  7-salvataggio vettore termini noti\n");
      printf("  8-salvataggio intensità\n");
      printf("  9-calcolo Cpr nei null point\n? ");

      scanf("%s",&c);
      r=c-48;

      if (r==2) { risolvi(); savedati(); }
      if (r==3) calvelnp(r);
      if ((r>3)&&(r<8)) savealtri(r);
      if (r==8) savedati();
      if (r==9) calvelnp(r);
   }
   while (!(r==1)) ;

   printf("Esecuzione terminata\n");
   exit(0);
}

/**********************************LOADDATI************************************/
/* routine di caricamento dati */
/* legge da disco i dati di input dal file nome.din */
/* i dati sono in formato binario */
/* primi 4 byte (intero) : n° elementi del corpo in esame */
/* successivi nelem*4*3*8 byte (TIPOD) : array punto */ 

void loaddati(void)
{
   FILE *fp ;
   size_t sod ;
   int j,k,b,c ;
   char r ;

/* caricamento dati da disco */

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

   printf("Tipo di simmetria applicabile :\n- 0 : nessuna\n");
   printf("- 1 : z -> -z\n");
   printf("- 2 : x -> -x\n- 3 : y -> -y");
   do
   {
      printf("\n? ");
      scanf("%s",&r) ;
      tiposimm=r-48 ;
   }
   while((tiposimm<0)||(tiposimm>3));

   if((fp=fopen(strcat(nome,".din") , "rb"))==NULL)
   {  printf("errore apertura file");
      exit(1);  }

   printf("Caricamento dati dal file %s\n",nome);
   nome[strlen(nome)-4]=0;                         /* toglie l'estensione */

   sod=sizeof(TIPOD);
   fread(&nelem,sizeof(int),1,fp);                 /* carica nelem */
   punto=(TIPOD *)calloc((size_t)(12*nelem),sod); /* alloca punto */
   fread(punto,sod,(size_t)12*nelem,fp);           /* carica array punto */
   fclose(fp);

   printf("N'elementi caricati :%d\n",nelem);

/* eliminazione elementi con area nulla */

   printf("Controllo elementi ad area nulla :");

   for(j=0;j<nelem;j++)
      if (  ((punto(j,0,0)==punto(j,2,0))&&
             (punto(j,0,1)==punto(j,2,1))&&
             (punto(j,0,2)==punto(j,2,2)) )
          ||((punto(j,1,0)==punto(j,3,0))&&
             (punto(j,1,1)==punto(j,3,1))&&
             (punto(j,1,2)==punto(j,3,2)) ) )
      {
         nelem--;
         for(k=j;k<nelem;k++)
            for(b=0;b<4;b++)
               for(c=0;c<3;c++)
                  punto(k,b,c)=punto((k+1),b,c);

         printf("\n Eliminato %d°elemento \n",j) ;
      }

   printf(" OK\n");

/* allocazioni degli array globali */

   pmedio=(TIPOD *)calloc((size_t)(3*nelem),sod);
   rotazione=(TIPOD *)calloc((size_t)(9*nelem),sod);
   np=(TIPOD *)calloc((size_t)(3*nelem),sod);

   minf=(TIPOD *)calloc((size_t)(nelem*nelem),sod);
   vnex=(TIPOD *)calloc((size_t)nelem,sod);
   intens=(TIPOD *)calloc((size_t)nelem,sod);

   prel=(TIPOD *)calloc((size_t)(nelem*8),sod);

   d12=(TIPOD *)calloc((size_t)nelem,sod);
   d23=(TIPOD *)calloc((size_t)nelem,sod);
   d34=(TIPOD *)calloc((size_t)nelem,sod);
   d41=(TIPOD *)calloc((size_t)nelem,sod);

   m12=(TIPOD *)calloc((size_t)nelem,sod);
   m23=(TIPOD *)calloc((size_t)nelem,sod);
   m34=(TIPOD *)calloc((size_t)nelem,sod);
   m41=(TIPOD *)calloc((size_t)nelem,sod);

   etasud12=(TIPOD *)calloc((size_t)nelem,sod);
   etasud23=(TIPOD *)calloc((size_t)nelem,sod);
   etasud34=(TIPOD *)calloc((size_t)nelem,sod);
   etasud41=(TIPOD *)calloc((size_t)nelem,sod);

   xisud12=(TIPOD *)calloc((size_t)nelem,sod);
   xisud23=(TIPOD *)calloc((size_t)nelem,sod);
   xisud34=(TIPOD *)calloc((size_t)nelem,sod);
   xisud41=(TIPOD *)calloc((size_t)nelem,sod);

   area=(TIPOD *)calloc((size_t)nelem,sod);

/* controllo allocazioni */

   printf("Controllo allocazioni :");
   if (!(punto||pmedio||rotazione||np||minf||vnex||intens||prel||d12||d23||d34||
        d41||m12||m23||m34||m41||etasud12||etasud23||etasud34||etasud41||
        xisud12||xisud23||xisud34||xisud41||area))
      {
       printf("MEMORIA INSUFFICENTE\n");
       exit(1);
      }
   printf("OK , memoria sufficiente\n");

   return;          
}

/*******************************CALGEO*****************************************/
/**********************Routine di calcolo geometrico***************************/

void calgeo(void)
{
   int k,j,i ;
   int n ;
   TIPOD d1[3];
   TIPOD d2[3];
   TIPOD ver[3];
   TIPOD modulo,t,ro ;

   printf("Calcoli geometrici :");

   for(n=0;n<nelem;n++)
   {
   /* calcolo diagonali */

      for(j=0;j<3;j++)
      {
         d1[j] = punto(n,3,j) - punto(n,1,j) ;
         d2[j] = punto(n,2,j) - punto(n,0,j) ;
      }

   /* calcolo versore */

      ver[0] = d1[1]*d2[2] - d1[2]*d2[1];
      ver[1] = d1[2]*d2[0] - d1[0]*d2[2];
      ver[2] = d1[0]*d2[1] - d1[1]*d2[0];

      modulo = sqrt(SQ(ver[0])+SQ(ver[1])+SQ(ver[2]));

      ver[0] /= modulo;
      ver[1] /= modulo;
      ver[2] /= modulo;

   /* calcolo del centroide */

      for(j=0;j<3;j++)
      {
         pmedio(n,j)=0 ;
         for(i=0;i<4;i++)
            pmedio(n,j) += punto(n,i,j) ;
         pmedio(n,j) /= 4;
      }

   /* proiezione vertici dell'elemento */

      for(k=0;k<4;k++)
      {
         t=0;
         for(i=0;i<3;i++)
            t += (ver[i] * ( pmedio(n,i) - punto(n,k,i) ));

         for(i=0;i<3;i++)
            punto(n,k,i) += (ver[i] * t);
      }

   /* diagonale elemento proiettato */

      for(j=0;j<3;j++)
         d2[j] = punto(n,2,j) - punto(n,0,j) ;

      ro=sqrt( SQ(d2[0]) + SQ(d2[1]) + SQ(d2[2]) );      /* = modulo di d2 */

   /* tensore di rotazione */

      for(k=0;k<3;k++)
      {
         rotazione(n,k,2) = ver[k] ;
         rotazione(n,k,0) = d2[k]/ro ;
      }
      rotazione(n,0,1)=ver[1]*d2[2]/ro - ver[2]*d2[1]/ro ;
      rotazione(n,1,1)=ver[2]*d2[0]/ro - ver[0]*d2[2]/ro ;
      rotazione(n,2,1)=ver[0]*d2[1]/ro - ver[1]*d2[0]/ro ;

   }
   printf("Eseguiti\n");

   return;
}

/**************************TRASF----ANTITRASF**********************************/
/*   funzioni di rototraslazione elemento  *
 *   passaggio dalle coordinate assolute   *
 *   a quelle relative all'elemento n°     *
 *  (elemento n°nel piano xy e y//diag12 ) */

void trasf(TIPOD *assc,int n,int trasla)
{
   int k;
   TIPOD relc[3];

   if (trasla)
      for(k=0;k<3;k++)
         assc(k) -= pmedio(n,k) ;

   for(k=0;k<3;k++)
      relc(k) = assc(0) * rotazione(n,0,k) +
                assc(1) * rotazione(n,1,k) +
                assc(2) * rotazione(n,2,k) ;

   for(k=0;k<3;k++)
      assc(k) = relc(k);

   return;
}

void antitrasf(TIPOD *relc,int n,int trasla)
{
   int k;
   TIPOD assc[3];

   for(k=0;k<3;k++)
      assc(k)=relc(0) * rotazione(n,k,0) +
              relc(1) * rotazione(n,k,1) +
              relc(2) * rotazione(n,k,2) ;

   if (trasla)
      for(k=0;k<3;k++)
         assc(k) += pmedio(n,k);

   for(k=0;k<3;k++)
      relc(k)=assc(k);

   return;
}

/**********************************CALUTILI************************************/
/* funzione di calcolo quantità utili    *
 * ripetitive nei calcoli di nullpoint   *
 * e matrice di influenza                */

void calutili(void)
{
   void trasf(TIPOD *abs_coord,int elem_rif,int flag_trasla);
   int n;
   int k,j;
   TIPOD tempor[3];

   printf("Calcolo costanti utili : ");

   for(n=0;n<nelem;n++)
   {
   /* coordinate relative dei vertici di ogni elemento */

      for(k=0;k<4;k++)
      {
         for(j=0;j<3;j++)
            tempor[j]=punto(n,k,j);     
         trasf(tempor,n,1);
         prel(n,k,0)=tempor[0];
         prel(n,k,1)=tempor[1];
      }

   /* altre quantità utili */

      d12(n)=sqrt(SQ(prel(n,1,0)-prel(n,0,0))+SQ(prel(n,1,1)-prel(n,0,1)));
      d23(n)=sqrt(SQ(prel(n,2,0)-prel(n,1,0))+SQ(prel(n,2,1)-prel(n,1,1)));
      d34(n)=sqrt(SQ(prel(n,3,0)-prel(n,2,0))+SQ(prel(n,3,1)-prel(n,2,1)));
      d41(n)=sqrt(SQ(prel(n,0,0)-prel(n,3,0))+SQ(prel(n,0,1)-prel(n,3,1)));

      if (prel(n,1,0)!=prel(n,0,0))
         m12(n)=(prel(n,1,1)-prel(n,0,1))/(prel(n,1,0)-prel(n,0,0));
      else m12(n)=MAXNUM*(prel(n,1,1)-prel(n,0,1));
      if (prel(n,2,0)-prel(n,1,0))
         m23(n)=(prel(n,2,1)-prel(n,1,1))/(prel(n,2,0)-prel(n,1,0));
      else m23(n)=MAXNUM*(prel(n,2,1)-prel(n,1,1));
      if (prel(n,3,0)-prel(n,2,0))
         m34(n)=(prel(n,3,1)-prel(n,2,1))/(prel(n,3,0)-prel(n,2,0));
      else m34(n)=MAXNUM*(prel(n,3,1)-prel(n,2,1));
      if (prel(n,0,0)-prel(n,3,0))
         m41(n)=(prel(n,0,1)-prel(n,3,1))/(prel(n,0,0)-prel(n,3,0));
      else m41(n)=MAXNUM*(prel(n,0,1)-prel(n,3,1));

      if (d12(n))
         etasud12(n)=(prel(n,1,1)-prel(n,0,1))/d12(n);
      else etasud12(n)=MAXNUM*(prel(n,1,1)-prel(n,0,1));
      if (d23(n))
         etasud23(n)=(prel(n,2,1)-prel(n,1,1))/d23(n);
      else etasud23(n)=MAXNUM*(prel(n,2,1)-prel(n,1,1));
      if (d34(n))
         etasud34(n)=(prel(n,3,1)-prel(n,2,1))/d34(n);
      else etasud34(n)=MAXNUM*(prel(n,3,1)-prel(n,2,1));
      if (d41(n))
         etasud41(n)=(prel(n,0,1)-prel(n,3,1))/d41(n);
      else etasud41(n)=MAXNUM*(prel(n,0,1)-prel(n,3,1));

      if (d12(n))
         xisud12(n)=(prel(n,0,0)-prel(n,1,0))/d12(n);
      else xisud12(n)=MAXNUM*(prel(n,0,0)-prel(n,1,0));
      if (d23(n))
         xisud23(n)=(prel(n,1,0)-prel(n,2,0))/d23(n);
      else xisud23(n)=MAXNUM*(prel(n,1,0)-prel(n,2,0));
      if (d34(n))
         xisud34(n)=(prel(n,2,0)-prel(n,3,0))/d34(n);
      else xisud34(n)=MAXNUM*(prel(n,2,0)-prel(n,3,0));
      if (d41(n))
         xisud41(n)=(prel(n,3,0)-prel(n,0,0))/d41(n);
      else xisud41(n)=MAXNUM*(prel(n,3,0)-prel(n,0,0));

   /* aree e momenti degli elementi */

      area(n)=0.5*(prel(n,2,0)-prel(n,0,0))*(prel(n,1,1)-prel(n,3,1));

   }
   printf("Eseguito \n");

   return;
}

/*********************************CALNP****************************************/
/*******************calcolo dei null points (metodo di Newton)*****************/

void calnp(void)
{
   void antitrasf(TIPOD *rel_coor,int elem_rif,int flag_trasla);
   int k;
   int n;
   TIPOD cnp[2];
   TIPOD cnpnext[3];
   TIPOD r1,r2,r3,r4;
   TIPOD r1dx,r1dy,r2dx,r2dy,r3dx,r3dy,r4dx,r4dy;
   TIPOD loguno,logdue,logtre,logqtr;
   TIPOD logunodx,logduedx,logtredx,logqtrdx;
   TIPOD logunody,logduedy,logtredy,logqtrdy;
   TIPOD vx,vy,vxdx,vxdy,vydx,vydy;
   TIPOD errmax,passox,passoy ;

   printf("Calcolo Null Points :");

   for(n=0;n<nelem;n++)
   {

      errmax=sqrt(area(n))/100;

      cnpnext[0]=( prel(n,0,0)+prel(n,1,0)+prel(n,2,0)+prel(n,3,0) ) / 4 ;
      cnpnext[1]=( prel(n,0,1)+prel(n,1,1)+prel(n,2,1)+prel(n,3,1) ) / 4 ;
      cnpnext[2]=0;

      do
      {
         cnp[0]=cnpnext[0];
         cnp[1]=cnpnext[1];

         r1=sqrt( SQ(cnp[0]-prel(n,0,0)) + SQ(cnp[1]-prel(n,0,1)) );
         r2=sqrt( SQ(cnp[0]-prel(n,1,0)) + SQ(cnp[1]-prel(n,1,1)) );
         r3=sqrt( SQ(cnp[0]-prel(n,2,0)) + SQ(cnp[1]-prel(n,2,1)) );
         r4=sqrt( SQ(cnp[0]-prel(n,3,0)) + SQ(cnp[1]-prel(n,3,1)) );

         if (!(( r1+r2-d12(n) )/( r1+r2+d12(n) )))
            loguno=-MAXNUM;
            else loguno=log( ( r1+r2-d12(n) )/( r1+r2+d12(n) ) );

         if (!(( r2+r3-d23(n) )/( r2+r3+d23(n) )))
            logdue=-MAXNUM;
            else logdue=log( ( r2+r3-d23(n) )/( r2+r3+d23(n) ) );

         if (!(( r3+r4-d34(n) )/( r3+r4+d34(n) )))
            logtre=-MAXNUM;
            else logtre=log( ( r3+r4-d34(n) )/( r3+r4+d34(n) ) );

         if (!(( r4+r1-d41(n) )/( r4+r1+d41(n) )))
            logqtr=-MAXNUM;
            else logqtr=log( ( r4+r1-d41(n) )/( r4+r1+d41(n) ) );

         vx= etasud12(n) * loguno + etasud23(n) * logdue +
             etasud34(n) * logtre + etasud41(n) * logqtr ;
         vy= xisud12(n) * loguno + xisud23(n) * logdue +
             xisud34(n) * logtre + xisud41(n) * logqtr ;

         if (!((vx==0)&(vy==0)))
         {
            r1dx=( cnp[0] - prel(n,0,0) ) / r1;
            r1dy=( cnp[1] - prel(n,0,1) ) / r1;
            r2dx=( cnp[0] - prel(n,1,0) ) / r2;
            r2dy=( cnp[1] - prel(n,1,1) ) / r2;
            r3dx=( cnp[0] - prel(n,2,0) ) / r3;
            r3dy=( cnp[1] - prel(n,2,1) ) / r3;
            r4dx=( cnp[0] - prel(n,3,0) ) / r4;
            r4dy=( cnp[1] - prel(n,3,1) ) / r4;

            logunodx=(2.0*d12(n)*(r1dx+r2dx))/((r1+r2-d12(n))*(r1+r2+d12(n)));
            logunody=(2.0*d12(n)*(r1dy+r2dy))/((r1+r2-d12(n))*(r1+r2+d12(n)));
            logduedx=(2.0*d23(n)*(r2dx+r3dx))/((r2+r3-d23(n))*(r2+r3+d23(n)));
            logduedy=(2.0*d23(n)*(r2dy+r3dy))/((r2+r3-d23(n))*(r2+r3+d23(n)));
            logtredx=(2.0*d34(n)*(r3dx+r4dx))/((r3+r4-d34(n))*(r3+r4+d34(n)));
            logtredy=(2.0*d34(n)*(r3dy+r4dy))/((r3+r4-d34(n))*(r3+r4+d34(n)));
            logqtrdx=(2.0*d41(n)*(r4dx+r1dx))/((r4+r1-d41(n))*(r4+r1+d41(n)));
            logqtrdy=(2.0*d41(n)*(r4dy+r1dy))/((r4+r1-d41(n))*(r4+r1+d41(n)));

            vxdx= etasud12(n) * logunodx + etasud23(n) * logduedx +
                  etasud34(n) * logtredx + etasud41(n) * logqtrdx ;
            vydx= xisud12(n) * logunodx + xisud23(n) * logduedx +
                  xisud34(n) * logtredx + xisud41(n) * logqtrdx ;
            vxdy= etasud12(n) * logunody + etasud23(n) * logduedy +
                  etasud34(n) * logtredy + etasud41(n) * logqtrdy ;
            vydy= xisud12(n) * logunody + xisud23(n) * logduedy +
                  xisud34(n) * logtredy + xisud41(n) * logqtrdy ;

            passox= (vx*vydy-vy*vxdy)/(vxdx*vydy-vxdy*vydx) ;
            passoy= (vy*vxdx-vx*vydx)/(vxdx*vydy-vxdy*vydx) ;

            cnpnext[0]=cnp[0]-passox ;
            cnpnext[1]=cnp[1]-passoy ;

            while ((cnpnext[0]>prel(n,2,0))||(cnpnext[0]<prel(n,0,0))||
                    (cnpnext[1]>prel(n,1,1))||(cnpnext[1]<prel(n,3,1)))
            {
               passox /= 2 ;
               passoy /= 2 ;
               cnpnext[0] = cnp [0] - passox ;
               cnpnext[1] = cnp [1] - passoy ;
            }
         }
      }
      while ((sqrt(SQ(cnp[0]-cnpnext[0])+SQ(cnp[1]-cnpnext[1]))>errmax)&&
             (!((vx==0)&(vy==0))));

      antitrasf(cnpnext,n,1);
      for(k=0;k<3;k++)
         np(n,k) = cnpnext[k];
   }
   printf("Eseguito\n");

   return;
}

/******************************CALMINF*****************************************/
/******************rutine di calcolo della matrice di influenza****************/

void calminf(void)
{
   void antitrasf(TIPOD *rel_coord,int elem_rif,int flag_trasla);
   void trasf(TIPOD *ass_coord,int elem_rif,int flag_trasla);
   void esatto(TIPOD *npr,int n,TIPOD *vel);

   int n,k,j;               /* contatori */
   TIPOD vel[3],npr[3];    /* velocità influente null point relativo */

   printf("Calcolo della matrice di influenza :     ");

   for(n=0;n<nelem;n++)
   {

      printf("\b\b\b\b%3.0f%",((TIPOD)100*n/nelem));   /* contatore */

      for(k=0;k<nelem;k++)
      {
         if (k==n)
            minf(k,n)=6.2831853;
         else
         {
         /* calcola coordinate relative a n del nullpoint k-esimo */

            for(j=0;j<3;j++)
               npr[j]=np(k,j);
            trasf(npr,n,1);

            esatto(npr,n,vel);

         /* calcolo dell'elemento */

            antitrasf(vel,n,0);
            minf(k,n)=vel[0]*rotazione(k,0,2)+vel[1]*rotazione(k,1,2)
                      +vel[2]*rotazione(k,2,2);
         }

         /* presenza di simmetria */

         if (tiposimm!=0)
         {
            if (tiposimm==2) npr[0]=-np(k,0);
              else npr[0]=np(k,0);
            if (tiposimm==3) npr[1]=-np(k,1);
              else npr[1]=np(k,1);
            if (tiposimm==1) npr[2]=-np(k,2);
              else npr[2]=np(k,2);

            trasf(npr,n,1);
            esatto(npr,n,vel);
            antitrasf(vel,n,0);

            if (tiposimm==2) minf(k,n)-=vel[0]*rotazione(k,0,2);
              else minf(k,n)+=vel[0]*rotazione(k,0,2);
            if (tiposimm==3) minf(k,n)-=vel[1]*rotazione(k,1,2);
              else minf(k,n)+=vel[1]*rotazione(k,1,2);
            if (tiposimm==1) minf(k,n)-=vel[2]*rotazione(k,2,2);
              else minf(k,n)+=vel[2]*rotazione(k,2,2);
         }
      }
   }
   printf("\b\b\b\b\bEseguito\n");

   return;
}

/**********************************ESATTO**************************************/
/********************calcolo della velocità relativa***************************/

void esatto(TIPOD *npr,int n,TIPOD *vel)
{
   TIPOD r[4];     /* qtà utili */
   TIPOD e[4];     /*  "    "   */
   TIPOD h[4];     /*  "    "   */
   TIPOD loguno,logdue,logtre,logqtr;
   int j;

/* quantità ripetute */

   for(j=0;j<4;j++)
   {
      r[j]=sqrt(SQ(npr(0)-prel(n,j,0))+SQ(npr(1)-prel(n,j,1))+SQ(npr(2)));
      e[j]= SQ( npr(2) ) + SQ( npr(0)-prel(n,j,0) );
      h[j]=( npr(1)-prel(n,j,1) ) * ( npr(0)-prel(n,j,0) );
   }

   if (!( ( r[0]+r[1]-d12(n) )/( r[0]+r[1]+d12(n) ) ))
      loguno=-MAXNUM;
      else loguno=log( ( r[0]+r[1]-d12(n) )/( r[0]+r[1]+d12(n) ) );
   if (!( ( r[1]+r[2]-d23(n) )/( r[1]+r[2]+d23(n) ) ))
      logdue=-MAXNUM;
      else logdue=log( ( r[1]+r[2]-d23(n) )/( r[1]+r[2]+d23(n) ) );
   if (!( ( r[2]+r[3]-d34(n) )/( r[2]+r[3]+d34(n) ) ))
      logtre=-MAXNUM;
      else logtre=log( ( r[2]+r[3]-d34(n) )/( r[2]+r[3]+d34(n) ) );
   if (!( ( r[3]+r[0]-d41(n) )/( r[3]+r[0]+d41(n) ) ))
      logqtr=-MAXNUM;
      else logqtr=log( ( r[3]+r[0]-d41(n) )/( r[3]+r[0]+d41(n) ) );


/* velocità */

   *(vel) = etasud12(n)*loguno + etasud23(n)*logdue +
            etasud34(n)*logtre + etasud41(n)*logqtr ;
   *(vel+1) = xisud12(n)*loguno + xisud23(n)*logdue +
              xisud34(n)*logtre + xisud41(n)*logqtr ;
   *(vel+2) = atan( (m12(n)*e[0]-h[0]) / (npr(2)*r[0]) )-
              atan( (m12(n)*e[1]-h[1]) / (npr(2)*r[1]) )+
              atan( (m23(n)*e[1]-h[1]) / (npr(2)*r[1]) )-
              atan( (m23(n)*e[2]-h[2]) / (npr(2)*r[2]) )+
              atan( (m34(n)*e[2]-h[2]) / (npr(2)*r[2]) )-
              atan( (m34(n)*e[3]-h[3]) / (npr(2)*r[3]) )+
              atan( (m41(n)*e[3]-h[3]) / (npr(2)*r[3]) )-
              atan( (m41(n)*e[0]-h[0]) / (npr(2)*r[0]) );
   return;
}

/*******************************RISOLVI****************************************/
/************************calcolo del sistema risolutivo************************/

void risolvi(void)
{
   TIPOD pivot,molt,quot;
   int n,k,j;

/* inserimento velocità esterna */

   vex[0]=0;vex[1]=0;vex[2]=0;
   printf("\nInserire velocita' esterna  \n");
   if (tiposimm!=2)
   {
      printf("componente x : ");
      scanf("%f",vex);
   }
   else printf("componente x : 0\n");
   if (tiposimm!=3)
   {
      printf("componente y : ");
      scanf("%f",(vex+1));
   }
   else printf("componente y : 0\n");
   if (tiposimm!=1)
   {
      printf("componente z : ");
      scanf("%f",(vex+2));
   }
   else printf("componente z : 0\n");

   printf("Soluzione sistema finale :      ");

/* calcolo componenti normali della velocità esterna (vettore termini noti)*/

   for( n=0 ; n<nelem ; n++ )
      vnex(n)=-vex[0]*rotazione(n,0,2)-vex[1]*rotazione(n,1,2)-
               vex[2]*rotazione(n,2,2);

/* eliminazione gaussiana */

   for( n=0 ; n<(nelem-1) ; n++ )
   {
      printf("\b\b\b\b%3.0f%",((TIPOD)100*n/nelem));    /* contatore */

      pivot = minf(n,n);

/* triangolarizzazione array minf */

      for( k=(n+1) ; k<nelem ; k++ )
      {
         molt=minf(k,n)/pivot;
         minf(k,n)=0;
         for(j=(n+1);j<nelem;j++)
            minf(k,j) -= minf(n,j) * molt ;
         vnex(k) -= vnex(n) * molt ;
      }
   }

/* sostituzioni all'indietro */

   intens(nelem-1)=vnex(nelem-1)/minf(nelem-1,nelem-1);

   for( k=(nelem-2) ; k>(-1) ; k-- )
   {
      quot = vnex(k);
      for( j=(k+1) ; j<nelem ; j++ )
         quot -= minf(k,j) * intens(j) ;
      intens(k) = quot / minf(k,k) ;
   }

   printf("\b\b\b\b\bEseguita\n");

   return;
}

/*********************************CALVELNP*************************************/
/****************calcola velocità o pressioni nei null points******************/

void calvelnp(int s)
{
   void trasf(TIPOD *ass_coord,int elem_rif,int flag_trasla);
   void antitrasf(TIPOD *rel_coord,int elem_rif,int flag_trasla);
   int n,k,j;
   TIPOD vel[3],vel2[3];
   TIPOD modv;
   TIPOD *velf,*cp;
   TIPOD npr[3];
   char nomei[20]=" ";
   FILE *fp;

   velf = (TIPOD*)calloc((size_t)(nelem*3),sizeof(TIPOD));
   cp = (TIPOD*)calloc((size_t)(nelem),sizeof(TIPOD));

   for(k=0;k<(nelem*3);k++)
      *(velf+k)=0;

   printf("Calcolo delle velocità nei null point :       ");

   for(n=0;n<nelem;n++)
   {

       printf("\b\b\b\b%3.0f%",((TIPOD)100*n/nelem));    /* contatore */

       for(k=0;k<nelem;k++)
       {
       /* calcola coordinate relative a n del null point k-esimo */

          if (k==n)
          {
             vel[0]=0;
             vel[1]=0;
             vel[2]=6.2831853;
          }
          else
          {
             for(j=0;j<3;j++)
                npr[j]=np(k,j);
             trasf(npr,n,1);

             esatto(npr,n,vel);
          }

          antitrasf(vel,n,0);

          /* presenza di simmetria */

          if (tiposimm!=0)
          {
             if (tiposimm==2) npr[0]=-np(k,0);
               else npr[0]=np(k,0);
             if (tiposimm==3) npr[1]=-np(k,1);
               else npr[1]=np(k,1);
             if (tiposimm==1) npr[2]=-np(k,2);
               else npr[2]=np(k,2);

             trasf(npr,n,1);
             esatto(npr,n,vel2);
             antitrasf(vel2,n,0);

             if (tiposimm==2) vel[0]-=vel2[0];
               else vel[0]+=vel[0];
             if (tiposimm==3) vel[1]-=vel2[1];
               else vel[1]+=vel[1];
             if (tiposimm==1) vel[2]-=vel2[2];
               else vel[2]+=vel[2];
          }
          for(j=0;j<3;j++)
             *(velf+k*3+j)+=vel[j]*intens(n);
       }
   }

   printf("\b\b\b\b\bEseguito\n");

   for ( n=0;n<nelem;n++)
      for (k=0;k<3;k++)
         *(velf+n*3+k)+=vex[k];          /* aggiunge velocità esterna */

   if (s==3)
   {
      printf("Nome per il file velocità nei NP (max 8 char): ");
      scanf("%s",nomei);
      nomei[8]=0;

      printf("Salvataggio risultati nel file %s.dvn\n",nomei);

      if((fp=fopen( strcat(nomei,".dvn") , "wb" ))==NULL)
      { printf("errore apertura file");
          exit(1);
      }
      fwrite(&nelem,sizeof(int),1,fp);
      fwrite(velf,sizeof(TIPOD),(size_t)(nelem*3),fp);
   }

   if (s==9)
   {
      modv=sqrt(SQ(vex[0])+SQ(vex[1])+SQ(vex[2]));  /* modulo velocità esterna*/

      for ( n=0;n<nelem;n++ )
         *(cp+n)=1.0-SQ(sqrt(SQ(*(velf+n*3))
                            +SQ(*(velf+n*3+1))+SQ(*(velf+n*3+2)))/modv);

      printf("Nome per il file Cp nei NP (max 8 char): ");
      scanf("%s",nomei);
      nomei[8]=0;

      printf("Salvataggio risultati nel file %s.dcp\n",nomei);

      if((fp=fopen( strcat(nomei,".dcp") , "wb" ))==NULL)
      { printf("errore apertura file");
          exit(1);
      }
      fwrite(&nelem,sizeof(int),1,fp);
      fwrite(cp,sizeof(TIPOD),(size_t)(nelem),fp);
   }

   fclose(fp);

   return;
}

/***********************************SAVEDATI**********************************/
/* funzione di salvataggio della souzione (array intens) *
 * salva in formato formato binario :                    *
 * 4_byte_int : n° di elementi (int nelem)               *
 * successivi (nelem * 8_byte_TIPOD) : array intens     */

void savedati(void)
{
   FILE *fp;
   char nomei[20]=" ";

   printf("Nome per il salvataggio delle intensita' (max 8 char): ");
   scanf("%s",nomei);
   nomei[8]=0;

   printf("Salvataggio risultati nel file %s.don\n",nomei);

   if((fp=fopen( strcat(nomei,".don") , "wb" ))==NULL)
   { printf("errore apertura file");
       exit(1);
   }

   fwrite(&nelem,sizeof(int),1,fp);

   fwrite(intens,sizeof(TIPOD),(size_t)nelem,fp);
   fclose(fp);
   return;
}
/*****************************SAVEALTRI****************************************/
/**********************routine di salvataggio dati vari************************/
void savealtri(int r)
{
   FILE *fp;
   char nomei[20]=" ";

   printf("\nNome per il file (max 8 char): ");
   scanf("%s",nomei);
   nomei[8]=0;

   printf("Salvataggio risultati nel file %s",nomei);

   switch (r)
   {
     case 4 :
        printf(".dan\n");
        if((fp=fopen( strcat(nomei,".dan") , "wb" ))==NULL)
        { printf("errore apertura file");
            exit(1);
        }
        fwrite(&nelem,sizeof(int),1,fp);
        fwrite(punto,sizeof(TIPOD),(size_t)(nelem*12),fp);
        break;
     case 5 :
        printf(".dnu\n");
        if((fp=fopen( strcat(nomei,".dnu") , "wb" ))==NULL)
        { printf("errore apertura file");
            exit(1);
        }
        fwrite(&nelem,sizeof(int),1,fp);
        fwrite(np,sizeof(TIPOD),(size_t)(nelem*3),fp);
        break;
     case 6 :
        printf(".dmt\n");
        if((fp=fopen( strcat(nomei,".dmt") , "wb" ))==NULL)
        { printf("errore apertura file");
            exit(1);
        }
        fwrite(&nelem,sizeof(int),1,fp);
        fwrite(minf,sizeof(TIPOD),(size_t)(nelem*nelem),fp);
        break;
     case 7 :
        printf(".dtr\n");
        if((fp=fopen( strcat(nomei,".dtr") , "wb" ))==NULL)
        { printf("errore apertura file");
            exit(1);
        }
        fwrite(&nelem,sizeof(int),1,fp);
        fwrite(vnex,sizeof(TIPOD),(size_t)nelem,fp);
        break;
   }
   fclose(fp);
   return;
}