//Libraries reading
# include <stdlib.h>
# include <stdio.h>
# include <math.h>
# include <complex.h>
# include <time.h>
#include <ctype.h>
/*
**********************************************************************************************************
* THIS PROGRAM SOLVES NUMERICALLY THE BURGERS EQUATION FOR TEMPERATURE-INHOMOGENEOUS FLUID AS HOMOGENEOUS*                                          *
**********************************************************************************************************
*/

/* Declaration of global variables and functions */
int M;
double st,alfa,beta,cA,cB,p0;
double w,g,R,TA,rA,L,Q,G,N,A;
double s,n_steps,ss;
double complex Pn[201]; /* Field for complex parts of the pressure - set by the norm*/
double complex TPn[201]; /* Temporary field for complex parts of the pressure - set by the norm*/

double abso(double q); /* Absolute value function */
double T0(double x); /* Temperature function */
double dT0(double x); /* Derivation of temperature function */
double c0(double x); /* Sound velocity function */
double r0(double x); /* Density function */
/* Complex right-hand side of differential equation- set by the norm*/
double complex Fn(int n, double x);



const float pi = 3.1415926537;

int main()
{ /* begin */
/* Declaration of local variables */
FILE *cR, *cI, *Data;
int i;
double complex k1[201];
double complex k2[201];
double complex k3[201];
double complex k4[201];

/* -------------------------------------------------------------------------*/
/* --------------------Program body starts here--------------------------- */
/* -------------------------------------------------------------------------*/
/* Set ZEROs to the fields */
for (i = 0; i <= 200; i++)
{
        Pn[i] = 0.+0.*I;
        TPn[i] = 0.+0.*I;
}

/* Setting of input parameters */
p0 = 101325;
R = 287.058;
g = 1.402;
L = 1.;
TA = 298.15;//357.78;
cA = 346.4;
rA = 1.193;
alfa = 0.00004578;
beta = 1.201;
w = 10000;
A = 0.5;

/* Boundary condition */
Pn[1] = -(1000/(rA*cA*cA))*I;
//Pn[1] = -0.5*I;

s = 0.; /* Setting of dimensionless distance */

//Setting of input parameters
printf("Set the whole distance ss=");
scanf("%lf", &ss);
printf("Set the whole number of Fourier series members (max. 200) M=");
scanf("%d", &M);
printf("Set the number of steps=");
scanf("%lf", &n_steps);
/* the distance of the numerical step */
st = ss/n_steps;

//Linear dependence of temperature on the distance: T0=TA*(1+a*x)

cB = c0(L);
//Parameters calculation
Q = w*L/(2*cA);
N = beta*w*L/(cA);
G = alfa*w*w*L/(2*rA*cA*cA*cA);


/* Pn = cR/2 - j*cI/2 */
Data = fopen("Data.txt", "wt"); /* Information about resonator and the calculation is been writing here */
fprintf(Data, "The number of Fourier series members M= %d \n",M);
fprintf(Data, "The number of steps= %lf \n",n_steps);
fprintf(Data, "step= %lf \n",st);
fprintf(Data, "The whole distance ss=%.8lf \n",ss);
fprintf(Data, "Q=%.8lf \n",Q);
fprintf(Data, "G=%.10lf \n",G);
fprintf(Data, "N=%.8lf \n",N);
fprintf(Data, "A=%.8lf \n",A);
fprintf(Data, "cA=%.8lf \n",cA);
fprintf(Data, "cB=%.8lf \n",cB);
fclose(Data);

//Eigenvalue calculation starts here
   for (i = 1; i<=M; i++)
	{
      TPn[i] = Pn[i];
	}

while (s<=ss)
{ /*begin while */

//RK4 method
	for (i = 1; i<=M; i++)
	{
      k1[i] = Fn(i,s);
	}
    for (i = 1; i<=M; i++)
	{
      TPn[i] = Pn[i]+k1[i]/2.;
	}
    for (i = 1; i<=M; i++)
	{
      k2[i] = Fn(i,s+st/2.);
	}
    for (i = 1; i<=M; i++)
	{
      TPn[i] = Pn[i]+k2[i]/2.;
	}
	for (i = 1; i<=M; i++)
	{
      k3[i] = Fn(i,s+st/2.);
	}
    for (i = 1; i<=M; i++)
	{
      TPn[i] = Pn[i]+k3[i];
	}
	for (i = 1; i<=M; i++)
	{
      k4[i] = Fn(i,s+st);
	}
    for (i = 1; i<=M; i++)
	{
      TPn[i] = Pn[i]+(1./6.)*(k1[i]+2.*k2[i]+2.*k3[i]+k4[i])*sin(i/25.)/(i/25.);
	}

for (i = 1; i<=M; i++)
	{
        Pn[i] = TPn[i];
	}
printf("  s= %lf \n",s);//writing the distance to the display

s = s+st; //distance incrementation
} /* end while */

cR = fopen("cR.txt", "wt"); /* Real parts of the pressure harm. components are been writing here*/
cI = fopen("cI.txt", "wt"); /* Imaginary parts of the pressure harm. components are been writing here */
for (i = 1; i<=M; i++)
if (i==M)
	{
		fprintf(cR, "%.8e \n",creal(Pn[i]));
		fprintf(cI, "%.8e \n",cimag(Pn[i]));
	}
else
	{
        fprintf(cR, "%.8e \t",creal(Pn[i]));
		fprintf(cI, "%.8e \t",cimag(Pn[i]));
	}
fclose(cR);/* finish of the text file for the real Fourier series members */
fclose(cI);/* finish of the text file for the imaginary Fourier series members */
} /* end */

/* Functions definitions */

double T0(double x)
{
    double TT;

    TT = TA*(1+A*x);
    return TT;
}

double dT0(double x)
{
    double ddT;

    ddT = TA*A*x;
    return ddT;
}

double c0(double x)
{
    double cc;
    cc = sqrt(g*R*T0(x));
    return cc;
}

double r0(double x)
{
    double rr;
    rr = p0/(R*T0(x));
    return rr;
}

double complex Fn(int n, double x)
{
int m;
double complex konvoluce,F;

konvoluce = 0.+0.*I;

    for (m = 1; m < n ; m++)
        konvoluce += TPn[m]*TPn[n-m];

    for (m = (n+1); m <= M; m++)
        konvoluce += 2*(TPn[m]*conj(TPn[m-n]));

//   konvoluce = 0.+0.*I; //Solve linear equation by comment out

        F = I*(n*N*TA/(2.*T0(x)))*konvoluce
        -(n*n*G*TA/T0(x))*TPn[n]
        -dT0(x)*TPn[n]/(2*T0(x))
        -I*n*Q*((TA-T0(x))/T0(x))*TPn[n];

        F = st*F;
return F;
}


double abso(double q)
{
double ab;
ab = q;
if (ab < 0)
ab = -ab;
return ab;
}