/*
 *  idc3d.c
 *  
 *
 *  Created by Pierre CHAINAIS on Thu Mar 10 2005.
 *  Copyright (c) 2005 __MyCompanyName__. All rights reserved.
 *
 */

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

#define PI 3.14159265

typedef struct
{
	double xi;
   double yi;
   double zi;
   double ri;
}Boule;



int main(int argc,char *argv[])
{
   if (argc==8 || argc==9 || argc ==10)
   {

      double Xmax,Ymax,Zmax,res,deltat,N2,C,mu,sigma2,temp,K;
      int nbpasX,nbpasY,nbpasZ,N;
 /*     double law; */
      int i,j,k,l;
      double umax,H1;
      double *u,*W,*tmp,***Qr,*xt,*yt,*zt;
      Boule *b;
      time_t t;
      FILE *f2;

      int **BX,**BY,**BZ;
      int law;

      int Arrondi(double x);    /*Arrondi un reel à l'entier le plus proches*/
      int TirageAleatoire(int n); /*Tirage uniforme entre 0 et n-1*/
      double NormalBoxMuller(double mean,double s); /*Tirage loi normale*/



/*Affectation des arguments de la ligne de commandes*/
      Xmax=(float)atof(argv[1]);
      Ymax=(float)atof(argv[2]);
      Zmax=(float)atof(argv[3]);
      res=(float)atof(argv[4]);
      deltat=(float)atof(argv[5]);
      law=(int)atoi(argv[6]);
      printf ("argc=%d\n",argc);

/* choix aleatoire d'un germe pour la fonction rand avec l'heure calendaire */      
      srand((unsigned) time(&t));

       switch(argc)
      {
          case 8 : switch(law)
                        {
                           case 2 : printf("Poisson composee log-exponentielle\n");
                                         temp=(float)atof(argv[7]);
                                         C=1;break;
                           case 3 : printf("Constante\n");
                                         K=(float)atof(argv[7]);
                                         C=1;break;
                           default : printf("Erreur\n");exit(0);
                        }
                        break;
          case 9 : switch (law)
          					{
                           case 1 : printf("Poisson composee Log Normale\n");
                                        mu=(float)atof(argv[7]);
                                        sigma2=(float)atof(argv[8]);
                                        C=1;
                                        break;
                           case 2 : printf("Poisson composee log-exponentielle\n");
                                        temp=(float)atof(argv[7]);
                                        C=(float)atof(argv[8]);
                                        break;
                           case 3 : printf("Constante\n");
                                        K=(float)atof(argv[7]);
                                        C=(float)atof(argv[8]);
                                        break;
                           default : printf("Erreur\n");exit(0);
                        }
                        break;
          case 10 : switch(law)
                        {
                           case 1 : printf("Poisson composee Log Normale\n");
                                        mu=(float)atof(argv[7]);
                                        sigma2=(float)atof(argv[8]);
                                        C=(float)atof(argv[9]);
                                        break;
                           default : printf("Erreur\n");exit(0);
                        }
                        break;
          default :   printf("usage : nom_executable <Xmax> <Ymax> <Zmax> <resolution> <pas_de_espace> <choix_loi> <param_loi1> (<param_loi2>) (<constante:densifier tirages> sinon constante=1) \n"); exit(0);
      }

/*Calcul du nb de pas en fonction de Xmax,Ymax,Zmax et de deltat*/
      nbpasX=Arrondi(Xmax/deltat);    /* arrondi a l'entier le plus proche */
      nbpasY=Arrondi(Ymax/deltat);
      nbpasZ=Arrondi(Zmax/deltat);

/*Definition umax*/
      umax=1/(res*res*res)-1;

/*Calcul N*/
      N2=((double)C*2*umax*(Xmax+1)*(Ymax+1)*(Zmax+1))/PI;
      printf("N2=%f\n",N2);
      N=Arrondi(N2);
      printf("N=%d\n",N);



   b=(Boule *)malloc(N*sizeof(Boule));
   u=(double *)malloc(N*sizeof(double));
   for (i=0;i<N;i++)
   {
      b[i].xi=-0.5+rand()/(double)RAND_MAX*(Xmax+1);         /* tirage aleatoire des xi dans [-0.5,Xmax+1] */
      b[i].yi=-0.5+rand()/(double)RAND_MAX*(Ymax+1);         /* tirage aleatoire des yi dans [-0.5,Ymax+1] */
      b[i].zi=-0.5+rand()/(double)RAND_MAX*(Zmax+1);         /* tirage aleatoire des zi dans [-0.5,Zmax+1] */
      u[i]=rand()/(double)RAND_MAX*umax;                         /* tirage aleatoire des ri */
      b[i].ri=pow((double)(1/(res*res*res)-u[i]),-(double)1/3);            /* calcul des ri */
   }

   W=(double *)malloc(N*sizeof(double));
   tmp=(double *)malloc(N*sizeof(double));
   switch(law)
   {
          case 1 :  for (i=0;i<N;i++)
                         {
                             tmp[i]=NormalBoxMuller(mu,sigma2);
                             W[i]=exp(tmp[i]);
                        }
                        H1=C*(1-exp(mu+sigma2/2));
                        break;

          case 2 : for (i=0;i<N;i++)
                        {
                            tmp[i]=rand()/(double)RAND_MAX;
                            W[i]=pow(pow(1+temp,(double)1/temp)*tmp[i],temp);
                        }
                        H1=0;
                        break;

          case 3 : for (i=0;i<N;i++)   W[i]=(double)K;
                        H1=C*(1-K);
                        break;

          default : printf("Erreur lors de l'initialisation de W\n");
   }

   free(tmp);

   Qr=(double ***)malloc(nbpasX*sizeof(double));
   for (i=0;i<nbpasX;i++)
   {
      Qr[i]=(double **)malloc(nbpasY*sizeof(double));
      for (j=0;j<nbpasY;j++)
         Qr[i][j]=(double *)malloc(nbpasZ*sizeof(double));
   }

/* Initialisation Qr, x, y, z */
   xt=(double *)malloc(nbpasX*sizeof(double));
   yt=(double *)malloc(nbpasY*sizeof(double));
   zt=(double *)malloc(nbpasZ*sizeof(double));
   for (i=0;i<nbpasX;i++)  xt[i]=i*deltat;
   for (i=0;i<nbpasY;i++)  yt[i]=i*deltat;
   for (i=0;i<nbpasZ;i++)  zt[i]=i*deltat;

   for(i=0;i<nbpasX;i++){
     for (j=0;j<nbpasY;j++){
       for (k=0;k<nbpasZ;k++){
            Qr[i][j][k]=1;
       }
     }
   }


  BX=(int **)malloc(N*sizeof(int));
  BY=(int **)malloc(N*sizeof(int));
  BZ=(int **)malloc(N*sizeof(int));
  for(i=0;i<N;i++)
  {
     BX[i]=(int *)malloc(nbpasX*sizeof(int));
     BY[i]=(int *)malloc(nbpasY*sizeof(int));
     BZ[i]=(int *)malloc(nbpasZ*sizeof(int));
  }


   for (l=0;l<N;l++)
   {
     for(i=0;i<nbpasX;i++)
  	 if (fabs(xt[i]-b[l].xi)<=b[l].ri/2)
  	 {
  	BX[l][i]=1;
  	 }
     for(j=0;j<nbpasY;j++)
            if (fabs(yt[j]-b[l].yi)<=b[l].ri/2)
  	 {
  	  BY[l][j]=1;
  	 }
     for(k=0;k<nbpasZ;k++)
            if (fabs(zt[k]-b[l].zi)<=b[l].ri/2)
  	 {
  	  BZ[l][k]=1;
  	 }

   }

     /* calcul des Qr  */
     for (l=0;l<N;l++)
      if (b[l].xi<=(Xmax+b[l].ri/2)&& b[l].yi<=(Ymax+b[l].ri/2)&& b[l].zi<=(Zmax+b[l].ri/2))
       {
        for(i=0;i<nbpasX;i++)
         if (BX[l][i]!=0)
         {
          for (j=0;j<nbpasY;j++)
           if (BY[l][j]!=0)
           {
             for (k=0;k<nbpasZ;k++)
  	           if (BZ[l][k]!=0)
                {
  	               if ((xt[i]-b[l].xi)*(xt[i]-b[l].xi)+ (yt[j]-b[l].yi)*(yt[j]-b[l].yi) + (zt[k]-b[l].zi)*(zt[k]-b[l].zi)<=(b[l].ri/2)*(b[l].ri/2))
                    {
                       Qr[i][j][k]*=W[l];
                    }
                }
           }
  	    }
       }


    /* normalisation de Qr   */
    /* et inscription dans un fichier   */
    f2=fopen("sortie.txt","wt");
    for(i=0;i<nbpasX;i++)
     for (j=0;j<nbpasY;j++)
       for (k=0;k<nbpasZ;k++)
       {
           Qr[i][j][k]=Qr[i][j][k]*pow(res,-H1);
           fprintf(f2,"%f	%f	%f	%f\n",xt[i],yt[j],zt[k],Qr[i][j][k]);
       }
    fclose(f2);

   free(u);

   free(W);

  }
  else
  {
            printf("usage : nom_executable <Xmax> <Ymax> <Zmax> <resolution> <pas_de_espace> <choix_loi> <param_loi1> (<param_loi2>) (<constante:densifier tirages> sinon constante=1) \n"); exit(0);
  }
  return 0;
}
/**********************Fin du main*************************************************/

/*********************************************************************************/
/**********************Fonctions**************************************************/
/*********************************************************************************/

int Arrondi(double x)
{
	return (((x-(int)x)<=0.5) ? (int)x : (int)x+1);
}


/***********************************
	Tirage aléatoire d'un entier entre 0 et n-1
***********************************/
int TirageAleatoire(int n)
{
	return ((n==0) ? 0 : rand()%n);
}

/***************************************************************************************
 	Simulation d'une loi normale par la methode de Box Muller : Si
    x1 et x2 suivent une loi uniforme dans]0,1[, y1=sqrt(-2ln(x1))cos(2M_PI.x2)
    et y2=sqrt(-2ln(x1))sin(2M_PI.x2) suivent une loi normale réduite. On en déduit que
    m+sy1 et m+sy2suivent une loi N(m,s)
****************************************************************************************/
double NormalBoxMuller(double mean,double s)
{
	 double 	x1;
    long    pas;

    /* tirage d'un reel uniforme dans [0,1]  */
    pas = 100;
    x1 = (double)TirageAleatoire(pas)/(double)pas;
    x1 = (x1 <1E-6 ? 1E-6 : x1);
    x1 = (x1 ==1 ? 1-1E-6 : x1);

    return (mean + s*sqrt(-2*log(x1))*cos(2*PI*x1));
}




