#include <fstream.h> 
#include <time.h> 
#include <math.h> 

int main(void) {
	time_t first, second;
	first = time(NULL);  /* Gets system time */
	ifstream f0("file.IN");
	ofstream f1("junk.dat");
	if (!f0) cerr << "Cannot open infile for input";
	if (!f1) cerr << "Cannot open outfile1 for output"; 
//***************** Euler Integration ************************ 
// this section for calculating a trajectory of a classical rubidium atom
leaving a MOT after 
// undergoing some inelastic process giving it some initial velocity. this
1-d cooling model 
// calculates the trap depth for a given zmax, the trap radius, a given
detuning and intensity 
// of the trap laser, etc. 
// using the Euler formula for integration
	double delta, I, detuning, depth;
	double a, v, z, eta, d, n, q, vinit, t, x, y;
	long i, j, k; 
//*********************INPUT PARAMS HERE*******************************
	const double M = 87/(6.022e23);
	const double isat = 1.62;
	const double zmax = 0.3;
	const double bgrad = 4.8;
        const double tstep = 1e-6;
	k = 1;
//********************************************************************
	while (k <= 4) {
	  cout << k << "\n";
          j = -30;
	  vinit = 25;
	  while (j <= 20)  {
		I = isat*exp(j * log(10)/10);
		detuning = -k;
                z = 0;
		while (z-zmax < -0.01 || z-zmax > 0.01) {
			v = vinit;
			z = 0;
			while (v > 0) { // do the integration (1st-order)
				z = z + v*tstep;
				eta = 2.175e-3 * v + 0.2374 * bgrad * z;
				delta = detuning*detuning + 0.25 + 
0.125*I/isat;
				q = (delta - eta*eta)*(delta - eta*eta) + 
eta*eta;
				d = (1 + 4*eta*eta)*(q +
0.125*I/isat*(4*delta + 3*I/(8*isat) + 4*eta*eta)) +
I*I/(16*isat*isat)*(delta + 3*I/(8*isat) - 3*eta*eta);
				n = detuning * eta*(1 + 4*eta*eta)*I/isat;
				v = v + 1.093e7*n/d*tstep;
			}
			// This line iteratively finds the vinit that
			// launches an atom from the trap which comes to 
			// rest at point z where 0.99*zmax < z < 1.01 zmax
			// this is very crude - if the routine does not
			// converge for your parameters try changing this
			// line.
			vinit = vinit*(1 - 0.1*(z-zmax)/zmax); // cout <<
vinit << " " << z << "\n";
		}
		depth = M*vinit*vinit/2;
        	cout << "depth =" << depth << " I=" << I << " detuning="
<< detuning << "zfinal=" << z << "\n";
		f1 << detuning << " " << I << " " << depth << "\n";
		++j;
          }
          ++k;
        }
	// this section for reporting on the total runtime
	second = time(NULL); /* Gets system time again */
	cout << "runtime: " << difftime(second,first) << "\n";
	return(0); }