# include <iostream.h>
# include <math.h>
# include "complex.h"
# include <conio.h>
# include <stdio.h>
# include "stdlib.h"

long double c0, bt, fr,alpha, Z0, I0, u0, a, k0, r0, zmax, rmax;
long double pi= atan(1);
FILE *fic1,*fic2;
long double ma, ld;
long double dr, dz;
complex *ui,*uf;

void init (int iz, int jr)
{
   ///////////////////////////////////////////////////////////////
   ///////////////  Initialisation des constantes  ///////////////
   ///////////////////////////////////////////////////////////////
   
   //Définition des constantes
   
   printf("Initialisation des constanes en cours\n");
   c0=1480;                            // La célérité du son dans le milieu
   bt=3.5;                             // Le paramètre de la nonlinéarité
   fr=2.2e6;                           // La fréquence d'émission
   alpha=0.025;                        // L'absorption acoustique du milieu
   Z0=1.48e6;                          // L'impédance acoustique du milieu
   //////////////////////////////////////////////////////////////////////////
   I0=1;                               // La puissance à la source en w/cm²
   u0=sqrt(2*I0*1e4/Z0);               // La velocité des particules à la surface de la source
   a=5e-3;                             // Le rayon de la source
   k0=2*pi*fr/c0;                      // Le nombre d'onde
   r0=k0*a*a/2;                        // La distance de Rayleigh
   //////////////////////////////////////////////////////////////////////////
   zmax=2.5*r0;                        // La distance maximal sur la direction axiale
   dz=zmax/iz;                         // Le pas d'échantillonnage sur la direction z
   rmax=10*a;                          // Le rayon maximal sur la direction radiale
   dr=2*rmax/jr;                       // Le pas d'échantillonnage sur la direction r
   ma=u0/c0;                           // Le nombre de Mach
   ld=c0/(bt*ma*fr);                   // La distance de shock
   
   if ((dz/(dr*dr)) >= (r0/(2*a*a)))   // Test pour la divergence
      cout << "ERREUR: Risque de divergence\n";

   printf("Initialisation des constante terminee avec succes\n");
}

void source(int jr, int N)
{
	  
   ///////////////////////////////////////////////////////////////
   ///////////////     Definition de la source     ///////////////
   /////////////////////////////////////////////////////////////// 
  int j,n,k;
  long double r;

  printf("Definition de la source\n");
  for (n=1;n<N+1;n++)
  {
	  for (j=1;j<jr+2;j++)
	  {
		  r=-rmax+(j-1)*dr;
                  k=(n-1)*jr+j;
		  if (n==1)
		  {
			 if ((r>-a-(dr/2))&&(r<a+(dr/2)))
				  ui[k] =complex(u0);
			 else
				  ui[k] =complex(0.0);
			 fprintf(fic1,"%lf\n",norm(ui[k]));
		  }
		  else
		  {
			 ui[k] =complex(0.0);
			 if (n==2)
				fprintf(fic2,"%lf\n",norm(ui[k]));
		  }
	   }
  }
  printf("Definition de la source est terminee avec succes\n");
}
void main ()
{
	
	int i,j;
	long double z, r;

	int iz=1000;
	int jr=200;
	int Nh=5;

	//clrscr();
	// Initialisation
	init (iz,jr);

	fic1= fopen("z.txt","w");
	fic2= fopen("r.txt","w");

	for (i=1;i<iz+1;i++)
	{
		z=(i-1)*dz;
		fprintf(fic1,"%lf\n",z);
	}

	for (j=1;j<jr+2;j++)
	{
		r=-rmax+(j-1)*dr;
		fprintf(fic2,"%lf\n",r);
	}

	fclose (fic1);
	fclose (fic2);

	fic1= fopen("harm1.txt","w");
	fic2= fopen("harm2.txt","w");
	
	source(jr,Nh); // Definition de la source
    getch();
}