#include <stdio.h>
#include <stdlib.h>
 
int main()
{
printf("********************************\n");
printf("********************************\n");
printf("*                              *\n");
printf("*           LW Solver          *\n");
printf("*                              *\n");
printf("********************************\n");
 
printf("This program solves the hyperbolic EDP of a bar vibration\n");
printf("The dimensions of the bar are L=1 m, D= 0,3 m. The material of the\n");
printf("bar is steel, with E=200 GPa and rho=7850 kg/m^3. dx=0.05 m\n");
printf("dt=0.05 s. The solution will be write in a .dat file\n\n");
 
/* Variables */
 
int i,j;
double u[20][20];
 
FILE *fp;
 
/* Escritura del fichero solucion */
fp = fopen("Solution_LW-Solver_005.dat","w");
 
/*----------------------------------------------------------------------------*/
 
fprintf(fp,"********************************\n");
fprintf(fp,"********************************\n");
fprintf(fp,"*                              *\n");
fprintf(fp,"*           LW Solver          *\n");
fprintf(fp,"*                              *\n");
fprintf(fp,"********************************\n");
 
fprintf(fp,"This program solves the hyperbolic EDP of a bar vibration\n");
fprintf(fp,"The dimensions of the bar are L=1 m, D= 0,3 m. The material of\n");
fprintf(fp,"bar is steel, with E=200 GPa and rho=7850 kg/m^3. dx=0.05 m,\n");
fprintf(fp,"dt=0.05 s.\n\n");
 
 
/*----------------------------------------------------------------------------*/
 
printf("Solving...\n");
fprintf(fp, "Solution for dt=0.05m\n\n");
fprintf(fp,"*************************************************************\n\n");
fprintf(fp,"Calculated values are:\n\n");
/* Condiciones de contorno e iniciales */
 
/* Condicion de desplazamiento nulo en el origen */
for (j=0; j<=20; j++) {
    u[0][j]=0;
}
printf("\n\n1....");
/* Condicion de desplazamiento inicial nulo */
for (i=0; i<=20; i++) {
    u[i][0]=0;
}
printf("\n\n2....");
/* Ecuaciones iniciales */
 
/* Ecuacion inicial para todos los nodos salvo para el primero y el ultimo */
for (i=1; i <=19; i++) {
    u[i][1]=1;
}
printf("\n\n3....");
/* Ecuacion del nodo del extremo */
u[20][1]=0;
printf("\n\n4....");
/* Resolucion iterativa */
/* Step 1 par; Step 2 impar */
 
for (j=2; j<=20; j++) {
    if (j%2 == 0)
    {
        for (i=1; i<=19; i++) {
          u[i][j]=0.5*(u[i+1][j-1]+u[i-1][j-1])-0.08*(u[i+1][j-1]-u[i-1][j-1]);
        }
            u[20][j]=u[19][j-1];        
    }
    else if (j%2 !=0) 
    {
        for (i=1; i<=19; i++) {
          u[i][j]=u[i][j-2]-0.16*(u[i+1][j-1]-u[i-1][j-1]);
        }
          u[20][j]=u[20][j-2];
    } 
}
 
printf("\n\n5....");
 
for (i=0; i<=20; i++) {
    for (j=0; j<=20; j++) {
        fprintf(fp,"U%d%d= %12.4e C\n", i, j, u[i][j]);
    }
}
 
printf("\n\nSolved!");
 
fprintf(fp,"\n\nEND DATA");
fclose(fp);
 
  printf("\n\nPress enter to continue ...");
  getchar();
  return 0;
}