Are my implicit type conversions correct?

I'm having some issues with my code and I'm trying to rule out any syntax or programming mistakes. One of the things I'm trying to rule out is the different variable types used in my equations. I'm mostly operating with doubles, but I did try to ensure that all the integers in the equations were correct. Thank you in advance.

#include <iostream>
#include <math.h>
#include <stdio.h>
#include <fstream>
#include <cmath>
#include <algorithm>

#define pi 3.14159265

using namespace std;

//Define Fluid Properties
double rho_L = 998; //Liquid Density
double rho_LG = 828.9; //Liquid-Gas Density Ratio
double mu_L = 0.000798; //Liquid Viscosity
double mu_LG = 40.24; //Liquid-Gas Viscosity Ratio
double sigma = 0.0712; //Surface Tension
double nu_G = (mu_L/mu_LG)/(rho_L/rho_LG);

//Define Injector Properties
double Uinj = 56.7; //Injection Velocity
double Dinj = 0.0998; //Injector Diameter
double theta = 45.0*pi/180.0; //Spray Cone Angle
double L = 100.0*Dinj; //Atomization Length
double Ls = L/Dinj; //Normalized Atomization Length

//Define Solver Parameters
double K = 26470; //Viscous Dissipation Coefficient
double Eps = pow(10,-5); //Residual Error
double dx = 0.01; //Step Size
double Ui = 15.0; //Initial Guess
//const int Z = static_cast<const int>((Ls-1)/dx + 1) + 1;
const int N = 9901;//Z;

double deriv (double U, double X, double delta, double m)
{
double dudx;
dudx = -(1.0/delta)*(1.0/U)*(U - sqrt(1.0 - U)/sqrt(m*X*X))*(U - sqrt(1.0 - U)/sqrt(m*X*X));
return (dudx);
}

int main()
{
//Declare Variables
int max_step;
double ERR;
int step;
double DEN;
double SMD;
double m;
double Ug;
double Re;
double Cd;
double delta;
double K1;
double K2;
double K3;
double K4;

//Allocate Memory From Heap
double *U = new double [N];
double *X = new double [N];

//Initialize Vectors and Variables
DEN = 0.5*rho_L - (4.0/3.0)*K*(mu_L)/(Uinj*Dinj*Dinj)*L;
m = 4.0/rho_LG*tan(theta)*tan(theta);

for (int i = 0; i < N; i++)
{
X[i] = 1.0 + dx*i;
}
U[0] = 1.0;

max_step = 10000;
ERR = 1;
step = 0;
while(abs(ERR) > Eps && step < max_step)
{

//Calculate Ug
Ug = sqrt(1.0 - (Ui/Uinj))/sqrt(m*Ls*Ls)*Uinj;

//Calculate SMD
SMD = 6.0*sigma/(DEN*(Uinj*Uinj - Ui*Ui));

//Calculate Re # and Drag Coefficient
Re = abs(Ui-Ug)*SMD/nu_G;

if(Re <= 0.01)
{
Cd = (0.1875) + (24.0/Re);
}
else if(Re > 0.01 && Re <= 260.0)
{
Cd = (24.0/Re)*(1.0 + 0.1315*pow(Re,(0.32 - 0.05*log10(Re))));
}
else
{
Cd = (24.0/Re)*(1.0 + 0.1935*pow(Re,0.6305));
}

//Determine New U
delta = (4.0/3.0)*(1.0/Cd)*(rho_LG)*(SMD/Dinj);

//RK4
for (int i = 0; i < N-1; i++)
{
K1 = deriv(U[i],X[i],delta,m);
K2 = deriv(U[i]+0.5*dx*K1,X[i]+0.5*dx,delta,m);
K3 = deriv(U[i]+0.5*dx*K2,X[i]+0.5*dx,delta,m);
K4 = deriv(U[i]+dx*K3,X[i+1],delta,m);
U[i+1] = U[i] + dx/6.0*(K1 + 2.0*K2 + 2.0*K3 + K4);
}

ERR = abs(U[N-1]*Uinj - Ui)/Ui;

Ui = U[N-1]*Uinj;

step = step + 1;
}

SMD = 6.0*sigma/(DEN*(Uinj*Uinj - Ui*Ui));

/*cout << "U = " << Ui << endl;
cout << "SMD = " << SMD << endl;
cout << "DEN = " << DEN << endl;
cout << "Ug = " << Ug << endl;
cout << "m = " << m << endl;
cout << "delta = " << delta << endl;
cout << "Re = " << Re << endl;
cout << "Cd = " << Cd << endl;
cout << "U* = " << U[N-1] << endl;
cout << "Error = " << ERR << endl;
cout << "step = " << step << endl;*/

//Output Data Into Text File
ofstream outputdata("result.txt");
for (int i = 0; i < N; i++)
{
outputdata << X[i] << " " << U[i] << '\n';
}
outputdata.close();

delete [] U;
delete [] X;

return 0;
}

-
This is likely off-topic as per the FAQ - 4. To the best of my knowledge, does the code work? –  Yuushi Feb 13 '13 at 1:43
The code is functioning, yes. The results are slightly off of what I expected, so I was hoping to localize the inaccuracies. –  user1562781 Feb 14 '13 at 3:41

There are no implicit type conversions here: the only variables that are not double are the step and max_step and they are not used with the other ones. Otherwise, you're using doubles everywhere so there is no issue.
Slight inaccuracies are possible when doing a lot of iterations, the most likely reason is that floating point numbers are simply not fully accurate. If you're seeing a difference when comparing this to another program, it's possible that the order of operations is the issue, too. On a computer, a * (b * c) is slightly different from (a * b) * c.
You probably meant a * (b * c) is different from (a * b) * c . –  Sylvain Defresne Feb 16 '13 at 20:54
@user1562781 Yes, indeed, 4/3 == 1, while 4.0/3.0 is very close to 1.25, which is what you want. Your question only contains doubles, either because the variable is double (double K = 25670), or because it is written as a double (4.0). –  Quentin Pradet Feb 18 '13 at 6:43