//--------------------------------------
//| Main program |
//--------------------------------------
main()
{
int num_points = 0, rho_initial_x_i;
LD x[max_points], rho[max_points], deltax, initial_condition,
kn1, kn2, kn3, kn4;
vector S[max_points], F[max_points], P[max_points], U[max_points],
temp_s, temp_f, temp_u, temp_dudx;
char flag, junk[10];
file_input (num_points, deltax, initial_condition, x, F, S, U);
rho_initial_x_i = (num_points - 1) / 2;
rho[rho_initial_x_i] = initial_condition;
P[rho_initial_x_i].a = 0;
P[rho_initial_x_i].b = initial_condition;
for (int i = rho_initial_x_i + 1; i < num_points; i++)
{
kn1 = drhodx(rho[i-1], S[i-1], U[i-1], dudx(U, i, deltax));
temp_s = multiply (0.5, add(S[i],S[i-1]));
temp_f = multiply (0.5, add(F[i],F[i-1]));
temp_u = multiply (-1/normal(temp_f), temp_f);
temp_dudx = dudx(U,i, deltax);
kn2 = drhodx(rho[i-1] + deltax * kn1 / 2, temp_s, temp_u, temp_dudx);
kn3 = drhodx(rho[i-1] + deltax * kn2 / 2, temp_s, temp_u, temp_dudx);
kn4 = drhodx(rho[i-1] + deltax * kn3, S[i], U[i], dudx(U, i, deltax));
rho[i] = rho[i-1] + deltax / 6 * (kn1 + 2 * kn2 + 2 * kn3 + kn4);
P[i] = multiply(rho[i], U[i]);
}
for (i = rho_initial_x_i - 1; i >=0; i--)
{
kn1 = drhodx(rho[i+1], S[i+1], U[i+1], dudx(U, i, deltax));
temp_s = multiply (0.5, add(S[i],S[i+1]));
temp_f = multiply (0.5, add(F[i],F[i+1]));
temp_u = multiply (-1/normal(temp_f), temp_f);
temp_dudx = dudx(U,i, deltax);
kn2 = drhodx(rho[i+1] + deltax * kn1 / 2, temp_s, temp_u, temp_dudx);
kn3 = drhodx(rho[i+1] + deltax * kn2 / 2, temp_s, temp_u, temp_dudx);
kn4 = drhodx(rho[i+1] + deltax * kn3, S[i], U[i], dudx(U, i, deltax));
rho[i] = rho[i+1] - deltax / 6 * (kn1 + 2 * kn2 + 2 * kn3 + kn4);
P[i] = multiply(rho[i], U[i]);
}
cout << "\n\n\n\n\n\nDo you want file output? [Y/N] ";
cin.get(flag);
cin.getline(junk,10);
if ((flag == 'Y') || (flag == 'y'))
file_output(P, num_points);
return (0); //end main
}
