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

int main(void) {

	time_t first, second;
	first = time(NULL);  /* Gets system time */

	char ch;

	//ifstream f0("FILE.IN");
	ofstream f1("2u_dt.dat");
	ofstream f2("2u_de.dat");
	ofstream f3("2u_dr.dat");
	ofstream f4("2u_tex.dat");
	ofstream f5("2u_db.dat");

	//if (!f0) cerr << "Cannot open infile for input";
 	if (!f1) cerr << "Cannot open outfile1 for output";
	if (!f2) cerr << "Cannot open outfile2 for output";
	if (!f3) cerr << "Cannot open outfile3 for output";
	if (!f4) cerr << "Cannot open outfile4 for output";
	if (!f5) cerr << "Cannot open outfile5 for output";

// ***************** Euler Integration ************************
// this section for calculating a trajectory of a classical rubidium collision on the
// 5s+5p potential starting at the Trap condon radius of 141nm with velocity vinit and
// impact parameter binit.
// using the Euler formula for integration
//
//*************** constants in MKS units *******************************
	const double PI = 3.1415927;
	const double HBAR = 1.054572E-34; 			//Planck constant 
	const double NA = 6.0221367E23;				//Avogadro's number
	const double K = 1.380658E-23;				//Boltzmann Constant
	const double M = 85e-3/NA;
	const double D2 = .6525072417e-47;			//matrix element for Rb (J-V paper)
	const double GA = 2*PI*5.89E6;				//atomic decay rate
        const double TAU = 1/GA;                //molecular lifetime
	const double C3 = D2;
	const double A0 = 0, A2 = 0.2, A4 = -0.0107429;
	double energy,c2,temp,integral;
	double t,v,vi,b,bi;
	double r,l,d,lz,gm,pop;
	double Pv,Pb,Pt,delay,dP,dPred,lztrap,lzprobe;
	long x,pts;
//************************** INPUT PARAMS HERE in MKS units ************
	const double DETUNET = -1;
	const double DETUNEP = -170;			// trap and probe detunings in units of GA
        const double RNU = pow(-C3/(-HBAR*GA),0.33333); // critical radius for a given potential (J-V paper)
	const double RT = pow(-C3/(HBAR*DETUNET*GA),0.33333);
	const double RP = pow(-C3/(HBAR*DETUNEP*GA),0.33333);
      	const double DB = 1e-9;          		// impact parameter steps in m
      	const double BMAX = RT;
	const double DT = 0.1e-9;
	const double TLIMIT = TAU*20;         		// time step in seconds, time limit
	const double DV = 0.01;                      	// velocity steps in m/s
	const double VMAX = 0.25;			// velocity cutoff
	const double IT = 5; 				// I/Isat trap laser
	const double IP = 1000;                      	// I/Isat probe laser
	const double MT = A0 + A2*RT*RT/(RNU*RNU) + A4*RT*RT*RT*RT/(RNU*RNU*RNU*RNU);
	const double MP = A0 + A2*RP*RP/(RNU*RNU) + A4*RP*RP*RP*RP/(RNU*RNU*RNU*RNU);  // molecular corrections for landau-zener factors
	double prob[300];				// array size
	double erob[300];				// array size
	double rrob[300];				// array size
	double trob[300];				// array size
	double brob[300];				// array size
	temp = 50e-6;					// temperature in K
        cout << RT << " " << C3 << " MT=" << MT << " MP=" << MP << "\n";
	pts=0;
	x=0;
	while (x < 300) {
		prob[x]=0;
		erob[x]=0;
		rrob[x]=0;
		trob[x]=0;
		brob[x]=0;
		++x;
	}
//***********************************************************************************************
	lz = PI*HBAR*GA*GA/(18*C3);                                       // used in calculating LZ params
	vi = DV;
        cout << "starting\n";
	while (vi < VMAX){
		cout << vi << "  ";
		energy = M*vi*vi/4 - C3/(RT*RT*RT);                       //use reduced mass
		bi = DB;
		while (bi < BMAX){
			lztrap = 1 - exp( - lz * MT * MT * RT*RT*RT*RT * IT / (vi * sqrt(RT*RT - bi*bi)/RT)); // landau-zener factor for trap laser
			Pv = (M/(2*K*temp)) * vi * exp(-vi*vi*M/(4*K*temp)) * DV; //probability for this velocity group
			Pb = 2*bi*DB/(RP*RP);                             //probability for this impact parameter
			dPred = (bi < RP ? exp ( - lz * MP * MP * RP*RP*RP*RP * IP / (vi * sqrt(RP*RP - bi*bi)/RP)) : 1); // landau-zener factor for probe laser alone - used to subtract off later
			l = M*vi*bi/2;
			c2 = l*l/M;
			r = RT;
                        t = 0;
                        pop = 1;			//reset ex. st. population
			while (t < TLIMIT && r > RP && energy + C3/(r*r*r) > c2/(r*r)){
				r += - DT * sqrt(4/M) * sqrt(fabs(energy + C3/(r*r*r) - c2/(r*r)));
				v = sqrt((4/M)*(C3/(r*r*r) + energy));
				b = vi*bi/v;
				gm = GA*(A0 + A2*(r*r/(RT*RT)) + A4*(r*r*r*r/(RT*RT*RT*RT)));
                                pop = pop * (1 - DT * gm);
				Pt = pop * gm * DT;
				if (b < RP){
					lzprobe = 1 - exp ( - lz * MP * MP * RP*RP*RP*RP * IP / (v * sqrt(RP*RP - b*b)/RP)); // landau-zener factor for probe laser
					dP = dPred*Pv*Pb*Pt*lztrap*lzprobe;
					delay = t + (sqrt(r*r - b*b) - sqrt(RP*RP - b*b))/v;

					x = delay * 1e8;
      					if (x < 300 && x > 0){
						prob[x] += dP;
					}
					x = 100*(M/4)*(v*v-vi*vi)/(HBAR*GA);
					if (x < 300 && x > 0){
						erob[x] += dP;
					}
					x = (RT-r)*1e9;
					if (x < 300 && x > 0){
						rrob[x] += dP;
					}
					x = t*1e9;
					if (x < 300 && x > 0){
						trob[x] += dP;
					}
					x = bi*1e9;
					if (x < 300 && x > 0){
						brob[x] += dP;
					}
				}
				t += DT;
			}
			bi += DB;
		}
		vi += DV;
	} 
	x=0;
	integral=0;	
	while (x < 300){
		f1 << x << " " << prob[x] << "\n";
		f2 << x << " " << erob[x] << "\n";
		f3 << x << " " << rrob[x] << "\n";
		f4 << x << " " << trob[x] << "\n";
		f5 << x << " " << brob[x] << "\n";
		integral += prob[x];
		++x;
      	}
	cout << "\n integral = " << integral << "\n";
	// this section for reporting on the total runtime
	second = time(NULL); /* Gets system time again */
	cout << "runtime: " << difftime(second,first) << "\n";
	cout << pts << " points\n";	
	return(0);
}