// This program solves the time dependent equations for
// a 2 state system density matrix.
//
// The calculations have the amplitude for the E-field be constant
//
// This program uses the leapfrog algorithm 2nd order Runge-Kutta algorithm
// F(t+dt/2) = F(t) + (dt/2) F'(t)
// F(t+dt) = F(t) + dt F'(t+dt/2)
//

#include <stdio.h>   // standard IO
#include <string.h>  // string manipulations
#include <math.h>    // alternate math functions
#include <complex>   // package for doing complex manipulations
#include <cmath>

// define the type called couble complex
using namespace std ;
typedef complex<double> dcx;

//main program
int main()
{
dcx c0[2], crotw[2] ;
long int numtim, nmod ;
double tau, rabomeg, detun, dt, tstr, tfin, tim, p1, p2,
       genrab, gamma, dmc, dms, rhotil[3], drhotil[3], rhowrk[3] ;

//define the complex number i
const complex <double> icom (0.0,1.0) ;
//define the number pi
const double pi = double(2)*asin(double(1)) ;

//the output will be in whatever is defined here
FILE *output ; output = fopen("densmat_const_2.0_1.0_2.0_cpp.txt","w") ;
FILE *output2; output2= fopen("blochsp_const_2.0_1.0_2.0_cpp.txt","w") ;

// the electric field is assumed proportional to exp(-t^2/tau^2) cos(omega t)
tau = 0.1e-6 ;
// define Omega
rabomeg = 2.0*2.0*pi/1.0e-6 ;
// define the detuning
detun = 1.0/tau ;
// define the general Rabi frequency
genrab = sqrt(detun*detun + rabomeg*rabomeg) ;
// define the decay rate in units of 10^6 s^-1
gamma = 2.0/1.0e-6 ;
      
// set the time step (you should check convergence with respect to dt)
dt = 0.0001e0*2.e0*pi/genrab ;
// define the starting (tstr) and final time (tfin)
tstr = 0.0 ;
tfin = 15.0*2.0*pi/sqrt(detun*detun + rabomeg*rabomeg) ;
// define the number of steps in the numerical integration
numtim = (((tfin-tstr)/dt)) ;
// initialize the coefficients in the two state calculation
c0[0] = dcx(1.0,0.0) ;
c0[1] = dcx(0.0,0.0) ;
// initialize the coefficients in the two state calculation
//    see page 8 of the Chap7.pdf for definition in terms of coefficients
//    rhotil[0] = rhotilde_x, rhotil[1] = rhotilde_y, rhotil[2] = rhotilde_z
rhotil[0] = 0.5*real(c0[0]*conj(c0[1]) + conj(c0[0])*c0[1]) ;
rhotil[1] = 0.5*real((c0[0]*conj(c0[1]) - conj(c0[0])*c0[1])*icom) ;
rhotil[2] = 0.5*(norm(c0[0]) - norm(c0[1])) ;
      
// nmod is roughly how many steps to skip so the number of output points
//    are the denominator
nmod = numtim/2000 ;
      
// the next loop is the time steps
//  because of how the leapfrog algorithm works, the times are stepped in 2's
 for(long int j = 0 ; j <= numtim ; j++)
 {
// define the time for step j
 tim = tstr + double(j)*dt ;
// compute the drhotilde/dt at time tim
 drhotil[0] =-detun*rhotil[1] - 0.5*gamma*rhotil[0] ;
 drhotil[1] = detun*rhotil[0] - rabomeg*rhotil[2] - 0.5*gamma*rhotil[1] ;
 drhotil[2] = rabomeg*rhotil[1] - gamma*(rhotil[2]-0.5) ;
// do the 1/2 time step and store in a work array
 rhowrk[0] = rhotil[0] + 0.5*dt*drhotil[0] ;
 rhowrk[1] = rhotil[1] + 0.5*dt*drhotil[1] ;
 rhowrk[2] = rhotil[2] + 0.5*dt*drhotil[2] ;
// compute the drhotilde/dt at time tim+dt/2
 drhotil[0] =-detun*rhowrk[1] - 0.5*gamma*rhowrk[0] ;
 drhotil[1] = detun*rhowrk[0] - rabomeg*rhowrk[2] - 0.5*gamma*rhowrk[1] ;
 drhotil[2] = rabomeg*rhowrk[1] - gamma*(rhowrk[2]-0.5) ;
// do the full time step
 rhotil[0] = rhotil[0] + dt*drhotil[0] ;
 rhotil[1] = rhotil[1] + dt*drhotil[1] ;
 rhotil[2] = rhotil[2] + dt*drhotil[2] ;

// rhotil is now at time given by tim+dt
 tim = tstr + double(j+1)*dt ;

// every nmod steps output the populations
  if(j%nmod == 0)
  {
// write the time in microsecs, P_1, P_2, P_1+P_2 to the screen
  p1 = 0.5+rhotil[2] ;
  p2 = 0.5-rhotil[2] ;
  printf("%15.7E %15.7E %15.7E %15.7E\n",tim*1.e6,p1,p2,p1+p2) ;
// crotw is the rotating wave approximation to the c
  dmc = cos(0.5*genrab*tim) ;
  dms = sin(0.5*genrab*tim) ;
  crotw[0] = c0[0]*dcx(dmc,-dms*detun/genrab) -
             c0[1]*icom*rabomeg*dms/genrab ;
  crotw[1] = c0[1]*dcx(dmc, dms*detun/genrab) -
             c0[0]*icom*rabomeg*dms/genrab ;
// write the timein microsecs, P_1, P_2, P_1+P_2, P_1RW, P_2RW to the output file
  fprintf(output,"%15.7E %15.7E %15.7E %15.7E %15.7E %15.7E\n",tim*1.e6,p1,p2,p1+p2,norm(crotw[0]),norm(crotw[1])) ;
  fprintf(output2,"%15.7E %15.7E %15.7E %15.7E\n",tim*1.e6,rhotil[0],rhotil[1],rhotil[2]) ;
  }
 }
// end of time loop
 fclose(output) ;
 fclose(output2) ;
return 0 ;
}
